Valerio Olevano

Benoît Sklénard

François Triozon

Valerio

Philippe Blaise

Di Pranzo Carlotta

Anna Lisa

Andrea Giulio

Giovanni

Fernando

Alessandro Sota

Paola

Sara Giulia

Aurelio

Cinzia

Davide Galli

Fabio Finocchi

Philippe Depondt

Cinzia, Ezgi Alessandro Aurelio

Matteo Silvia

Deborah

Bodil, Katya Andria Federico

Camille Augustin

Giovanna Vorrei

Manfredo

Compa

Alberto Charlotte

Martina Cesco

Emi- Lie Fabri

Roberta E Tantissimi

Mia Genitori Massimo E Marina

Sorella Anita

Mie Nonni Mario

Paola E Dina

Mio Zio Marco E Mio Cugino Artur

Tessie

Lilla

Kira Lapo

stati fondamentali sia per lo svolgimento della mia ricerca che sul piano personale.

Introduction

Resistive non-volatile random-access memories based on oxides (OxRAM) have been proposed for their performances as a replacement technology in flash memories and are nowadays developed by several companies. Their operation is based on the reversible modification of intrinsic characteristics of the oxide allowing to change its electrical resistivity. The process is achieved by applying a voltage to form and to break a conductive filament in the oxide, resulting in a low and a high resistive state. OxRAM memories maintain some good properties of flash memories, like non-volatility and rather good scalability, adding two important advantages: the access and programming speed, which is several orders of magnitude higher, and the much lower power consumption. On the other hand, they exhibit a significant variability of resistivity both cycle-to-cycle and cell-to-cell, which is an important drawback for the development of a robust technology. The microscopic mechanisms underlying their working principle, which are responsible for this variability, are not yet fully understood. The object of this thesis is the theoretical study of electronic transport mechanisms in hafnium oxide OxRAM memories. Since these mechanisms are deeply influenced by the device microscopic structure, ab-initio simulation techniques turn out to be suitable for this kind of study, because they permit to precisely control the microscopic structure, and see the effect of small changes on the electronic conductance. The standard procedure to calculate electronic transport properties, which is based on "DFT + Landauer" quantum transport theory, has some disadvantages. Firstly DFT is a ground state theory, which systematically underestimates the band gap in insulators, and seems to overestimate the conductance of bulk junctions, even in the 0-bias limit. Moreover, there is no rigorously exact formulation of DFT for out-of-equilibrium and open systems (i.e. varying number of electrons), preventing calculations at high voltages. Then a practical implementation of the Landauer theory requires a localized functions basis-set, which is suitable for molecules, but not very common and adapted for bulk systems. Finally, the Landauer formalism is rigorously exact Introduction only in the ballistic regime of conduction, losing reliability when scattering events, which are crucial for electronic transport, occur. Therefore the purpose of this thesis is to overcome these problems by performing many-body corrections, like GW calculations, to open the DFT gap, and by applying the non-equilibrium Green function (NEGF) theory to metal-insulator-metal (MIM) systems, and specifically to OxRAM. This PhD work requires a theoretical study and the development of new computational methods, since theories such as GW or NEGF have never been applied to MIM systems electronic transport in literature. Then the application of these new methods to systems of technological interest, like OxRAM, PCRAM, CBRAM, 2D materials or metal-semiconductor contacts, will be useful to improve the understanding and to drive the technological development.

Chapter 1

Memory devices background

Electronic devices are nowadays widespread and fundamental in our everyday life. The performance of all their components is increasing at a high-speed rate to satisfy the demand for more and more challenging tasks. Among these components, memory devices should also evolve not to prevent the development of the final electronic product. The features that a good memory device must posses are the following:

• Data retention, which should be as long as possible to not lose data.

If the data retention is too short the memory device requires a power supply to maintain the information.

• Density of memory cells. It depends on the dimension and shape of the single cell. A high density allows storing more data at a given volume, saving space and reducing costs.

• Writing/reading speed. It can deeply affect the overall performance of the electronic device.

• Endurance, that is the number of cycles the memory device is capable to perform before failure. It should be as high as possible since high endurance translates in a long operating life.

• Production cost, which clearly should be as low as possible.

• Power consumption, which should be minimized.

Real memory devices never possess all these properties together. On the contrary, they always are the result of a compromise. Memory devices are usually classified into two main different categories: volatile and non-volatile, which often posses several complementary features.

Memory devices background

Volatile memory devices

Volatile memory devices are characterized by very short data retention. Therefore they have the disadvantage to require a power supply to periodically refresh the information. On the contrary, in general, the main advantage is that the writing/reading speed is very high. Other advantages are that they have a very good endurance, which permits to overwrite the information for a high number of times, and low power consumption to write and read the memory cells. That given, they are not suitable for big data storage, which would require considerable and continuous power consumption to maintain the information, even when the device is turned off. Moreover some volatile memories, like for instance the SRAM (Static Random Access Memory), also have a very low density, which would imply a waste of space for big data storage. On the contrary, volatile memories are perfectly suitable to quickly store temporary information. Typical applications are for instance CPU cache memory or RAM (Random Access Memory) in the computer architecture.

The CPU cache is usually made of SRAM memories. SRAM uses six transistors per memory cell to store information. The transistors create a bistable latching circuitry (flip-flop) to store each bit. The SRAMs reading/writing speed is very high, in the order of the nanosecond, and the endurance is excellent [START_REF] Natarajan | Searching for the dream embedded memory[END_REF]. However, the complexity of the single memory cell makes them poorly optimized in term of density, and thus more expensive at given storage capacity. The computer RAM is instead usually made of DRAM memories (Dynamic Random Access Memory). DRAM has a simpler structure, made of a transistor in series with a capacitor, which naturally gives rise to two memory states when the capacitor is charged or empty. The main advantage of this technology lies in the simplicity of its structure, which makes it cheap and easy to integrate due to the high density [START_REF] Natarajan | Searching for the dream embedded memory[END_REF]. Other advantages are the very high writing/reading speed, even if lower than SRAM, and excellent endurance. On the contrary, the disadvantage is that the leakage current flowing from the capacitor allows the DRAM to retain the information for a time in the order of milliseconds. Therefore the system has to periodically refresh the memory cell to maintain the information, even when the device is turned on. This is the reason why DRAM is called "dynamic" and SRAM is called "static". Due to the absence of capacitors, SRAM does not need to refresh the information, which is retained as long as the memory is powered. Therefore SRAM consumes much more power when they are operational than when they are idle. Whereas DRAM power consumption is continuous and constant due to the presence of the capacitor and the necessity to refresh information.

1.2. Non-volatile memory devices 5

Non-volatile memory devices

Non-volatile memories are instead characterized by long data retention. The stored information is maintained for long time even when the device is turned off. This advantage, together with a very high density (for some kind of technology, e.g. flash), makes them suitable to store big data. Nevertheless, contrary to volatile devices, they are limited by very low writing/reading speed, worse endurance, and high power consumption. The latter, however, only occurs during the writing/reading operations. For these reasons non-volatile memories are not suitable for rapid and temporary information storage (RAM or CPU cache), because of the great speed difference with modern CPUs and the limited number of cycles before failure.

The most classical and widespread (nowadays together with flash memories) non-volatile memory for big data storage is the hard disk drive (HDD). Invented in 1954 by an IBM team, the HDD consists of a film of ferromagnetic applied on a rotating disk, whose magnetic state can be used as a logical bit. The magnetic state is read and written by two mobile heads, one per side of the disk. Because of its structure, which is composed of several mobile parts, the HDD has notable dimensions, and a very much limited writing/reading speed. The latter, especially due to the mobile write/read heads, cannot go under some tens of milliseconds, which is seven order of magnitude slower than CPUs working frequency. Mainly for this reason, but also to reduce dimensions, in the last years several industries were led to develop a faster, more compact, non-volatile technology.

Flash

Flash memory technology was originally developed by Toshiba in the '80s and commercialized the first time by Intel in 1988. In the last ten years flash memory devices have invaded the market, largely replacing HDD, which remains a good competitor due to its much lower price. Flash technology belongs to the solid state drive (SSD) class of memory devices. The main advantage is that all mobile parts are removed, which entails a remarkable increase in speed performance. The flash writing/reading speed is in the order of tens of microseconds, one thousand times faster than HDD. The structure is also very simple and compact, simplifying the integration of high capacity memories into electronic devices thanks to their high density. The structure consists of a floating-gate MOS transistor, where a floating gate is inserted between transistor control gate and channel [START_REF] Bez | Introduction to flash memory[END_REF] (see fig. 1.1 a)). The latter is isolated from the channel by a tunnel dielectric (usually an oxide), and from the gate by a control dielectric. When a positive voltage is applied to the control gate, The two apparent voltage states used as logical bit. Image adopted from ref. [START_REF] Hubert | Optimisation de mémoires PCRAM pour générations sub-40 nm: intégration de matériaux alternatifs et structures innovantes[END_REF].

electrons are injected from the substrate the channel is made of to the floating gate through the tunnel oxide, via tunnel effect. The electrons trapped into the floating gate shift the apparent threshold voltage seen from the control gate. When a negative voltage is applied electrons are injected back into the substrate via tunnel effect, and the apparent voltage at the control gate shifts also back to its original value. The two apparent voltage states are used as a logical bit (see fig. 1.1 b)).

Besides speed, compared to volatile memories flash devices have much lower endurance, and higher operational power consumption, but also a lower price. Despite the improvement compared to HDD, flash memories have reached a scaling limit [START_REF] Lai | A multilayer stackable thin-film transistor (TFT) NAND-type flash memory[END_REF][START_REF] Tanaka | Bit cost scalable technology with punch and plug process for ultra high density flash memory[END_REF], which prevents to further reduce the memory device dimensions, and the speed performance of volatile memories are still far to be reached. If CPUs work at frequencies in the gigahertz order of magnitude, which means one operation per nanosecond, flash memories can respond to the CPU only after ten thousand nanoseconds. Therefore non-volatile memory devices for big data storage constitute a speed bottleneck in modern electronic devices. For this reason, several companies, universities, and research centers, in the last years are looking for alternatives to flash memories. Figure 1.2 shows the hierarchy of the present memory technologies, pointing out how a large gap of performance and price separates volatile from non-volatile devices. Therefore there is a large space for the development of new devices with in-between characteristics. The ambitious challenge is to develop a non-volatile memory, which may have speed performance close to modern CPUs, that is at least lower than one hundred nanoseconds, maintaining and possibly increasing the flash technology density and endurance.

Figure 1.2: Hierarchy of memories. On top very quick and expensive volatile memories for CPU applications. On the bottom slow, cheap, non-volatile memories for data storage. Adopted from ref. [START_REF] Nail | Investigation of hybrid CBRAM/OXRAM non-volatile memories for low consumption and high reliability[END_REF] Good candidates, which are arising to satisfy these characteristics, are the resistive state based memory devices (RRAM), which use their resistive state to define two logical different states. Their performances promise to be close to those of volatile devices, with the important advantage of being non-volatile (see fig. 1.3). The RRAM category usually refers to two different emerging technologies: OxRAM and CBRAM. Other emerging memories, which in some cases also exploit their intrinsic resistive state as a logical bit, are FeRAM, MRAM, and PCRAM (presented in section 1.4). according to materials employed and especially resistive switching mechanisms. A memory device that uses the resistive switching characteristic of a transition metal oxide is called OxRAM. In the next section, we describe the main features of OxRAMs, which are the memory devices of main interest in this work.

OxRAM

Figure 1.4: OxRAM structure and physics of resistive switching. Image adopted from ref. [START_REF]Functional Metal Oxide Nanostructures[END_REF].

An OxRAM has a very simple structure, made by a transition metal oxide sandwiched between a top (TE) and a bottom (BE) metallic electrode (see fig. 1.4). At pristine state, the oxide has a very high resistance value. By applying a voltage a conductive filament (CF) can be formed or broken, defining a low (LRS) and a high (HRS) resistive state. The most accepted hypothesis in literature is that the CF is made of oxygen vacancies (V O ). Essentially under the applied electric field oxygen atoms can migrate from the CF region towards oxide/electrode interfaces or into the oxide around the CF region. The resulting V O region constitutes the CF which carries current from one electrode to the other. The CF probably connects the two electrodes carrying current directly in the LRS, while it has some ruptures which limit the conduction in the HRS. The first CF forming from the pristine state to the LRS is called "forming", the break of the CF from LRS to HRS is called "reset", while the reforming of the CF from HRS to LRS is called "set" (see [START_REF] Traoré | Étude de cellules mémoires résistives RRAM à base de HfO 2 par caractérisation électrique et simulations atomistiques[END_REF] symmetric (unipolar), or asymmetric (bipolar) (see fig. 1.6). For unipolar devices, both high and low resistive states are programmed by using the same voltage polarity. Whereas bipolar devices require different voltage polarities to program the HRS or the LRS. The polar characteristic mainly depends on the metals used as electrodes. Several binary metal-oxides have been found to exhibit resistive switching effect, where most of the metals are transition metals. Among them TaO x , TiO x , AlO x , NiO and HfO x have raised more research interest [START_REF] Wong | Metal-Oxide RRAM[END_REF]. The met-Figure 1.6: Typical I/V characteristics of unipolar and bipolar devices. (a) Unipolar devices, which have a characteristic that does not depend on the polarisation imposed. (b) Bipolar devices, which have a characteristic that depends on the polarisation imposed. Image taken from [START_REF] Waser | Redox-based resistive switching memories -nanoionic mechanisms, prospects, and challenges[END_REF]. als most widely used in research as electrodes are: Ti, TiN, Zr, Hf, Pt, Ni, Cu, Ir, etc. [START_REF] Wong | Metal-Oxide RRAM[END_REF]. The main interest of this work is on HfO 2 based OxRAMs, which have acquired great interest due to their characteristic: switching speed below 1.3. Resistive memory devices (RRAM) 11 10 ns [START_REF] Govoreanu | 10×10 nm 2 Hf/HfO x crossbar resistive RAM with excellent performance, reliability and low-energy operation[END_REF], very small memory cell dimensions ( 10 × 10 nm 2 [START_REF] Govoreanu | 10×10 nm 2 Hf/HfO x crossbar resistive RAM with excellent performance, reliability and low-energy operation[END_REF][START_REF] Vianello | Resistive memories for ultra-low-power embedded computing design[END_REF]), which translate in very high density and good scalability, high endurance (> 10 8 cycles) [START_REF] Vianello | Resistive memories for ultra-low-power embedded computing design[END_REF], low operational voltage [START_REF] Govoreanu | 10×10 nm 2 Hf/HfO x crossbar resistive RAM with excellent performance, reliability and low-energy operation[END_REF][START_REF] Vianello | Resistive memories for ultra-low-power embedded computing design[END_REF], and very good data retention [START_REF] Zhao | Ultrathin (∼2nm) HfO xas the fundamental resistive switching element: Thickness scaling limit, stack engineering and 3D integration[END_REF][START_REF] Chen | Improvement of data retention in HfO 2 /Hf 1T1R RRAM cell under low operating current[END_REF].

Understanding of OxRAM working principles

Despite all the good qualities exhibited by OxRAMs the understanding of their working principles is still limited. This mainly concerns the origin of the polar characteristic, the CF forming and the switching mechanism, the CF structure, and the conduction mechanism. Surprisingly, the origin of the polar characteristic seems to be independent of the oxide, which means that by using the same oxide different memory devices can observe different polar switching. The origin seems more related to the metals the electrodes are made of. A general trend that is observed is that when the TE is an oxygen reactive metal like Ti or Hf, the switching type is mostly bipolar, while when both TE and BE are symmetrical and less reactive with oxygen the switching type is mostly unipolar. Nevertheless, on one hand, there are non-symmetrical devices (e.g. TiN/HfO 2 /Pt) which exhibit unipolar behavior in some works [START_REF] Wang | Effect of anodic interface layers on the unipolar switching of HfO 2 -based resistive RAM[END_REF], whereas bipolar in other works [START_REF] Zhao | Ultrathin (∼2nm) HfO xas the fundamental resistive switching element: Thickness scaling limit, stack engineering and 3D integration[END_REF]. On the other hand, some symmetric devices (e.g. TiN/HfO 2 /TiN) exhibit bipolar behavior in some works [START_REF] Traoré | Étude de cellules mémoires résistives RRAM à base de HfO 2 par caractérisation électrique et simulations atomistiques[END_REF], whereas unipolar in others [START_REF] Stefano | Semiconducting-like filament formation in TiN/HfO 2 /TiN resistive switching random access memories[END_REF]. Therefore the simple symmetry and oxygen active electrodes argument is probably insufficient to characterize the polarity of the device, which is more complicated and still not fully understood. If the CF hypothesis is widely accepted, the forming and switching mechanisms and the structure of the CF are still under debate. Some works point out the crucial role of grain boundaries in the oxide [START_REF] Bersuker | Metal oxide RRAM switching mechanism based on conductive filament microscopic properties[END_REF][START_REF] Lanza | Grain boundaries as preferential sites for resistive switching in the HfO 2 resistive random access memory structures[END_REF], which constitute leakage paths that ultimately transform into CF. Under voltage and temperature V 2+ O migrate and segregate at the grain boundaries, where they turn into V 1+ O by capturing electrons from the cathode. The released O 2-ions migrate towards the anode, and the conduction through the V 1+

O vacancies is achieved via trap assisted tunneling according to:

V 1+ O + e -→ V 0 O → V 1+ O + e -
This model is also similar to the one proposed by Vandelli et al. [START_REF] Vandelli | Comprehensive physical modeling of forming and switching operations in HfO 2 RRAM devices[END_REF]. Zhao et al. [START_REF] Zhao | Dynamic modeling and atomistic simulations of SET and RESET operations in TiO 2 -based unipolar resistive memory[END_REF] and Jeong et al. [START_REF] Jeong | Characteristic electroforming behavior in Pt/TiO 2 /Pt resistive switching cells depending on atmosphere[END_REF], in their work on Pt/TiO 2 /Pt, propose another model. V O migrate towards the cathode and grow into a cone-shaped CF made of Ti 4 O 7 with the larger side at the cathode. The oxygen ions are removed from the oxide and absorbed near the anode interface at the grain boundaries of the Pt electrode. The Ti ions in the oxygen deficient region might capture electrons from the cathode to fulfill the local charge neutrality:

T i 4+ + ne -→ T i (4-n)+
These reduced valence Ti states can form a metallic phase in the form TiO 2-n/2 when n > 1.5 approximately. According to Waser and Aono [START_REF] Waser | Nanoionics-based resistive switching memories[END_REF], the conductive TiO 2-n/2 acts as a "virtual cathode" closer to the anode. The only widely accepted point about the CF, common to these models, is that it is rich in oxygen vacancies. Their charge state and type are still under debate. Moreover, morphology and composition of the CF are still unknown.

The hypothesis according to which the metal electrode could migrate and diffuse into the oxide during the CF formation cannot be put aside. Therefore it remains very unclear whether the CF is a rich metal cation precipitate, a sub-oxide, a mixture of sub-oxide and electrode atoms, etc.

Regarding the conduction mechanism, in general, the LRS observes an Ohmic type of conduction, while there are controversies about the HRS. Among the 8) Tunneling from trap to anode. Adopted from ref. [START_REF] Wong | Metal-Oxide RRAM[END_REF].

conduction mechanisms proposed for the HRS there are Poole Frenkel [START_REF] Walczyk | Pulse-induced low-power resistive switching in HfO 2 metalinsulator-metal diodes for nonvolatile memory applications[END_REF], space charge limited current [START_REF] Lee | Evidence and solution of over-RESET problem for HfO X based resistive memory with sub-ns switching speed and high endurance[END_REF], trap assisted tunneling, Schottky emission [START_REF] Wei | Highly reliable TaO x ReRAM and direct evidence of redox reaction mechanism[END_REF]. Figure 1.7 shows the different possible conduction mechanisms in OxRAMs. The predominance of one type compared to another may be due to several parameters, like electrode/oxide interface properties, fabrication process, defect density, trap energy levels, etc. [START_REF] Wong | Metal-Oxide RRAM[END_REF].

Resistive memory devices (RRAM)

Present drawbacks

Besides the limited knowledge about the OxRAMs working principles, the main operational drawback, which is delaying a mass fabrication for commercial purposes, is the non-uniformity of switching parameters. Devices show an exponential dependence of forming time on the applied voltage, which indicates that the forming process is not a single simultaneous event, but a process resulting from an upsurge of defects [START_REF] Buh | Time-dependent electroforming in NiO resistive switching devices[END_REF][START_REF] Wong | Metal-Oxide RRAM[END_REF]. Therefore, the forming process is naturally stochastic, and the same macroscopic parameters of current and voltage can lead to different CF structures, and thus different resistance values, both device-to-device and cycle-to-cycle. Indeed, for instance, any small variation of diameter in the CF can change the electronic transport in the LRS, while any small variation in the thickness of the CF ruptures can influence the HRS. Another typical drawback, possibly related to variability, is the LRS failure, which gets stuck to the HRS. This can often happen after a rather limited number of cycles, which is 10 4 or 10 5 . According to Chen et al. [START_REF] Chen | Understanding of the endurance failure in scaled HfO 2 -based 1T1R RRAM through vacancy mobility degradation[END_REF] this failure could be attributed to a progressive degradation of V O diffusion barrier in the CF due to high local temperature. However, the failure could also be due to the reinforcement of the interfacial layer between electrode and oxide, making more difficult the "set" operation. Sometimes the failure occurs also for the HRS, which remains stuck in the LRS. According to Wong et al. [START_REF] Wong | Metal-Oxide RRAM[END_REF] the reason could be the generation of too many V O at or near the electrode/oxide interface, near the CF, or in the oxide in general. Nevertheless, a full comprehension of the failure phenomena is still far to be reached. Nowadays also price, which is very high, constitutes an important drawback for commercial purposes. However, once the technology will be mature, a mass fabrication could strongly reduce the prices.

CBRAM

Conductive bridge RAM (CBRAM) memories have a structure similar to OxRAMs: the cell is composed of a solid electrolyte (e.g. GeS 2 [START_REF] Gopalan | Demonstration of conductive bridging random access memory (CBRAM) in logic CMOS process[END_REF]) sandwhiched between two metal electrodes, one of which is electrochemically active. If the solid electrolyte is in its original initial state the CBRAM is in the HRS. Applying a positive voltage to the active electrode a CF is formed essentially made of active metal electrode atoms in the solid electrolyte, connecting the two electrodes. Thus the device turns into its LRS. To bring it back to the HRS a negative voltage is applied to the active electrode. The result is the dissolution of the electrochemical active metal in the CF. Similarly to OxRAM memories, low voltage, high speed, and high density are some of [START_REF] Celano | Three-dimensional observation of the conductive filament in nanoscaled resistive memory devices[END_REF].

the performance characteristics demonstrated by CBRAM devices. However, they observe high temperature instability [START_REF] Vianello | Sb-doped GeS 2 as performance and reliability booster in conductive bridge RAM[END_REF], which requires a deeper understanding of the working principles to pursue technological development. Figure 1.8 shows the CBRAM devices' structure, switching mechanisms, and operational current and voltage values.

Other emerging memories

Among the emerging non-volatile memory technologies we briefly present the phase change memory (PCM), which seems to be the most mature technology nowadays, the Ferroelectric RAM (FeRAM), the Magnetoresistive RAM (MRAM), and the Mott insulator memory. .9: Sketch of structure (left) and programming electrical power (right) of a PCM. Adopted respectively from [START_REF] Pershin | Memory effects in complex materials and nanoscale systems[END_REF] and [START_REF] Burr | Phase change memory technology[END_REF].

A PCM uses chalcogenide based materials, such as GeSbTe, which can exist at both crystalline and amorphous phases. The two phases have different resistance values, which can be used as a logical bit. Resistance is low in the crystalline phase and increases even of many orders of magnitude in the amorphous phase [START_REF] Hubert | Optimisation de mémoires PCRAM pour générations sub-40 nm: intégration de matériaux alternatifs et structures innovantes[END_REF]. The structure of a PCM is sketched in figure 1.9 (left). The materials PCMs are made of are very well characterized since they were largely adopted for the development of CDs and DVDs. This fact has significantly sped up the development of PCM, which is probably the most mature RRAM technology nowadays. The switching between LRS and HRS is driven by heat, which is induced by applying an electrical pulse on the active region (see fig. 1.9 left). The electrical pulse is shorter and at higher voltages to pass from crystalline to amorphous phase, that is from LRS to HRS, than vice versa (see fig. 1.9 right). This is required to melt and quench the active region of the memory cell. The main limitation of PCM cells is the high current required to program the HRS. Although some improvements have been achieved [START_REF] Hubert | Carbon-doped Ge 2 Sb 2 Te 5 phase-change memory devices featuring reduced reset current and power consumption[END_REF][START_REF] Lacaita | The race of phase change memories to nanoscale storage and applications[END_REF], the high HRS switching current translates into larger memory selector regions, which occupies the most layout area of the cell [START_REF] Wong | Phase change memory[END_REF], limiting the integration density of the devices. PCM also observe some reliability issues to solve [START_REF] Ielmini | Phase change materials in non-volatile storage[END_REF].

FeRAM

FeRAM technology is based on two stable polarization states proper to ferromagnetic materials, generally made of zirconate titanate [START_REF] Meena | Overview of emerging nonvolatile memory technologies[END_REF]. FeRAM devices possess an intrinsic dipole, which switches when applying an electric field resulting in two distinguished resistance states (see fig. is a rather old technology, possessing high read/write speed and low power consumption, but a very low density and high price.

MRAM

The storage in MRAM devices is based on the change of the magnetization state of a free layer (FL) associated with another reference layer (RL) [START_REF] Julliere | Tunneling between ferromagnetic films[END_REF]. The FL and the RL are separated by an insulating barrier in a magnetic tunnel junction (MTJ), whose dimensions are below 100 nm. When the FL magnetization state (majority electron spins) is parallel to the RL magnetization spin the resistance of the MTJ is low, whereas when the magnetization states are anti-parallel the resistance is high (see fig. 1.11). MRAM is very Figure 1.11: MTJ junction showing a) parallel and b) antiparallel spins corresponding to low and high resistance states of an MRAM, respectively. Adopted from ref. [START_REF] Traoré | Étude de cellules mémoires résistives RRAM à base de HfO 2 par caractérisation électrique et simulations atomistiques[END_REF].

fast, similarly to DRAM and RRAM, and has very good endurance, com-1.5. Conclusions about OxRAM devices parable with DRAM or SRAM. Early MRAM suffered from very high write current. However, the recent development of Spin Torque Transfer MRAM (STT-MRAM [START_REF] Khvalkovskiy | Basic principles of STT-MRAM cell operation in memory arrays[END_REF]) has largely improved the current levels. Nowadays the main drawback remains the integration density, which is much lower than DRAM or flash devices [START_REF] Gallagher | Emerging nonvolatile magnetic memory technologies[END_REF].

Mott insulator memory

Mott insulators are a class of materials that are predicted to behave as conductors according to conventional band theories, but in fact behave as insulators (in particular at low temperatures). This effect is due to electron-electron interactions, which are not accounted for by conventional band theory. By changing some parameters, which may be composition, pressure, strain, voltage, or magnetic field, Mott insulators can turn to conductors. This effect, known as a Mott transition, can be exploited to engineer memory devices, which take advantage of the two different conductive states as logical bit [START_REF] Son | A heterojunction modulation-doped mott transistor[END_REF][START_REF] Zhou | Correlated electron materials and field effect transistors for logic: A review[END_REF]. In this chapter we have outlined the emerging OxRAM memory technology, detailing its several interesting qualities, which make them a promising candidate to replace flash memories in massive data storage applications. Recently OxRAM devices have also been proposed as embedded memories mimicking synapses at the nanometer scale in convolutional neural networks Memory devices background [START_REF] Garbin | HfO 2 -based OxRAM devices as synapses for convolutional neural networks[END_REF]. Indeed they seem to possess very interesting qualities for this kind of application: they offer the possibility to modulate their conductance by applying low biases, they can be easily integrated with CMOS-based neuron circuits, they have good data retention and low power consumption. Another related, active, and very recent field of research tries to develop multi-level OxRAM cells for fast training in neuromorphic applications [START_REF] Ramadan | Adaptive programming in multi-level cell reram[END_REF][START_REF] García-Redondo | Self-controlled multilevel writing architecture for fast training in neuromorphic RRAM applications[END_REF]. However, this chapter has also pointed out the present drawbacks of OxRAM technology, which is mainly resistive state variability and failure. To continue the development in order to engineer and fabricate OxRAM devices for a wide market, the present drawbacks have to be solved. This purpose requires making intensive efforts to achieve a full understanding of the working principles in OxRAM devices. Deep researches have still to be done about the properties of the CF: structure and mechanisms of conduction, forming and switching above all. The focus of this work is mainly on the conduction mechanism and properties, which we choose to study taking advantage of ab initio simulations. Ab initio simulations make us able to perfectly control the CF structure and morphology, which are unknown, and deeply influence the electrons conduction in OxRAMs. Therefore, by ab initio simulations, we can study the impact of different guessed forms of CF on the conduction. The next chapter will give a complete overview of the present "state of the art" theories to compute the electronic conduction properties of bulk OxRAM-like structures by ab initio simulations.

Conclusions about OxRAM devices

Chapter 2

Electronic transport state of the art

In this chapter we present the state of the art theoretical context for the calculation of electronic transport properties, and specifically the conductance, in systems of mesoscopic and nanoscopic dimension, characterizing the most used approach based on "DFT + Landauer", and outlining its limits of applicability. The first section describes the Landauer formalism, while the second shows how DFT, combined with the use of the Green function formalism, can provide the required input terms to apply the Landauer formula. All the concepts presented in this chapter are well described in two books of reference by S. Datta [START_REF] Datta | Electronic Transport in Mesoscopic Systems[END_REF] and M. Di Ventra [START_REF] Ventra | Electronic Transport in Nanoscale Systems[END_REF].

Landauer formalism

Consider a system composed of a conductor placed between two larger contacts, as in figure 2.1. If the dimensions of the conductor were macroscopic then its conductance G, which is defined as the inverse of the resistance G = R -1 , would be given by an Ohmic linear equation:

G = σ S L (2.1)
where σ is the conductivity of the material the conductor is made of, which is independent of the sample dimensions, L is the length of the conductor along the transport direction, which goes from one contact to the other, and S is the conductor surface perpendicular to the transport direction.

When the conductor goes to mesoscopic or smaller dimensions there are two experimental observations which are not predicted by law 2.1: 1. The measured conductance G approaches a limiting value when the length of the conductor L becomes shorter than the electrons mean free path. On the contrary equation 2.1 predicts the conductance to grow indefinitely as L is shortened.

2. As the surface S of the conductor is reduced the conductance goes down in discrete steps, each of height 2e 2 /h. Whereas, according to equation (2.1), the relation between G and S is strictly linear.

Rolf Landauer, mainly in three works published in 1957 [START_REF] Landauer | Spatial variation of currents and fields due to localized scatterers in metallic conduction[END_REF], 1970 [START_REF] Landauer | Electrical resistance of disordered one-dimensional lattices[END_REF], and 1985 [START_REF] Büttiker | Generalized manychannel conductance formula with application to small rings[END_REF], proposed a simple and very general formula which incorporates both of these features, based on an elegant interpretation of the two experimental phenomena.

Firstly he relates the conductance to the number of conductive modes into the conductor, each carrying a quantum of conductance G 0 = 2e 2 /h = 7.75 • 10 -5 Ω -1 . The conductive modes are the electronic states where electrons can propagate along the transport direction, and contribute to the conductance, or the current if a bias is applied. A general result of solid state physics is that only partially filled states carry conductance, whereas filled states do not conduct, because electrons cannot propagate (Ashcroft and Mermin, 1975 [61]). Therefore, considering for simplicity zero temperature, so that the electrons occupation is given by a step function, a conductive mode is a band of the dispersion relation E( k) which crosses the Fermi energy, where k indicates a wave vector along the transport direction (figure 2.2). Secondly Landauer considers that the total conductance should be proportional to the probability of an electron in a conductive mode to be transmitted from one contact to the other. This probability can be taken into account by simple transmission and reflection coefficients for electrons, which can be elastically reflected by the contacts, but propagate in ballistic regime into the conductor. On the contrary, if electrons are subjected to inelastic scattering into the conductor, a more complicated microscopic theory is required to account for the transmission probability, because scattering events would imply complicated quantum mechanics transition rules governed by the Pauli principle. For most of the mesoscopic and nanoscopic samples, the electrons mean free path is longer than the conductor length. Hence Landauer assumes that the electronic transport into the conductor is ballistic, and accounts for the transmission probability by a simple coefficient. Therefore, leaving the demonstration in the next section, the Landauer formula at zero temperature is given by:

G = G 0 M (E f ) i T i (E f ) = G 0 M (E f ) T (E f ); G 0 = 2e 2 h (2.2) Equation (2.
2) perfectly represents the Landauer interpretation of conductance. There are M (E f ) conductive modes at the Fermi energy, each contributing to the conductance with a quantum of conductance G 0 . This value Electronic transport state of the art is then normalized by the average probability T = T i /M for an electron in a conductive mode to be transmitted from one contact to the other. Remark that for each conductive mode the correspondent transmission probability T i can, in general, be different from the others.

From this new perspective the interpretation of the two experimental phenomena not predicted by the Ohm law 2.1 becomes clearer:

1. Since electrons are fermions, so that they cannot occupy the same quantum state of another electron, and since the number of states carrying conductance is limited, the conductance observes a maximum intrinsic value equal to G 0 • M . This intrinsic limit is high and practically never reached for macroscopic conductors, which have a very high number of conductive modes M, and low for mesoscopic or smaller conductors, which have a very low M value. As long as the conductor length L shortens the scattering events become less frequent, the resistance encountered by electrons drops, and the conductance increases. When L becomes shorter than the electron phase-breaking mean free path scattering events stop to occur, the conductive regime becomes ballistic, and equation 2.2 rigorously holds. In ballistic regime T only depends on the interface geometry between conductor and contacts, which in general can reflect electrons back into the conductor. Therefore T becomes independent of L, and the conductance doesn't increase any more as L is further shortened.

2.

Varying the conductor surface corresponds to vary M. If the surface has macroscopic dimensions, M is a high number (often in the order of 10 3 or 10 4 ), and a variation of one single mode over thousands of modes is seen as a continuous change in the conductance. On the contrary, if the surface has mesoscopic or nanoscopic dimensions, M is a low number (often in the order of unity), and a variation of one single mode over few modes is seen as a discrete change in the conductance.

Equation 2.2 can be generalized to non-zero temperatures. In this case the electronic states occupation is not given by a simple step function, but by an occupation function η(E) whose step is smoothed by the thermal energy k B θ. Due to thermal energy the states which carry conductance, where electrons can propagate, are not only the states crossing the Fermi energy, but all the states in an energy range

E f -few k B θ < E < E f + few k B θ.
All the states of energy E such that η(E) = 0 and η(E) = 1, that is such that ∂ η(E)/∂E = 0, are partially filled, and contribute to the conductance increasing the number of conductive modes (figure 2.3). Therefore the gen-eralized formula turns out to be:

G = G 0 +∞ -∞ - ∂ η(E) ∂E M (E) T (E) dE (2.3)
As shown in figure 2.3, the term -∂η(E)/∂E selects the conductive modes and weighs their contribution to the conductance. In practice, for the purpose of this work, the Fermi distribution is always a reasonable choice for the electron occupation function:

η f (E) = 1 exp E-E f k B θ + 1 (2.4)
Therefore it is easy to prove that the weighting factor of equation (2.3) is given by: Equation 2.5 is plotted in figure 2.4. In the limit of θ → 0 it reduces to a Dirac delta centered on E f , and equation 2.2 is recovered. Increasing the temperature, it broadens, thus selecting a wider range of energies around E f . Equation 2.5 is always normalized, indeed:

- ∂ η f (E) ∂E = n f (E)(1 -n f (E)) k B θ (2.5)
+∞ -∞ - ∂ η f (E) ∂E dE = -η f (E) +∞ -∞ = 1
Thus it can be seen as a probability distribution which gives the probability for a state of energy between E and E +dE to contribute to the conductance. It is interesting to remark that, independently of the θ value, equation 2.5 is very similar to a Gaussian distribution of mean value equal to E f , and width given by a standard deviation σ 1.6 • k B θ (figure 2.4):

g(E) = 1 σ √ 2π exp - (E -E f ) 2 2σ 2 ; σ 1.6 k B θ (2.6)
Hence, since for a Gaussian distribution the range of 3 σ around the mean value represents about the 99% of the total probability, the states in a range of approximately

E f ± 3 • 1.6 k B θ E f ± 5 k B θ
constitute the total number of states contributing to the conductance. This is confirmed by the fact that:

+5 k B θ -5 k B θ - ∂ η f (E) ∂E dE = -η f (E) +5 k B θ -5 k B θ

∀θ

That given, if M (E)T (E) is constant and independent on energy on a range

E f ± 5 k B θ, equation 2.3 becomes: G = G 0 M (E f )T (E f ) +∞ -∞ - ∂ η f (E) ∂E dE = G 0 M (E f )T (E f )
which is the Landauer equation at zero temperature. Hence the conductance becomes independent of temperature. For low temperatures 5 k B θ is small and function 2.5 is very narrow around E f . Thus M (E)T (E) turns out to be almost constant on the small range E f ± 5 k B θ, and the conductance can be approximated with its zero-temperature value. At room temperature (θ 300 K) 5 k B θ 0.13 eV , which is a rather narrow energy range. Therefore, except for very rapid variations around E f of M (E)T (E), the room temperature conductance can be approximated by the zero temperature conductance: The current I → flowing through the conductor in a given direction, from left to right or from right to left, is carried by N electrons provided by one of the reservoirs, with average velocity v. These electrons have a mean free path which is assumed to be longer than the conductor length L, and have average probability 1 -T to be reflected by the lead back into the conductor. Therefore, by definition of electric current, and considering that reflected electrons do not contribute to the current:

G(E; θ = 300 K) G(E; θ = 0 K)

Current and proof of Landauer equation

I → = ∂Q ∂t T = ∂Q ∂L ∂L ∂t T = N e L v T (2.7)
where Q is the charge. From equation 2.7, to compute the current, the required terms are the transmission coefficient, the total number of electrons, and the velocity. The latter, according to the Bloch theory of electrons in solids [START_REF] Ashcroft | Solid State Physics[END_REF], for an electron of momentum h k, is simply:

v k = 1 h ∂E( k) ∂ k (2.8)
where E( k) is a Bloch band. The electrons into the conductor are organized in a band structure, which provides the permitted states for the electrons.

Thus the total number of electrons can be calculated, per unit of energy, as the product of density of states D(E) and occupation function η(E) times 2 (for spin degeneracy). Since the problem is one-dimensional, the density of states for a single momentum h k is given by:

D k (E) dE = L 2π d k → D k (E) = L 2π ∂ k ∂E = L h v k (2.9)
while the occupation function η(E) can be left general, without specifying its form. Finally every electron of momentum h k is characterized by a transmission coefficient T k (E), whose calculation will be the object of the next section of this chapter. Now that all the needed terms are specified, to calculate the current we have to sum over all the available momentum h k per each energy (equal to M (E)), and then integrate over energy:

I → = 2e L M (E) k D k (E) v k T k (E)η(E -µ)dE = 2e h M (E)T (E)η(E -µ)dE (2.10)
where T (E) = T k (E)/M (E) is the average transmission probability as function of energy, µ is the chemical potential of the reservoir providing the electrons, and the second equality is obtained by using equation 2.9. The net current flowing through the conductor is the difference between the current I LR flowing from left to right and the current I RL flowing from right to left. If the two leads are at a different chemical potential µ L and µ R , with for instance µ L > µ R , from equation 2.10 the net current is:

I = I LR -I RL = 2e h M L (E)T L (E)η(E -µ L ) -M R (E)T R (E)η(E -µ R )dE
(2.11) Note that in equation 2.11 we have defined two different function M L and M R for the number of modes close to the left and the right contact. This remarks the fact that the number of modes is evaluated at different energies, that is µ L and µ R , close to the two different contacts. Therefore they can in general be different. Nevertheless, this is not strictly necessary, since into the conductor there is a single band structure, and the number of modes at a given energy is given by a unique function M (E). In absence of bias, with µ L = µ R , the current from left should perfectly balance the current from right, resulting in zero net current. This can be assured only if:

M L (E)T L (E) = M R (E)T R (E) = τ (E) (2.12)
As a consequence of the continuity equation for the density charge flux through a surface, which is always conserved, equation 2.12 always holds, also under bias or out of thermal equilibrium. This property also rigorously follows from the fact that the scattering matrix, which relates the outgoing wave amplitudes to the incoming wave amplitudes at different leads, is unitary (see section 2.2.1). Remark that equation 2.12 is in general not true if the electrons into the conductor undergo inelastic scattering events , but this possibility is here assumed to be prevented. As consequence of equation 2.12 equation 2.11 can be further simplified: the Landauer equation. If the applied bias is small, that is |µ L -µ R | 1, the difference between the occupation functions may be expanded by means of the following Taylor expansion:

I = 2e h [η(E -µ L ) -η(E -µ R )]τ (E)dE (2.
η(E -µ L ) -η(E -µ R ) - ∂η(E -E f ) ∂E [µ L -µ R ] (2.14)
where:

E f µ L + µ R 2 (2.15)
By means of equation 2.14, equation 2.13 can be rewritten as:

I = 2e 2 h - ∂η(E) ∂E τ (E) dE [µ L -µ R ] e (2.16)

Landauer formalism

29

which is in the Ohmic linear form I = GV , where V = [µ L -µ R ]/e, and:

G = 2e 2 h - ∂ η(E -E f ) ∂E τ (E) dE (2.17)
which is the Landauer equation 2.3 for non-zero temperature. If θ = 0 then η(E) becomes a step function, and equation 2.13 turns out to be:

I = 2e h µ L µ R τ (E) dE 2e 2 h τ (E f ) µ L -µ R e (2.18)
where the second equality again holds under small bias hypothesis. By linearity it immediately follows that G = G 0 τ (E f ), which is the Landauer equation 2.2 at zero-temperature.

Limits of applicability

The derivation of the Landauer equation in the previous section has required two approximations. The first assumes that the electrons mean free path is longer than the conductor length, which allows to describe the transmission probability as a simple coefficient. The second supposes that the applied voltage is small, which allows to linearize the relation between current, conductance and voltage.

For macroscopic samples the electrons mean free path is always much shorter than the conductor length, and many inelastic scattering events can occur, for instance among electrons or electrons and phonons. These events give rise to the electric resistance, which can be accounted for by an average coefficient, without caring about the microscopic details. Nevertheless, in nanoscopic samples an average description becomes inadequate (see for instance the predictive failures of law 2.1), and the microscopic behavior of each electron carrying current has to be characterized. This task becomes hard in presence of inelastic scattering events, which are extremely complicated to describe singularly mainly because of the Pauli exclusion principle. For this reason Landauer introduces the hypothesis, which in many cases is justified, of ballistic conduction. Indeed the form of equation 2.

f (x + h) -f (x -h) = f (x) 2h + O(h 2 )
Consider also that the exponential argument in the Fermi distribution can be rewritten as:

E -µ L,R k B θ = E k B θ - µ L,R k B θ
Therefore, for the difference η(E -µ L ) -η(E -µ R ) in equation 2.13, the increment h can be identified by [µ L -µ R ]/2k B θ, and the truncation error Err from equation 2.14 turns out to be:

Err = O µ L -µ R 2 k B θ 2 (2.19)
Hence a good criterion to assure the linearity of I(V ) is:

|µ L -µ R | 2 k B θ (2.20)
There are also situations where the relation could be linear also for large bias and for |∆µ| 2 k B θ. This case occurs when τ (E) is constant over an energy range from µ R -few k B θ to µ L + few k B θ. Consider firstly zero temperature, if τ is constant equation 2.13 becomes:

I = 2e h τ (E f ) η(E -µ L ) -η(E -µ R ) dE = 2e 2 h τ (E f ) µ L -µ R e (2.
η f (E -µ L ) -η f (E -µ R )
shown for θ = 0 (red line) and θ > 0 (blue line). The areas under the two functions, that is the integrals, are the same equal to µ L -µ R . Thus equation 2.21 holds for every temperature.

the current by equation 2.13. This correction should in general accounted for by taking into account the electron density inside the conductor and obtaining the electric field from the Poisson equation (see section 2.2.4 for more details). On the contrary, in the linear response regime, this aspect of the problem can be neglected, and τ can be computed at 0-bias.

To conclude, the Landauer formalism rigorously holds as long as the conductor size along transport remains smaller than the electron mean free path, or, in other words, as long as the electrons transport is ballistic. This is a mandatory precondition to guarantee the validity of equations 2.13 and 2.3 for current and conductance derived in this section. In presence of small bias (|∆µ| 2k B θ), or peculiar slowly varying transmission function τ (E) in the range ∆µ, the current is given by the linear Ohmic relation I = GV , where G is given by the Landauer equation 2. be well approximated in the linear response regime by using τ (E) computed at equilibrium. When the response regime stops to be linear τ (E) has to be computed out of equilibrium. Since so far there is no rigorously exact formulation of DFT out of equilibrium, the calculation of τ (E) out of equilibrium involves an iterative procedure by means of the Poisson equation. This procedure will be detailed in section 2.2.4. We consider a conductor placed between two semi-infinite leads, as in figure 2.10. Such an open-system is closer to the usual experimental situation in which the leads are much bigger than the conductor, and can simulate the effect of the reservoirs providing electrons.

Transmission function

This section is about how to compute the transmission function τ (E) through the conductor under the assumption of coherent transport, that is without inelastic scattering events in the conductor. Firstly we will introduce the concept of scattering matrix (S-matrix), which is strictly related to the transmission function. Then, we will relate the S-matrix to the Green function, whose definition and properties will be briefly reminded. Finally we will show how, by DFT on a localized basis-set, we can compute an effective Hamiltonian, corrected by a term called self-energy to simulate the open-system, to obtain the Green function, and then the transmission function.

S-matrix

The S-matrix was introduced by John Archibald Wheeler in 1937 to relate at each energy the outgoing wave amplitudes to the incoming wave amplitudes of a physical system undergoing a scattering process. To characterize the S-matrix consider the simple one-dimensional case of a free electrons beam of energy E incident on a localized squared potential barrier U (x) (see fig. The solution of the Schrödinger equation outside the barrier are two plane waves, for the left and the right side of the barrier:

ψ L (x) = L in e ikx + L out e -ikx ψ R (x) = R out e ikx + R in e -ikx
where k = 2mE/h 2 . The S-matrix relates the outgoing amplitudes [L, R] out Electronic transport state of the art to the incoming amplitudes [L, R] in in a linear relation:

L out R out = s ll s lr s rl s rr L in R in (2.22)
that is, in a compact form, A out = S A in . Equation 2.22, which completely characterizes the scattering properties of the barrier U (x), constitutes a set of two equations in the four unknown variables s mn . Nevertheless, the property of density charge conservation, which implies density current conservation through the continuity equation, can provide the missing equations to infer the S-matrix. The density current J of the wave-function ψ(x) is defined as:

J = h 2mi ψ * ∂ψ ∂x -ψ ∂ψ * ∂x
Thus:

J L = hk m (L 2 in -L 2 out ) J R = hk m (R 2 out -R 2 in )
which should be equal as consequence of current conservation. Therefore:

L 2 in + R 2 in = L 2 out + R 2 out (2.23)
This means that the total incoming squared amplitude must be conserved in the total outgoing squared amplitude. Equation 2.23 can be rewritten as:

{A in } † {A in } = {A out } † {A out }
And since A out = S A in : Equation 2.24 expresses the fact that the S-matrix is unitary, and provide additional equations to derive all the coefficient of the S-matrix. The transmission and reflection coefficients for a beam from left to right are defined as:

{A in } † {A in } = {SA in } † {SA in } = {A in } † [S] † [S]{A in } Hence: [S] † [S] = s 2 ll +
T L = R 2 out L 2 in = s 2 rl R L = L 2 out L 2 in = s 2
where the second equality in both equations is obtained putting R in = 0 in equation 2.22. As a consequence of unitarity of the S-matrix, from equation 2.24:

s 2 rl + s 2 ll = T L + R L = 1 (2.25)
which means that the total probability is normalized to one, that is the electron must go somewhere, that is current must be conserved. These considerations about the S-matrix of a potential barrier can be ap- The transmission probability from one lead to another is the sum of all the transmission probability between conductive modes belonging to different leads:

T RL = m∈L n∈R T nm
The unitarity of the S-matrix implies:

M T m=1 s 2 mn = M T m=1 s 2 nm = 1
which is the equivalent of equation 2.25 for the potential barrier, and constitute a demonstration of equation 2.12, which is fundamental in simplifying the current expression and make it linear, as shown in the previous section. The S-matrix, and then the transmission function, can always be computed at given energy when having complete knowledge of the wave-function, as done for the simple case of a rectangular potential barrier. However, the elements of the S-matrix can be expressed in terms of the conductor Green function which, in practice, can be sometimes simpler to compute. If the S-matrix represents the response in a lead to a wave incident on the conductor in the same or the other lead, the Green function G(r, r ) represents the response in any point r of the system to an excitation in another point r , even into the conductor. The Green function is a more powerful object, which contains more information than the S-matrix. For this reason it is also suitable to extend the theory to inelastic scattering transport through the conductor. Indeed while the S-matrix is completely unaware of the dynamics into the conductor, which hence should be ballistic, the Green function can describe the inelastic interactions into the conductor. This is the basis of the NEGF theory (Non Equilibrium Green Function), which allows to account for scattering events into the conductor (e.g. electron-electron or electron-phonon scattering) and finite bias calculations (out of the linear response regime).

Green function

The Green function was introduced in the 1830s by the British mathematician George Green, who developed a method to solve non-homogeneous linear differential equations. Given an arbitrary linear operator L r , acting on a function space in the variable r, a differential equation is written as:

L r u( r) = f ( r) (2.26)
The Green function G( r, r ) associated to the L r operator is defined as:

L r G( r, r ) = δ( r -r )
Thanks to the Dirac delta properties:

f ( r) = f ( r ) δ( r -r ) d r = f ( r ) L r G( r, r ) d r
Moving L r , which only acts on r, outside the integral on the right hand side of equation, and substituting the left hand side by equation 2.26 gives:

L r u( r) = L r f ( r ) G( r, r ) d r Therefore: u( r) = f ( r ) G( r, r ) d r + q( r) (2.27)
where q( r) is a solution of the associated homogeneous differential equation L r q( r) = 0 completely defined by the boundary conditions. Equation 2.27 states that the convolution of G( r, r ) and f ( r) is the solution to the inhomogeneous differential equation for f ( r). In other words, since f ( r) is a sum of delta functions, the solution u( r) is also a sum of Green functions, by linearity of L r . From a physical point of view, the Green function in equation 2.27 plays the role of an impulse response function, which describes the reaction of any dynamic system at r in response to some external change at r . This makes it suitable to treat the problem of computing a scattering transmission function, which is the purpose here pursued. Considering the time independent Schrödinger equation:

[E -H( r)] ψ( r) = 0
the Green function associated to the H( r) operator can be found by solving:

[E -H( r)] G( r, r ; E) = δ( r -r ) (2.28)
Assume for simplicity that H( r) = h2 2m

∂ 2
∂ r 2 -U 0 , which represents a single free particle propagating in a constant potential U 0 , then both the following solutions satisfy equation 2.28:

G R ( r, r ; E) = -A exp i k • | r -r | G A ( r, r ; E) = A exp -i k • | r -r | where | k| ≡ 2m(E -U 0 ) h
Nevertheless, they correspond to different boundary conditions: G R , called retarded Green function, corresponds to outgoing waves that originate at the point of excitation r , while G A , called advanced Green function, corresponds to incoming waves that disappear at the point of excitation. A way to incorporate boundary conditions into the equation itself is to add an infinitesimal imaginary part to energy. Equation 2.28 becomes

[ -H( r)] G( r, r ; ) = δ( r -r ); = E ± iη; η > 0 (2.29)
If = E + iη the wave-vector k turns out to be:

| k | = 2m(E + iη -U 0 ) h = 2m(E -U 0 ) h 1 + iη E -U 0 ≈ 2m(E -U 0 ) h 1 + i η 2(E -U 0 ) = | k|(1 + i∆)

Electronic transport state of the art

Substituting this expression into the advanced Green function gives:

A exp -i k(1 + i∆) | r -r | = A exp -i k • | r -r | exp[∆ | r -r |]
Consequently the advanced Green function grows indefinitely when moving far away from the point of excitation r . This makes the retarded Green function the only acceptable solution, since a proper solution must be bounded.

-H] G( ) = I
This equation, which is equivalent to 2.29, implies that the only needed knowledge to compute the Green function of a given system is the Hamiltonian matrix H. Indeed the Green function is given by a standard operation of matrix inversion:

G( ) = [ I -H] -1 ; = E ± iη; η > 0 (2.33)
Once G( ) is known it can be used to recover the wavefunction ψ or to extract directly some property of the system, like for instance the transmission function.

As already pointed out, on one hand the Green function is defined by adding an infinitesimal imaginary part to energy to automatically take into account boundary conditions. On the other hand this is also necessary to well define the Green function by a limiting procedure in correspondence with the Hamiltonian eigenvalues, where the Green function diverges. Indeed, rewriting equation 2.33 using the second identity 2.32 gives:

G( ) = 1 -H n |φ n φ n |
1 -H can be moved inside the sum, since it is a constant. Then using equation 2.31 to substitute H with its eigenvalues gives the following relationship:

G( ) = 1 -H n |φ n φ n | = n 1 -H |φ n φ n | = n |φ n φ n | -E n ; = E n
This result shows that G( ) has poles in correspondence with the discrete eigenvalues E n of H. Nevertheless, G( ) can be defined at = E n by a limiting procedure, which is not possible without adding the imaginary part to energy:

G ± ( = E n ) ≡ lim η→0 + G( = E n ± iη) G + is
defined on Im( ) > 0 and corresponds to the retarded Green function, while G -is defined on Im( ) < 0 and corresponds to the advanced Green function.

Hamiltonian and self-energy by DFT

An effective and standard method to compute the Hamiltonian matrix H of a given system is the density functional theory (DFT) [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF][START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF], developed by Walter Kohn, Pierre Hohenberg and Lu Jeu Sham in the first half of the sixties. DFT is mainly based on two different concepts: Electronic transport state of the art 1. Electronic density is sufficient to completely characterize a quantum system, and can therefore be used as basic variable.

2. The problem of computing the ground state properties of a quantum system, composed by interacting electrons, with a given electronic density, can be replaced by the problem of computing the ground state properties of a quantum system composed by non-interacting electrons, with the same electronic density of the interacting system.

The celebrated Kohn-Sham (KS) equations write:

[ H H ( r; [n( r)]) + V xc ( r; [n( r)])]ψ j ( r) = E j ψ j ( r) n( r) = j |ψ j ( r)| 2 (2.34)
H H is the Hartree Hamiltonian operator. It consists in a semi-classical approximation which accounts for the electrostatic repulsion felt by one single electron as a mean field produced by the others:

H H ( r; [n( r)]) = - h2 ∇ 2 2m + U ext ( r) + e 2 n( r ) | r -r | d r (2.35)
V xc is called exchange correlation functional, and contains all the many-body corrections beyond the Hartree mean filed theory. V xc is in practice never known for real systems, thus its approximation deeply influences the DFT accuracy. The most popular approximation is certainly the local density approximation (LDA), which defines V xc as a local contribution taken from the homogeneous electron gas theory:

E LDA xc = n( r) hom xc (n( r)) d r
where hom xc (n( r)) is the exchange and correlation energy per electron for a homogeneous electron gas of density n, and is known. A refinement of the LDA is the generalized gradient approximation (GGA) [START_REF] Perdew | Generalized gradient approximation made simple[END_REF], where the exchange correlation functional also depends on the density gradient:

E GGA xc = f [n( r), ∇n( r)] d r
Several forms of f [n( r), ∇n( r)], which is a general function, exist and are used in literature. It is important to understand that, even if V xc were exact, the KS eigenvalues and eigenfunctions would not be those of the true, interacting system, but those of the auxiliary non-interacting system. Hence they should be considered just as mathematical objects without any physical meaning. We now show briefly how the KS equations are solved in practice, because this is useful to understand how H, and G in turn, are computed. Since both H 0 and V xc are a functional of the electronic density n( r), which in turn is defined by the orbitals ψ j ( r), the KS equations must be solved selfconsistently starting from an initial guess. To solve the KS equations at each self-consistent cycle consider a basis-set of known functions f α ( r) such that the KS orbitals can be expanded as follows:

ψ j ( r) = α c j,α f α ( r) (2.36)
Substituting this expression into the first KS equation 2.34, multiplying both sides by a function f * β ( r) belonging to the basis-set, and integrating over space gives:

α c j,α f * β ( r) H( r)f α ( r)d r = α c j,α E j f * β ( r)f α ( r)d r α c j,α [H α,β -E j S α,β ] = 0 (2.37)
where H = H H + V xc , H α,β and S α,β are respectively the elements of the H matrix and the overlap matrix S, expressed in the basis-set f α ( r), and defined by the integrals in the first line of equation 2.37. The homogeneous linear system 2.37 admits non zero solutions for the unknown coefficients c j,α only if:

Det[H -E j S] = 0 (2.38)
Therefore the problem of solving the KS equations is transformed into a problem of linear algebra, which can be solved by standard diagonalization techniques. The solution of equation 2.38 provides the eigenvalues E j , which then are substituted into equation 2.37 to obtain the coefficients c j,α . The orbitals ψ j ( r) are finally available through equation 2.36. In the end the orbitals ψ j ( r) constitute a complete basis-set on which the matrix H is diagonal, filled with the eigenvalues E j . If the considered system is a bulk material with some symmetry and periodicity the Bloch theorem states that the orbitals ψ j ( r) should satisfy:

ψ n k ( r + r l ) = e i k• r l ψ n k ( r) (2.39)
where r l is a lattice translation defined by r l = 3 i l i a i , a i are the primitive lattice vectors and l i are arbitrary integers. n and k are band and crystal vector indexes, respectively. That given the orbitals can be decomposed in a product of a phase factor and a function u n k having the crystal periodicity, that is a Bloch state:

ψ n k ( r) = 1 √ Ω e i k• r u n k ( r); u n k ( r) = u n k ( r + r l )
where Ω = a 1 • ( a 2 × a 3 ) is the volume of the primitive cell.

It is then a natural choice to expand the periodic function u n k ( r) in a discrete sum over the reciprocal lattice vectors:

u n k ( r) = G e i G• r u n( k+ G) (2.40)
where G = 3 i l i b i are lattice translations in the reciprocal space, b i = 2π/ a i are the reciprocal vectors, l i are integers, and u n( k+ G) are Fourier coefficients in the reciprocal space. Therefore DFT codes to simulate bulk materials are often implemented by using a plane-waves basis-set of the form:

ψ n k ( r) = 1 √ Ω G e i( k+ G)• r u n( k+ G)
The problem of solving the KS equations is then reduced to the linear algebra problem of finding the coefficients u n( k+ G) . In principle, to represent exactly any function u n k ( r), the sum 2.40 should run on an infinite number of plane-waves, that is the basis-set should be complete. However, this is not feasible for numerical applications on a computer. From the basic principles of Fourier analysis, the need of including wave vectors of increasing magnitude is connected to the smoothness of the function to be expanded: the more smooth the function, the faster the expansion converges. Functions that vary rapidly in space need a very large number of wave vectors when expanded in plane-waves. Therefore, the sum in 2.40 is generally rearranged summing up over wave vectors of increasing magnitude. Then, it is cut off by excluding wave vectors such that

h2 | k+ G| 2 2m > E cut : u n k (r) = G: h2 | k+ G| 2 2m ≤Ecut e i G• r u n( k+ G) (2.41)
Once the KS equations are solved, and the electronic density is converged, the Hamiltonian matrix H per each k-point k can be equally expressed in the plane-waves basis-set, or in the orbitals ψ n k basis-set. In the first case the H matrix elements are given by:

k + G H k + G = e -i( k+ G )• r H( r)e i( k+ G)• r d r = e i( G-G )• r H( r) d r
where H( r) in the last equality contains the kinetic term

h2 ( k+ G) 2 2m
, which is the result of applying the kinetic operator contained in H( r) (first term of eq. 2.35) on the plane wave e i( k+ G)• r . In the second case the H elements are instead:

ψ n k H ψ n k = E n k ψ n k ψ n k = E n k δ n n
where here δ n n is the Kronecker delta. H is thus diagonal. Both of them can be inverted by applying equation 2.33 to find the Green function matrix G( ) in the correspondent basis-set. Then, as we will see, the transmission function τ ( ) is basically given by the trace of G( ), which is an invariant property under basis-set transformation. Hence to express H in different basis-sets does not change the final transmission function.

This procedure to obtain τ , that is to solve the KS equations to obtain H, to invert H to obtain G, and to compute the trace of G to finally compute τ , has one main disadvantage. As a consequence of equation 2.39, the KS equations can be solved for the primitive cell applying periodic boundary conditions (PBC). For transport calculations, the primitive cell becomes a supercell that contains the conductor and the leads. This permits to well describe the conductor-leads interface, which is critical to correctly account for the transmission coefficient. The KS equations are thus solved for this supercell replicated in three dimensions, and H represents a system with a succession of replicated leads and conductors in the transport direction. On the contrary, we are interested in computing the transmission function of a system where the periodicity along the transport direction is broken, as in figure 2.10. For this purpose, the Bloch states should be modified to describe the semi-infinite leads before computing H. Nevertheless, this is not easy in practice. To overcome this problem it is worth considering the use of a space-localized functions instead of a plane-waves basis-set. In this case, the H elements are given by:

f β |H|f α = f * β ( r) H( r)f α ( r)d r (2.42)
where f α ( r) and f β ( r) are localized functions of the basis-set. The advantage of writing the Hamiltonian matrix as in equation 2.42 is that the matrix elements can assume a real spatial interpretation, representing the coupling between two localized functions in the basis-set. This allows extracting Hamiltonian sub-matrices of localized regions belonging to the supercell: for instance the region of the conductor, or the region of a lead. These different sub-matrices can then be combined to build an effective Hamiltonian which represents the open-system with semi-infinite leads. Localized functions basis-sets have also other advantages: firstly they make it possible to perform an efficient k-point interpolation [START_REF] Slater | Simplified LCAO method for the periodic potential problem[END_REF] to refine the DFT grid for the following transport calculation, which often requires more k-points than DFT to converge; and secondly, they are suitable to extend the theory to full NEGF calculations. On the contrary the main disadvantage of localized functions basis-sets is that they are not as easy as plane-waves basis-sets to optimize, and convergence is often non-trivial to achieve. Indeed for localized functions, there is not a simple rule as 2.41, which easily permits to perform convergence tests.

As described in several works [START_REF] Datta | Electronic Transport in Mesoscopic Systems[END_REF][START_REF] Nardelli | Electronic transport in extended systems: Application to carbon nanotubes[END_REF][START_REF] Papior | Improvements on non-equilibrium and transport green function techniques: The next-generation transiesta[END_REF][START_REF] Calzolari | Ab initio transport properties of nanostructures from maximally localized wannier functions[END_REF][START_REF] Ventra | Electronic Transport in Nanoscale Systems[END_REF], for the open-system of figure 2.10, the Green function can be partitioned into sub-matrices that correspond to the conductor and leads individual subsystems:

  g L g LC g LCR g CL G C g CR g LRC g RC g R   =   -h L -h LC 0 -h † LC -H C -h CR 0 -h † CR -h R   -1 (2.43)
where the matrix [ -H C ] represents the finite "isolated" conductor with no coupling elements to the leads, [ -h R,L ] represent the semi-infinite leads, and h CR and h LC are the coupling matrices between the conductor and the leads. Direct couplings between the two leads are assumed to be negligible. As a convention, lower case letters are used for semi-infinite matrices and upper case for finite dimension matrices. Writing explicit equations involving the second column of the left hand side matrix in equation 2.43 gives:

     g LC ( -h L ) -G C h LC = 0 -g LC h † LC + G C ( -H C ) -g RC h CR = I -G C h † CR + g RC ( -h R ) = 0
Then solving the first and third equations for g LC and g RC , and substituting the resulting expressions in the second equation to solve for G C leads to:

G C ( ) = [ -H C -Σ L ( ) -Σ R ( )] -1 (2.44)
where:

Σ L ( ) = h † LC [ -h L ] -1 h LC Σ R ( ) = h † CR [ -h R ] -1 h CR Σ R,L
are defined as the self energy functions due to the semi-infinite leads. These terms can be viewed as effective Hamiltonian matrices which arise from the coupling of the conductor with the leads. Indeed, the Schrödinger equation corresponding to 2.44 turns out to be:

(H C -Σ L -Σ R )ψ C = Eψ C
where the self-energy functions clearly act as an effective potential on the central region due to the leads. If G C ( ) were known by equation 2.44, the correspondent transmission function through the open-system would be given by the trace formula, derived by Meir and Wingreen [START_REF] Meir | Landauer formula for the current through an interacting electron region[END_REF]:

τ ( ) = Tr[Γ L ( ) G r C ( ) Γ R ( ) G a C ( )] (2.45)
where Γ L,R are functions which describe the coupling of the conductor to the leads. These functions are strictly related to the self energy functions Σ R,L :

Γ R,L ( ) = i Σ R,L ( ) -Σ † R,L ( ) (2.46)
As already said equation 2.45 is invariant under transformations of the basisset, which can be chosen to express the Hamiltonian matrices in the most convenient form. However, it seems that the formalism until here described by equation 2.44 and 2.45 would bring no gain because of two reasons. Firstly DFT cannot simulate an infinite system, but only a supercell with finite leads replicated along the transport direction. By using a localized basis-set the sub-matrix H C can be extracted from the supercell, but the infinite matrices h L and h R remain unknown. Secondly, even if h L and h R were known, it would be unfeasible to invert them since they are infinite. To overcome these difficulties consider that any solid can be viewed as an infinite stack of principal layers with nearest-neighbor interactions. This corresponds to transform the original system into a linear chain of principal layers. The Hamiltonian matrix of a principal layer can then be extracted from the supercell by using the localized basis-set properties. For simplicity consider a system where conductor and leads are made of the same material. Thus the conductor can be viewed as one principal layer sandwiched between two semi-infinite stacks of principal layers. That given the Hamiltonian matrix can be divided into submatrices for the different layers and for the coupling between first neighbors. Moreover in a bulk system is assumed that all the Hamiltonian sub-matrices corresponding to the different layers are equal, as well as the first neighbors coupling sub-matrices. Therefore in the principal layer approach equation 2.33 can be written as:

  -H 00 -H 01 0 -H † 01 -H 00 -H 01 0 -H † 01 • • •   •   G 00 G 01 G 02 G 10 G 11 G 12 G 20 G 21 • • •   =   I 0 0 0 I 0 0 0 • • •   (2.47)
where all the sub-matrices are of finite dimensions. It is then straightforward to express the Green function of an individual layer, that is the diagonal Electronic transport state of the art elements, in terms of the coupling Green functions with the preceding and following one. For instance, from 2.47, G 11 obeys the following expression:

-H † 01 G 01 + ( -H 00 )G 11 -H 01 G 21 = I
By introducing the transfer matrices T and T , defined such that:

G n+1,n = T G nn , G nn = T G n+1,n
G 11 becomes:

G 11 = -H 00 -H † 01 T -H 01 T -1
(2.48) Equation 2.48 is equivalent to equation 2.44 by the following identification:

Σ L = H † 01 T , Σ R = H 01 T
Thus the knowledge of the transfer matrices T and T is sufficient to calculate the self energy functions without inverting infinite matrices, which would be computationally not feasible. Lopez-Sancho et al. describe a method to compute the transfer matrices T and T by an iterative procedure [START_REF] Lopez-Sancho | Quick iterative scheme for the calculation of transfer matrices: application to Mo (100)[END_REF][START_REF] Lopez-Sancho | Highly convergent schemes for the calculation of bulk and surface green functions[END_REF]. From equation 2.47 it is straightforward to obtain a series of matrix equations for the Green function:

( -H 00 )G 00 = I + H 01 G 10 ( -H 00 )G 10 = H † 01 G 00 + H 01 G 20 • • • ( -H 00 )G n0 = H † 01 G n-1,0 + H 01 G n+1,0 (2.49)
The last expression in 2.49 is general, valid for the G n0 matrix, with arbitrary values of n, and for the G n+1,0 and G n-1,0 matrices as well. Therefore, introducing:

t 0 = [ -H 00 ] -1 H † 01 t0 = [ -H 00 ] -1 H 01 (2.50)
it is correct to write:

G n0 = t 0 G n-1,0 + t0 G n+1,0 , n ≥ 1 G n-1,0 = t 0 G n-2,0 + t0 G n0 , n ≥ 2 G n+1,0 = t 0 G n0 + t0 G n+2,0 , n ≥ 0 (2.51)
Then substituting the second two lines into the first leads to:

G n0 = t 0 (t 0 G n-2,0 + t0 G n0 ) + t0 (t 0 G n0 + t0 G n+2,0 )
and solving for G n0 writes:

G n0 = t 1 G n-2,0 + t1 G n+2,0 n ≥ 2 (2.52)
where:

t 1 = I -t 0 t0 -t0 t 0 -1 t 2 0 t1 = I -t 0 t0 -t0 t 0 -1 t2 0
Since the first line of equation 2.51 is isomorphic to equation 2.52, the process can be repeated iteratively. After i iterations:

G n0 = t i G n-2 i ,0 + ti G n+2 i ,0 n ≥ 2 i (2.53)
where:

t i = I -t i-1 ti-1 -ti-1 t i-1 -1 t 2 i-1 ti = I -t i-1 ti-1 -ti-1 t i-1 -1 t2 i-1
(2.54) Equation 2.53 defines the following chain:

G 10 = t 0 G 00 + t0 G 20 G 20 = t 1 G 00 + t1 G 40 • • • G 2 n 0 = t n G 00 + tn G 2 n+1 0
Then substituting the second line into the first:

G 10 = (t 0 + t0 t 1 )G 00 + t0 t1 G 40
And iterating the substitution gives:

G 10 = (t 0 + t0 t 1 + • • • + t0 • • • tn-1 t n )G 00 + t0 • • • tn G 2 n+1 0 (2.55)
The iterations have to be repeated until t n+1 , tn+1 < ξ as small as wished.

Then the second term in equation 2.55 vanishes reducing the expression to:

G 10 = T G 00 where T = t 0 + t0 t 1 + t0 t1 t 2 + • • • + t0 • • • tn-1 t n (2.56)
which is the desired expression for the transfer matrix T. t 0 is given by equation 2.50, and t i by equation 2.54. Similarly T is given by:

T = t0 + t 0 t1 + t 0 t 1 t2 + • • • + t 0 • • • t n-1 tn (2.57)
The procedure outlined above can also be applied to the case of electron transmission through one or more interfaces, between different media. This case requires a more advanced technique, known as the surface Green function matching (SGFM) theory, pioneered by [START_REF] Garcia-Moliner | Theory of Single and Multiple Interfaces[END_REF]. The scattering region, defined by an Hamiltonian H C , is coupled to the leads by the matrices H LC and H CR . The leads are represented in the principal layer approximation by a set of sub-matrices H {L,R}

nm . That given it can be shown that equation 2.33 is transformed into the following secular equation [START_REF] Garcia-Moliner | Theory of Single and Multiple Interfaces[END_REF]:

  G L G LC G LR G CL G C G CR G RL G RC G R   =   -H L 00 -(H L 01 ) † T -H LC 0 -H CL -H C -H CR 0 -H RC -H R 00 -H R 01 T   -1
Similarly to what done to derive equation 2.44 from equation 2.43, it is straightforward to write an expression for the central region Green function:

G C ( ) = [ -H C -Σ L ( ) -Σ R ( )] -1 (2.58)
where:

Σ L ( ) = H † LC -H L 00 -(H L 01 ) † T L -1 H LC Σ R ( ) = H CR -H R 00 -H R 01 T R -1 H † CR (2.59)
A one dimensional chain of principal layer, which is physical only for systems like nano-tubes or quamtum wires, is here assumed. The extension to a three dimensional case is straightforward using the Bloch theorem.

Non-linear response to bias: Poisson equation

So far DFT is in principle rigorously formulated only for equilibrium and closed systems (i.e. fixed number of electrons). However, some codes which implement the Landauer formalism on top of DFT, such as TranSiesta [START_REF] Papior | Improvements on non-equilibrium and transport green function techniques: The next-generation transiesta[END_REF], allow applying a finite bias voltage between the electrodes, still using the equilibrium density functionals. Despite the incompleteness of such approach [START_REF] Vignale | Incompleteness of the Landauer formula for electronic transport[END_REF], it allows addressing realistic devices within a reasonable computing time and provides valuable information about the transmission of electrons through interfaces.

The procedure, which is iterative, approximately works as follow: in a system at equilibrium the electronic density ρ eq can be computed as function of the Green function 2.58:

ρ eq = i 2π BZ [G k ( ) -G † k ( )] η f (E) e -i k• r d k dE (2.60)
where η f is the Fermi distribution of both electrodes, which are at equilibrium. Non-equilibrium arises when the two electrodes have different Fermi distributions η f,L and η f,R due to a temperature gradient or to an applied bias. In this case the total density ρ can be decomposed in the equilibrium term 2.60 plus a non-equilibrium correction:

ρ = ρ eq + ∆ L R ≡ ρ neq ∆ L
R represents the non-equilibrium density correction of the left electrode L due to the right electrode R, and it is given by:

∆ L R = 1 2π BZ A R,k ( ) [η f,R (E) -η f,L (E)] e -i k• r d k dE
where A R,k ( ) is the right electrode spectral function given by:

A R,k ( ) = G k ( ) Γ R,k ( ) G † k ( )
The non-equilibrium density can then be injected into the Poisson equation to obtain a corrected value for the Hartree potential φ H ( r):

∇ 2 φ H ( r) + ρ( r) 0 = 0
The Hartree potential is finally directly linked to the electric potential V ( r) = -e φ H ( r), which can be used to directly correct the Hamiltonian matrix elements. The on-site energies, that is the Hamiltonian diagonal matrix elements, become:

H ii → H ii + φ * i ( r)V ( r)φ i ( r) d r + δV xc (ρ)
where φ i ( r) constitute the complete basis-set of localized functions to expand the wavefunctions. The new non-equilibrium density also changes the exchange correlation functional, which is directly dependent on density. A correction δV xc (ρ) can thus be included in the Hamiltonian updating. If V ( r) does not vary rapidly close to the on-site position r i , given by the localized function φ i , it exits the integral which reduces to one. In this case V ( r i ) directly corrects the Hamiltonian:

H ii → H ii + V ( r i ) + δV xc (ρ)
A similar discussion can be done for the non-diagonal matrix elements. An important difference is that, if V ( r) has no rapid oscillations, and if the basisset is orthogonal, the contribution due to the electric potential is zero. The reason is that φ * i ( r)φ j ( r) d r = 0. A new Green function is then computed using the updated Hamiltonian and a new iteration begins. The procedure is iterated until the desired convergence is reached.

Chapter 3 The many body Landauer transport

The "DFT + Landauer" formalism, described in chapter 2, is implemented in several widely used transport codes (Transiesta [START_REF] Papior | Improvements on non-equilibrium and transport green function techniques: The next-generation transiesta[END_REF], based on the Siesta code [START_REF] Ordejón | Self-consistent order-N density-functional calculations for very large systems[END_REF][START_REF] Soler | The SIESTA method for ab initio order-N materials simulation[END_REF]; Wannier90 [START_REF] Calzolari | Ab initio transport properties of nanostructures from maximally localized wannier functions[END_REF][START_REF] Mostofi | Wannier90: A tool for obtaining maximally-localised Wannier functions[END_REF][START_REF] Mostofi | An updated version of Wannier90: A tool for obtaining maximally-localised wannier functions[END_REF]), which are able to compute the transmission function, and then the conductance, of a bulk system of nanoscopic dimensions. On one hand this approach has some limitations and constraints to be applied, and on the other hand, the obtained results are qualitatively good but quantitatively imprecise. In this chapter we present a development of the "DFT + Landauer" formalism . This improved methodology, elaborated within the frame of this doctoral work, can significantly improve the reliability of the results, while the limitations remain unvaried. To introduce the new approach, which takes the name of "Many Body Landauer Transport", we briefly characterize the limitations and the reliability of the "DFT + Landauer" formalism:

Limitations Some limitations are due to the Landauer theory which, as pointed out in section 2.1.2, requires the electron mean free path to be longer than the length of the active channel of conduction. As a consequence the Landauer theory cannot afford to account for events like electron-phonon or electron-electron scattering, which in general can always occur even in nanoscopic systems. Another limitation is due to DFT. So far there exists no exact formulation of DFT for out of equilibrium or open systems (i.e. varying number of electrons). Therefore DFT is in principle only suitable to calculate the current at voltages in the linear response regime, where it could be well approximated by its equilibrium value. Even by the Poisson iterative procedure of section 2.2.4 the current calculation is in theory not rigorous out of the linear response The many body Landauer transport regime, which is very narrow at room temperature.

Reliability In principle, Kohn-Sham DFT is an exact many body (MB) theory exclusively to calculate the ground state total energy and electronic density. Kohn-Sham DFT has no access to excited state energies and wavefunctions. Indeed, the Kohn-Sham energies and wavefunctions are just only mathematical objects referring to the fictitious Kohn-Sham auxiliary system, and they have nothing to do with the physical quasiparticle energies and wavefunctions of the real system. This is the case also for the exact Kohn-Sham energies and wavefunctions solution of the non-approximated exact exchange-correlation potential V xc : this has been explicitly checked in systems where the exact V xc is mathematically known or inferred by reverse engineering [START_REF] Li | Comparing manybody approaches against the helium atom exact solution[END_REF]. Therefore, the use of the Kohn-Sham electronic structure to calculate the Green function, instead of a more accurate quasi-particle (QP) electronic structure, has to be considered as just only an approximation which strongly depends on the chosen xc functional, and which turns out to be in many cases unreliable. An example of unreliability is the systematical underestimation of the band-gap [START_REF] Perdew | Density functional theory and the band gap problem[END_REF] in insulators and semiconductors, or of the HOMO-LUMO gap in molecules, by the most popular and physical functionals in the LDA or GGA/PBE approximations. Another example, closer to the topic of this thesis, is the systematic overestimation of the 0bias conductance by even several order of magnitudes in tunneling junctions [START_REF] Gangineni | Interfacial electronic transport phenomena in single crystalline Fe-MgO-Fe thin barrier junctions[END_REF][START_REF] Vignale | Incompleteness of the Landauer formula for electronic transport[END_REF] by DFT-LDA/GGA Landauer calculations. This work is focused on the reliability of the transmission function calculation, leaving the problem to overcome the limitations as perspective. Treating scattering events and out of equilibrium systems would require the development of a more general theory, which takes the name of Non Equilibrium Green Function (NEGF). The development of this theory, which is not described in details here, constitute the natural continuation of this work in perspective. To improve the reliability the "Many Body Landauer" technique goes beyond DFT towards more appropriate MB, QP approaches. Between these approaches, we particularly consider two many body perturbation theories (MBPT): GW [START_REF] Hedin | New method for calculating the one-particle Green's function with application to the electron-gas problem[END_REF] and COHSEX (Coulomb hole screened exchange) [START_REF] Hedin | New method for calculating the one-particle Green's function with application to the electron-gas problem[END_REF]. We also consider the HSE [START_REF] Becke | A new mixing of Hartree-Fock and local density functional theories[END_REF] hybrid functional, which cannot rigorously be considered an approximation within DFT like LDA or GGA functionals, but rather a semi-empirical MB approximation, with a non-local self-energy half-way between the Hartree-Fock (HF) [START_REF] Slater | A simplification of the Hartree-Fock method[END_REF] exchange Fock operator and the COHSEX screened-exchange static self-energy. Compared to DFT, these QP approaches certainly provide a more physical electronic structure, and this can translate in an improved description also of quantum transport properties, that is transmission function and conductance. Conductance calculations by Green function formalism and Landauer formula on top of QP calculations are applied in literature only to molecular systems [START_REF] Rangel | Transport properties of molecular junctions from manybody perturbation theory[END_REF][START_REF] Rangel | Many-body correlations and coupling in benzene-dithiol junctions[END_REF]. Localized functions basis-sets, which are required to compute the transmission function, are particularly suitable for molecular systems, which hence are somehow easier to treat. There are no existing DFT codes which can perform QP calculations on bulk material systems by using a localized functions basis-set since this task is much easier by using a plane-waves basis-set. The challenge of this work is thus to combine the features of well known plane-waves codes (VASP [START_REF] Kresse | Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set[END_REF][START_REF] Kresse | From ultrasoft pseudopotentials to the projector augmented-wave method[END_REF], Quantum Espresso [START_REF] Giannozzi | QUANTUM ESPRESSO: a modular and opensource software project for quantum simulations of materials[END_REF][START_REF] Giannozzi | Advanced capabilities for materials modelling with quantum ESPRESSO[END_REF], Abinit [START_REF] Gonze | First-principles computation of material properties: The ABINIT software project[END_REF][START_REF] Gonze | A brief introduction to the ABINIT software package[END_REF][START_REF] Gonze | ABINIT: First-principles approach to material and nanosystem properties[END_REF][START_REF] Gonze | Recent developments in the ABINIT software package[END_REF]), which are highly developed for QP calculations on bulk materials, but are not suitable for transport calculations, and the features of localized functions basis-sets, which allow to compute the Green function 2.58, but are not suitable for bulk materials. Even if, at first sight, the features of two different basis-sets may appear incompatible, a methodology based on a basis-set transformation can fulfill this task. This idea, even if not free from difficulties and challenges, constitutes the basis to build the "Many Body Landauer" method. This method is a puzzle of well known techniques, which combined constitute an original approach, at our knowledge so far never implemented in any distributed code. The procedure is divided into four steps, each of which detailed in the next sections:

1. Building of a supercell which represents the bulk system that one desires to study, and which is suitable for the next step of the procedure. The Hamiltonian sub-matrix of the bulk electrode principal layer must correspond as much as possible to the isolated bulk structure. The best strategy is thus to extract it from the cubic cell closest to the border of the supercell, that is far from the insulator. The electrode should thus be long enough to avoid the insulator to interact with the metal at the border of the supercell. Moreover, since the supercell will be the object of a simulation which makes use of periodic boundary conditions (PBC), we need to isolate spurious replicas of the insulator. Thus the length of the electrodes in the supercell should also be sufficient to avoid direct interaction between the insulator and its PBC replicas along the transport direction. At the interface between electrodes and insulator, a relaxation always occurs of both the atomic and electronic structures (e.g. transfers of charge, which can always occur even at equilibrium and 0 bias). Therefore the Hamiltonian matrix of the central part, which represents the conductor in the Landauer formalism, should not contain the bare insulator alone, but also the interface with some layers of metal. As we will see in section 3.4.2, we choose to consider the whole ab initio supercell as the central conductor. Figure 3.1 shows an example of supercell built for electronic transport calculations. Building the supercell can be a tricky task since it is the result of a compromise between the features to satisfy and the size of the supercell. The latter should indeed be reduced as much as possible to save computational resources, which can badly scale when performing heavy QP calculations in step 2.

Step 2: Electronic structure correction

The purpose of this step is to correct the DFT electronic structure, that is eigenvalues and optionally also eigenfunctions, by means of a QP calculation. The QP electronic structure, calculated on a plane-waves basis and then transformed on a Wannier functions basis, turns out to be a more coherent object to build the Green function and the transmission function, which assume a more relevant physical meaning. An important point to underline is that if this step is skipped, that is Bloch states and eigenvalues are computed by DFT without applying any correction, the expected conductance at the end of the procedure should be very close to existing transport codes based on localized functions basis-sets DFT, like for instance Transiesta or Wannier90. This is particularly useful to test the reliability of the method, and especially the quality of the basisset transformation (step 3), which should not modify the physics. In other words, a plane-waves DFT calculation followed by a basis-set transformation to Wannier functions is supposed to give results very similar to a direct localized functions DFT calculation. We now briefly introduce the second quantization, where the MB Green function can be more easily expressed. Then we introduce the QP equation and the Hedin equations, which are the basis for the GW and COHSEX approximations. Finally, we shortly describe the HSE hybrid functional. To perform this step many different plane-waves DFT codes, like Quantum Espresso, Abinit, or VASP, may be used.

The many body Landauer transport

Second quantization and many body Green function

Consider the time dependent Schrodinger equation for a single particle, or many non-interacting particles:

ih ∂ ∂t ψ( r, t) = H( r)ψ( r, t) (3.1)
where the Hamiltonian is in the general form:

H( r) = - h2 2m ∇ r -U ( r)
and U ( r) is an external local potential. The Green function associated to equation 3.1 is:

ih ∂ ∂t -H( r) G(r, t; r , t ) = δ(t, t ) δ(r, r ) (3.2)
which, Fourier transformed, turns out to be: The impossibility to treat wavefunctions in the form ψ( r 1 , . . . , r N ; t) was one of the motivations which pushed Paul Dirac to develop the second quantization formalism [START_REF] Mahan | Many-Particle Physics[END_REF][START_REF] Stefanucci | Nonequilibrium Many-Body Theory of Quantum Systems, A Modern Introduction[END_REF] in the end of the twenties. In standard quantum mechanics, also called first quantization, the quantum many-body state (QMBS) is represented on a basis-set of single particle wave functions ψ 1 , . . . , ψ N in the Hilbert space. The QMBS is described by answering a series of questions like "which particle is on which state". However, these are not physical questions because the particles are identical, and it is impossible to tell which particle is which in the first place. The seemingly different states ψ 1 , . . . , ψ N are actually redundant names of the same quantum many-body state. Thus the symmetrization, or anti-symmetrization, must be introduced to eliminate this redundancy in the first quantization description.

[hω -H( r)]G(r, r , ω) = δ(r, r ) (3.3)
In the second quantization formalism, the QMBS is represented on the occupation number basis. A basis state |n α denotes n particles in the singleparticle state ψ α . Basis states are also known as Fock states, which form a complete basis-set for the many-body Hilbert space, named the Fock space. In the second quantization language, instead of asking "each particle on which state", the question becomes "how many particles are there on each state". Because this description does not refer to the labeling of particles, it contains no redundant information, and hence leads to a precise and simpler description of the QMBS.

Even if first and second quantization are formally equivalent, treating manybody interacting systems by using the second quantization presents great advantages:

• Indexes on the total, often very big, number of particles are avoided.

• The bosonic symmetrization and the fermionic anti-symmetrization are automatically taken into account.

• Systems with a variable number of particles are treatable.

• The Green function can be generalized to a correlation function of field operators. In the new definition the Green function is a function of only two space-time variables, and can be used as fundamental degree of freedom of the theory instead of the more complicated many body wavefunction ψ( r 1 , . . . , r N , t), function of 3N spatial coordinates and one time variable. Therefore a much simpler object contains the whole physics of the system.

The definition of Green function generalized to a correlation function is [START_REF] Mahan | Many-Particle Physics[END_REF]:

G( r, t; r , t ) = -i Ψ 0 |T C( r, t) C † ( r , t )|Ψ 0 (3.4)
Equation 3.4 holds in particular for an Hamiltonian of the type:

H = H 0 + V
where H 0 is the non-interacting Hamiltonian, whose eigenfunction and eigenvalues are known, V contains the interaction, and H is the interacting Hamiltonian. In equation 3.4 Ψ 0 represents the ground state of H, which is unknown, and is actually what we want to compute by means of the Green function. T is the time ordering operator, which has the effect to produce two different solutions for t > t and t < t . The first case represents the propagation of a particle from the point (r , t ) to the point (r, t), while the second represents the propagation of a hole. The two solutions correspond respectively to the retarded and the advanced Green function introduced in section 2.2.2. Finally the C and C † are operators in the Heisenberg representation, that is dependent on time, defined by:

C( r, t) = e i Ht k a k ψ k ( r) e -i Ht C † ( r, t) = e i Ht k a † k ψ * k ( r) e -i Ht (3.5)
where ψ k are the known eigenfunctions of H 0 , while a † k and a k are respectively the creation and annihilation operators, which creates and destroys a particle in a particular state k. These operators constitute the core of the second quantization. They set the main difference between classical fields and quantum fields, that is the appearance of particles. In other words, a quantum field is no longer thought as a continuous field but as an entity that is composed of discrete particles. Remark that the part in square brackets on the right hand side of equation 3.5 is the definition of C and C † in the Schrödinger picture. Thus equation 3.5 also constitutes the relation between Heisenberg and Schrödinger pictures. It can be shown that the non-interacting "first quantization" Green function 3.3 is equivalent to the generalized "second quantization" Green function 3.4 when the latter is meant as a two-point correlation function of a noninteracting system. Omitting the formal mathematical details, to show the equivalence C and C † in equation 3.4 have firstly to be rewritten in the Schrödinger picture (eq. 3.5). Then separating C and C † in two different integrals and Fourier transforming to frequencies equation 3.4 becomes:

G( r, r , ω) = s L s ( r) L * s ( r ) hω -s ; s = E s ± iη (3.6)
where, as done in section 2.2.2, an arbitrary small imaginary part ±iη is added to the eigenvalues of H 0 to take into account boundary conditions and well define the Green function. L s ( r) and L * s ( r ) are called Lehman amplitudes, defined by:

L s ( r) = Ψ N 0 C( r) Ψ N +1 s ; if E N +1 s -E N 0 ≥ µ Ψ N -1 s C( r) Ψ N 0 ; if E N 0 -E N -1 s < µ
where µ is the chemical potential of the system. The Lehman amplitudes connect the Ψ N 0 state to the excited states Ψ N +1 s and Ψ N -1 s states by creation and annihilation of one electron, respectively. That given it is clear that the poles of the Green function 3.6 are the addition and removal energies of the system. In the case of a non-interacting system the Lehman amplitudes reduce to the eigenfunctions ψ s of the non-interacting Hamiltonian H 0 , and equation 3.6 reduces to:

G( r, r , ω) = s ψ s ( r) ψ * s ( r ) hω -s ; s = E s ± iη (3.7)
which is formally equivalent to equation 3.3. The difference between the interacting and the non-interacting case, that is between equation 3.6 and 3.7, can be better appreciated when looking at the spectral function, which is nothing but the imaginary part of the Green function:

A( r, r , ω) = 1 π Im [G( r, r , ω)] = s L s ( r) L * s ( r )δ(hω -s ); interacting s ψ s ( r) ψ *
s ( r )δ(hωs ); non-interacting (3.8) Figure 3.2, taken from ref. [START_REF]Exchange and Correlation in the Electronic Structure of Solids, from Silicon to Cuprous Oxide: GW Approximation and beyond[END_REF], shows A ii (ω) = ψ i |A|ψ i for the interacting and non-interacting cases, that is for the two different lines of equation 3.8. In the non-interacting case A ii (ω) is a delta peak, which indicates that this state is stable. In the interacting case, A ii (ω) presents a finite width, which indicates that the state presents a life-time. This means that the excited state oscillates for some time before relaxing back to the ground state. These states are known as quasiparticles. A second peak, called satellite plasmon peak, may appear in the interacting case. Therefore, not only the excitations of the non-interacting system, when going to the interacting system, get renormalized as by a shift of their energy (Re[Σ]) and acquire a broadening/lifetime (Im[Σ]); but also the excitation spectrum of the interacting system is enriched of a new sector (non-coehrent part of the excitation spectrum), constituted by excitations with a character beyond the simple single-particle electronhole type and including collective character such as plasmons (quasiparti-cle+plasmon). In principle, to compute the interacting Green function 3.4 or the equivalent 3.6, a complete knowledge of the |Ψ interacting eigenstates is necessary, Figure 3.2: Schematic representation of the spectral function A ii (ω). In the non-interacting case A ii (ω) is a delta peak. In the interacting case, the matrix elements have a width and a satellite peak may appear. Taken from ref. [START_REF]Exchange and Correlation in the Electronic Structure of Solids, from Silicon to Cuprous Oxide: GW Approximation and beyond[END_REF].

which nevertheless are what we actually want to compute in the end. This situation may look like an impasse. However, as for DFT, where the KS equations separate the interacting and non-interacting part, including the whole latter in the exchange correlation term V xc , in MBPT a quasiparticle equation can be set where a self energy term Σ( r, r , ω) contains the particle correlation. This object is more complicated than the local V xc of the KS theory, giving rise to complex poles and particle life-times. The quasiparticle equation [START_REF] Onida | Electronic excitations: densityfunctional versus many-body Green's-function approaches[END_REF] can be written as:

H 0 ( r) + V H ( r) ψ QP s ( r) + Σ( r, r , ω) ψ QP s ( r ) d r = QP s ψ QP s ( r) (3.9)
where QP s are the poles of the Green function, which give the addition and removal energies of the system. Due to the complex form of Σ, the poles QP s are now complex energies. Σ here act as a non-hermitian, non-local, and frequency dependent operator. V H ( r) is instead the Hartree potential. As in the DFT case, the crucial point is to find a suitable approximation for Σ.

Hedin equations and GW approximation

In 1965 Lars Hedin proposed a set of equations [START_REF] Hedin | New method for calculating the one-particle Green's function with application to the electron-gas problem[END_REF][START_REF] Onida | Electronic excitations: densityfunctional versus many-body Green's-function approaches[END_REF] which implicitly define the exact self energy Σ:

G(1, 2) = G 0 (1, 2) + G 0 (1, 3) Σ(3, 4) G(4, 2) d(34) (3.10) Γ(1, 2; 3) = δ(1, 2) δ(1, 3) + ∂Σ(1, 2) ∂G(4, 5) G(4, 6) G(7, 5) Γ(6, 7; 3) d(4567) (3.11) χ(1, 2) = -i G(1, 3) G(4, 1) Γ(3, 4; 2) d(34) (3.12) W (1, 2) = V C (1, 2) + V C (1, 3) χ(3, 4) W (4, 2) d(34) (3.13) Σ(1, 2) = i G(1, 3) Γ(3, 2; 4) W (4 , 1)d(34) (3.14)
To write the Hedin equations we have adopted a different notation where integer numbers denote time and position, that is for instance G(1, 2) = G( r, t; r , t ). The first equation (eq. 3.10) is called the Dyson equation [START_REF] Stefanucci | Nonequilibrium Many-Body Theory of Quantum Systems, A Modern Introduction[END_REF], which implicitly connect the non-interacting and interacting Green functions.

We have called G 0 the non-interacting Green function 3.3, which is known, and G the unknown interacting Green function 3.4. The two are connected by means of the self energy Σ, which contains the particle correlation neglected in G 0 , but taken into account in G.

When a particle (electron or hole) is added the electronic structure is modified and the the system is polarized. The polarizability χ depends on the propagation dynamics of the electron-hole pairs. If the propagation is not characterized by collisions χ turns out to be the product of the propagators G(1, 2) G(2, 1). In the general case χ is given by equation 3.12, where the vertex function Γ (eq. 3.11) takes into account the effect of collisions and interaction between electron-hole pairs. The polarizability defines how the bare Coulomb interaction is modified and transformed into a screened interaction W (eq. 3.13), which is the interaction that the particles effectively perceive. Equation 3.14 finally provides the self energy as convolution of G, W , and Γ.

In principle starting from a first guess for Σ the Hedin equations permit to

The many body Landauer transport compute G, then Γ and χ, and finally W and Σ. The loop can then be iterated to obtain Σ with arbitrary high precision. This procedure is shown in figure 3.3 a). In practice solving the complete set of Hedin equations is unfeasible. Taken from [START_REF] Darancet | Théorie et simulation du transport quantique dans les nanostructures[END_REF].

The vertex function Γ especially is inaccessible, since it contains the greatest unknown part of the problem. The calculation of the functional derivative Σ(1, 2)/∂G(4, 5) is the main difficulty. Therefore in the GW approximation the vertex function is extremely simplified, neglecting the integral which contains the second order interaction terms between pairs. that is:

Γ GW (1, 2; 3) = δ(1, 2) δ(1, 3)
Actually this also corresponds to a development at first order of the screening

W : W (1, 2) = V C (1, 2) + W (1, 3) χ(3, 4) V C (4, 2) d(34) (3.15)
The polarizability reduces to:

χ(12) = -i G(1, 2) G(2, 1+) ≡ χ GW -RP A (3.16)
and finally the self energy becomes:

Σ(1, 2) = i G(1, 3) W (2, 3) ≡ Σ GW (3.17)
which gives the name to the approach. The GW approximation is justified if the sum over the neglected higher terms in the W development are small with respect to the first order terms of Σ GW and χ GW . This new set of equations, which constitute the GW approximation of the Hedin equations, can be solved iteratively in a self-consistent loop (named GW , see fig. In practice both approaches require an initial guess for Σ 0 and G 0 . The choice of the first guess is particularly delicate for G 0 W 0 , while it is supposed to become irrelevant after few iterations in the self-consistent GW . It is today well known that the exchange correlation potential LDA V LDA xc is a very good local and static approximation for the exact many body Σ. Even if there is no chance to obtain the true electronic structure by using a local and static approximation, the latter can constitute a good starting point for the QP calculation. Therefore the initial non-interacting propagator G 0 is in practice always built by equation 3.7 by using the KS eigenstates ψ KS i ( r), that is G KS 0 . On the contrary HF is proved to be a much better initial guess for molecules and isolated systems, but a worse guess for solids. Once G KS 0 is computed the polarizability is easily obtained by equation 3.16, while the screening calculation requires huge computational resources, and in practical implementations it is usually computed by a different version of equation 3.15:

W ( r 1 , r 2 , ω) = -1 ( r 1 , r 3 , ω) V C ( r 3 , r 2 ) d r 3 = 1 (2π) 3 BZ G 1 , G 2 e i( q+ G 1 ) r 1 -1 G 1 , G 2 ( q, ω) 4π ( q + G 2 ) 2 e -i( q+ G 2 ) r 2 d q (3.18)
where is the dielectric function of the system. The second equality is performed transforming the Coulomb interaction into the Fourier space. q are vectors in the first Brillouin zone (BZ), while G 1 and G 2 are the reciprocal lattice vectors. Computing -1

G 1 , G 2
( q, ω) is extremely onerous because the sum over the reciprocal lattice vectors of equation 3.18 requires a high number of terms to converge. This number is often very far from some empty states close to the Fermi energy, and makes badly scale the resources needed for the calculation. Since this step is very heavy, in self-consistent calculations usually it is performed only in the first iteration and then kept fixed to the W 0 value, while G is updated at each iteration. The resulting approach is often called self-consistent GW 0 . Once W 0 is computed the last step of the first iteration is to compute Σ GW . This task requires in principle to integrate over an infinite number of fre-The many body Landauer transport quencies:

Σ GW ( r, r ; ω) = i 2π W 0 ( r, r ; ω ) G 0 ( r, r ; ω -ω ) e iηω dω (3.19)
In practical implementations this integral is approximated mainly by two different methods: The plasmon-pole model, and the contour deformation method. The plasmon pole model is based on the observation that the main contribution to Im [ -1 ( q, ω)] comes from the plasmon excitation peak, which is modeled by a function describing its width and a delta defining its position. The real part of -1 ( q, ω) is then obtained via Kramers-Kronig relation. Some plasmon pole models are: Hybertsen-Louie [START_REF] Hybertsen | Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies[END_REF], Godby-Needs [START_REF] Godby | Metal-insulator transition in Kohn-Sham theory and quasiparticle theory[END_REF], Von der Linden-Horsch [START_REF] Der Linden | Precise quasiparticle energies and Hartree-Fock bands of semiconductors and insulators[END_REF] and Engel-Farid [START_REF] Engel | Generalized plasmon-pole model and plasmon band structures of crystals[END_REF]. The contour deformation method (details in [START_REF] Giantomassi | Core-electrons and self-consistency in the GW approximation from a PAW perspective[END_REF]) is a more robust approach which transforms the real axis frequency integral into an integral over a contour in the complex plane (residue theorem [START_REF] Rudin | Real and Complex Analysis[END_REF]) plus contributions from outside this contour.

Once the first iteration is completed Σ GW can be used to solve the QP equation 3.9 in the G 0 W 0 approximation in a perturbative way starting from the KS equations. This is feasible by observing that the QP equation is very similar to the KS equation 2.34. Moreover it is nowadays well known that ψ QP ψ KS is a good and justified assumption. Therefore, neglecting the non-Hermitian part of Σ GW which gives rise to QP life times and imaginary eigenvalues, that is only considering the diagonal part of Σ GW , we can write:

ψ KS i Σ GW ( QP i ) ψ KS i -ψ KS i V xc ψ KS i = QP i -KS i (3.20)
Since Σ GW has to be evaluated at the QP energies, which are unknown, it is convenient to expand it around the KS energies. Stopping at the first order gives:

Σ GW ( QP i ) Σ GW ( KS i ) + ( QP i -KS i ) ∂Σ GW (ω) ∂ω hω= KS i (3.21)
By substituting equation 3.21 into equation 3.20 [START_REF] Hybertsen | Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies[END_REF][START_REF] Hybertsen | First-principles theory of quasiparticles: Calculation of band gaps in semiconductors and insulators[END_REF]:

QP i = KS i + Z ψ KS i Σ GW ( KS i ) ψ KS i -ψ KS i V xc ψ KS i Z = 1 -ψ KS i ∂Σ GW (ω) ∂ω hω= KS i ψ KS i -1 (3.22)
where Z is a normalization factor. Equation 3.22 constitutes the G 0 W 0 approximation, which only corrects the KS energies. On the contrary the eigenfunctions are not updated, and remain equal to the KS eigenstates. Differently, to go on with a self-consistent GW , the self energy Σ GW can 3.2.

Step 2: Electronic structure correction 67 be used to compute a new V xc exchange correlation functional, to solve selfconsistently the KS equations by using this new V xc , and then to recompute a new self energy by the GW approximated Hedin equations until convergence is reached. The new V xc can be written as function of Σ GW as [START_REF] Kotani | Quasiparticle selfconsistent gw method: A basis for the independent-particle approximation[END_REF]:

V xc = 1 2 ij |ψ i {Re[ ψ i |Σ GW ( i )|ψ j ] + Re[ ψ i |Σ GW ( j )|ψ j ]} ψ j | (3.23)
which represents the Hermitian part of Σ GW . At each step, the wavefunctions |ψ i are updated via the non-diagonal elements of Σ GW . Finally the initial KS states evolve into QP states. The number of iterations required to reach a good convergence is usually very low (often < 5). This is particularly due to the Hedin idea to use W instead of the bare V C , which effectively helps the perturbation theory in converging quickly.

Applying the GW approximation to real system is a formidable task, which requires huge computational resources. Nevertheless, nowadays GW is a standard technique to calculate more accurate gaps and electronic structures in solids, implemented in several plane-waves codes (Quantum Espresso, Abinit, VASP). with experiments, whereas DFT-LDA systematically underestimates and HF systematically overestimates the band gap values. The important point to stress for transport is that the QP eigenstates, together with the QP eigenvalues, once expressed on a Wannier functions basisset in step 3, are formally suitable to build the Green function by equation The many body Landauer transport 2.33, and the electrodes self energies by equation 2.59, which assume a direct physical meaning. On the contrary in the "DFT+Landauer" approach, where the electronic structure has no real physical meaning, also the Green function and the electrodes self energies does not have a physical meaning. The subtle point is that equation 2.33 holds in principle for non-interacting particles. Nevertheless, equation 2.33 is now applied by using the QP eigenstates and values, corrected by the self energy Σ GW , instead of by using the KS ones. Consequently equation 2.33 becomes formally equivalent to the interacting Green function 3.4, and suitable to compute a physical, quantitatively reliable transmission function and conductance value in the end.

The COHSEX approximation

COHSEX (Coulomb hole screened exchange) is an approach based on HF which uses a static screening, modified into a sum of an instantaneous Coulomb term and a polarization term of response, to better describe the Coulomb interaction in bulk materials. Due to the approximation on the screening, which requires to sum only on the occupied states, COHSEX needs less computational resources than GW , but is also less precise in reproducing experiments. Indeed, as HF, it systematically overestimates the band gap. However, the error is always smaller than HF. Following [START_REF] Gordillo | Many-body perturbation theory and maximally-localized Wannier functions: a combined tool for first-principles electronic structure and quantum transport calculations[END_REF], by Fourier transforming from energy to time W and G in equation 3.19, the self energy Σ GW can be rewritten as:

Σ( r, r ; ω) = i 2π W ( r, r ; τ ) G( r, r ; τ ) e iωτ dτ (3.24)
Introducing the approximation of static screening, that is:

W ( r, r ; τ ) = 1 2π W ( r, r ; ω = 0) δ(τ ) (3.25)
substituting equation 3.25 into equation 3.24:

Σ( r, r ; ω) = i 2π W ( r, r ; ω = 0) G( r, r ; τ ) δ(τ + η) e iωτ dτ (3.26)
where η is a small positive number. Finally the screened exchange (SX) self energy is obtained by introducing and rearranging the MB definition of Green function 3.4 in equation 3.26 [START_REF] Gordillo | Many-body perturbation theory and maximally-localized Wannier functions: a combined tool for first-principles electronic structure and quantum transport calculations[END_REF]:

Σ SX ( r, r ) = -W ( r, r ; ω = 0) Ψ 0 | C † ( r ) C( r)|Ψ 0 (3.27)
where the frequency dependency was dropped due to the fact that all of the terms e i(hω-s)η → 1, when η → 0. The SX approximation 3.27 contains only a screened-exchange term. This approximation is similar to HF, but the bare Coulomb interaction V C is replaced by a static screened Coulomb interaction W = V C -1 . COHSEX goes beyond the SX approximation by splitting the screening in two terms: W = V C + W p , where V C is the Coulomb part, which is instantaneous, while W p is the polarizable part, which has some frequency dependency. This modification makes an additional Coulomb hole term [START_REF] Hedin | Solid State Physics[END_REF] appear in the self energy expression. Substituting W = V C + W p into equation 3.24:

Σ( r, r ; ω) = i 2π [V C ( r, r ) δ(τ +η)+W p ( r, r ; ω = 0) δ(τ )] G( r, r ; τ ) e iωτ dτ
(3.28) The first term is the HF term, which contains contributions only from the occupied states. The second term also contains contributions only from the occupied states, since η is removed from its delta function. As done for the SX approximation the MB Green function 3.4 can be substituted into equation 3.28, giving (details in [START_REF] Gordillo | Many-body perturbation theory and maximally-localized Wannier functions: a combined tool for first-principles electronic structure and quantum transport calculations[END_REF]):

Σ CHSX ( r, r ) = -Ψ 0 | C † ( r ) C( r)|Ψ 0 W ( r, r ; ω = 0) + 1 2 δ( r -r ) W p ( r, r ; ω = 0) (3.29)
The first term is the SX term and the second one is known as the Coulomb hole (CH) term:

Σ CH ( r, r ) = 1 2 δ( r -r ) W p ( r, r ; ω = 0) (3.30)
The CH term 3.30 accounts for a classical response of the medium after the addition of a point charge (electron or hole). This classical term is still a crude approximation from the quantum point of view. Nevertheless, it greatly simplify the problem with respect to the GW approximation, where the Coulomb-hole term is treated in a full quantum mechanics frame.

As for GW , Σ CHSX can be evaluated one-shot, or self consistently until the desired convergence is reached. In the first case COHSEX updates only eigenvalues, while in the second both eigenvalues and eigenfunctions.

The HSE hybrid functional

DFT-LDA and DFT-GGA have the well known inconvenient to systematically underestimate the band gap of insulators. On the contrary, Hartree-Fock is well known to strongly overestimate the same physical quantity. The use of hybrid functionals aim to improve the calculation of the band gap,

The many body Landauer transport of the wavefunctions, and of other physical properties by combining DFT-LDA/GGA with HF. The method is quite empirical, and requires to adjust a mixing parameter, but can provide better results than simple DFT or HF in many cases.

On one hand the typical problems of "classical" DFT functionals are mainly two: they do not cancel the electron self-interaction error of the Coulomb term, and they give a rough representation of non-local correlation phenomena, due to the local and semi-local character of LDA and GGA approximations respectively. On the other hand the HF approximation only contains an exact non-local exchange term which cancels the Coulomb self-interaction, but no electronic correlation term. That given, hybrid functionals incorporate a portion of exact HF exchange term with the complementary portion of LDA/GGA exchange-correlation term. The HF exchange term writes:

V HF x = - 1 2 ψ * i ( r)ψ * j ( r ) 1 r -r ψ j ( r)ψ i ( r ) d r d r (3.31)
In this work we adopt the HSE hybrid functional [START_REF] Heyd | Hybrid functionals based on a screened coulomb potential[END_REF], which uses an error function screened Coulomb potential to calculate the exchange portion of the energy, combined with the PBE exchange-correlation DFT functional [START_REF] Perdew | Generalized gradient approximation made simple[END_REF], which is an improved and simplified version of GGA functionals:

V xc = α V HF,SR x (w) + (1 -α)V P BE,SR x (w) + V P BE,LR x (w) + V P BE c (3.32)
where SR and LR refers to short and long range, α is the mixing parameter, and w is a parameter controlling the border between short-range and longrange regions of interaction. A standard value for the screening parameter is w = 0.2, which is usually referred as HSE06 [START_REF] Krukau | Influence of the exchange screening parameter on the performance of screened hybrid functionals[END_REF]. The screening cancels the singularity of the Coulomb potential at r = r , significantly improving the quality of the functional. A typical value for the mixing parameter is α = 0.25. Nevertheless, it must be always optimized to obtain the best result according to experiments or more advanced QP calculations.

In the HSE approach the KS equations are solved by using the hybrid functional defined by equations 3.32 and 3.31. Usually the KS wavefunctions are used as starting point for the self-consistent calculation. Once convergence is reached the resulting wavefunctions and eigenvalues constitute an improved electronic structure, which is in principle closer to the true, physical one. The values of several physical quantities, like for instance band gap, total energy difference with chemical accuracy, atomization and formation energy, Van der Waals interaction, and others, turn out to be more accurate.

Step 3: Basis-set transformation

The basis-set transformation from plane-waves to WF is probably the most delicate part of the whole procedure. If it turns out to be successful the electronic structure corrected in step 2 on the plane waves basis-set remains well represented also in the WF basis-set, which is then suitable to compute a more reliable MB transmission function (step 4). To perform the transformation we use the Marzari and Vanderbilt formulation [START_REF] Marzari | Maximally localized generalized Wannier functions for composite energy bands[END_REF] of WF, implemented in the Wannier90 code. We now briefly recap the properties of WF in their original formulation, underlying the most important ones for transport calculations. Then we will characterize the Marzari and Vanderbilt formulation (section 3.3.3), which maintains all the properties of the original definition, and is more suitable for practical applications.

Original continuous definition of Wannier functions

The transformation, in the original formulation of WF by Gregory H. Wannier in 1937 [START_REF] Wannier | The structure of electronic excitation levels in insulating crystals[END_REF], is defined for any band n by:

w Rn ( r) = 1 V BZ BZ e -i k • R ψ n k ( r) d k (3.33)
where V BZ = (2π) 3 /V cell is the volume of the first Brillouin zone, ψ n k ( r) are Bloch states, and the integral runs over the first Brillouin zone. R is a vector of the real lattice, defined by: R = p 1 a 1 + p 2 a 2 + p 3 a 3 (3.34) where p i are arbitrary integers and a i are the real lattice basis vectors of the unit cell. WF are not uniquely defined, since even in the absence of degeneracies, each Bloch state is defined up to an arbitrary phase factor. The variation of this phase factor with k impacts the Brillouin zone integral of equation 3.33, hence the shape of the Wannier function. The orthogonality of Bloch states and their normalization in one unit cell imply that WF are orthonormal:

Rn R n = δ R, R δ n,n (3.35)
where Rn is the WF w Rn ( r) in bra-ket notation defined by w Rn ( r) ≡ r Rn . Note that while Bloch states have infinite norm in the infinite crystal, WF are normalized, which suggests that they are localized in some region of space.

Discrete formulation

In practice, in step 2, only a finite number of Bloch states is computed on a finite Monkhorst-Pack (MP) k-points grid of dimensions N 1 × N 2 × N 3 . The number of Bloch states is equal to N × N b , where N = N 1 × N 2 × N 3 is the total number of k-points, and N b is the number of ab initio bands. That given the integral 3.33 is approximated by a discrete sum over k-points:

Rn = 1 N k e -i k• R ψ n k (3.36)
R is always given by equation 3.34, where now a i denote the basis vectors of the ab initio supercell defined in step 1, which in the following will be referred to as the "reference supercell" (RS), and the p i integers assume values from 0 to N i -1. As consequence of the finite MP grid, R does not span over an infinite crystal, like for equation 3.33, but only over an "extended supercell" (ES) which is made by N i replicas of the RS along each direction a i (see fig.

T = p 1 N 1 a 1 + p 2 N 2 a 2 + p 3 N 3 a 3 (3.37)
Indeed the k-points of the MP grid are of the form:

k = k 1 N 1 b 1 + k 2 N 2 b 2 + k 3 N 3 b 3
where k i are integers, and b i are the reciprocal lattice vectors, defined by:

a i • b j = 2πδ i,j
Hence k • T is a multiple of 2π, and equation 3.36 yields:

w Rn ( r + T ) = w Rn ( r)
When the k-points density increases, the ES size increases as well, WF replicas are pushed further, and the infinite crystal limit of equation 3.33 is approached. As we will see in chapter 5 the dimensions of the ES, and thus the distance between closest PBC Wannier replicas, can in some cases be crucial in the calculation of the Wannier Hamiltonian matrix. This turns out to be

w Rn ( r) = 1 N k e -i k• R ψ n k ( r) = 1 N k e i k•( r-R) u n k ( r -R) = w 0n ( r -R) (3.
38) where we have used the definition of Bloch state ψ n k ( r) = e i k• r u n k ( r), and the fact that u n k ( r) = u n k ( r + R) is periodic along the RS. In practice, this means that a given WF |0n is a linear combination of the Bloch states belonging to the band n, which interfere destructively almost everywhere into the ES except in a localized area into the RS. Any WF Rn is also a linear combination of the same Bloch states, but where each Bloch state is shifted by a different phase factor e i k• R given by k which changes in the sum 3.36. Since the product k • R is always a fraction of 2π the non-trivial result is that the constructive interference of Bloch states maintains the same shape and translates of R outside the RS. On the contrary into the RS the interference becomes destructive and Rn goes to zero (see figure 3.5). Equation 3.38 constitutes the basis for an efficient interpolation at arbitrary k-points, which is particularly useful to compute the band structure or to converge a transmission function calculation with respect to the k-points grid transverse to transport. To understand how interpolation works consider the inverse transformation of equation 3.36, which permits to express the Bloch states in terms of WF. It can be found considering that:

Rn ψ n k = 1 N k e i k • R ψ n k ψ n k = e i k• R
where the last equality follows from the orthogonality of Bloch states and their normalization to N into the ES. Also remark that, in this discrete formulation, the scalar products are integrals over the ES. The inverse WF transformation thus reads:

ψ n k = R e i k• R Rn (3.39)
The sum in equation 3.39 is in principle restricted over the WF into the ES, which corresponds to the Bloch states computed on the original MP grid. However, equation 3.38 easily permits to build WF outside the ES to satisfy Bloch states PBC over a larger ES. In practice, this is obtained by applying a phase factor, and thus by simple translation of the WF belonging to the RS outside the ES. These new WF can extend the sum 3.39 to obtain Bloch states on a finer MP grid. This is similar in the spirit to the interpolation method of Slater and Koster [START_REF] Slater | Simplified LCAO method for the periodic potential problem[END_REF], based on linear combinations of atomic orbitals. When computing the transmission function the periodicity is broken along the transport direction, whereas the surface transverse to transport is enlarged with respect to the RS according to the number of transverse k-points. Thus the transmission function computed on top of the basis-set transformation is actually the one through this enlarged lateral surface of the ES along the transport direction. To go back to the RS transverse surface an average over the transverse k-points has to be done. The obtained transmission function should then be interpreted as an average transmission per unit of surface. This average procedure often requires a denser k-points grid to converge with respect to the grid used for DFT. Therefore the interpolation property of WF is very important to save computational resources, which can badly scale with the number of k-points in many body calculations like for instance GW in step 2. One can perform the many body calculation by using the minimum necessary number of k-points, and then refine the grid for the transport calculation. As briefly discussed above WF are normalized, and thus localized. Furthermore, the variation of the Bloch states arbitrary phase factors affects the WF shape, and thus their localization. The WF property of being localized is very important for electronic transport since it allows to compute Hamiltonian sub-matrices representing localized areas of the RS. In general, the more the WF are localized, the more the transmission function calculation becomes qualitatively good and reliable. It is hence important to choose the Bloch states arbitrary phase factors in order to maximize the WF localization. It can be shown that if the Bloch states are infinite differentiable with respect to k, and if they satisfy boundary conditions at the edges of the first Brillouin zone ψ n, k+ G = ψ n, k , the WF decay "faster than power law", that is in fact exponentially. Bloch states are eigenvectors of the Kohn-Sham Hamiltonian operator. This operator varies analytically with respect to each of the 3 components of k considered as complex variables, and hence it is infinitely differentiable. Forgetting mathematical rigor also its eigenvectors are expected to be infinitely differentiable. However, the eigenvectors obtained from numerical solvers have an arbitrary phase factor which may be discontinuous with k. Obtaining well-localized WF requires changing the phase factors in order to restore the smoothness of the Bloch states.

Marzari and Vanderbilt formulation

The original definition of WF described up to now is in practice suitable only in the absence of bands degeneracy, that is if all the bands and consequently

The many body Landauer transport the Bloch states are well isolated in energy. Indeed when the sum 3.36 arrives at a degeneracy point while running over k at fixed band index n, a numerical ambiguity occurs regarding which k path to follow to not change the band index n. To solve this ambiguity consider that any linear combination of degenerate Bloch states is also a Bloch eigenstate of the Hamiltonian. Thus a possible solution is to choose, for each band, at each degeneracy point, the linear combination which gives the best matching (i.e. continuity) with the Bloch state at the neighboring, non-degenerate k-points. This operation becomes difficult if there are many bands and degeneracies, as happens in the large supercells here considered with interfaces, defects, and different materials. The more general approach developed by Marzari and Vanderbilt [START_REF] Marzari | Maximally localized generalized Wannier functions for composite energy bands[END_REF] is to allow any linear combination of bands at every k-point, even for non-degenerate eigenstates, and to apply a minimization algorithm to find the WF set with the best localization. In this formulation, the transformation between Bloch states and WF becomes:

| Rn = 1 N k e -i k• R N b m=1 U ( k) mn |ψ m k (3.40)
While equation 3.36 directly transforms Bloch states into WF, equation 3.40 builds WF from the states:

ψ n k = N b m=1 U ( k) mn |ψ m k (3.41)
which are linear combinations of Bloch states defined by the matrices U ( k) . Since the mixing between non-degenerate k-points is allowed, ψ n k are in general not eigenstates of the Hamiltonian. Nevertheless, the transformation U ( k) is unitary and preserves all the properties valid for the WF in the original formulation 3.36. Indeed these properties in practice only follow by the orthogonality condition 3.35, which is preserved also for ψ n k . The matrices U ( k) are optimized to minimize the sum of the quadratic spreads of the WF, which in the continuous limit (infinite number of k-points) reads:

Ω = n 0n|r 2 |0n -| 0n|r|0n | 2 (3.42)
The minimization of the spread assures to obtain maximally localized WF. This is also related to finding the best continuity of the Bloch states over k, at both degenerate and non-degenerate points, and thus to conveniently Mathematically speaking the space spanned by the Bloch states contains the space spanned by the WF, which in turn contains the space spanned by the Bloch states in the frozen window. N w can be in principle always chosen equal to N b , so that U ( k) become squared matrices and all the "ab initio" Bloch states are reproduced. Nevertheless, this considerably increases the constraints in the linear combination of eigenstates and can often lead to unreliable results. The solution is thus to compute more "ab initio" Bloch states than the desired frozen window to reproduce by WF, which is usually quite easy and fast by plane-waves codes.

For transport, according to voltage and temperature, only a tight energy range around the Fermi energy is relevant. Hence the frozen window can be adapted to well describe only states in an appropriate range around E f . This is useful for several reasons: the optimization of the WF basis-set becomes easier if less Bloch states are required to be exactly reproduced; high energy Bloch states in the conduction band, which can often lead to less localized WF, can be excluded; and the size of the WF basis-set can be reduced, lowering the computational cost of the transmission function calculation (matrix operations like inversion scale as the cube of the matrix size). As all minimization algorithms, the optimization of U ( k) to minimize the spread 3.42 requires an initial guess as a starting point. The latter is built by projecting a set of N w initial trial orbitals |Φ n in the RS onto the subspace spanned by the N b × N Bloch states, obtaining some rough initial trial WF:

Φ n = N b m, k ψ m k ψ m k |Φ n (3.43)
The many body Landauer transport Φ n are then orthonormalized via the Löwdin symmetric orthogonalization procedure [START_REF] Löwdin | On the non orthogonality problem connected with the use of atomic wave functions in the theory of molecules and crystals[END_REF]:

Φ ⊥ n = N b m (S -1/2 ) mn Φ n (3.44)
where

S mn = Φ m | Φ n . Φ ⊥ n
are the trial WF which, substituted in the left hand side of equation 3.40 at the place of |0n , define the initial matrices U ( k) . The choice of the initial trial orbitals |Φ n , which is often non-trivial, is crucial to obtain a good WF basis-set. However, if the choice is good the WF obtained as the projection of the trial orbitals, without performing any minimization step, can constitute a very good basis-set. This property can become very useful for transport calculation. Indeed in the present methodology, as we will see in step 4, to build the Hamiltonian of the periodic electrodes the WF basis must satisfy the periodicity of the electrodes. The WF of a cubic cell of the electrode must be as similar as possible to the ones of another cubic cell and must be indexed in the same order. The easiest way to ensure this periodicity is to directly use Φ ⊥

Step 4: Wannier Hamiltonian construction and transmission function calculation

In this step the WF basis-set of step 3, on which the QP electronic structure computed in step 2 is well represented, is used to extract the Hamiltonian matrices required to build the electrodes self energies 2.59 and the Green function of the open system 2.58. To fulfill this task we have written an original code inspired to the transport module of Wannier90, which builds the Hamiltonian matrices and implements the Green function formalism of section 2.2.3. As mentioned in the previous section, the transmission function and conductance calculation is very sensitive to the number of k-points transverse to transport. Therefore, contrary to Wannier90, which only permits to compute the transmission function at the Γ-point, our new implementation allows to treat and interpolate transverse k-points. In the next sections we firstly introduce how we build a Wannier Hamiltonian, starting from the WF basis-set, which guarantees a good k-points interpolation and a coherent physical meaning. Then we explain how the electrodes and central region Hamiltonian matrices are chosen and computed.

Tight-binding-like Hamiltionian on the Wannier functions basis-set

From equation (3.40), the expression of the QP Hamiltonian in the Wannier basis, which represents the periodic ES, is written as:

Rn H R n ES = 1 N k e i k•( R-R ) N b m=1 U ( k) * mn m k U ( k) mn (3.45)
where m k are the Bloch states energies. We instead prefer a different tightbinding-like Hamiltonian construction, proposed in the Wannier90 user guide, which represents the infinite crystal with finite range coupling instead of the periodic ES. For this purpose localized orbitals |φ Rn are defined as the restrictions of each WF | Rn to the Wigner-Seitz (WS) cell centered around R + r n , where r n is the average position of each WF |0n . The WS considered here is the one associated to the lattice spanned by the lattice vectors N i a i , which has the same volume of the ES. In other words, each Wannier function is cut at the border of the WS cell, and cleared from its periodic replicas (see fig. ). In practice, the matrix elements are tabulated with R running over the ES. For each one, we search for the vector T such that ( R + r n + T ) is closest to ¯ r n , i.e. within the WS cell centered around r n . This yields the following matrix element between orbitals:

φ 0n H φ ( R + T )n ≡ 0n H R n ES (3.46)
In the left-hand side, the matrix element is meant as an integral over the infinite crystal (which restricts to the finite volume where the orbitals are non-

φ Rn |H|φ R n = φ 0n H φ ( R -R)n (3.47)
These conventions for defining the Hamiltonian set a cutoff distance: each orbital is only coupled to neighbors within its WS cell. Couplings within a limited region of space are important for transport, since we should be able to extract Hamiltonian matrices of partial regions only interacting with first neighbor regions. However, equation 3.46 is inaccurate for couplings at distance in the order of the WS cell size. If more accuracy is needed, the MP grid must be refined, which increases the WS cell size. When refining the MP grid in the ab initio simulation, the size of the WS cell is increased, hence more distant couplings are included. As will be shown in the results chapter, the Hamiltonian matrix elements decrease rapidly with the distance, due to the localization of the orbitals, and distant couplings may be neglected. Moreover the tight-binding-like Hamiltonian (eq.s 3.46 and 3.47) yields a more accurate interpolation of the band structure on k-points not belonging to the MP grid of the ab initio calculation [START_REF] Souza | Maximally localized Wannier functions for entangled energy bands[END_REF][START_REF] Marzari | Maximally localized Wannier functions: Theory and applications[END_REF] (similar to Slater and Koster [START_REF] Slater | Simplified LCAO method for the periodic potential problem[END_REF]). The main reason is that, while the phase factors of all the k-points belonging to the ab initio MP grid equal 0 modulo 2π, the phase factors of the interpolated k-points can in general assume any value. Thus, while the coupling matrix elements do not change for MP k-points if considering original WF or replicas, they can considerably change for interpolated k-points. The tight-binding-like Hamiltonian, which is built by localized orbitals with couplings in a finite region (the WS cell), has smaller phase factors than the "direct" Hamiltonian matrix 3.45, and is proved to guarantee a better stability. The interpolation on the Hamiltonian 3.45 can indeed introduce non-physical oscillations due to interference. To practically implement the interpolation by using the tight-binding-like Hamiltonian we consider k-points associated to a finer

N f 1 × N f 2 × N f 3 MP grid.
Bloch states are expressed on the following basis set:

|n k = 1 √ N f R e i k•( R+ rn) |φ Rn (3.48)
where

N f = N f 1 × N f 2 × N f 3
, and R spans the infinite crystal. This basis-set is orthonormal in the N f 1 × N f 2 × N f 3 extended supercell ES f , and for each k-point the Hamiltonian matrix elements read:

H ( k) nn ≡ n k H n k ES f = 1 N f R, R e i k•( R + r n -R-rn) φ Rn |H|φ R n ES f = R e i k•( R + r n -rn) φ 0n |H|φ R n (3.49)
In the last equality, obtained by translation symmetry, the matrix elements are meant as integrals over the infinite crystal, as defined in equation 3.46. R runs over all coupled orbitals. The diagonalization of H ( k) for each k yields the interpolated band structure, which is often very accurate as long as the ab initio MP grid is fine enough to obtain a WS cell much larger than the WF mean spread. It can be shown that if k belongs to the ab initio MP grid, one recovers exactly the ab initio Bloch states in the frozen energy window. As already mentioned in the previous section we are interested in computing the conductance per unit area through the considered junction. The reason is that the transmission function strongly depends on the wavevector parallel to the interfaces, k = k 2 b 2 + k 3 b 3 (note that here we are assuming a 1 as transport direction). Therefore an average over the transmission functions at each transverse k-points becomes necessary. However, the ab initio MP grid is rarely fine enough to resolve this dependence, and should be refined by interpolation. Consider thus an extended supercell ES made of N 2 × N 3 repetitions of the ab initio supercell along the transverse directions a 2 and a 3 . The many-channel Landauer transmission is computed through this ES with periodic boundary conditions along a 2 and a 3 . If the surface is large enough, the conductance per unit area becomes nearly insensitive to the boundary conditions and convergence is reached. The transverse k-points compatible with these conditions are:

k = q 2 N 2 b 2 + q 3 N 3 b 3 (3.50)
where q i are integers. The wavefunctions are developed on the following basis, indexed by the lattice vector R ⊥ = p 1 a 1 along the transport direction, the orbital index n, and the transverse wavevector k :

| R ⊥ n k = 1 √ N R e i k •( R + rn) |φ ( R ⊥ + R )n (3.51)

Step 4: Wannier Hamiltonian construction and transmission function calculation 83

where N = N 2 × N 3 , and R runs over the infinity of transverse lattice vectors: R = p 2 a 2 + p 3 a 3 (3.52)

These basis states are localized along the transport direction a 1 , around position R ⊥ + r n , and extended along the transverse directions. They form an orthonormal basis when integrated over the ES along the transverse directions and over infinity along the transport direction. The Hamiltonian only couples states with the same wavevector, and its matrix elements read:

H ( k ) R ⊥ n R ⊥ n = R e i k •( R + r n -rn) φ R ⊥ n Ĥ φ ( R ⊥ + R )n (3.53)
where R is a transverse lattice vector and runs over all coupled orbitals. Equation 3.53 is exactly how we compute the Wannier Hamiltonian for transport calculations.

In the code we also introduced the possibility to discard the matrix elements beyond a cut-off distance L cut along the transport direction. The discarding condition reads:

R ⊥ + r n -R ⊥ -r n • â1 > L cut (3.54)
where â1 is the unit vector along the transport direction. As we will see in chapter 5, when the dimensions of the ES are too small along the transport direction, the tails of the orbitals |φ R ⊥ n raise. This is due to the fact that WF are localized, but also periodic along the ES. The tails, which well describe the periodic system, but are overestimated for the semi-infinite one (see section 3.1), produce anomalous high non-diagonal terms in the Hamiltonian matrix when applying eq. 3.49, invalidating the transmission function calculation. The extreme situation is when the ES coincides with the RS, that is N 1 = 1. The best solution, from a formal point of view, would be to refine the ab initio grid along transport (i.e. to increase N 1 ), to increase the size of the ES. However, this solution often entails a great increase of computational resources. Since the coupling between orbitals should quickly decay with distance, the corresponding Hamiltonian matrix elements should also quickly tend to zero. Thus another solution, which does not entail any increase of computational resources, is to discard the non-diagonal elements beyond a cut-off distance, to remove the anomalous terms.

Wannier Hamiltonian matrices for transport

As pointed out in section 2.2.3, to compute the transmission function by the Green function -Landauer method, we must be able to extract sub-matrices

The many body Landauer transport representing the central region, the principal layers of the electrodes, and the couplings between these matrices. In this section we explain how we choose and how we compute all these matrices. As pointed out in section 3.1 (step 1), the central region should contain the insulator plus some cells of metal. Hence the central region is chosen as the whole RS. For each transverse wavevector (omitted to lighten notations), the N w × N w Hamiltonian matrix H C of the central region reads:

H C nn = H ( k ) 0n 0n (3.55)
where the left hand side is given by equation 3.53. This means that the couplings between the RS ( R ⊥ = 0) and its replicas along the transport direction are removed. The Hamiltonian of the periodic electrodes, as described in 3.1, has to be extracted far enough from the insulator such that its influence is screened. We also need to define a cutoff distance beyond which the couplings between orbitals in the metal electrode are negligible. For the structures studied in this work the couplings between one cubic unit cell of metal and its 4 th neighbors and beyond always turned out to be negligible. Anyway a convergence study should always be performed when considering a new system. We define ( R ⊥ = a 1 ) are made of the relevant matrix elements:

H ( k ) 0n a 1 n (3.57)
Note that the 0i indexes of H 0i are related to electrode unit cells, and not directly related to the RS replica indexes of equations 3.56 and 3.57. All the blocks H 0i have dimension N cc w × N cc w , where N cc w is the number of Wannier functions in a single metal cubic cell. They are used to build the Hamiltonians of the electrodes, which are block-tridiagonal, each block corresponding to 3 cubic cells. The right electrode Hamiltonian reads:

H R =    D τ τ † D . . . . . . . . .    (3.58)
with the diagonal block:

D =   H 00 H 01 H 02 H 01 † H 00 H 01 H 02 † H 01 † H 00   (3.59)
and the coupling block:

τ =   H 03 0 0 H 02 H 03 0 H 01 H 02 H 03   (3.60)
The Hamiltonian of the left electrode is defined with the same blocks. τ also defines the coupling between the electrodes and the central region. Note that this construction of the electrode Hamiltonian is valid as long as the WF satisfy the periodicity of the metal cubic cell. As mentioned in section 3.3 (step 3) this is assured by a compromise in the WF spread minimization procedure, that is in the optimization of the WF basis-set. More iterations are used in the minimization, more the spread turns out to be low, but less the periodicity is assured, since WF may be more hybridized between the starting guess orbitals. On the contrary less iterations translate in higher spreads, but more periodic WF positions.

Finally the transmission function is obtained by the Green function method like in the standard "DFT + Landauer" approach. Following section 2.2.3 the decimation algorithm of Lopez-Sancho et al. [START_REF] Lopez-Sancho | Quick iterative scheme for the calculation of transfer matrices: application to Mo (100)[END_REF][START_REF] Lopez-Sancho | Highly convergent schemes for the calculation of bulk and surface green functions[END_REF] is used to compute the surface retarded Green's functions at energy E of the left and right electrodes, g L (E) and g R (E). This yields the left and right electrodes self-energies:

Σ L (E) = τ † g L (E)τ (3.61) Σ R (E) = τ g R (E)τ † (3.62)
The many body Landauer transport

These blocks are added respectively to the upper diagonal block (first 3 metal cubic cells) and lower diagonal block (last 3 metal cubic cells) of H C to obtain the retarded Green's function in the central region:

G(E) = E -H C -Σ L (E) -Σ R (E) -1 (3.63)
In the end, the sum of transmission coefficients is obtained from the trace formula [START_REF] Meir | Landauer formula for the current through an interacting electron region[END_REF]:

T ( k , E) = Tr Γ L (E)G(E)Γ R (E)G † (E) (3.64) where Γ L,R (E) = i Σ L,R (E) -Σ † L,R (E) (3.65)
The sum of transmission coefficients per metal cubic cell area is:

T (E) = 1 N k T ( k , E) (3.66)
where k runs over the refined N 2 ×N 3 MP grid. The N 2 and N 3 values must be adjusted to obtain good convergence of T (E). The zero-bias conductance per metal cubic cell area is given by the Landauer formula:

G = G 0 T (E) - ∂η f ∂E dE (3.67)
where η f is the Fermi distribution 2.4. Since the entire procedure formally holds only in equilibrium conditions, the actual implementation does not give the possibility to the user to compute the transmission function at non-zero voltages. The Poisson method (section 2.2.4) is usually not implemented in plane-waves DFT codes, and would be anyway not correct to apply. The conductance value provided by equation 3.67 is at zero bias and temperature θ given by the Fermi distribution η f , but anyway suitable to approximate the current as I = G∆V in the linear response regime (∆V 2k B θ/e, eq. 2.20). For larger voltages a full NEGF frame would be necessary.

Chapter 4 HfO 2 GW study

This chapter is the first of three results chapters, which constitute a summary of the main studies and researches carried out during the three years of PhD. This first chapter of results is a deep study of the electronic properties of hafnium oxide (HfO 2 ), which is among the most promising materials for real technology applications such as OxRAMs. The study is carried out by advanced GW ab initio simulations, where convergence and theoretical framework are pushed to the present possible limits. The second chapter of results has the purpose to verify the robustness and to analyze preliminary results of the "many body Landauer " method for bulk materials developed in chapter 3. Here we are especially interested in the difference between different QP methods used to compute the transmission function. We consider DFT, which should match state of the art calculations (e.g. Transiesta), G 0 W 0 , HSE and self consistent COHSEX. For this testing purpose, we consider a junction made of silver and magnesium oxide (Ag/MgO/Ag), which is easier to simulate than more realistic OxRAM devices (e.g. TiN/HfO 2 /Ti, Pt/HfO 2 /Pt, or TiN/HfO 2 /Pt). The more realistic OxRAM devices are considered in the third chapter of results, which is dedicated to transport calculations on a HfO 2 junction with platinum electrodes, that is Pt/HfO 2 /Pt. The Pt/HfO 2 /Pt supercell built in this section is still a rough model of real OxRAM, which is the result of a compromise to save computational resources. However, simulations provide interesting results and information. A form of conductive filament (CF) is hypothesized for the high and the low resistive states, and the resulting transmission function, computed by DFT and HSE, is analyzed compared to the pristine state.

Hafnium oxide properties

Hafnium oxide is a transition metal oxide which possesses very interesting optical and electrical insulating properties for several technological applications and research fields. The main interest in this work is for the increasing interest of HfO 2 as the oxide insulating layer in OxRAM devices [START_REF] Wong | Metal-Oxide RRAM[END_REF][START_REF] Bersuker | Metal oxide RRAM switching mechanism based on conductive filament microscopic properties[END_REF]. Other typical applications are optical coatings in the near-ultraviolet to infrared wavelength range [START_REF] Gilo | Study of HfO 2 films prepared by ion-assisted deposition using a gridless end-hall ion source[END_REF], or high permittivity dielectric in submicrometer silicon-based technologies [START_REF] Robertson | High dielectric constant gate oxides for metal oxide si transistors[END_REF]. At ambient pressure bulk HfO 2 presents three stable crystalline structures according to thermodynamical conditions. At low temperature, under 2000 K, the energetically most favorable structure is the monoclinic with P 2 1 /c symmetry [START_REF] Wolten | Diffusionless phase transformations in zirconia and hafnia[END_REF]. Above temperatures around 2000 K HfO 2 observes a phase transition to the tetragonal structure with P 4 2 /nmc. Finally, at very high temperatures the most stable phase becomes the cubic one with F m3m symmetry. Typically the deposited films of HfO 2 are initially amorphous and crystallize after annealing [START_REF] Ho | Morphology and crystallization kinetics in HfO 2 thin films grown by atomic layer deposition[END_REF]. The crystal phase in films is usually composed by a prevalence of the monoclinic phase, with the presence of also tetragonal and metastable orthorhombic phases [START_REF] Ho | Morphology and crystallization kinetics in HfO 2 thin films grown by atomic layer deposition[END_REF][START_REF] Nguyen | Sub-bandgap defect states in polycrystalline hafnium oxide and their suppression by admixture of silicon[END_REF][START_REF] Park | Optical properties of thermally annealed hafnium oxide and their correlation with structural change[END_REF][START_REF] Hill | The relationship between local order, long range order, and sub-band-gap defects in hafnium oxide and hafnium silicate films[END_REF]. The purpose of this chapter is to perform an accurate comparison between experimental and highly converged simulated values of the optical and electronic band gaps of HfO 2 . Therefore it is worth to give an overview of the most used characterization experimental techniques. The electronic structure is generally investigated by means of both X-ray (XPS) and ultraviolet (UPS) photoemission spectroscopy [START_REF] Sayan | Band alignment issues related to HfO 2 /SiO 2 /p-Si gate stacks[END_REF][START_REF] Bersch | Band offsets of ultrathin high-κ oxide films with Si[END_REF]. The optical and dielectric properties are instead usually investigated by using X-ray or optical absorption [START_REF] Balog | Chemical vapor deposition and characterization of HfO 2 films from organo-hafnium compounds[END_REF][START_REF] Toyoda | Chemistry and band offsets of HfO 2 thin films on Si revealed by photoelectron spectroscopy and x-ray absorption spectroscopy[END_REF], spectroscopic ellipsometry (SE) [START_REF] Zhu | Effect of Al inclusion in HfO 2 on the physical and electrical properties of the dielectrics[END_REF][START_REF] Edwards | Status and prospects for VUV ellipsometry (applied to high K and low K materials)[END_REF][START_REF] Modreanu | Optical characterization of high-k dielectrics HfO 2 thin films obtained by MOCVD[END_REF][START_REF] Nguyen | Sub-bandgap defect states in polycrystalline hafnium oxide and their suppression by admixture of silicon[END_REF][START_REF] Park | Optical properties of thermally annealed hafnium oxide and their correlation with structural change[END_REF][START_REF] Hill | The relationship between local order, long range order, and sub-band-gap defects in hafnium oxide and hafnium silicate films[END_REF], and energy-loss spectroscopy (EELS) [START_REF] Yu | Energy gap and band alignment for (HfO 2 ) x (Al 2 O 3 ) 1-x on (100) Si[END_REF][START_REF] Puthenkovilakam | An accurate determination of barrier heights at the HfO 2 /Si interfaces[END_REF][START_REF] Cheynet | Crystal structure and band gap determination of HfO 2 thin films[END_REF][START_REF] Ikarashi | Electronic structure analysis of Zr silicate and Hf silicate films by using spatially resolved valence electron energyloss spectroscopy[END_REF][START_REF] Guedj | Evidence for anisotropic dielectric properties of monoclinic hafnia using valence electron energy-loss spectroscopy in high-resolution transmission electron microscopy and ab initio time-dependent densityfunctional theory[END_REF][START_REF] Hung | Interpretation of monoclinic hafnia valence electron energy-loss spectra by time-dependent density functional theory[END_REF]. Band gap values are then extracted from experimental data by using linearization and extrapolation techniques over measured spectra. The optical values range from 5.1 eV [START_REF] Toyoda | Chemistry and band offsets of HfO 2 thin films on Si revealed by photoelectron spectroscopy and x-ray absorption spectroscopy[END_REF] to 5.95 eV [START_REF] Modreanu | Optical characterization of high-k dielectrics HfO 2 thin films obtained by MOCVD[END_REF], while the electronic values range between 5.7 eV [START_REF] Bersch | Band offsets of ultrathin high-κ oxide films with Si[END_REF] and 5.86 eV [START_REF] Sayan | Band alignment issues related to HfO 2 /SiO 2 /p-Si gate stacks[END_REF]. Surprisingly, the two ranges overlap, making unclear the distinction between the optical and electronic gap. For a complete comparison analysis, we are also interested in theoretical studies of the HfO 2 band gap. LDA or GGA calculations underestimate the experimental band gap by about 30%, giving values from 3.8 eV [START_REF] Grüning | Quasiparticle calculations of the electronic properties of ZrO 2 and HfO 2 polymorphs and their interface with Si[END_REF] to 4.08 eV [START_REF] Ondračka | Accurate prediction of band gaps and optical properties of HfO 2[END_REF]. By using advanced semi-empirical functionals like TbmBJ [START_REF] Ondračka | Accurate prediction of band gaps and optical properties of HfO 2[END_REF], or many body perturbation theory in the GW approximation [START_REF] Grüning | Quasiparticle calculations of the electronic properties of ZrO 2 and HfO 2 polymorphs and their interface with Si[END_REF][START_REF] Jiang | Electronic band structure of zirconia and hafnia polymorphs from the gw perspective[END_REF][START_REF] Van Der Geest | Ab initio study of the electrostatic dipole modulation due to cation substitution in HfO 2 /SiO 2 interfaces[END_REF], some recent works have reconciled a good agreement with experimental data, giving band gap values in the range 5.7-5.9 eV. An important point to re-mark is that these studies do not account for excitonic effects. Nevertheless, they are consistent with optical SE and EELS gaps. This is the first clue that the good agreement between theory and experiments may be caused by some kind of error compensation. Thus there are two small inconsistencies in this general HfO 2 characterization framework: the difficulty in distinguishing optical and electronic band gap values in experiments, and the good agreement between experimental optical measurements with theoretical calculations which do not take into account excitonic effects. To clarify the situation in this work we compute the many body electronic structure of HfO 2 by using the GW approximation [START_REF] Hedin | New method for calculating the one-particle Green's function with application to the electron-gas problem[END_REF][START_REF] Strinati | Dynamical correlation effects on the quasiparticle bloch states of a covalent crystal[END_REF][START_REF] Strinati | Dynamical aspects of correlation corrections in a covalent crystal[END_REF], applying self-consistency on wavefunctions (QSGW ) [START_REF] Faleev | All-electron selfconsistent GW approximation: Application to Si, MnO, and NiO[END_REF], and pushing to the limit the convergence of the calculation. The core electrons are also included in the calculation by using the projector-augmented wave (PAW) method. On top of QSGW we add a Bethe-Salpeter equation calculation (BSE) [START_REF] Salpeter | A relativistic equation for bound-state problems[END_REF][START_REF] Hanke | Dielectric response in the Wannier representation: Application to the optical spectrum of diamond[END_REF][START_REF] Hanke | Many-particle effects in the optical excitations of a semiconductor[END_REF] of the optical band gap, which correctly takes into account electron-hole interaction, that is excitonic effects. Ab initio calculations based on density functional theory [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF][START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF] in the LDA or PBE [START_REF] Perdew | Generalized gradient approximation made simple[END_REF] approximations are carried out using the VASP code [START_REF] Kresse | Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set[END_REF][START_REF] Kresse | From ultrasoft pseudopotentials to the projector augmented-wave method[END_REF]. The core-valence interaction are described with PAW datasets including the semicore 5s and 5p states for Hf. Electron wavefunctions are expanded in a plane waves basis set with kinetic energy cutoff of 500 eV, and the Brillouin zone is sampled using 4 × 4 × 4, 6 × 6 × 4 and 6 × 6 × 6 Γ-centered Monkhorst-Pack meshes for monoclinic, tetragonal and cubic HfO 2 phases, respectively. We study cubic (spacegroup F m3m), tetragonal (spacegroup P 4 2 /nmc) and monoclinic (spacegroup P 2 1 /c) phases of HfO 2 . Each structure is relaxed until the maximum residual forces are less than 10 -3 eV/ Å. For the relaxation we use a fine k-mesh of 12×12×12, 12×12×8 and 8×8×8 for cubic, tetragonal and monoclinic phases, respectively. Calculated structural parameters are summarized in table 4.1.

HfO

Convergence of QP and BSE calculations

Many body effects are accounted for by three different levels of quasi-particle (QP) calculations on top of DFT: G 0 W 0 , GW 0 , which perform self-consistency only on eigenvalues, and QSGW , which updates also the wavefunctions, but fixing W at W 0 . Indeed, fixing W has been shown to improve the agreement with experimental band gaps than using self-consistent W [START_REF] Shishkin | Self-consistent GW calculations for semiconductors and insulators[END_REF]. Excitonic effects are accounted for by solving the BSE equation on top of the GW QP band structure. The convergence of all this calculations is particularly cared in this work. We carefully examine the convergence of our GW calculations in order to achieve QP band gap values converged within 0.1 eV. To calculate the response function and the correlation part of self-energy, a summation over empty states is required and quasi particle (QP) energies exhibit a very slow convergence with respect to the number of virtual orbitals. According to ref. [START_REF] Klimeš | Predictive GW calculations using plane waves and pseudopotentials[END_REF], for a correct convergence study, the number of empty bands N b , the corresponding orbital basis-set N pw (controlled by the plane waves cutoff E pw ) and the size of the response function basis-set N χ pw (controlled by the plane waves cutoff E χ pw ) have to be increased simultaneously. Therefore we fix E χ pw = 2/3E pw and we use the complete plane waves basis set (i.e. N b = N pw ) for each considered E pw . We consider the monoclinic phase of HfO 2 with E pw ranging from 350 eV to 700 eV with a step of 50 eV. A frequency grid with 100 frequency points is used to represent the polarizability and a 4 × 4 × 4 Γ-centered k-mesh is used to sample the Brillouin zone. Figure 4.1 shows the QP corrections (∆ ) of the Kohn-Sham eigenvalues as a function of 1/N b (which equals 1/N pw ) for a G 0 W 0 calculation on top of PBE (a similar behavior is observed for an LDA starting point). The dotted lines show the linear extrapolation of QP corrections (∆ (N b ) = A/N b + ∆ (∞)) where only the values corresponding to E pw ≥ 500 eV are included in the fit. The extrapolated QP corrections of the valence band maximum (VBM) at Γ and conduction band minimum (CBM) at B, that is the ∆ (∞) extracted by the two fits, are respectively -1.33 eV and 0.51 eV, giving an extrapolated band gap of 5.85 eV. By using a 500 eV plane waves cutoff (3603 plane waves and E χ pw = 333 eV) the band gap values turns out to be around 5.77 eV, which is converged within 80 meV. For this basis-set, the number of empty bands can be decreased to 2000 without deteriorating the convergence (see figure 4.2). We also check the influence of the number of frequency points and k-mesh. Self-consistent GW calculations are carried out with the same parameters as for G 0 W 0 . For GW 0 and QSGW calculations, we include respectively 256 and 512 QP energies in the self-consistent procedure for monoclinic HfO 2 . In the case of GW 0 , 4 self-consistent iterations are enough to get a converged band gap. For QSGW, 8 self-consistent GW iterations allowed to get converged QP energies and band gaps within 5 meV as shown in Fig. 4.3.

BSE

The excitonic properties are determined by solving the Bethe-Salpeter equation (BSE) within the Tamm-Dancoff approximation [START_REF] Sander | Beyond the Tamm-Dancoff approximation for extended systems using exact diagonalization[END_REF] on top of the GW quasiparticle band structure. BSE calculations usually require fine k-point sampling to converge exciton spectra. However, such calculations are computationally very expensive for both GW and BSE. In the case of monoclinic HfO 2 we test the k-point convergence of the BSE calculation on top of G 0 W 0 using 4 × 4 × 4 and 6 × 6 × 6 grids. Fig. 4.4 shows the imaginary part of the dielectric function ε 2 (we ignore the polarization dependence and assume ε = (ε xx + ε yy + ε zz ) /3). A simple visual inspection suggests that the two k-meshes give very similar spectra. More quantitatively, the exciton binding energies of 4×4×4 and 6×6×6 meshes are 0.56 eV and 0.52 eV, respectively. A more complete convergence study with respect to k-mesh could be carried out using a model BSE (mBSE) scheme or interpolation techniques but have not been considered in this work.

Electronic and optical band gap results

In this section we analyze the results obtained in this work about the electronic and the optical band gap of HfO 2 , trying to clarify the two inconsistencies underlined in section 4.1, that is the confusion between experimental optical and electronic gap and the agreement between optical experimental measurements and theoretical simulations which does not account for exciton HfO 2 GW study energies. Between the DFT and the GW approximation, there is some rearrangement of the bands, but the most important effect is a shift of both valence and conduction bands which increases band gaps. The same behavior is observed for the cubic and the tetragonal phases. The shift is similar but smaller for the GW 0 , and especially the G 0 W 0 , method. We pass now to analyze the fundamental band gap values of HfO 2 . A first comment is worth about the effect of our convergence parameters with respect to previous studies. The comparison can be done on G 0 W 0 and GW 0 calculations, whereas it is difficult to find previous works about QSGW . ) from different GW calculations. The calculations differ by the level of selfconsistency, starting mean-field theory and GW convergence parameters. The convergence parameters include: the k-mesh, the number of empty bands (N b ) and the method to describe dynamical screening (Dyn. Scr.). ble 4.2 shows that the fundamental band gap of monoclinic HfO 2 turns out to be higher when computed by using the convergence parameters of our study. The GW 0 values around 6.2 ∼ 6.3 eV are about 0.4 ∼ 0.5 eV higher than previous theoretical and experimental studies in literature. This is a first clue that the true electronic band gap of HfO 2 could be higher than so far believed. Table 4.3 gives a more complete overview of the fundamental minimum band gaps for all phases, and for all the approximations considered, including also QSGW (some values are reported from table 4.2). At the DFT level both the LDA and PBE results of this work, which use PAW datasets, are more in agreement with the FP-LAPW (full potential linearized augmented plane wave) LDA and PBE gaps of, respectively, ref. [START_REF] Jiang | Electronic band structure of zirconia and hafnia polymorphs from the gw perspective[END_REF] and [START_REF] Ondračka | Accurate prediction of band gaps and optical properties of HfO 2[END_REF] (both are all-electron calculations). The agreement is instead less good with the normconserving pseudopotentials plane waves (NCPP PW) LDA gap of ref. [START_REF] Grüning | Quasiparticle calculations of the electronic properties of ZrO 2 and HfO 2 polymorphs and their interface with Si[END_REF]. At the QP level the G 0 W 0 band gaps of this work is systematically larger than in the other considered works (ref.s [START_REF] Grüning | Quasiparticle calculations of the electronic properties of ZrO 2 and HfO 2 polymorphs and their interface with Si[END_REF][START_REF] Jiang | Electronic band structure of zirconia and hafnia polymorphs from the gw perspective[END_REF]), but closer to the NCPP PW calculations performed in ref. [START_REF] Grüning | Quasiparticle calculations of the electronic properties of ZrO 2 and HfO 2 polymorphs and their interface with Si[END_REF]. The bare PAW

Electronic band gap

G 0 W 0 correction (∆E GW g = E GW g -E DFT g
) is in good agreement with the one of ref. [START_REF] Grüning | Quasiparticle calculations of the electronic properties of ZrO 2 and HfO 2 polymorphs and their interface with Si[END_REF] (see also table 4.2), and the final difference is explained by the different starting DFT gap. On the contrary, the starting gap of ref. [START_REF] Jiang | Electronic band structure of zirconia and hafnia polymorphs from the gw perspective[END_REF], computed by FP-LAPW, is closer to the one of this work, but the correction turns out to be much lower (see tab. 4.2), giving a final band gap 0. [START_REF] Grüning | Quasiparticle calculations of the electronic properties of ZrO 2 and HfO 2 polymorphs and their interface with Si[END_REF]149] 5.2 5.8 5.7 LDA+GW 0 [START_REF] Grüning | Quasiparticle calculations of the electronic properties of ZrO 2 and HfO 2 polymorphs and their interface with Si[END_REF] 5.5 6.0 5.9 DFT-LDA [START_REF] Jiang | Electronic band structure of zirconia and hafnia polymorphs from the gw perspective[END_REF] 3.55 4.36 3.95 LDA+G 0 W 0 [START_REF] Jiang | Electronic band structure of zirconia and hafnia polymorphs from the gw perspective[END_REF] 4.91 5.78 5.45 LDA+GW 0 [START_REF] Jiang | Electronic band structure of zirconia and hafnia polymorphs from the gw perspective[END_REF] 5.20 6.11 5.78 DFT-PBE [START_REF] Ondračka | Accurate prediction of band gaps and optical properties of HfO 2[END_REF] 3.77 4.79 4.08 TB-mBJ orig [START_REF] Ondračka | Accurate prediction of band gaps and optical properties of HfO 2[END_REF] 5.88 6.54 5.76 TB-mBJ imp [START_REF] Ondračka | Accurate prediction of band gaps and optical properties of HfO 2[END_REF] 6.17 6.81 6.01 TB-mBJ semi [START_REF] Ondračka | Accurate prediction of band gaps and optical properties of HfO 2[END_REF] 6.74 7.35 6.54 Experimental works UPS+IPS (straight line extrapolation) [START_REF] Bersch | Band offsets of ultrathin high-κ oxide films with Si[END_REF] 5.7 XPS+IPS (straight line extrapolation) [START_REF] Sayan | Band alignment issues related to HfO 2 /SiO 2 /p-Si gate stacks[END_REF] 5.86 XPS+IPS (comparison with shifted DFT) [START_REF] Sayan | Band alignment issues related to HfO 2 /SiO 2 /p-Si gate stacks[END_REF] 6.7

Table 4.3: Fundamental minimum band gap for the cubic (X → X), tetragonal (Z → Γ), and monoclinic (Γ → B) HfO 2 phases, as calculated in various approximations. Ref. [START_REF] Grüning | Quasiparticle calculations of the electronic properties of ZrO 2 and HfO 2 polymorphs and their interface with Si[END_REF] uses a norm-conserving pseudopotentials plane waves approach, whereas ref. [START_REF] Jiang | Electronic band structure of zirconia and hafnia polymorphs from the gw perspective[END_REF] and ref. [START_REF] Ondračka | Accurate prediction of band gaps and optical properties of HfO 2[END_REF] use FP-LAPW. The determination of the band gap from PES/IPS experiments relies on post-processing data analysis to remove tails due to impurities and fits: for the XPS+IPS experiment of ref. [START_REF] Sayan | Band alignment issues related to HfO 2 /SiO 2 /p-Si gate stacks[END_REF] we quote two estimates removing or not the effect of tails.

consistency on eigenvalues only, using the GW 0 approach, further increases the gap by ∼ 0.4 eV in our study and by 0.2 ∼ 0.3 eV in ref.s [START_REF] Grüning | Quasiparticle calculations of the electronic properties of ZrO 2 and HfO 2 polymorphs and their interface with Si[END_REF][START_REF] Jiang | Electronic band structure of zirconia and hafnia polymorphs from the gw perspective[END_REF]. Full self-consistency also on the wavefunctions [START_REF] Rangel | Transport properties of molecular junctions from many-body perturbation theory[END_REF], using the QSGW approach [START_REF] Faleev | All-electron selfconsistent GW approximation: Application to Si, MnO, and NiO[END_REF], removes any influence of the LDA or PBE starting point in our study, reducing the gap difference to less than 0.1 eV, a residual due to the use of

Electronic and optical band gap results

97

different relaxed LDA and PBE atomic structures. Table 4.3 also provides the results obtained in two different works [START_REF] Bersch | Band offsets of ultrathin high-κ oxide films with Si[END_REF][START_REF] Sayan | Band alignment issues related to HfO 2 /SiO 2 /p-Si gate stacks[END_REF] by photoemission experiments, which are planned to extract an estimate for the fundamental electronic band gap of solids. For the monoclinic phase, the GW 0 estimates of ref.s [START_REF] Jiang | Electronic band structure of zirconia and hafnia polymorphs from the gw perspective[END_REF][START_REF] Grüning | Quasiparticle calculations of the electronic properties of ZrO 2 and HfO 2 polymorphs and their interface with Si[END_REF][START_REF] Ondračka | Accurate prediction of band gaps and optical properties of HfO 2[END_REF] are of 5.9 eV, 5.78 eV and 5.76 eV respectively, which are in very good agreement with the 5.7 eV and 5.86 eV electronic band gap values determined by photoemission experiments. On the contrary, the results of this work seem to indicate the fundamental electronic gap to be higher than what generally accepted in literature. Indeed the PBE+QSGW calculation, which is the most advanced here considered, clearly provides a fundamental band gap of 6.33 eV. If we could find some arguments, mainly related to convergence and theoretical framework of our calculations, to claim that our result may be more trustworthy than previous theoretical works, how can it be consistent with photoemission experiments?

The answer may be that the determination of electronic band gaps from photoemission experiments requires a linear extrapolation of the photoemission band edges to the background intensity. As discussed in sections III.C and F of ref. [START_REF] Sayan | Band alignment issues related to HfO 2 /SiO 2 /p-Si gate stacks[END_REF], this procedure has some uncertainty due to the presence of band-tail and defects in the vicinity of the valence band maximum and conduction band minimum. For this reason, we prefer to directly compare the experimental PES+IPS spectra to our calculated DOS (fig. 4.6). The same ref. [START_REF] Sayan | Band alignment issues related to HfO 2 /SiO 2 /p-Si gate stacks[END_REF] suggests this comparison to provide a safer estimate of the real electronic band gap. By using a DFT-LDA DOS and by evaluating the scissor operator shift to make theoretical and experimental DOS coincide, they estimate a band gap of 6.7 eV for monoclinic HfO 2 [START_REF] Sayan | Band alignment issues related to HfO 2 /SiO 2 /p-Si gate stacks[END_REF] (last line of tab. 4.3). The QSGW DOS computed in this work favorably compares with photoemission spectra, especially on the shape, even though we have not taken into account extrinsic and finite state effects which are evident when comparing XPS with UPS shapes. As it can be estimated by the deviation of theory and experiment in the conduction edge (fig. 4.6), the QSGW band gap of 6.33 eV is still an underestimation of about 0.2 eV of the real band gap. GW 0 provides a similar, but slightly greater, underestimation, whereas G 0 W 0 strongly underestimates the band gap about more than 0.7 eV. Thus we prudently conclude that the real electronic band gap of monoclinic HfO 2 would be E g > 6.33 eV, and probably E g 6.5 eV. This is consistent with the result of Ondračka et al. [START_REF] Ondračka | Accurate prediction of band gaps and optical properties of HfO 2[END_REF] obtained by modifying the TB-mBJ functional ("semi" version in table 4.3) to target the experimental DOS. The QSGW approximation has been reported to systematically overestimate band gaps in all studied materials [START_REF] Bruneval | Quasiparticle Self-Consistent GW Method for the Spectral Properties of Complex Materials[END_REF]. In this work, for monoclinic HfO 2 , the close agreement between QSGW and experimental spectra may be due to fortuitous error cancellation with other effects not taken into account, such in the G 0 W 0 , GW 0 , and QSGW approximations compared to UPS+IPS [START_REF] Bersch | Band offsets of ultrathin high-κ oxide films with Si[END_REF] and XPS+IPS [START_REF] Sayan | Band alignment issues related to HfO 2 /SiO 2 /p-Si gate stacks[END_REF] spectra. The theoretical DOS has been interpolated using MLWF on a 40 × 40 × 40 k-point grid and convoluted with a gaussian broadening of 0.7 eV. Experimental spectra have been aligned at the valence band-edge (the zero of the energy for all theoretical DOS), and the conduction IPS and valence PS independent measurements have been rescaled separately to match the height of the theoretical DOS. The error in this procedure can be estimated by the maximum deviation among valence band-edges in the steepest rising linear region: 50 meV, which is less than our theoretical error bar. as electron-phonon [START_REF] Poncé | Verification of first-principles codes: Comparison of total energies, phonon frequencies, electron-phonon coupling and zeropoint motion correction to the gap between ABINIT and QE/Yambo[END_REF], and both single-particle (e.g. spin-orbit) or manybody (e.g. Breit interaction) relativistic corrections. It may also be that, after applying these corrections, our QSGW calculations would provide an overestimation of the electronic band gap, whereas previous works based on other approaches (e.g. G 0 W 0 or GW 0 of ref. [START_REF] Grüning | Quasiparticle calculations of the electronic properties of ZrO 2 and HfO 2 polymorphs and their interface with Si[END_REF] or [START_REF] Jiang | Electronic band structure of zirconia and hafnia polymorphs from the gw perspective[END_REF]) could provide a better agreement with the DOS spectra of figure 4.6. [START_REF] Nguyen | Sub-bandgap defect states in polycrystalline hafnium oxide and their suppression by admixture of silicon[END_REF] 5.6-5.8 SE [START_REF] Hill | The relationship between local order, long range order, and sub-band-gap defects in hafnium oxide and hafnium silicate films[END_REF] 5.5-6.0 Table 4.4: Direct electronic band gap (DFT or GW ) and optical gap (BSE first exciton eigenvalue energy) for the cubic (X → X), tetragonal (Γ → Γ), and monoclinic (B → B) HfO 2 phases. We report also the exciton binding energy equal to the difference between the QSGW direct band gap and the BSE optical gap,

Optical band gap

E exc b = E d g -E o g .
For monoclinic HfO 2 , BSE gap agrees well with SE (spectroscopic ellipsometry) optical onset. Energies are in eV.

surements. The optical gap in solids is typically measured by performing experiments of optical or X-ray absorption, spectroscopic ellipsometry (SE), or energy-loss (EELS). Ref.s [START_REF] Nguyen | Sub-bandgap defect states in polycrystalline hafnium oxide and their suppression by admixture of silicon[END_REF][START_REF] Hill | The relationship between local order, long range order, and sub-band-gap defects in hafnium oxide and hafnium silicate films[END_REF] report values derived from SE spectra in the ranges 5.6 -5.8 eV and 5.5 -6.0 eV, respectively. The uncertainty is due, as for the band edges in the density of states, to the method (e.g. Tauc-Lorentz) used to linearly extrapolate to the background. However,compared to photoemission, optical experiments are less affected by defects, surface, interface or substrate effects, and are more sensitive to the bulk. Contrary to what observed for the electronic gap, the QSGW +BSE optical gap computed in this work lies in the range of the experimental reports. Nevertheless, it is again worth [START_REF] Weissker | Dynamic structure factor and dielectric function of silicon for finite momentum transfer: Inelastic xray scattering experiments and ab initio calculations[END_REF] to compare the raw SE spectra to the calculated optical absorption ε 2 (ω). Fig. 4.7 shows that the QSGW +BSE dielectric function significantly improves the lower level of approximation QSGW +RPA (Ran- dom Phase Approximation), and achieves a very good agreement with SE spectra. BSE introduces electron-hole interaction effects and gives rise to the exciton, which is represented by the first peak of the BSE spectrum, absent in the RPA. However, we remark a 0.1 ∼ 0.2 eV red shift of the exciton peak compared to its position in SE spectra. Part of this red shift could be corrected by k-point sampling extrapolation to zero [START_REF] Fuchs | Efficient O(N 2 ) approach to solve the Bethe-Salpeter equation for excitonic bound states[END_REF]. Nevertheless, as for the electronic band gap, the simulation of the optical gap probably suffers an underestimation. Therefore they have both to be regarded as lower bounds of the real values. The nature of the small peak observed at the onset of SE optical spectra [START_REF] Nguyen | Sub-bandgap defect states in polycrystalline hafnium oxide and their suppression by admixture of silicon[END_REF][START_REF] Park | Optical properties of thermally annealed hafnium oxide and their correlation with structural change[END_REF][START_REF] Hill | The relationship between local order, long range order, and sub-band-gap defects in hafnium oxide and hafnium silicate films[END_REF] of crystalline samples has been attributed to different causes. By combining X-ray absorption (XAS), X-ray diffraction (XRD) and SE techniques, Hill et al. [START_REF] Hill | The relationship between local order, long range order, and sub-band-gap defects in hafnium oxide and hafnium silicate films[END_REF] conclude that this feature could be intrinsic to the monoclinic phase. According to the analysis in this work, it constitutes a real bound [START_REF] Sottile | Efficient ab initio calculations of bound and continuum excitons in the absorption spectra of semiconductors and insulators[END_REF] exciton peak, as also interpreted in ref.s [START_REF] Edwards | Status and prospects for VUV ellipsometry (applied to high K and low K materials)[END_REF][START_REF] Ondračka | Accurate prediction of band gaps and optical properties of HfO 2[END_REF], and not a defect state, as instead interpreted by Nguyen et al. [START_REF] Nguyen | Sub-bandgap defect states in polycrystalline hafnium oxide and their suppression by admixture of silicon[END_REF]. The exciton only appears in the y polarization (see Fig. in the dielectric properties of monoclinic HfO 2 [START_REF] Guedj | Evidence for anisotropic dielectric properties of monoclinic hafnia using valence electron energy-loss spectroscopy in high-resolution transmission electron microscopy and ab initio time-dependent densityfunctional theory[END_REF][START_REF] Hung | Interpretation of monoclinic hafnia valence electron energy-loss spectra by time-dependent density functional theory[END_REF][START_REF] Botti | Ab initio calculations of the anisotropic dielectric tensor of GaAs/AlAs superlattices[END_REF]. Our calculations show an exciton also in the cubic and tetragonal phases of HfO 2 , but their oscillator strength tends to zero. Hence they are dark excitons, not detectable in SE spectra. This property can provide an easy and safe criterion to characterize the HfO 2 monoclinic phase with respect to all other phases, by simply detecting the presence or absence of the 5.85 eV exciton peak in SE spectra.

Conclusions about HfO 2 properties

We conclude that the confusion between the electronic and optical gap in monoclinic HfO 2 is probably caused by two main reasons. On one hand, previous theoretical calculations probably underestimate the electronic band gap. This may be due to some corrections not taken into account (e.g. electronphonon or relativistic), and/or to a possibly imperfect choice of convergence parameters. On the other hand, there is an uncertainty in the extrapolation of the same electronic gap by photoemission measurements, which lead to similar underestimations. This is inferred by observing the good agreement between the photoemission spectra and the DOS computed by a well converged PBE+QSGW calculation, which indicates 6.33 eV as lower bound for the electronic gap. Therefore the electronic gap extrapolations from photoemission measurements, which lie well below 6 eV, should probably be reviewed to higher values. The electronic and optical gaps of monoclinic HfO 2 should thus be well distinct, separated by at least 0.4 ∼ 0.6 eV. We also revisit the strange agreement between experimental measurements of the optical gap of HfO 2 and theoretical simulations which do not take into account the creation of exciton energy levels due to electron-hole interaction. In this case, the experimental measurements are more trustworthy, whereas the theoretical result is the product of error compensations: an underestimation due to the same reasons as for the electronic gap compensates the overestimation produced by non-subtracting the exciton energy. A well converged PBE+QSGW calculation of the direct B → B gap of monoclinic HfO 2 gives a value of 6.41 eV, which decreases to 5.85 eV after applying the Bethe-Salpeter equation restoring the agreement with experiments. The exciton binding energy is thus equal to 0.56 eV, which is also consistent with SE measurements. A direct comparison between SE and BSE spectra reinforces the whole analysis.

Chapter 5 Many Body Landauer test on Ag/MgO/Ag

Before studying the conduction properties of HfO 2 based OxRAM-like structures by means of the more physically reliable "many body Landauer " method, the latter needs to be tested and validated. HfO 2 , which is the material of main interest in this work, is not easy nor fast to simulate (see chapter 4). Moreover the interfaces between HfO 2 and the most used metals as electrodes, that is for instance titanium (Ti), aluminum (Al), or Platinum (Pt), are often really disordered and difficult to relax to the global energy minimum. Therefore, for the testing phase, we chose simpler materials and structures.

Structure

We consider an Ag/MgO/Ag junction. In theory, it qualitatively possesses all the resistive characteristics of OxRAM HfO 2 based devices but is easier and faster to treat. Indeed MgO has a much simpler electronic structure than HfO 2 , without d or f orbitals, and the Ag/MgO interface is very simple. Both Ag and MgO belong to space group Fm-3m (225), which corresponds to a face centered cubic (FCC) Bravais lattice. The atomic positions in the cubic cells are visualized in figure 5.1. Moreover the experimental cubic lattice constants at room temperature are close: 4.212 Å for MgO [START_REF] Wyckoff | Crystal Structure[END_REF], and 4.079 Å for Ag [START_REF] Davey | Precision measurements of the lattice constants of twelve common metals[END_REF]. These features lead to a clean and ordered interface, which is easy and fast to relax. To build the supercell firstly a cell volume relaxation is performed on bulk Ag and bulk MgO. The atomic forces are computed by DFT using a GGA functional (details in section 5.2). The lattice constants turn out to be slightly modified compared to the experimental values: 4.20 Å for MgO, and 4.14 Å for Ag. Then 4 different supercells are built by hand as a stack of relaxed cubic cells in the [START_REF] Onida | Electronic excitations: densityfunctional versus many-body Green's-function approaches[END_REF] direction, which is chosen as transport direction. The supercells are made by 4 Ag cubic cells at each side of n MgO cubic cells, with n ranging from 1 to 4. Hence the 4 supercells contain 40, 48, 56 and 64 atoms respectively. At the interface, O atoms are positioned in front of Ag atoms, which is known to be the most favorable configuration [START_REF] Schintke | Insulator at the ultrathin limit: MgO on Ag(001)[END_REF], at a guessed distance around 2 Å, which is coherent with the lattice constants. The lattice constants along directions [010] and [001], parallel to the interfaces, are fixed to the lattice constant of bulk Ag (4.14 Å) in order to mimic thin layers of MgO epitaxially sandwiched between thick unstrained Ag layers. Hence the MgO layer is compressively strained in the in-plane directions. The length along direction [START_REF] Onida | Electronic excitations: densityfunctional versus many-body Green's-function approaches[END_REF] of these rough supercells built by hand is optimized by performing relaxations of the atomic positions with different supercell lengths, and by doing a parabolic fit of the total energy vs the length. This method permits to avoid Pulay stress, which can affect a direct optimization within a single DFT calculation. The last relaxation of atomic positions is performed with the supercell length fixed at its optimal value. The final supercell length turns out to be extended compared to the starting point. The interfacial distance becomes d 2.69 Å between the Ag and the MgO atomic planes, in good agreement with previous studies [START_REF] Schintke | Insulator at the ultrathin limit: MgO on Ag(001)[END_REF], while the MgO lattice constant turns out to be about 4.33 Å perpendicularly to the interfaces.

At this point, the Ag/MgO interface is well constructed and relaxed, but we have to ensure the perfect periodicity of the first 3 and last 3 Ag unit cells. At a larger distance than 3 Ag cubic cells, the couplings between Wannier orbitals can be neglected, since they have been checked to have no impact on the interpolated band structure and the transmission of the pristine electrode. Thus 3 Ag periodic cells are suitable to define the Hamiltonian of the semi-infinite electrodes. To ensure the periodicity the first and last 3 Ag cubic cells of the relaxed supercells are substituted with 3 Ag cells with the bulk lattice constant of 4.14 Å. The central part of the relaxed supercells (1Ag/nMgO/1Ag), which represents the metal/oxide interface, is kept fixed. This method is justified by the fact that, after enough distance from MgO, Ag should recover its perfect bulk structure. Due to the very similar structure Ag is not very much modified by MgO, and one single Ag cell is sufficient to connect the interface to the bulk part. Moreover, since we are limited by computational resources, the length of the electrode, which should be as long as possible from a formal point of view, is always the result of a compromise. The final supercells contain 4 Ag cubic cells, n MgO cubic cells, and 4 Ag cubic cells (n from 1 to 4), that is 4Ag/nMgO/4Ag (see figure 5

Computational details

The transmission function through the 4Ag/nMgO/4Ag supercells is computed on top of DFT, to test the reliability of the "many body Landauer" method comparing it to state of the art calculations, and on top of HSE, G 0 W 0 and self-consistent COHSEX (simply indicated as COHSEX from now on), to study the impact of different QP calculations. Ab initio calculations based on DFT are carried out using the VASP code [START_REF] Kresse | Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set[END_REF][START_REF] Kresse | From ultrasoft pseudopotentials to the projector augmented-wave method[END_REF] with the PBE [START_REF] Perdew | Generalized gradient approximation made simple[END_REF] functional. The core-valence interaction is described by Projector Augmented Waves (PAW) datasets including 4d and 5s states for Ag, 2s and 2p for O, and 3s for Mg. Electron wavefunctions are expanded in a plane-waves basis set with kinetic energy cut-off of 415 eV. The k-points Monkhorst pack (MP) grid to sample the Brillouin zone is chosen equal to 1 × 4 × 4. Calculations are also performed by an increased 2 × 4 × 4 grid for accuracy reasons (details in section 5.3). The chosen state of the art code for comparison purpose is Transiesta [START_REF] Papior | Improvements on non-equilibrium and transport green function techniques: The next-generation transiesta[END_REF]. The calculations by Transiesta are carried out on the same 4Ag/nMgO/4Ag supercells. As well as VASP the functional is PBE, while the pseudopotential is of the type norm-conserving Troullier-Martins [START_REF] Troullier | Efficient pseudopotentials for planewave calculations[END_REF], including the same states as the VASP PAW datasets. The pseudo-atomic-orbitals which compose the basis-set are defined by a double zeta polarized size, and a cut-off radius obtained by an EnergyShift parameter of 25 meV. The plane-waves cut-off for the real space grid is taken equal to 450 Ry, while the Brillouin zone is sampled using a 1 × 6 × 6 MP grid. HSE calculations are carried out using VASP with the HSE06 hybrid functional [START_REF] Krukau | Influence of the exchange screening parameter on the performance of screened hybrid functionals[END_REF]. According to Viñes et al. [START_REF] Viñes | Systematic study of the effect of hse functional internal parameters on the electronic structure and band gap of a representative set of metal oxides[END_REF] HSE06 requires an α mixing parameter higher than 0.4 to reproduce the MgO experimental gap. Thus α is fixed to 0.43, which gives an electronic band gap at Γ of 7.77 eV for the MgO primitive cell, in agreement with the experimental value which lies between 7.67 and 7.83 eV [START_REF] Viñes | Systematic study of the effect of hse functional internal parameters on the electronic structure and band gap of a representative set of metal oxides[END_REF]. The screening parameter is fixed to 0.2 Å-1 to select the HSE06 functional. For G 0 W 0 calculations on top of DFT, we use a cut-off energy of 210 eV for the response function and 64 frequency grid points. The number of empty bands used in the dielectric function and for the Green's function is 6144 for the structures with 1 and 2 MgO cubic cells, and 6656 for the structures with 3 and 4 MgO cubic cells. We remind that the one shot G 0 W 0 method only updates eigenvalues. To save computational resources only the first 768 bands are actually updated, which are sufficient for the following Wannierisation procedure, while all the other bands are only used in the linear combination to converge the screening W 0 (eq. 3.18). Eigenfunctions are instead kept unmodified to the Kohn-Sham (KS) ones. For COHSEX the number of bands is reduced to 2048, except for the biggest structure (4 MgO), which uses 3072 bands. The number of self-consistent cycles is taken equal to 3, which is a good compromise between the great computational cost and precision. Differently from G 0 W 0 , all energy bands are updated at each self-consistent cycle, to be coherent when updating the eigenfunctions. We also give some computational details about the basis-set transformation from Bloch states to maximally localized Wannier functions (WF), carried out by Wannier90 [START_REF] Mostofi | An updated version of Wannier90: A tool for obtaining maximally-localised wannier functions[END_REF], which implements the Marzari and Vanderbilt method [START_REF] Souza | Maximally localized Wannier functions for entangled energy bands[END_REF], and its interface with VASP. The frozen energy window is chosen between the bottom of the valence band up to about 3 eV above the Fermi level. This value always guarantees a good description of several conduction states of MgO. For Ag, we choose the set of trial orbitals proposed by Souza et al. [START_REF] Souza | Maximally localized Wannier functions for entangled energy bands[END_REF] for copper: 5 d orbitals centered on each atom, and interstitial s orbitals. For MgO, we choose 1 s and 3 p orbitals on oxygen atoms, and also interstitial s orbitals. Hence there are 28 Wannier functions per Ag cubic cell and 24 Wannier functions per MgO cubic cell, which yields, for instance, a total of 320 Wannier functions for the largest structure (4Ag/4MgO/4Ag). The number N b of ab initio bands included in the optimization of the Wannier subspace is chosen at least twice the number N w of Wannier functions, which is always smaller than the 768 bands corrected by G 0 W 0 . Using 100 iterations for optimizing the subspace ("disentanglement" step) and 100 iterations for optimizing the Wannier functions, we obtain well-localized functions with quadratic spreads not exceeding 2 Å2 . Finally, for the transmission function calculation, the 4 × 4 ab initio grid of k-points transverse to the transport direction is refined to 10 × 10 by WF interpolation 3.53. This grid is sufficient to reach convergence.

Wannier Hamiltonian accuracy: k-points and cut-off

For test purposes, the preliminary calculations were carried out at the DFT level on the four supercells considered. Their lengths range from 38.5 Å for 4Ag/1MgO/4Ag to 51.5 Å for 4Ag/4MgO/4Ag along the transport direction [START_REF] Onida | Electronic excitations: densityfunctional versus many-body Green's-function approaches[END_REF]. That given, for plane-waves DFT calculations one single k-point sampling along the transport direction is largely sufficient to obtain a converged result. Thus calculations were initially carried out by using a 1 × 4 × 4 ab initio MP grid. However, the final transmission functions were very often spiky, unstable and unreliable (fig. 5.3, top plain black line), especially for structures with thick oxides (3 and 4 MgO), which yield low conductance. The debugging, understanding, and solution of this problem was quite challenging and took a very long period of the second year of PhD. Firstly we discovered that the introduction of a cut-off (eq. 3.54) to discard the long range non-diagonal term of the Wannier tight-binding-like Hamiltonian could significantly increase the stability. The more the cut-off applied was short, the more the transmission function became lower and smoother until it stabilized (fig. 5.3, the two dashed black lines). However, the result was cut-off dependent, there was not a clear criterion neither to determine a correct cut-off value, nor to determine if the discarded matrix elements were physical or not, and thus there was no clear physical justification to discard them. Then we found that increasing the k-point sampling along the trans- Average transmission (E)

2 1 0 1 2 3 E E f (eV)
1 k-point, no cut-off 1 k-point, L cut = 20 Å 1 k-point, L cut = 13 Å 2 k-points, no cut-off 2 k-points, L cut = 13 Å Figure 5
.3: Average transmission function (eq. 3.66) for the 4Ag/4MgO/4Ag supercell, at the DFT level for different k-points and cut-offs. If a single k-point is used the spectrum strongly depends on the cut-off distance L cut (black lines). If 2 k-points are used, the spectrum becomes nearly independent on L cut down to 13 Å (red lines). Moreover, the spectra with 1 or 2 k-points are nearly identical at L cut = 13 Å. Note that the highest black line corresponds to a "natural cut-off" of 26 Å (half of the ab initio supercell, which is about 52 Å). By Hamiltonian construction, more distant couplings are not accessible when using 1 single k-point (see also fig. 5.6). port direction from 1 to 2 points, which may appear redundant, could solve the instability without cutting any matrix element, matching the previous 1 k-point stabilized cut-off results (fig. 5.3, red "dots-lines" line). This fact opened the way to the interpretation of the problem. On one hand, 2 k-points calculations are more trustworthy and reliable: as we will see in the next section they agree with other codes results, and they are independent on the choice of a cut-off distance (fig. 5.3, red dashed lines), solving the ambiguity of choosing which matrix elements to discard. This means that the physics mainly lies very close to the Hamiltonian diagonal, and the non-diagonal terms tend rapidly to zero when increasing the distance (see fig. 5.5 and fig. 5.12 of section 5.5.1). On the other hand, 1 k-point calculations can reproduce the same spectrum by applying a given cut-off to discard some matrix elements. Hence these discarded elements must contain a negligible physical contribution and a majority of non-physical noise. The questions now become: Which is the source of this non-physical noise in the non-diagonal Hamiltonian matrix elements? And how using 2 k-points in the DFT calculation removes the noise? The answer to these questions has to be researched in the construction of the tight-binding-like Wannier Hamiltonian. Since WF are localized, but also periodic, they do not tend to zero, but instead raise again at half distance from its replicas, that is at the border of its Wigner-Seitz (WS) cell (fig. 5.4 b)), where they assume a minimum value. Sampling the transport direction by 1 single k-point entails that the reference ab initio supercell (RS) and the extended Wannier supercell (ES) coincide along this direction. Therefore, as visible in figure 5.4, the WF replicas are in the first RS neighbor cell, at a distance L from the original WF (L being the RS length along transport). The border of the WS cell is thus at distance L/2. Using 2 k-points doubles these distances, L becoming the distance of the WS cell border and 2L the one of the replicas. Figure 5.4 b) also shows that, if the WF shape at the center of the WS cell does not change when using 1 or 2 k-points, the decay Many Body Landauer test on Ag/MgO/Ag becomes different with distance from the center. The WF corresponding to 1 k-point, close to the border of its WS cell, feels its replica and, as consequence, decreases slower than the WF corresponding to 2 k-points, which assumes lower values. The Hamiltonian coupling matrix element between two WF Rn H R n ES is defined as an integral over the ES. Assuming for instance Rn to be the WF of figure 5.4 b), in the 1 k-point case the calculation of the coupling with another WF R n would require to integrate the black plain line over the blue area, or equivalently over the WS cell, in fig. 5.4 b). In the 2 k-points case, the calculation would instead require to integrate the red dashed line over the whole x-axis range, which is equivalent to the ES in the 2 k-point case. Clearly, the coupling will turn out to be higher when using 1 k-point than 2 k-point due to the higher values of the WF tails. This overestimation produces the noise in the transmission spectra. The WF basis-set computed by using 1 k-point is more sensitive to the periodicity of the ab initio simulation, and thus suitable and physical to describe the periodic system. However, our purpose is to break the periodicity to simulate an infinite system, and the use of 2 k-points helps in approaching this limit, whereas the use of 1 k-point produces a Hamiltonian matrix which is somehow less physical to describe the infinite system. The overestimation of the matrix elements causing the noisy spectra mainly concerns couplings between distant WF. Indeed the couplings between close WF are always rather strong, and the overestimated tails turn out to be negligible in the integration to compute the corresponding matrix element, which is close to the diagonal. On the contrary couplings between distant WF, corresponding to far non-diagonal elements, are supposed to be weak, and the overestimated tails become dominant in the integration, which in turn becomes overall overestimated. This analysis is reinforced by observing that the transmission function usually stabilizes when the cut-off distance goes under the MgO thickness. This corresponds to discard the Ag-Ag tunneling couplings through MgO, which are supposed to be very weak and tend to zero, especially for thick oxides, but instead assume non-negligible overestimated values. Figure 5.5 confirms that the Ag-Ag tunneling couplings through MgO are much higher when computed by using 1 k-point and thus must be discarded. The matrix elements computed by using 2 k-points are instead much lower. The fact to keep or discard them becomes irrelevant, and the final shape of the transmission function turns out to be cut-off independent (fig. 5.3, red lines). By construction, the Hamiltonian has a "natural cut-off" beyond which the supercell as a function of distance along the transport direction, defined by the left-hand side of eq. 3.54. Before a distance of approximately 17 Å, which approaches the MgO thickness, the Ag-Ag couplings are direct, strong and decay exponentially with distance. For these couplings, there is no relevant difference between 1 and 2 k-points cases. Beyond 17 Å only Ag-Ag tunneling couplings through MgO are possible. In the 2 k-points case these couplings go to negligible values (under 10 -4 eV), in the 1 k-point case they are 1, 2 order of magnitude higher because of the overestimated WF tails values (see also fig. 5.4 b)). Note that, due to the Hamiltonian construction, for 1 k-point the maximum coupling distance is L/2 (the ES length is L). For 2 k-points the maximum coupling distance is doubled up to L (the ES length is 2L).

couplings cannot be defined. This cut-off is equal to L/2 • N 1 , where N 1 is the number of k-points along the transport direction. The reason is that the couplings are meant at the closest distance by translation 3.37 (see also fig. 3.7). Therefore, by using 1 k-point, the couplings are only defined over half supercell, while by 2 k-points the couplings can be evaluated over the whole supercell (see fig. 5.6, left and center). However, including or not these additional long distance couplings has a negligible influence on the transmission spectrum (fig. 5.6, right). The accuracy of the obtained Hamiltonian can be checked by interpolating the band structure of the ab initio supercell on a fine grid of k-points along the transport direction, and by comparing it to the band structure computed b): Same plot with the energy range restricted to the band close to E F . The interpolation from 1 k-point strongly overestimates the band width. Note that k 1 = 0.5 belongs to the 2 k-points ab initio MP grid, but not to the 1 k-point grid. Thus the band structure at this point is replicated exactly by the 2 k-points Wannier basis, but is interpolated when using the 1 k-point Wannier basis.

This confirms the fact that some direct couplings between Ag electrodes are overestimated in the transport simulations. Conversely, the better accuracy obtained with the 2 × 4 × 4 interpolation of the band structure gives confidence that the transport properties will be accurately calculated, since the Hamiltonian describes well the coupling of Ag states through MgO. The main conclusion of this section is that, when the transport direction is sampled by 1 single k-point, Wannier replicas inherent to the ab initio periodic supercell lead to an overestimate of some distant couplings. This error becomes dominant in the transmission spectrum when the tunneling current is very small, that is for thick oxides. This error is nearly completely removed when using 2 k-points along the transport direction. Residual spurious spikes in the transmission spectra can be removed by imposing a cut-off distance L cut .

Comparison with Transiesta

In this section, we verify the agreement of our original method with Transiesta, which is among the most widely used codes for electronic transport calculations. At the DFT level, our approach is based on a basis-set transformation from Kohn-Sham Bloch states to Wannier functions, while Transiesta takes advantage of DFT calculations directly performed on a localized orbitals basis-set. In the end, both methods compute transmission function and conductance by means of the Landauer-Green function approach. De- Considering the analysis of section 5.3 about k-points and cut-off, from now on all the results shown in this chapter are computed by using 2 k-points along the transport direction. However, note that the transmission function of 4Ag/1MgO/4Ag and 4Ag/2MgO/4Ag at the DFT level would be acceptable even if computed by using 1 single k-point. The reason is that the high transmission through the thin MgO layers, which is also systematically overestimated by DFT, hides the noise produced by the overestimated WF tails. Nevertheless, when passing to the QP level, the transmission function is lowered (see section 5.5), and also 1 and 2 MgO structures have to be treated by 2 k-points.

To compare the spectra to Transiesta we firstly studied the impact of the refinement of the interpolated k-points grid transverse to the transport direction. Figure 5.8 shows the transmission function of the 4Ag/3MgO/4Ag structure as computed by different k-points interpolated grids using our plane-waves+Wannier method. The analysis shows that the ab initio grid is largely insufficient to compute a converged spectrum for the transmission function and that a 10 × 10 grid is at least required. The analysis is similar and remains valid for the other structure with different MgO thickness. Remark that the interpolation and the final transmission function calculation has a very low computational cost compared to DFT or QP calculations, or to the basis-set transformation. Therefore there is no reason to not take advantage of a highly converged grid. Figure 5.9 shows the average transmission function (eq. 3.66) for the four supercells considered, computed by our plane-waves+Wannier method at DFT-PBE level and by Transiesta. The first observation is that the agreement is very good, especially for the two structures with thinnest oxides. The spectra of the two structures with thickest oxides appear a little higher than the Transiesta ones. There are many differences in the implementation which may explain this small discrepancy. Among them, the most important probably concerns the basis-set. We use from 248 (for 1MgO) to 320 (for 4MgO) WF to exactly reproduce the wavefunctions computed on the plane-waves basis-set. The latter is highly converged and rather complete, counting around 10 5 plane-waves. Siesta instead uses a few thousands of orbitals, which are directly optimized and converged at the DFT level. For localized functions DFT codes there is no clear criterion to perform a convergence study of the basis-set, like by using an energy cut-off (eq. 2.41) in plane-waves codes. Thus the Siesta basis-set is certainly less complete and converged than the plane-waves one. That given, on one hand, if the basisset transformation from plane-waves to WF is successful, the WF basis-set should be more complete than the Siesta one. On the other hand, the transformation is not uniquely defined and rather delicate, possibly turning out to be imperfect. As already pointed out, WF may be sensible to their periodicity which may produce an overestimation of their tails and some Hamiltonian matrix elements. Therefore this could also produce an overestimation of the transmission when the latter is small. For high transmission functions, these issues are instead hidden at lower orders of magnitude. The tails of the Siesta orbitals are instead put to zero beyond a given cut-off radius, to avoid any numerical problem. The difference between our plane-waves+Wannier method and Transiesta is anyway very small, and of second order compared to QP corrections, as we will point out in the next section. We can now start to give an interpretation of the transmission curves of figure 5.9. Below around -2.2 eV, the spectra have a maximum of conductance and the shape is independent on the MgO thickness. This region corresponds to the MgO valence band, where both MgO and Ag states are present, and the density of states is maximum. At -2.2 eV there is a drastic drop of con-ductance for all thicknesses which corresponds to the valence band edge of MgO. Beyond this edge, in the MgO gap, the transport enters the tunneling regime and the transmission decreases exponentially with the thickness of MgO. This is confirmed by looking at the zero-bias conductance values, which correspond to τ (E f ). The values are given in table 5.1 and their exponential decrease with the MgO thickness is shown in figure 5.14 (top blue line). This behavior is similar to what is obtained by Butler et al. in DFT+Landauer simulations of Fe/MgO/Fe magnetic tunnel junctions with parallel magnetizations [START_REF] Butler | Spin-dependent tunneling conductance of Fe/MgO/Fe sandwiches[END_REF]. The decay rate of the conductance with the MgO thickness found by Butler et al. is consistent with the one extracetd by figure 5.9 (see also fig. 5.14 and table 5.1). This decay is associated to evanescent states in the gap of MgO (also discussed in appendix B), which have ∆ 1 symmetry at k = 0 [START_REF] Butler | Spin-dependent tunneling conductance of Fe/MgO/Fe sandwiches[END_REF]. Indeed, the analysis of the k -resolved transmission T ( k , E f ), shown in figure 5.10, proves that the main contribution comes from transverse wavevectors close to zero. Observing figure 5.9 it is not clear at which energy MgO conduction states start again to contribute to transport. We could expect to observe a behavior similar to the valence edge, with a drastic rise of conduction independent on thickness. However, contrary to the valence edge, which is thickness independent, the conduction edge seems to be dependent on the number of MgO layers. We could also think that the edge of the conduction band lies where the smooth part of the spectrum reaches a maximum of around 1.7 eV and then starts to decrease irregularly. Nevertheless, the gap value computed by DFT-PBE on bulk MgO is equal to 4.72 eV (see table 5.2). Thus the MgO conduction edge should lie around 2.5 eV, in the middle of the irregular decrease. This discrepancy is explained by symmetry considerations. As the evanescent states, the first conduction band of MgO has ∆ 1 symmetry. However, at E > E f + 1.6 eV, Ag observes a transition to a band with a different symmetry, which cannot be transmitted through the ∆ 1 states of MgO. This fact prevents the transmission spectrum to raise again to high values at the MgO conduction edge. This symmetry analysis is detailed in appendix C.

4 3 2 1 0 1 2 3 E E f (eV)

Influence of QP calculations on the transmission functions

Before showing and interpreting the transmission functions computed at different levels of QP calculation, we analyze the stability of the transmission calculation as done at the DFT level in section 5.3.

Stability of the transmission calculation by different QP corrections

As already said at this point all calculations are performed by using a 2×4×4 ab initio MP grid, which guarantees smooth and cut-off independent transmission spectrum at the DFT level. QP calculations are expected to lower the spectra of 1 or 2 order of magnitude, which is the typical gap between DFT based transport calculations and experiments. Therefore some numerical instability coming from the overestimation of the WF tails could in principle become comparable with the transmission spectrum values, introducing some non-physical noise. For HSE and COHSEX this does not occur for the first three structures (from 1 to 3 MgO cubic cells), whose spectrum remains smooth and cut-off independent even when corrected by HSE or COHSEX. Note that all these transmission functions never assume values lower than the 4 MgO DFT one, whose minimum is around 3 This scenario completely changes when considering G 0 W 0 . Unexpectedly, figure 5.11 (bottom left) shows that already for the 2 MgO structure, at values between 10 -3 and 10 -2 , the spectrum turns out to be extremely noisy. The clue for interpreting this unexpected behavior is the observation that G 0 W 0 is the only considered QP calculation which only updates eigenvalues, keeping the DFT Kohn-Sham wavefunctions unmodified. At first sight, this observation can be more confusing and surprising, since if wavefunctions are the same as DFT the spectrum should also look similar. However, the modification of the eigenvalues by the G 0 W 0 correction is equivalent to a modification of the Hamiltonian which is local in energy, hence non-local in space. Thus it induces artificial high long-range couplings between Wannier functions. This was already noticed by Ferretti et al. [START_REF] Ferretti | Ab initio complex band structure of conjugated polymers: Effects of hydrid density functional theory and GW schemes[END_REF] who studied the impact of G 0 W 0 corrections in conjugated polymers. This hypothesis is supported by figure 5.12, which shows how the maximum modulus of the Wannier Hamiltonian matrix elements decreases with the distance between Wannier centers along the transport direction, defined by eV for all approximations except G 0 W 0 , which saturates before the others, around 13 Å, to higher values of around 10 -2 eV. This effect is also proven in appendix A for a simple one-dimensional, two levels, tight-binding model. Since the long-range couplings should always rapidly go to zero, and rapidly become negligible in the transmission spectrum calculation, we conclude that the important noise observed in the G 0 W 0 spectra (e.g. figure 5.11 bottom left) is due to the artificially raised long-range couplings produced when updating only the eigenvalues, keeping the wavefunctions unchanged. Thus, since the relevant physics for transport phenomena mainly lies close to the Wannier Hamiltonian diagonal, we decide to apply a cut-off L cut = 13 Å to discard the non-physical long-range couplings. The choice of this cut-off value is motivated by the fact that for both the graphs of figure 5.12, and thus for two different structures, G 0 W 0 starts to separate from the other approaches around 13 Å. The G 0 W 0 cut-off transmission spectra are smooth and regular (orange line of figure 5.11 bottom left), but subject to caution. From this analysis we conclude that probably G 0 W 0 should not be considered as a reliable method to compute the transport properties of tunnel junctions. This conclusion will be reinforced in the next section. We expect to observe a reduction of the transmission spectra compared to the DFT ones. This is confirmed for COHSEX and HSE, which for the 3 and 4 MgO structures reduce the T (E f ) value of more than 1 order of magnitude.

QP transmission functions analysis

On the contrary G 0 W 0 does not provide a significant reduction. For a more complete analysis, it is interesting to compare our results to the experimental results of Gangineni et al. [START_REF] Gangineni | Interfacial electronic transport phenomena in single crystalline Fe-MgO-Fe thin barrier junctions[END_REF]. In this work, measurements of the tunneling conductance have been performed on Fe/MgO/Fe junctions at different MgO thicknesses, up to 4 cubic cells. The experimental conductances for parallel magnetizations of the Fe electrodes have been then compared to the theoretical values obtained by Butler et al. [START_REF] Butler | Spin-dependent tunneling conductance of Fe/MgO/Fe sandwiches[END_REF] with the DFT+Landauer approach. For 4 MgO cubic cells, the deviation between experiment and theory turns out to be between 1 and 2 orders of magnitude (fig. 3 of ref. [START_REF] Gangineni | Interfacial electronic transport phenomena in single crystalline Fe-MgO-Fe thin barrier junctions[END_REF]). This is comparable to the deviation observed in our Ag/4MgO The values of zero-bias conductance (at zero temperature) are shown in table 5.1 for the four structures and four approximations. The values are converged with respect to the k-points interpolated grid, that is they are converged on a unit of lateral surface. Remember indeed that increasing the k-points grid corresponds to increase the number of replicas which compose the extended supercell, whose conductance is computed. The results of the simulations can thus be linearly reported to a given lateral surface, to allow the comparison with experimental samples of given size. A more reliable comparison with experiments should involve the conductance evaluated at room temperature. This evaluation can be done by applying eq. 3.67, weighing the transmission function by the -∂η f /∂E function. As discussed in chapter 2 the effect of this weighing function is that the transmission function gives 5.5. Influence of QP calculations on the transmission functions 123 a net contribute to the conductance only in a range E f ± 5k B θ, which is quite narrow at room temperature (5k B θ 0.13 eV ). Since all the transmission functions of figure 5.13 does not show rapid variations very close to the Fermi level, the conductance values at room temperature does not change significantly with respect to the zero temperature values. Therefore the values of table 5.1 are also a very good approximation for room temperature conductance values, as long as electron-phonon scattering can be neglected. Figure 5.14 shows that the decay with the MgO thickness is exponential, with a net separation between DFT/G 0 W 0 and HSE/COHSEX. For DFT the decay rate is very close to what obtained by Butler et al. [START_REF] Butler | Spin-dependent tunneling conductance of Fe/MgO/Fe sandwiches[END_REF]. We obtain G ∝ exp(-κ × n M gO ), where n M gO is the number of MgO cubic cells and κ 3.3. κ increases for HSE and COHSEX to values around 4.2. This entails that the overestimation error of the DFT+Landauer approach becomes more important in higher tunneling regimes. The decay rate obtained in HSE and COHSEX (κ 4.2) is also in good agreement with the rate obtained for the A sample in ref. [START_REF] Gangineni | Interfacial electronic transport phenomena in single crystalline Fe-MgO-Fe thin barrier junctions[END_REF]. On the contrary, the decay rate obtained in DFT (κ 3.3) is closer to the B sample (fig. 3 of ref. [START_REF] Gangineni | Interfacial electronic transport phenomena in single crystalline Fe-MgO-Fe thin barrier junctions[END_REF]). A comparison of the decay rates with the ones obtained from the complex band structure of bulk MgO is performed in appendix B.

The fact that G 0 W 0 does not provide a significant reduction of transmission and conductance (fig. 5.13 and 5.14) leads to another interesting analysis. The DFT systematical underestimation of the insulators band-gap [START_REF] Perdew | Density functional theory and the band gap problem[END_REF] is suspected to be the main cause of the overestimation of the tunneling transmission. This is inferred by assuming that the band offsets between the Fermi level E f of the electrodes and the valence and conduction band edges of MgO are both underestimated by DFT. This underestimation of the energy barriers should lead to an overestimation of tunneling. MB corrections on eigenvalues are thus expected to give more accurate conductance values. Table 5.2 shows the electronic band-gap of bulk MgO calculated in different approximations: DFT-PBE, G 0 W 0 , COHSEX, and HSE. G 0 W 0 gives a bandgap close to the experimental values, whereas COHSEX, as expected, overestimates the band-gap. Our observation is that the decrease of the tunneling transmission is not simply correlated to the opening of the MgO band-gap. In the G 0 W 0 approximation, the zero-bias conductances remain close to the DFT values, despite the band-gap opening from 4.72 to 7.44 eV. The HSE approximation gives similar gap opening but a smaller transmission. The COHSEX approximation gives conductances very close to the HSE ones, while the band-gap is much larger. The band offset between E f and the VB edge of MgO is also not simply correlated to the band-gap opening. In particular, the offset is larger in HSE than in COHSEX. The HSE offset is around 4 eV, which is close to the experimental value obtained by scanning tunneling spectroscopy on thin MgO layers deposited on an Ag substrate [START_REF] Schintke | Insulator at the ultrathin limit: MgO on Ag(001)[END_REF]. These results point to the fact that it is not the correction of the band-gap and energies which drives the reduction of the zero-bias conductance, but rather the correction of wavefunctions. Indeed, G 0 W 0 , which keeps the Kohn-Sham eigenfunctions unmodified, does not provide any significant correction to the DFT zero-bias conductance, while HSE and COHSEX, which update also eigenfunctions, provide significant and similar corrections to the zero-bias conductance, despite the different MgO band-gaps. Figure 5 around the Fermi energy is necessarily caused by the modification of eigenfunctions. This follows by the fact that the trace formula used to compute the transmission is local in energy. In other words, the transmission at energy E can be obtained by building diffusion eigenstates made of incident, reflected, and transmitted waves. These states are linear combinations of the ab initio Bloch eigenfunctions at energy E. Hence a modification of these wavefunctions is needed to modify the DFT transmission close to E f since eigenvalues remain basically unmodified. This point seems to have some analogies with the conclusions of Ferretti et al. [START_REF] Ferretti | Ab initio complex band structure of conjugated polymers: Effects of hydrid density functional theory and GW schemes[END_REF], who computed the decay of evanescent states in conjugated polymers under various DFT+MB approximations. In section 5.5.1 G 0 W 0 was pointed out as possibly not suitable as a method to compute transport properties in metal/oxide junctions. Indeed the updating of only eigenvalues produces artificially raised long range couplings, which must always be discarded by a cut-off to obtain a reliable result. Nevertheless, the use of such a cut-off is not completely justified. The analysis of Wannier Hamiltonian, Wannier function tails, and the number of k-points along the transport direction is not valid for G 0 W 0 . The spurious long range interaction is intrinsic to the method, and cannot be removed by using 2 k-points. Moreover, in this section, G 0 W 0 is proved to provide very small changes of the transmission spectrum compared to DFT, still far from the expected one order of magnitude lowering. The latter is indeed more related to the modification of wavefunctions than eigenvalues. Therefore we reinforce the conclusion that G 0 W 0 should not be considered as a reliable method for electronic transport calculations in metal/oxide bulk junctions.

Conclusion about the Many Body Landauer method

In this chapter we have successfully tested the many body Landauer method developed within this work. However, this method should be used with particular attention to some critical aspects.

The method proved to be very versatile, being able to compute transmission function and conductance on top of many different QP calculations. HSE and COHSEX revealed to provide corrections that go in the direction of a better agreement with experiments, whereas G 0 W 0 turned out to be not suitable for this kind of calculations. The reason is due to the interesting argument that the lowering of the transmission function is mainly driven by the modification of wavefunctions, which remain unmodified in G 0 W 0 . A full self-consistent GW calculation would be in principle more accurate than HSE or COH-SEX. Nevertheless, due to the long supercells needed for these calculations, the computational resources revealed to be a parameter to particularly care about when using the many body Landauer method.

The QP calculations on plane-waves, implemented in highly developed codes, are very stable, even if computationally demanding. The transmission function calculation by means of the Green function formalism on the Wannier functions Hamiltonian, carried out by our original code, proved also to be very stable, and actually rather fast. The critical part of the whole procedure turned out to be the basis-set transformation. Here a good choice of orbitals 5.6. Conclusion about the Many Body Landauer method 127 as initial Wannier functions guess, an appropriate choice of transformation parameters of the Wannier90 code (e.g. the number of optimization steps or the frozen energy window), and most meaningfully the use of 2 k-points along the transport direction for low tunneling regimes, revealed to be crucial to obtain a good and trustworthy result.

Remark anyway that here we are measuring a zero-bias conductance, that is an equilibrium property which can approximate a non-equilibrium property, like the current, in a narrow range of linear response to the applied voltage. In the tunneling regime, we obtain very small conductance values (< 10 -11 Ω -1 per cubic cell for the thickest structure). The developing of a full NEGF formalism, which also includes MB corrections, would allow the simulations of samples under bias, where the current is much higher. In this situation, closer to true experimental ones, all the issues about the Wannier Hamiltonian would probably be negligible. Hence these issues probably become important only in calculations under extreme conditions.

Chapter 6

Transport calculation on an OxRAM model: Pt/HfO 2 /Pt

After the tests and the interesting information about the robustness and the usage of the "many body Landauer " method, the latter can be applied to systems more similar to true OxRAM devices, to study the behavior of its pristine, high resistive (HRS), and low resistive (LRS) states.

In this chapter, we consider a junction made of platinum and hafnium oxide, which is proposed in some works as a possible structure of OxRAM device [START_REF] Cagli | Experimental and theoretical study of electrode effects in HfO 2 based RRAM[END_REF][START_REF] Cabout | Role of Ti and Pt electrodes on resistance switching variability of HfO 2 -based resistive random access memory[END_REF]. We build three different structures, corresponding to the pristine state, and two guesses of HRS and LRS. The transmission function and conductance of these three structures are then computed at the DFT and HSE levels. The latter was demonstrated to be the best compromise between accuracy, robustness, and computational cost in the tests of chapter 5 on the Ag/MgO junction. HSE turned out to be more suitable than G 0 W 0 to compute electronic transport properties, as well as self-consistent COHSEX, which however requires more computational resources.

Structure

The construction of the Pt/HfO 2 /Pt junction looks similar to what has been done for the Ag/MgO one. We consider the monoclinic phase of HfO 2 (see fig. 6.1), which is the most stable at low temperature (< 2000 K [START_REF] Wolten | Diffusionless phase transformations in zirconia and hafnia[END_REF]). The monoclinic phase is proved to constitute a stable interface with Pt, and is widely studied in literature to investigate the properties of the Pt/HfO 2 interface by first principles (for instance ref.s [START_REF] Gavrikov | Oxidation of the Pt/HfO 2 interface: The role of the oxygen chemical potential[END_REF][START_REF] Cho | Electronic structure of Pt/HfO 2 interface with oxygen vacancy[END_REF][START_REF] Zhu | Phase equilibria at Si-HfO 2 and Pt-HfO 2 interfaces from first principles thermodynamics[END_REF]). The monoclinic HfO 2 structure is relaxed by using the GGA functional. The resulting lattice parameters are 5.15 Å, 5.19 Å, and 5.33 Å, with an angle between the first and the third primitive vector of 99.71 • . Pt has an FCC structure, whose lattice parameter turns out to be 3.957 Å after a GGA relaxation. The interface is built taking a (111) Pt surface and a (001) HfO 2 surface, as also proposed in several works [START_REF] Gavrikov | Oxidation of the Pt/HfO 2 interface: The role of the oxygen chemical potential[END_REF][START_REF] Cho | Electronic structure of Pt/HfO 2 interface with oxygen vacancy[END_REF][START_REF] Zhu | Phase equilibria at Si-HfO 2 and Pt-HfO 2 interfaces from first principles thermodynamics[END_REF]. The additional difficulty compared to Ag/MgO is that the HfO 2 primitive cell is not symmetric along the transport direction. Therefore the two interfaces relax to very different ions positions, causing differences in the two electrodes. The two electrodes interact through PBC in the plane-waves calculation, breaking the perfect periodicity of bulk Pt, which is required to build the semi-infinite electrodes for the transport calculation. To avoid this effect we build an orthogonal supercell made of 3 Pt cubic cells, 2 HfO 2 monoclinic cells, and again 3 Pt cubic cells, whose lateral size is determined by the Pt lattice parameter. The whole supercell is relaxed keeping the atoms in the two external Pt cubic cells fixed, to guarantee the perfect periodicity of electrodes. The cubic cells closest to HfO 2 are relaxed together with HfO 2 to do the connection between the interface and the electrodes. Ideally, this buffer zone would work better with more than one single Pt cubic cell. However, the Pt cubic cell along the <111> direction contains 12 atoms and 3 mono-atomic planes, contrary to 4 atoms and 2 planes along the <100> direction. Thus we reduce as possible the number of cells to limit the supercell dimension and to save computational resources. The final structure of the pristine state is shown in figure 6.2. It contains 96 atoms and its length is about 52.9 Å. The final atomic forces are very low (< 0.05 eV/ Å), only the forces at the interface between bulk periodic Pt and buffer zone are higher, around 0.1 eV/ Å. To reduce this value we should add two Pt cubic layers (one per side), that is 24 additional atoms, to increase the buffer zone. Nevertheless, this would highly increase the computational cost.

To build the structures corresponding to LRS and the HRS states, we choose an idealized filament model. We identify a straight line of 4 oxygen atoms directly connecting the two electrodes. This straight line is visible in figure 6.3 (top) for the pristine state at the center of the supercell. The LRS is obtained by eliminating the 4 oxygen atoms of the straight line, obtaining a straight filament of oxygen vacancies. The structure is relaxed by using a GGA functional, keeping fixed the atomic positions of the periodic bulk the supercell, whereas when a filament is introduced PBC can produce an interaction with filament replicas. To avoid this interaction the lateral size Transport calculation on an OxRAM model: Pt/HfO 2 /Pt should be increased. However, the computational cost would hugely increase. Moreover, since the structure of the filament is unknown, an interaction between different regions of oxygen vacancies cannot be a priori excluded. For these reasons we keep the lateral size equal to the one of a single Pt cubic cell.

The supercell here considered is a rough model of OxRAM device: the oxide thickness is smaller than in true devices (1 nm versus 5 ∼ 10 nm), the oxygen atoms are not displaced but eliminated to create the filament, and the form of the filament is very simplified. However, a transport calculation on these structures is anyway interesting as a preliminary test to go towards more realistic structures, to check the behavior of the conductance with respect to these filament models, and to check the behavior of the many-body corrections on the pristine state, LRS, and HRS.

Computational details

The calculations of this chapter are mainly carried out by using the Quantum Espresso code (QE) [START_REF] Giannozzi | QUANTUM ESPRESSO: a modular and opensource software project for quantum simulations of materials[END_REF][START_REF] Giannozzi | Advanced capabilities for materials modelling with quantum ESPRESSO[END_REF], which has a more developed interface with Wannier90 than VASP, and a new HSE implementation based on a more efficient adaptively compressed formulation of the Fock exchange operator, which greatly reduces the computational cost without loss of accuracy [START_REF] Lin | Adaptively compressed exchange operator[END_REF]. QE is used to perform both PBE and HSE calculations. For comparison purposes, the transmission function of the pristine structure at the PBE level are also computed by VASP and TranSiesta.

For the QE/PBE calculations, we have expanded the wavefunctions on a plane-waves basis-set defined by a cut-off energy equal to 90 Ry. The pseudopotentials are of the norm-conserving type, with 12 valence electrons for Hf, 6 for O and 18 for Pt. The reciprocal space is sampled by a 4 × 4 × 1 kpoints grid, which in this case is actually sufficient to obtain a stable and converged transmission function at the PBE level. The VASP/PBE calculation on the pristine structure uses a plane-waves cut-off energy of 450 eV; PAW pseudo-potentials with 4 valence electrons for Hf, 6 for O, and 10 for Pt; and a 4×4×1 k-point grid. The same calculation has been carried out by TranSiesta using a polarized double-zeta basis-set of atomic orbitals, with a cut-off radius yielding an energy shift parameter of 40 meV. The planewaves cut-off to define the real space grid is taken equal to 500 Ry and the k-points grid is again 4 × 4 × 1.

For HSE calculations the attention to computational resources has been crucial. The k-points grid is reduced to a 2×2×1 sampling, used in combination with a Wannier Hamiltonian cut-off of 13 Å. To validate the use of the cut-off the result was tested with a single calculation performed by a 2×2×2 k-point grid on the pristine structure. The plane-waves energy cut-off is reduced to 80 Ry and the kinetic energy cutoff for the exact exchange operator is taken equal to 120 Ry. The fraction of exact HF exchange is α = 0.25, which gives a band gap of E g = 5.68 eV for monoclinic bulk HfO 2 . As pointed out in chapter 4 this E g value is consistent with some theoretical and experimental works, but not with our accurate highly converged QSGW calculations, which give E g = 6.33 eV (see chapter 4 and ref. [START_REF] Sklénard | Optical vs electronic gap of hafnia by ab initio bethe-salpeter equation[END_REF]). A value of α > 0.34 would be necessary to obtain E g = 6.33 eV. However, in this chapter, we choose to keep α = 0.25 to be consistent with the majority of the publications about this topic and to save computational resources, since the use of a higher α parameter tends to worsen the convergence. The screening parameter is fixed to 0.2 Å-1 , to select the HSE06 method.

It should be noted that the pseudo-potentials used by QE include more valence electrons than in the VASP PAW datasets. The 5 s and 5 p orbitals of Hf (8 atoms) and Pt (72 atoms) give 320 additional core energy bands. These core low energy bands create some instability in the optimization of the Wannier functions basis-set. Thus they were excluded in the basis-set transformation taking advantage of an option of the Wannier90 code. Specifically, they are excluded in the sum 3.41. The frozen window goes from the bottom of the valence bands, excluding the first 320, up to 3 eV above the Fermi level. The set of trial orbitals chosen for the transformation is s, p and d orbitals for Hf; s and p orbitals for O; sp3 and d orbitals for Pt. This gives a total of 768, 776 and 784 Wannier functions for the LRS, HRS, and pristine state, respectively. Contrary to the Ag/MgO junction the best result is here obtained by directing using the projections of the trial orbitals on the Bloch state sub-space as Wannier functions to build the Hamiltonian (see eq.s 3.43 and 3.44). This corresponds to set zero iterations of Wannier functions optimization to minimize the quadratic spread (eq. 3.42). Probably the Wannier functions considerably move after a few optimization steps, breaking the electrodes periodicity. The Wannier functions anyway remain well localized (Ω < 3 Å2 ). On the contrary, the optimization of the Bloch state basis-set ("disentanglement step") remains fundamental to obtain a good result, and is performed by iterating 1000 times.

For the transmission function calculation, as for the Ag/MgO junction, a refined 10 × 10 k-points grid is sufficient to reach convergence.

Transport calculation on an OxRAM model: Pt/HfO 2 /Pt

Reliability tests on the pristine structure

In this section, we present some tests to check the reliability of QE for electronic transport calculations, and to check if the use of the Hamiltonian cutoff when using one single k-point along the transport direction is confirmed as a good approach for HfO 2 structures. Figure 6.4 shows the comparison very clean, regular and symmetric interface like Ag/MgO yields to a transition which strongly blocks the transmission, whereas a rough, irregular and asymmetric interface like Pt/HfO 2 yields a high transmission of electrons. HSE calculations on the Pt/HfO 2 structures are critical for the computational resources. Contrary to PBE, the HSE level requires the use of 2 k-points along the transport direction to obtain a reliable transmission function. Nevertheless, this greatly increases the need for computational resources. Hence it can be interesting to consider the use of one single k-point in combination with a Hamiltonian cut-off. Figure 6.5 shows that this is a valid strategy. The result Transport calculation on an OxRAM model: Pt/HfO 2 /Pt obtained by one single k-point and a cut-off is in good agreement with the 2 k-points and no cut-off one, but computationally less expensive. Therefore the transmission functions of the LRS and the HRS at the HSE level can safely be computed by one k-point and a cut-off to save computational resources. Firstly the thickness of HfO 2 is different (1 nm against 10 nm), and more important the form of the filament is unknown and certainly very different in the two cases. Remark that, due to the lateral 4 × 4 k-points sampling, our filament is composed of 16 mono-atomic straight lines of oxygen vacancies, which probably interact with each other, whose conduction is averaged to obtain the conduction of a single line. In real devices, the conductance in the HRS and LRS is mainly given by the conduction through the filament, and thus it is quite insensitive to the lateral area of the oxide (e.g. figure 7 of Cagli et al.). Thus the LRS and HRS values of table 6.1 represent the conductance of the single filament and are indeed given in Ω -1 . It would be nonsense to give a linear dependence on the oxide lateral area for the LRS and the HRS since this would correspond to multiply the number of filaments. Such linear dependence is appropriate only for the pristine state. We can thus try to compare the values directly extracted from our calculations to experiments, without scaling the oxide lateral area. This comparison represents the capacity to carry the current of the simulated guessed filament compared to the true one in the experimental sample. We could expect the simulated conductance to be much higher than the experimental one since the oxide is thinner and the tunneling regime should be stronger. The work of Cagli et al. gives values around 1.17 × 10 -3 Ω -1 for the LRS and 1.25 × 10 -5 Ω -1 for the HRS, calculated on a sample with HfO 2 thickness around 10 nm. The difference between the two is in very good agreement with the values of table 6.1, but the values are unexpectedly much higher. This indicates that the simple modeling of the conductive filament as a mono-atomic line of oxygen vacancies is naive and not suitable. Probably the true filament is composed of a bigger area of oxygen vacancies which can carry much more current, and probably the size of the filament plays a more important role than the thickness of the oxide. On the contrary, the fact that the difference between HRS and LRS conductance is well reproduced may prudently indicate that the model of HRS as an interrupted filament, with the re-insertion of some oxygen atoms, is rather acceptable independently on the true filament shape. The comparison between theoretical and experimental results should also take into account the fact that we are computing a zero-bias conductance, while to measure the resistance, and thus the conductance, of a true sample a small bias has always to be applied. However, read bias values are usually lower than 0.1 V [START_REF] Traoré | Étude de cellules mémoires résistives RRAM à base de HfO 2 par caractérisation électrique et simulations atomistiques[END_REF][START_REF] Nail | Investigation of hybrid CBRAM/OXRAM non-volatile memories for low consumption and high reliability[END_REF], which is close to the linear response regime at room temperature about 0.05 V (see eq. 2.20). Therefore, from this point of view, the comparison is almost acceptable. From the comparison with experiments, at this point, we have no track to claim if the HSE correction to eigenvalues and wavefunctions should be a fundamental ingredient for these kinds of structures. However, figures 6.8, 6.9 and 6.10 show a systematic reduction of the tunneling conduction always about 1 order of magnitude. Since in literature, from theoretical and empirical arguments, it is generally accepted that the transport calculation at the LDA or GGA level is systematically overestimated for basically all systems, we conclude that HSE corrections improve the transport properties calculation of HfO 2 based OxRAMs. Therefore, to go towards an improvement of

LRS and HRS by PBE and HSE calculations

Conclusion about OxRAM calculations

141 the agreement between theory and experiments on these particular systems, which is crucial to continue the development, we consider that the use of QP calculations is a fundamental step.

Conclusion about OxRAM calculations

In this chapter, we tried to build a supercell that is the closest possible to a true OxRAM device, in terms of materials and structure. By computing a more accurate electronic structure by HSE simulations the purpose was to go towards a better agreement between theory and experiments, and improve our knowledge about OxRAM devices transport mechanisms. The task was made rather hard by the compromise with the high computational cost, and by the fact that the structure of the conductive filaments is basically unknown. Hence the purpose was also to improve the knowledge of the conductive filament structure by matching experimental measurements. We built a very trivial model of filament, composed by a mono-atomic wire of oxygen vacancies, and we concluded that it cannot constitute a good model for a conductive filament in real devices. This result was expected and not surprising. On the contrary, we found that the HRS model as an interrupted filament may be consistent.

Since the PhD required a considerable amount of time to develop the many body Landauer method, there was less space for applications. Thus the results and the analysis of this chapter are preliminary but open to wide and interesting perspectives. To move forward the size of the supercell should be increased along the transport direction and laterally, to replicate the size of true OxRAM devices and avoid filament interaction due to PBC. Then different forms of filament should be guessed to improve the agreement with experiments. In any case, the analysis of this chapter and the whole work of this doctoral thesis seem to indicate QP calculations (or semi-empirical many body calculations as HSE) to play a fundamental role for this kind of study.

Conclusions

In the next decades, electronic devices will acquire more and more relevance in our everyday life. The development of new technologies will revolutionize our world in a way we cannot foresee today. This revolution will require all the components of electronic devices to increase their performances, to fulfill more and more demanding tasks. Within this background nowadays several small companies, multinational corporations, research centers, and universities are focusing their efforts on the development of more advanced memory devices for data storage. The technology currently used to store high volumes of data in electronic devices is the flash technology, which is a solid state disk with excellent properties of scalability, a low price, and the fundamental characteristic of being non-volatile. Nevertheless, they pay a gap of speed performances compared to modern CPUs. Only volatile memory devices can work at speeds comparable with modern CPUs, but they cannot be used to store high volumes of data permanently without a continuous power supply. Therefore, to improve the overall performances of electronic devices, the development of a new memory technology is necessary. This technology should ideally maintain, and possibly increase, all the good characteristic of flash memories, reaching the working speed of non-volatile memories. Good candidates for this new memory technology are the non-volatile resistive RAM (RRAM), whose working principle is based on the modification of their intrinsic resistive state. OxRAM devices, which exploit the resistive state of a metal oxide as logical bit, are among the most promising RRAMs. However, they are not the most developed and mature technology. The nature of the conduction mechanisms in the low and high resistive states (LRS and HRS), the structure and nature of the conductive filament (CF) in the LRS, the I/V characteristic, and other properties of OxRAM devices, are still under debate in literature, thus preventing to solve the drawbacks which still affect this technology. Among them, the most penalizing are the non-uniformity of switching parameters between the two resistive states, and the resistive state failure which still occur too often after a high, but limited, number of cycles. The initial purpose of this doctoral work was very ambitious and perspec-

Conclusions

tive. Try to reach a better theoretical comprehension of the OxRAMs working principle in order to lead the development of the present technology and overcome its limitation. This better understanding should include the conduction mechanism of HRS and LRS, the structure of the CF, the theoretical calculation of the I/V characteristic, and the dynamics of the CF forming and breaking. To fulfill this purpose we chose to use ab initio simulations to compute transmission function and conductance of OxRAM-like structures. The advantage of this approach is that the microscopic structure, which is proved to deeply influence the transport properties, can be perfectly controlled. By guessing different forms of CF we are able to study the impact to the conductance, trying to match experimental measurements to improve our knowledge about the CF structure in real devices, and about CF conduction mechanisms. The final ambitious purpose is to compute the I/V characteristic of OxRAM devices. This strategy requires a robust and reliable theory to compute transmission function and conductance of bulk junctions from first principles. The most common theory to fulfill this task is the Landauer theory. The input required by the Landauer formalism, that is basically the system Hamiltonian, is usually computed by density functional theory (DFT). The resulting "DFT + Landauer " method has some limitations and reliability issues. Limitations concern non-equilibrium calculations and scattering events (e.g. electron-electron or electron-phonon scattering), which cannot be taken into account. Reliability issues are due to the local, mean-field character of DFT, which systematically underestimates the band gap, and provides an electronic structure which corresponds to an auxiliary non-interacting system. This PhD work firstly aim to overcome these reliability issues. Therefore we developed a method, named "many body Landauer ", which uses quasiparticle (QP) calculations to compute a more accurate electronic structure as input for the Landauer formalism. The main hurdle to overcome is that highly tested QP calculations on bulk systems are practically only implemented in plane-waves codes, while the Landauer formalism requires an implementations based on a localized functions basis-set. The solution adopted in this work is to build a procedure based on a basis-set transformation from planewaves to Wannier functions, which are directly, but not uniquely, linked to Bloch states. The first phase of this doctoral work was focused on QP calculations. We considered Hafnium oxide (HfO 2 ), which is among the most promising materials for OxRAM devices, and we performed a deep study of its electronic and optical properties by means of highly converged GW calculations. In literature experimental and theoretical works generally agree about optical and electronic gap values of HfO 2 , which are both believed to lie between 5.5 and 6 eV for the monoclinic phase, the most stable at room temperature and pressure. Nevertheless this produces some confusion between electronic and optical gap. Moreover there are some doubts about the convergence of the simulations, and about the procedure to extract values for the electronic gap from photoemission measurements. To make clear the situation we pushed to the limits convergence and theoretical framework of our simulations, by using a QSGW approach for the electronic gap and solving the Bethe-Salpeter equation (BSE) on top for the optical gap. Our study lead to the conclusion that the electronic gap (indirect Γ → B) is probably higher than what believed so far, with a lower bound around 6.33 eV. On the contrary the optical gap (direct Γ → Γ) is lowered to 5.85 eV by the presence of an exciton with binding energy of 0.57 eV, in agreement with spectroscopic ellipsometry measurements. Optical and electronic gaps are thus well separated, with the electronic higher than generally accepted. This characterization of HfO 2 can be an important indication for all HfO 2 based technologies, and also for OxRAM applications. This argument was object of publication during the second year of PhD [START_REF] Sklénard | Optical vs electronic gap of hafnia by ab initio bethe-salpeter equation[END_REF].

The longest and hardest part of this doctoral work was the development and the test of the many body Landauer method. The tests were performed on a junction made of Ag/MgO/Ag. The method proved to be able to replicate the results of state of the art codes at the DFT level. The method revealed also to be very versatile, able to compute transmission function and conductance beyond DFT, on top of many different QP calculations: HSE, G 0 W 0 and selfconsistent COHSEX. Among them G 0 W 0 turned out to be not suitable for electronic transport calculations since it only updates eigenvalues leaving the Kohn-Sham wavefunctions unmodified. This affects the non-diagonal terms of the Hamiltonian, which artificially raise producing noise in the transmission function. This fact also led to the interesting conclusion that it is the updating of wavefunctions to mainly drive the lowering of the transmission function towards more accurate values, whereas eigenvalues play a less crucial role. Differently from G 0 W 0 , HSE and COHSEX revealed to significantly improve the agreement with experiments, improving the reliability compared to the standard DFT + Landauer method. The basis-set transformation turned out to be the most delicate part of the procedure. The transformation is not unique and should be chosen carefully. Moreover, since Wannier functions are localized but also periodic, it can be sensible to the k-point sampling, which defines the period of Wannier functions. If the period is too short the Wannier functions tails are artificially raised, causing anomalous high non-diagonal terms in the Wannier Hamiltonian, invalidating the transmission function calculation especially when the conductance is weak. For this reason, the use of 2 k-points along the transport direction is often manda-Conclusions tory. A paper about the application of the many body Landauer method on Ag/MgO/Ag junctions was published in Physical Review B at the end of the third year of PhD [START_REF] Dragoni | Plane-wave many-body corrections to the conductance in bulk tunnel junctions[END_REF]. This method improves the reliability of transport calculations, but its limitations remain the same as the DFT + Landauer approach. The many body Landauer method constitutes a fundamental step to go towards a better agreement with experiments, but it still requires a great work of development to include calculations under bias and scattering events. Its inclusion in a full iterative NEGF scheme is the natural continuation and the challenging perspective of this work. In a NEGF scheme, the effect of the finite bias and the non-equilibrium situation is taken into account by introducing, beyond the retarded or advanced Green function, two further correlation functions, G < and G > , describing the out-of-equilibrium distributions of electrons and holes, respectively. The G < and G > then enter in all further steps of the calculation, for example in the calculation of the non-equilibrium electronic density which then affects (via the Poisson equation) the classical Hartree potential, but also directly in the exchange-correlation part of the self-energy and scattering functions. The scattering events are accounted for by using G < and G > to compute additional self-energy terms, which act as effective Hamiltonian. By combining these elements a new Green function is computed, and the procedure is iterated until convergence is reached. The many body Landauer method, which is based on the Green function formalism and uses a localized function basis-set, is in principle naturally suitable for an extension to NEGF. Nevertheless, the hurdles to reach a similar result are multiple, and several compromises about computational cost and implementation have certainly to be done.

In the last part of PhD, we tried to build structures more similar to real OxRAM devices, made of Platinum (Pt) and HfO 2 , trying to guess a form for the CF in the LRS and HRS. Applying the many body Landauer method, at the HSE level, to these structures gave interesting, but preliminary, results. The analysis showed that QP calculations seem to be fundamental also on OxRAM-like structures. The trivial model of CF as a mono-atomic oxygen vacancies wire revealed to be unrealistic. However, the hypothesis of HRS as corresponding to a broken filament of oxygen vacancies found some confirmations. Due to lack of time, the study of OxRAM structures performed in this work is only at the very beginning. Many guesses of CF have to be simulated, increasing the supercell size to more realistic values, but always caring about computational resources. In this doctoral work, we have reached partial, but very interesting and perspective results. An accurate characterization of the electronic and optical properties of HfO 2 has been performed, providing important information for any kind of HfO 2 based technology. A more reliable method to compute the electronic transport properties of bulk junctions has been developed. This could become in the close future a fundamental step to improve the agreement between theory and experimental measurements of conduction properties. In perspective, the insertion of the many body scheme in a full NEGF formalism will be fundamental to further improve the theoretical framework of transport calculations. Finally, we have started the study of OxRAM-like structures by means of the many body Landauer approach. However, a full comprehension of the OxRAM working principles will require in future many other studies. The final message we would like to underline in this work is that the correction of the DFT electronic structure by means of QP calculations before applying the Landauer formalism is a fundamental step to obtain reliable results of transmission function and conductance, and to reduce the distance with experimental measurements. This scissor operator is analogous to a G 0 W 0 correction, where only the eigenvalues are corrected. This simple model illustrates how an operator which is diagonal in Bloch space can increase the range of the Hamiltonian in a basis of localized orbitals. This may explain the long range couplings of the G 0 W 0 Hamiltonian in figure 5.12. Ferretti et al. found a similar behavior on a different tight-binding model [START_REF] Ferretti | Ab initio complex band structure of conjugated polymers: Effects of hydrid density functional theory and GW schemes[END_REF].

Appendix B

Conductance decay rate from bulk MgO complex band structure [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF][START_REF] Rungger | Algorithm for the construction of selfenergies for electronic transport calculations based on singularity elimination and singular value decomposition[END_REF]. The CBS can be used to infer the decay rate of the tunneling conductance of a junction composed by n MgO cubic cells between two electrodes, when n increases: G ∝ exp(-2κan). The factor 2 comes from the fact that the tunnelling probability is proportional to the square modulus of the evanescent state. In the MgO band gap, the evanescent states with the slowest decay (i.e. smaller κ) have the same symmetry as the states of the Ag electrodes at the Fermi energy (see appendix C for details). Hence the decay rate of the conductance is related to the smallest value of κ at the Fermi energy E f . However, the position of E f depends on the band alignment between Ag and MgO. In figure B.1, we have defined the Fermi level from the band alignement of the thickest structure (Ag/4MgO/Ag): for each approximation, we have extracted the difference between E f and the valence band edge from Fig. The separation between DFT and G 0 W 0 on one hand, and COHSEX and HSE on the other hand, is confirmed. Some differences may arise from the presence of the Ag/MgO interface when computing κ from the transmission spectra. It is not surprising that the decay factor of a junction is dependent on the material the electrodes are made of. Thus we don't expect to find exactly the same values, especially for thin structures as the ones considered here, where the bulk properties of MgO are surely influenced by the presence of the Ag electrodes. Moreover, the offset between the Fermi level of Ag and Resistive non-volatile memories based on oxides (OxRAM) are recently acquiring a wide interest for their performances, which make them promising candidates as storage memories to replace flash technology, and as embedded memories for neural network applications. Nevertheless, emerging OxRAM devices still present some drawbacks, like non-uniformity of switching parameters and switching failures. Overcoming these drawbacks requires a deeper comprehension of the OxRAM working principles, so far not completely understood. This can be achieved by means of ab initio simulations. Hence this work presents a careful characterization of HfO 2 , which is within the most promising materials to build OxRAM devices, by means of accurate quasiparticle (QP) calculations. A study of the electronic transport properties in OxRAM devices is also of primary importance. However, this requires a robust and reliable theoretical framework to compute the conductance of bulk metal/insulator junctions. The standard approach, based on density functional theory, Green function formalism, and Landauer formula, has some limitations and reliability issues. This work proposes a more reliable approach based on QP calculations, which provide a more accurate electronic structure to compute the conductance, and largely tests this new method on different junctions mimicking OxRAM devices.

Les mémoires résistives non volatiles basées sur les oxydes (OxRAM) acquièrent récemment un grand intérêt pour leurs performances, ce qui en fait des candidats prometteurs comme mémoire de stockage pour remplacer la technologie flash, et comme mémoires intégrées pour les applications réseau de neurones. Néanmoins, les dispositifs OxRAM émergents présentent encore certains inconvénients, comme la non-uniformité des paramètres de commutation et les défaillances de commutation. Surmonter ces inconvénients exige une compréhension plus profonde des principes de fonctionnement de l'OxRAM, jusqu'à présent pas complètement compris. Ceci peut être réalisé au moyen de simulations ab initio. Ce travail présente donc une étude approfondie de HfO 2 , qui fait partie des matériaux les plus prometteurs pour la construction de dispositifs OxRAM, au moyen de calculs précis des états de quasi-particules (QP). Une étude des propriétés du transport électronique dans les dispositifs OxRAM est également de première importance. Toutefois, cela nécessite un cadre théorique solide et fiable afin de calculer la conductance des jonctions métal/isolant. L'approche standard, basée sur la théorie fonctionnelle de la densité, le formalisme de la fonction de Green et la formule de Landauer, a quelques limites et soucis de fiabilité. Ce travail propose une approche plus fiable basée sur les calculs QP, qui fournissent une structure électronique plus précise pour calculer la conductance, et teste en grande partie cette nouvelle méthode sur différentes jonctions imitant les dispositifs OxRAM.

Figure 1 . 1 :

 11 Figure 1.1: a) Flash memory device structure. b) The two apparent voltage states used as logical bit. Image adopted from ref. [2].

Figure 1 . 3 :

 13 Figure 1.3: Non-volatile vs volatile memory devices. Their features are in general complementary and well separated. Volatile memories have better performances, but the great disadvantage of requiring a constant power supply. RRAM is proposed as a non-volatile candidate, which should posses features close to volatile memories.

  fig. 1.4). The current and voltage characteristics of these three operations are shown in fig. 1.5. The electrical I/V characteristic of OxRAMs can be

Figure 1 . 7 :

 17 Figure 1.7: Sketch of possible electron conduction paths in OxRAMs. (1) Schottky emission. (2) Fowler-Nordheim (FN) tunneling from cathode into conduction band (CB). (3) Direct tunneling. (4) Tunneling from cathode to trap. (5) Poole-Frenkel emission from trap to CB. (6) FN from trap CB. (7) Trap to trap Mott hopping. (8) Tunneling from trap to anode. Adopted from ref.[START_REF] Wong | Metal-Oxide RRAM[END_REF].

Figure 1 . 8 :

 18 Figure 1.8: Electrochemical processes during resistive switching of a CBRAM. The first three steps correspond to the set operation. (1) oxidation, (2) migration, (3) reduction. The fourth step is the reset operation. The corresponding values of current and voltage are shown in the center. Image taken from [38].

1. 4 .

 4 Figure1.9: Sketch of structure (left) and programming electrical power (right) of a PCM. Adopted respectively from[START_REF] Pershin | Memory effects in complex materials and nanoscale systems[END_REF] and[START_REF] Burr | Phase change memory technology[END_REF].

Figure 1 . 10 :

 110 Figure 1.10: Sketch of FeRAM polarisation states.

Figure 1 .

 1 Figure 1.12 summarizes the memory technologies described in this chapter, classifying them according to volatile/non-volatile characteristic, working principle, and state of development (emerging or current).

Figure 1 . 12 :

 112 Figure 1.12: Classification diagram of current memory technologies. Adopted from ref. [6].

Figure 2 . 1 :

 21 Figure 2.1: A conductor of length L and surface S connected to two larger contacts. In general conductor and contacts can be made of different materials.

Figure 2 . 2 :

 22 Figure 2.2: Example of dispersion relation E( k) into the conductor. k runs along the transport direction, while it is fixed in the two directions perpendicular to transport. At zero temperature the electrons occupation η is a function with a step located at the Fermi energy E = E f . In this situation, two bands cross E f which correspond to two conductive modes.

Figure 2 . 3 :

 23 Figure 2.3: Conductive modes for the same dispersion relation E( k) of figure 2.2 at a temperature θ = 0. The occupation function, which here is assumed equal to the Fermi distribution η f (E), and its negative derivative are also plotted. This latter weighs the contribution to the conductance of the states at a given energy in the integral 2.3. As consequence only states in a range of energy from E f -few k B θ to E f + few k B θ significantly contribute to the conductance. As effect of thermal energy the conductive modes increase from 2 to 3 compared to the zero temperature case of figure 2.2.

Figure 2 . 4 :

 24 Figure 2.4: Equation2.5 is plotted for different values of temperature θ (left y-axis and legend). The area under the function is equal to one and always conserved, but distributed on a wider energy range while increasing θ (remark that the blue arrow indicates a Dirac delta). Equation2.5 can be approximated by a Gaussian distribution (eq. 2.6) centered at E f , with standard deviation σ 1.6 • k B θ (right y-axis and legend).

Figure 2 . 5 :

 25 Figure 2.5: A conductor stretched between two leads, whose chemical potentials µ L and µ R are imposed by two reservoirs, which produce a voltage on the conductor.

13 )Figure 2 . 6 :

 1326 Figure 2.6: Conductive modes for the same dispersion relation E( k) of figure 2.2 and 2.3 at a temperature θ = 0, and with an applied voltage equal to (µ L -µ R )/e. Only states in a range of energy from µ R -few k B θ to µ L + few k B θ contribute to the conductance. As effect of the applied bias the conductive modes are here equal to 4.

Figure 2 . 9 :

 29 Figure2.9: Gold mono-atomic wire band structure at zero bias and correspondent τ (E), which is equal to the number of bands at each energy (be careful to degenerate levels). The k-points path goes from Γ (0,0,0) to Z (0,0,0.5) along the wire direction in the first Brillouin zone. The DFT cell is composed by three gold atoms at a distance of 2.56 Å (relaxed in GGA) surrounded by 5 Å of vacuum.

Figure 2 .

 2 Figure 2.10: A conductor connected to two semi-infinite leads.

Figure 2 . 11 :

 211 Figure 2.11: Free electrons waves incident on a one-dimensional squared potential barrier U (x) from left, with amplitude L in , and right, with amplitude R in . The electrons scatter through the barrier resulting in two outgoing waves: on the left, with amplitude L out , and on the right, with amplitude R out .

Figure 2 . 12 :

 212 Figure 2.12: Example of a conductor with three total conductive modes at a given energy, two for the left lead and one for the right one. The S-matrix is of dimension 3 × 3, relating the outgoing amplitudes b n to the incoming amplitudes a n .

The many body Landauer transport 3 . 1

 31 Step 1: Building the supercell In this work, the focus is on metal-insulator-metal (MIM) junctions, which are the fundamentals of many technological applications: 2D materials, metalsemiconductor contacts, PCRAMs, CBRAMs and especially OXRAMs. The supercell should be suitable to extract the Hamiltonian sub-matrices of the central conductor, of the bulk periodic electrode principal layer, and the couplings between different electrode layers and between central part and electrodes. They are necessary to build the Green function 2.58, which represents the open-system with semi-infinite leads. Therefore the supercell must contain not only the truly central conductor under study, but also some layers of the metal electrodes are made of.

Figure 3 . 1 :

 31 Figure 3.1: Example of Platinum (Pt), Hafnium oxide (HfO 2 ) junction. The supercell, shown along the transport direction, is composed by 3 cubic cells of Pt along the 111 direction on the two sides, and by two primitive cells of monoclinic HfO 2 along the 001 direction in the middle. The first two cubic cells of Pt on both sides are fixed to the bulk structure and identical, to guarantee the periodicity. The cubic cell close to HfO 2 on both sides is instead relaxed together with HfO 2 to create a good interface. The length along transport of one electrode is about 21 Å, which means a distance of 42 Å between two PBC replicas of HfO 2 . This distance assure a very weak coupling between the two replicas.

3. 2 .

 2 Step 2: Electronic structure correction 63 Notice that the quasiparticle equation 3.9 reduces to the Hartree equation when Σ = 0, to the HF one when Σ = iGV C (being V C the bare Coulomb potential), and to the KS equation when Σ = V xc .

Figure 3 . 3 :

 33 Figure 3.3: The Hedin pentagon in its complete and approximated versions. Every corner correspond to a term of the Hedin equations 3.10 -3.14. Arrows indicate the calculation path. a) Complete Hedin equations. b) Self-consistent GW approximation. c) Oneshot G 0 W 0 approximation. Taken from [101].

  3.3 b)), or one-shot to obtain Σ GW in one iteration (named G 0 W 0 , see fig. 3.3 c)).

Figure 3 .Figure 3 . 4 :

 334 Figure 3.4: Band gaps in solids as calculated by DFT-LDA, HF and G 0 W 0 in comparison with experimental measurements.

  3.5). Outside the ES PBC apply for the Bloch states, and for WF as consequence since they are linear combinations of Bloch states (see fig.3.5). The translation vectors between a given Wannier function and its PBC replicas outside the ES are given by:

Figure 3 . 5 :

 35 Figure 3.5: Sketch of some WF shown in the RS (black) and in the ES (grey) along a direction sampled by two k-points. Hence the ES is twice longer (2L) than the RS (L). ψ n k = e i k• r u n k (brown, purple and orange line, respectively) are also sketched for k = 0. Plain lines indicate the original objects (cells and functions) and dashed lines indicate the replicas. u n k has the periodicity of the RS, while the Bloch state ψ n k has the same periodicity of the phase factor e i k• r , which coincide with the ES length. The WF |0n (red line) belongs to the RS, as |0n (green line), which is built from a set of Bloch states belonging to another band n . Outside the RS, but inside the ES, the interference between the Bloch states is destructive and these WF approach zero. Outside the ES they are instead replicated with a period equal to 2L. Rn (blue line) is built from the same set n of Bloch state of |0n , but shifted in the RS replica into the ES.

  choose the arbitrary Bloch states phase factors. Equation 3.41 also combine the arbitrary phase factors of |ψ m k , which thus are automatically optimized in the minimization process to give the best continuity of ψ n k , and thus the best localization of the WF Rn . For metallic systems or in the conduction band of insulators, the Bloch bands are usually not isolated in energy. It is then necessary to use the extension of the Marzari and Vanderbilt formulation to non-isolated bands, developed by Souza et al.[START_REF] Souza | Maximally localized Wannier functions for entangled energy bands[END_REF]. Then the N w WF are built from a larger number N b of Bloch bands. Hence U ( k) are N unitary complex matrices with orthonormal columns of dimension N b ×N w , and the indexes m and n run over N b and N w , respectively, in equation 3.40. N w corresponds to the number of Bloch states in a frozen energy window which must be exactly replicated by the WF. The Bloch states outside the window are instead only used to extend the linear combination 3.41 to improve the optimization, but are not well reproduced in the WF basis-set.

n

  as WF basis-set. The trial orbitals |Φ n satisfy by construction the periodicity of the crystalline regions, and so do the WF Φ ⊥ n obtained by equations 3.43 and 3.44. Further optimization of the basis-set by minimization of the spread, if pushed too far, may break the periodicity by mixing orbitals differently on each electrode cubic cell.

Figure 3 .

 3 [START_REF] Nail | Investigation of hybrid CBRAM/OXRAM non-volatile memories for low consumption and high reliability[END_REF] shows an example of WF obtained by the method of Marzari and Vaderbilt for monoclinic HfO 2 .

Figure 3 . 6 : 3 . 4 .

 3634 Figure 3.6: Example of Wannier functions on the primitive cell of monoclinic HfO 2 . The one centered on the oxygen is obtained by a 2 s and a 2 p orbital initial guess, and exhibits sp shape and character. The one on the Hafnium comes from a 5 d orbital, and exhibits d character.

  Figure 3.7: a) Representation of a localized orbital |φ 0n centered into the RS. Outside its corresponding WS cell |φ 0n is set to zero. As consequence replicas disappear. b) Representation of how the coupling between two different localized orbital |φ 0n and |φ Rn has to be interpreted. The distance within the two is L = | R + r Rn -r 0n |, while L = | r 0n -R -r Rn + T | is the distance between |φ 0n and the |φ Rn replica translated by T . Since L < L, the coupling should be meant at a distance equal to L between the two localized orbitals, which is also physically more coherent.

3. 4 .

 4 Step 4: Wannier Hamiltonian construction and transmission function calculation 81 zero), and not as an integral over the ES. Note that replacing R by R + T in the right-hand-side would not change the matrix element, since the corresponding phase shift -k • T in equation 3.45 equals 0 modulo 2π for all k in the MP grid. If we find several equidistant nearest orbitals, which occurs only if they are at the surface of the WS cell, we consider that |φ 0n is coupled to each of these N r orbitals by a reduced matrix element 0n H R n ES /N r . Finally, other matrix elements are deduced by translation symmetry:

Figure 3 . 8 :

 38 Figure 3.8: Example of tunneling junction composed by 4 silver cubic cells electrodes and 2 cubic cells of magnesium oxide in the middle. The Hamiltonian blocks used for transport calculations are schematized: H C is the Hamiltonian of the central region and the blocks H 0i define the semi-infinite periodic Ag electrodes.

Figure 4 . 1 :

 41 Figure 4.1: QP corrections of the Kohn-Sham eigenvalues as a function of the inverse of the number of empty bands N b (or number of plane waves N pw ). The dotted lines show the linear extrapolation of QP corrections (see text for details).

Figure 4 . 2 :

 42 Figure 4.2: Convergence study of the monoclinic HfO 2 band gap at Γ with respect to the number of empty bands used in the G 0 W 0 calculation. The number of plane-waves is fixed to the value given by an energy cutoff of 500 eV.

Figure 4 . 3 :

 43 Figure 4.3: Convergence of QP band gap as a function of self-consistent GW iterations at GW 0 (top panel) and QSGW (bottom panel) level.

Figure 4 . 4 :

 44 Figure 4.4: Imaginary part of the dielectric function ε 2 for the monoclinic phase of HfO 2 calculated on top of PBE+G 0 W 0 for 4 × 4 × 4 and 6 × 6 × 6 k-mesh. A Gaussian broadening of 0.1 eV is used.

Figure 4 .

 4 Figure 4.5 shows the electronic band structure of monoclinic HfO 2 , as calculated by the PBE functional and by QSGW , which is the most advanced theoretical framework here considered, on top of it. In the monoclinic crystal

Figure 4 . 5 :

 45 Figure 4.5: Calculated band structure of monoclinic HfO 2 within the DFT-PBE (dotted line) and QSGW (plain red line) approximations. The QSGW band structure has been interpolated using MLWF. PBE and GW Fermi energies are aligned at zero.

Figure 4 . 6 :

 46 Figure 4.6: Density of states (DOS) of monoclinic HfO 2 in the G 0 W 0 , GW 0 , and QSGW approximations compared to UPS+IPS[START_REF] Bersch | Band offsets of ultrathin high-κ oxide films with Si[END_REF] and XPS+IPS[START_REF] Sayan | Band alignment issues related to HfO 2 /SiO 2 /p-Si gate stacks[END_REF] spectra. The theoretical DOS has been interpolated using MLWF on a 40 × 40 × 40 k-point grid and convoluted with a gaussian broadening of 0.7 eV. Experimental spectra have been aligned at the valence band-edge (the zero of the energy for all theoretical DOS), and the conduction IPS and valence PS independent measurements have been rescaled separately to match the height of the theoretical DOS. The error in this procedure can be estimated by the maximum deviation among valence band-edges in the steepest rising linear region: 50 meV, which is less than our theoretical error bar.

2

 2 Figure 4.6: Density of states (DOS) of monoclinic HfO 2 in the G 0 W 0 , GW 0 , and QSGW approximations compared to UPS+IPS[START_REF] Bersch | Band offsets of ultrathin high-κ oxide films with Si[END_REF] and XPS+IPS[START_REF] Sayan | Band alignment issues related to HfO 2 /SiO 2 /p-Si gate stacks[END_REF] spectra. The theoretical DOS has been interpolated using MLWF on a 40 × 40 × 40 k-point grid and convoluted with a gaussian broadening of 0.7 eV. Experimental spectra have been aligned at the valence band-edge (the zero of the energy for all theoretical DOS), and the conduction IPS and valence PS independent measurements have been rescaled separately to match the height of the theoretical DOS. The error in this procedure can be estimated by the maximum deviation among valence band-edges in the steepest rising linear region: 50 meV, which is less than our theoretical error bar.

Figure 4 . 7 :

 47 Figure 4.7: Imaginary part of the dielectric function ε 2 averaged over the three monoclinic HfO 2 lattice directions and convoluted with a gaussian broadening of 0.2 eV, compared to ellipsometry spectra (SE) [122].

  4.8), prooving the unusual anisotropy 4.4. Conclusions about HfO 2 properties

Figure 4 . 8 :

 48 Figure 4.8: Imaginary part of the dielectric function ε 2 for the three HfO 2 phases and all, (100), (010), and (001) polarization directions. The arrow indicates the only non-dark bound exciton in the (010) polarization of monoclinic HfO 2 .

Figure 5

 5 Figure 5.1: a) Ag atomic positions in the cubic cell. It is the standard FCC structure, atoms lie on the corners and the centers of the surfaces. The cell contains a total of 4 Ag atoms. b) MgO atomic positions in the cubic cell. Mg atoms have FCC structure, O atoms lie on a translated FCC lattice. Hence they lie in the middle of the cube sides and one oxygen is placed in the center of the cubic cell. There are a total of 4 Mg and 4 O. The Mg-O bonds are also shown.

  .

  2). Due to periodic boundary conditions, such a system can also be viewed as a superlattice of Ag layers with 8 cubic cells separated by MgO layers. The 8 Ag cubic cells are sufficient to screen the MgO/MgO interaction.

Figure 5 . 2 :

 52 Figure 5.2: 4Ag/2MgO/4Ag supercell for transport calculation.

Figure 5 . 4 :

 54 Figure 5.4: 4Ag/4MgO/4Ag structure along the transport direction. a): The RS (black dashed rectangle), and the superlattice formed by PBC. b): Sketch of a WF, shown as corresponding to 1 and 2 k-points. The first has the periodicity of the RS (blue area), and one replica is visible around 65 Å. Grey vertical lines indicate the borders of its Wigner-Seitz cell. The second has the periodicity of its ES, which is twice longer than the RS and coincides with the whole considered x -axis range. Replicas are thus not visible.

5. 3 .Figure 5 . 5 :

 355 Figure 5.5: Modulus of Ag-Ag maximum couplings for the 4Ag/4MgO/4Ag supercell as a function of distance along the transport direction, defined by the left-hand side of eq. 3.54. Before a distance of approximately 17 Å, which approaches the MgO thickness, the Ag-Ag couplings are direct, strong and decay exponentially with distance. For these couplings, there is no relevant difference between 1 and 2 k-points cases. Beyond 17 Å only Ag-Ag tunneling couplings through MgO are possible. In the 2 k-points case these couplings go to negligible values (under 10 -4 eV), in the 1 k-point case they are 1, 2 order of magnitude higher because of the overestimated WF tails values (see also fig.5.4 b)). Note that, due to the Hamiltonian construction, for 1 k-point the maximum coupling distance is L/2 (the ES length is L). For 2 k-points the maximum coupling distance is doubled up to L (the ES length is 2L).

Figure 5 . 6 :

 56 Figure 5.6: Hamitonian matrix at k = 0.0 b 1 -0.4 b 2 -0.4 b 3 representing the WF couplings in the transport direction for the 4Ag/4MgO/4Ag supercell. By using 1 ab initio k-point (left) the couplings are defined until L/2 26 Å. By using 2 k-points the couplings are defined on the whole supercell length (L 52 Å). On the right, the couplings are cut beyond 13 Å. The 1 k-point Hamiltonian is affected by tails overestimation, which yields an overestimation of some distant couplings (Ag-Ag into the black triangle), and a spiky overestimated transmission spectrum (top line of fig. 5.3). Surprisingly the 2 k-points and the 13 Å cut Hamiltonian matrices, which appear very different, give rise to the same transmission spectrum (bottom lines of fig. 5.3).This proves that the physics lies in a narrow region close to the Hamiltonian matrix diagonal, and distant non-diagonal terms can be neglected.

Figure 5 . 7 :

 57 Figure 5.7: a): Band structure of the 4Ag/4MgO/4Ag supercell along a k-line parallel to the transport direction: k = k 1 b 1 + 0.25 b 2 + 0.25 b 3 . The reference ab initio band structure is compared to the Wannier band structures interpolated from 1 or 2 k-points along the transport direction.b): Same plot with the energy range restricted to the band close to E F . The interpolation from 1 k-point strongly overestimates the band width. Note that k 1 = 0.5 belongs to the 2 k-points ab initio MP grid, but not to the 1 k-point grid. Thus the band structure at this point is replicated exactly by the 2 k-points Wannier basis, but is interpolated when using the 1 k-point Wannier basis.

Figure 5 . 8 :

 58 Figure 5.8: Average transmission function (eq. 3.66) computed by our plane-waves+Wannier method for the 4Ag/3MgO/4Ag supercell, shown on an energy range close to the Fermi energy E f , and for different k-points grids transverse to the transport direction. The top blue line corresponds to the ab initio MP grid (4 × 4), and thus to the spectrum computed by means of no interpolation. Increasing the grid significantly lowers the spectrum, which converges around a 10×10 grid. If 10×10 and 13×13 are almost overlapping, 7 × 7 is still overestimated at energies lower than E f .

Figure 5 . 9 :

 59 Figure 5.9: Average transmission function through MgO barriers of different thickness, from 1 cubic layer (1L) to 4 cubic layers (4L). Electrodes are as usual made by 4 Ag cubic layers. The plane-waves/Wannier spectra (plain lines labeled "Plane-waves") are computed at the DFT level and compared to those obtained by Transiesta (dotted lines).

Figure 5 .

 5 Figure 5.10: Conductance (eq. 3.67) of the 3MgO structure as function of reciprocal space transverse to the transport direction, at room temperature (26 meV). The k 2 and k 3 directions are both sampled by 50 interpolated k-points, yielding a 50 × 50 resolution of the conductance map.

2Figure 5 .

 5 Figure 5.11: Various transmission spectra obtained by applying the natural 2 k-points cut-off equal to L, and a shorter cut-off L cut = 13 Å: COHSEX and HSE for the Ag/4MgO/Ag structure (top), and G 0 W 0 for the Ag/2MgO/Ag structure (bottom left).

Figure 5 . 12 :

 512 Figure 5.12: Maximum Hamiltonian matrix element as a function of distance along the transport direction, shown for the 4 approximations considered and the Ag/2MgO/Ag and Ag/4MgO/Ag structures.

Figure 5 . 13 :

 513 Figure 5.13: Transmission function as calculated by various approximations, here shown for a tunneling barrier of different MgO thickness (1 to 4 cubic layers) with Ag electrodes. All the G 0 W 0 transmission functions, and the HSE and COHSEX ones for the thickest structure (4MgO), are calculated by applying a Hamiltonian cut-off of 13 Å (see section 5.5.1).

Figure 5 .

 5 Figure 5.[START_REF] Jameson | Conductive-bridge memory (CBRAM) with excellent high-temperature retention[END_REF] shows how G 0 W 0 , self consistent COHSEX, and HSE modify the DFT transmission function, for the 4 structures with increasing MgO thickness, from 1 to 4 cubic layers. We expect to observe a reduction of the transmission spectra compared to the DFT ones. This is confirmed for COHSEX and HSE, which for the 3 and 4 MgO structures reduce the T (E f ) value of more than 1 order of magnitude. On the contrary G 0 W 0 does not provide a significant reduction. For a more complete analysis, it is interesting to compare our results to the experimental results of Gangineni et al.[START_REF] Gangineni | Interfacial electronic transport phenomena in single crystalline Fe-MgO-Fe thin barrier junctions[END_REF]. In this work, measurements of

16 DFT 2 )Figure 5 . 14 :

 162514 Figure 5.14: Plot of the values of table 5.1. They are given in the same unit of the table (Ω -1 .µm -2 , right y-axis), and in unit of G 0 (left y-axis), corresponding to τ (E f ). This permits a direct visual comparison with the transmission functions of fig. 5.13. A linear regression is also performed on the conductance values logarithm for each approximation as function of the MgO thickness. This corresponds to an exponantial decay A • exp(-κ × n M gO ). Values of A and κ are given in the figure.

Figure 5 . 15 :

 515 Figure 5.15: Quasiparticle energies as function of the DFT-PBE eigenvalues shown for the Ag/1MgO/Ag structure.

Figure 6 . 1 :

 61 Figure 6.1: Monoclinic primitive cell of hafnium oxide.

Figure 6 . 2 :Figure 6 . 3 :

 6263 Figure 6.2: Pt/HfO 2 /Pt pristine state. Vertical lines separate different Pt cubic cells. The most external (2 per side) are kept fixed to bulk atomic positions. The central region (1Pt/HfO 2 /1Pt) is relaxed.

Figure 6 . 5 :

 65 Figure 6.5: Transmission function of the pristine structure computed at the HSE level by 2 k-points and 1 k-point in combination with a Hamiltonian cut-off L cut = 13 Å. The agreement between the two curves is very good.

Figure 6 .Figure 6 . 6 :

 666 Figure 6.6 shows the transmission function computed at the PBE level for the pristine state, the HRS and the LRS. The introduction of a mono-atomic

Figure 6 . 8 :

 68 Figure 6.8: Transmission function the pristine state, at the PBE and HSE levels.

Figure 6 . 9 :

 69 Figure 6.9: Transmission function the LRS, at the PBE and HSE levels.

Figure 6 . 10 :

 610 Figure 6.10: Transmission function the HRS, at the PBE and HSE levels.

Figure A. 2 :

 2 Figure A.2: Black lines: band structure of the 1D tight-binding model (see text). Red dotted line: rigid shift of the conduction band after application of a scissor operator of +2 eV.

Figure A. 3 :

 3 Figure A.3: Maximum modulus of the Hamiltonian matrix element vs. distance between orbitals, after application of the scissor operator in Bloch space. See model parameters in the main text.

Figure B. 1 :

 1 Figure B.1: Complex band structure of bulk MgO calculated at k = 0 within the 4 considered approximations.κ is the decay rate of the evanescent wavefunctions and a is the lattice constant. For each approximation, the difference between the Fermi level and the valence band edge is determined from the transmission spectrum of the Ag/4MgO/Ag structure (Fig.5.13).

Figure B. 1

 1 Figure B.1 shows the complex band structure (CBS) of bulk MgO within

Figure C. 2 :

 2 Figure C.2: Color maps of some Bloch wavefunctions at k = 0 in a (xy) slice of the cubic cell. State labelled 1 on the band structure plot (Fig. C.1) has symmetry ∆ 1 while state labelled 2 has symmetry ∆ 2 (it is doubly degenerate, the other one having mirror symmetry along the horizontal axis). The highest valence states and lowest conduction states of MgO both have ∆ 1 symmetry.

Figure C. 2

 2 shows a symmetry mismatch between the upper band of Ag and the conduction band of MgO. Color maps of the Bloch wavefunctions show that the lower band has symmetry ∆ 1 , the same as the states of MgO, while the upper band has symmetry ∆ 2 . Hence at the conduction band edge of MgO, which is well above E f + 1.6 eV even in the DFT case, the Ag states cannot be transmitted through the conduction states of MgO. They are probably transmitted by tunneling through evanescent states visible in the complex band structure of MgO (figure B.1).

  Ballistic conduction: There is a single scattering state of momentum k which carries current from |k, L to |k, R as long as it is occupied, without necessity of any state transition. The current is only limited by the interface barrier which electrons can tunnel with probability T .

	Electronic transport state of the art
	is small enough that the error resulting from applying equation 2.14, which
	approximates the difference of leads occupation functions by the occupation
	function derivative, is negligible. To find a reasonable criterion consider that,
	for an arbitrary function f , by definition of derivative f :
	Figure 2.7:
	a) can it be correct? The explanation is that viewing the current as arising from
	an electronic transition is only accurate if the leads are weakly coupled to
	the conductor. As the coupling becomes stronger the lead eigenstates evolve
	into scattering states consisting of an incident wave in lead L, together with
	a scattered wave in lead R (see figure 2.7 a). It carries current from L to
	R as long as it is occupied, without need of any state transition. Therefore
	equation 2.11 is correct in absence of inelastic scattering events which involve
	transition rules. On the contrary, if an inelastic scattering event produce a
	transition between two different energy states, equation 2.11 becomes inad-
	equate, and in general the simple inclusion of a factor [1 -η(E -µ R )] is
	not sufficient to correct the expression and compute a reliable current value.
	11 does not take into State transitions are affected by the exclusion principle in a complicated way,
	account the Pauli exclusion principle. We could instinctively think that the and it is difficult to do justice to this without a detailed microscopic theory.
	current from one lead to another arises from an electronic transition between If the ballistic conduction assumption drops also equation 2.12 loses validity,
	an eigenstate |k, L localized in lead L and an eigenstate |k , R localized in and equation 2.11 cannot be simplified into equation 2.13. This would make
	lead R (see figure 2.7 b). This transition should be forbidden if the final state complicated the linearization of I(V ), since it is hard to imagine a linear
	|k , R is occupied, and equation 2.11 should take into account this fact, for relationship from equation 2.11, while it is straightforward from equation
	example by including a factor [1 -η(E -µ R )]. Equation 2.11 is actually com-2.13, assuming a small voltage. However, how small should it be to assure
	pletely unaware of transition rules governed by the Pauli principle, thus how the linearity of I(V ) and the validity of the Landauer equation? The answer

b) Inelastic scattering conduction: The conduction implies an energy state transition from an initial state |k, L to a different final one |k , R . The current is always limited by the interface transmission probability, but also by the Pauli exclusion principle, which blocks the transition if |k , R is occupied.

  s 2 rl s ll s lr + s rr s rl s ll s lr + s rr s rl

	s 2 rr + s 2 lr	= I	(2.24)

Table 4 .

 4 2 structure and computational details 1: Structural parameters of cubic (spacegroup F m3m), tetragonal (spacegroup P 4 2 /nmc) and monoclinic (spacegroup P 2 1 /c) phases of HfO 2 .

	4.2.1 LDA and GGA calculations
		cubic		tetragonal	monoclinic
		PBE	LDA	PBE	LDA	PBE	LDA
	a ( Å)	5.082	4.994	3.594	3.533	5.145	5.048
	b ( Å)	-	-	-	-	5.206	5.142
	c ( Å)	-	-	5.225	5.076	5.326	5.206
	β ( • )	-	-	-	-	99.63	99.53

Table 4 .

 4 2: Comparison of QP minimum fundamental band gaps (E g ) of monoclinic HfO 2 (Γ → B) and GW corrections (∆E GW

	Ta-

g = E GW g -E DFT g

Table 4 .

 4 [START_REF] Lai | A multilayer stackable thin-film transistor (TFT) NAND-type flash memory[END_REF] shows all the DFT and GW direct band gaps, with the addition of the optical gaps calculated by solving the BSE on top of QSGW . The table also shows the first exciton binding energy, defined as the difference between the direct band gap energy and the optical gap energy, E exc

	b	= E d g -E o g . The

  • 10 -6 (see fig.5.13 and 5.14). We can thus conclude that over values of around 10 -6 the eventual noise due to the limit in the Wannier Hamiltonian accuracy is negligible. On the contrary, figure5.11 (top two graphs) shows that some moderate residual noise arises in the HSE and COHSEX calculations for the 4 MgO structure, which has the thickest oxide and the lowest transmission and conductance. The minima of the 4 MgO HSE and COHSEX spectra lie under 10 -7 , which probably is around the accuracy threshold of the Wannier Hamiltonian. The noise can be removed by applying a cut-off shorter than the MgO thickness (e.g. 13 Å), preserving reasonable reliability.

  /Ag system between DFT and COHSEX or HSE. It indicates that the correction introduced by COHSEX and HSE goes in the right direction to improve the agreement with experiments. This comparison is qualitative since the electrodes considered in the experiment are made of Fe. But it remains physically relevant since, for both Ag electrodes and Fe electrodes with parallel magnetizations, the exponential decay of the conductance with the MgO thickness is related to evanescent states of symmetry ∆ 1 in MgO. × 10-1 3.3 × 10 -2 1.5 × 10 -3 PBE + G 0 W 0 27 5.8 × 10 -1 2.3 × 10 -2 1.1 × 10 -3 COHSEX141.6 × 10 -1 2.6 × 10 -3 3.9 × 10 -5 HSE 13 1.6 × 10 -1 2.7 × 10 -3 4.7 × 10 -5

	Method	1 MgO	2 MgO	3 MgO	4 MgO
	DFT-PBE	31	8.9

Table 5 .

 5 1: Zero-bias tunneling conductance (in Ω -1 .µm -2 ) through a MgO barrier of different thickness, from 1 to 4 cubic cells, with Ag electrodes, as calculated by various approximations. The conductances are evaluated at zero temperature as G 0 • τ (E f ) / S, where τ (E f ) is given by eq. 3.66, and S = (4.14 × 10 -4 ) 2 is the surface of the ab initio supercell in µm2 . Their values at room temperature (eq. 3.67), not given here, do not differ by more than a few percents.

Table 5 .

 5 2: Electronic band gap of bulk MgO, calculated by various techniques. The lattice parameter is equal to 4.212 Å, which is the experimental value at room temperature[START_REF] Wyckoff | Crystal Structure[END_REF]. Values are in eV.

	.15, which shows

  Average transmission function (eq. 3.66) computed by different codes at the PBE level. The agreement turns out to be very good. between the average transmission functions computed by QE, VASP and TranSiesta. The very good agreement between the three codes is proof of the robustness of the approach developed in this work. A first remark is that the transmission spectra are much more irregular and rough than the spectra of the previously studied Ag/MgO junction. The greater irregularity of structure and interface leads also to more irregular transmission functions, especially in the tunneling regime close to the Fermi energy. It is also interesting to notice that, contrary to the Ag/MgO junction, for the Pt/HfO 2 junction the energy at which the conduction band of HfO 2 starts to give a contribution is easily recognizable in figure6.4. Around 2.2 eV the transmission function raises back almost to its values corresponding to the valence band of HfO 2 , that is almost to one. Even though it is difficult to give a precise estimate the valence band of HfO 2 could give a contribution until around -2 eV. This would give a bulk HfO 2 band gap of 4.2 eV, which is higher, but close, to the PBE band gap value (4.01 eV[START_REF] Sklénard | Optical vs electronic gap of hafnia by ab initio bethe-salpeter equation[END_REF]). This kind of analysis remains qualitative since in the considered supercell Pt and HfO 2 states are rearranged and hybridized. Contrary to Ag/MgO, the transmission of electrons with energy corresponding to the conduction band of HfO 2 is almost totally not blocked. This probably occurs because the Pt states have the good symmetry to be transmitted by the evanescent states in the gap of HfO 2 . Nevertheless, the rigorous symmetry analysis of the evanescent and Pt states is not performed here. It is also interesting to remark how a

	Average transmission (E)	10 4 10 3 10 2 10 1 10 0 10 1		Pristine Pt/HfO 2 by PBE		Q. Espresso TranSiesta VASP
		10 5	4	3		2	1 E -E f (eV) 0	1	2	3
	Figure 6.4: 5 4	3	2	1 0 E -E f (eV)	1	2	3	4

  The shape of the HRS transmission function is instead more similar to the pristine state one, but shifted towards lower energies similarly to the LRS one. As expected τ (E f ) raises with respect to the pristine state. The increase is almost one order of magnitude, while the difference between HRS and LRS is about two orders of magnitude. The latter result, even if the supercell considered in this work is a rough model of OxRAM, is in good qualitative agreement with what found in some works on true Pt/HfO 2 /Pt devices, like for instance in the work of Cagli et al.[START_REF] Cagli | Experimental and theoretical study of electrode effects in HfO 2 based RRAM[END_REF]. The same behavior is observed when the transmission functions are com-Transmission function at the HSE level for the pristine state, the HRS, and the LRS.puted at the HSE level (see figure6.7). As expected, all the curves are lowered by around one order of magnitude compared to DFT-PBE. Their shape remains similar to the DFT-PBE for the pristine state and HRS, whereas the LRS shows the opening of a conduction channel between -2 eV and 0 eV which was not clearly visible at the DFT-PBE level. To better appreciate the differences between PBE and HSE table6.1 shows the values of zero-bias conductance at zero temperature (2 methods applied to 3 resistive states), and to figures 6.8, 6.9 and 6.10, which show the comparison between PBE and HSE levels for the three resistive states. A direct comparison between the values of table 6.1 and experimental values, Transport calculation on an OxRAM model: Pt/HfO 2 /Pt like for instance the results of Cagli et al.[START_REF] Cagli | Experimental and theoretical study of electrode effects in HfO 2 based RRAM[END_REF], presents some difficulties.

	6 Figure 6.7: Method 4 10 5 10 4 10 3 10 2 10 1 10 0 10 1 Average transmission (E) Pt-HfO 2 by HSE PBE HSE	2 E E f (eV) Pristine 0 HRS LRS PS HRS 6.50 × 10 -9 4.04 × 10 -8 4.13 × 10 -6 2 LRS 1.29 × 10 -9 1.06 × 10 -8 6.76 × 10 -7	4

Table 6 .

 6 1: Zero-bias tunneling conductance (in Ω -1 ) through a HfO 2 barrier in its pristine state (PS), HRS and LRS, with Pt electrodes, as calculated at PBE and HSE levels. The conductances are evaluated at zero temperature by eq. 3.67. As for Ag/MgO (see section 5.5.2), their values at room temperature do not differ by more than a few percents. To infer the conductance of the pristine state as a function of the lateral supercell surface, the latter is equal to 27.1 Å2 . For the HRS and the LRS this operation would be incorrect.

Table B .

 B 5.13, and reported it in figure B.1. Table B.1 shows the values of 2κa extracted from the CBS compared to the ones extracted from the transmission spectra of the Ag/MgO/Ag junctions at k = 0. The values obtained by the two methods are reasonably consistent. 1: 2κa values extracted from the CBS of bulk MgO, and from the transmission functions of Ag/MgO/Ag junctions of different thickness, at k = 0. The values are computed within the 4 approximations considered in this work.

	Method	DFT G 0 W 0 COHSEX HSE
	CBS	2.87	3.09	3.68	3.85
	Transmission 2.73	2.84	3.78	3.62

Appendix A

Effect of a G 0 W 0 -like correction on a 1D tight-binding model [START_REF] Tomfohr | Complex band structure, decay lengths, and fermi level alignment in simple molecular electronic systems[END_REF].

We consider the 1D, two-bands, tight-binding model considered by Tomfohr and Sankey [START_REF] Tomfohr | Complex band structure, decay lengths, and fermi level alignment in simple molecular electronic systems[END_REF]. It is an infinite 1D periodic chain with alternating ψ A and ψ B orbitals obeying the Bloch theorem, onsite energies E A and E B , and nearest neighbor coupling -t between ψ A and ψ B . The band structure is given by:

, and a is the lattice constant.

It also shows the band structure after application of a "scissor operator" which rigidly shifts the conduction band by +2 eV. Keeping the Bloch wavefunctions unchanged, and just modifying the energies by the scissor operator, we have recomputd the Hamitonian in the tight-binding basis. Figure A.3 shows the decay of the new Hamiltonian matrix elements with the distance between orbitals. We see couplings above 10 -2 eV up to the 10 th neighbor, while the original Hamiltonian was limited to nearest neighbors couplings.

Appendix C

Symmetry mismatch between Ag and MgO Bloch functions at the conduction band edge This appendix details why the conduction band edge cannot be determined from the transmission spectra in figures 5.9 and 5.13. The DFT-PBE band structure of bulk Ag is shown in figure C.1. It is computed either in the primitive cell, or in the cubic cell at k = 0 and with the longitudinal