Je suis tout d'abord reconnaissant envers mon directeur de thèse Alexander Smogunov, un excellent chercheur et un aussi bon professeur qui m'a donné l'opportunité de travailler avec lui et de faire cette thèse. Il a toujours été patient, positif et encourageant à mon égard. Je te remercie du fond du coeur pour tout ce que tu m'as appris aussi bien scientifiquement que dans la vie de tous les jours. Je prends toujours du plaisir à apprendre de toi.

Je remercie aussi chaleureusement Fabienne Michelini, une super enseignantechercheuse avec qui j'ai eu l'honneur de travailler et qui m'a pris en stage de Master 2 et qui m'a soutenue à faire cette thèse.

I want to thank Jelena Sjakste (Ecole polytechnique), Stefano Sanvito (Trinity college of Dublin, Ireland) and Remi Avriller (Université de Bordeaux) to be the examinators of my jury; Thanks to Samuel Poncé (UC Louvain) and Andrea Ferretti (CNR, Itally) to give some of your time to read my Thesis; and finally I would like to thank particularly Jérôme Saint Martin (Université Paris-Saclay) for presiding my jury.

Je tiens aussi à remercier, Yannick, Binh, Mahé, Poonam, lilian et Ian-

Nous concevons une méthodologie générale et un code pour la simulation dans un cadre de mécanique quantique de différents porteurs, en premier lieu des électrons et des phonons, pour des systèmes à grande taille. Différents types d'interactions entre porteurs ont également été envisagés et implémentés dans notre code de transport, comme les interactions électron-phonon ou les interactions entre spin électronique et spin local (résidant sur une molécule magnétique par exemple), responsables d'effets inélastiques.

La méthode est basée sur des calculs réalistes de structure électronique effectués avec le package Quantum-ESPRESSO et le code Wannier90. Le code Wannier90 permet ensuite de construire des hamiltoniens électroniques sur une base localisée d'orbitales Wannier, tandis que le package PHonon calcule la matrice dynamique régissant la dynamique des phonons. Nous avons travaillé selon deux approches : i) les fonctions de Green hors équilibre (NEGF), permettant de calculer les spectres de transmission ainsi que les courants de charge, d'énergie ou de chaleur, et ii) la propagation en temps réel des paquets d'ondes en utilisant soit l'équation de Schrödinger, soit via les polynômes de Chebyshev pour développer l'opérateur d'évolution.

Ces méthodes ont été testées une première fois sur deux modèles 1D de deux chaînes Ag reliées par une molécule de benzène ou de vanadium. Enfin, nous avons appliqué la méthode à plusieurs matériaux 2D tels qu'une seule couche de phosphore noir ou de graphène avec des atomes de Co déposés sur ces matériaux.

Title: Ab initio modeling of quantum transport and dynamics of interacting carriers Keywords: Quantum-ESPRESSO, Quantum transport, spin-flip, NEGF, electronphonon, wave packet

Abstract :

We design a general methodology and a code for simulation in a quantummechanical framework different transport carriers, first of all electrons and phonons, for large scale systems. Different kinds of interactions between carriers have been also considered and implemented in our transport code, such as electron-phonon interactions or exchange coupling between electronic and local (residing on a magnetic molecule, for example) spins, responsible for inelastic effects.

The method is based on realistic electronic structure calculations carried out with Quantum-ESPRESSO package and Wannier90 code. We worked within two approaches: i) non-equilibrium Green's functions (NEGF), allowing to calculate the transmission spectra as well as the charge, energy, or heat currents and ii) real-time propagation of wave packets using either the Schrödinger equation or the Chebyshev polynomials for expanding the evolution operator.

The method was first applied to 1D model case of two Ag chains connected by a benzene and by a vanadocene molecule. Finally, we applied the method to several 2D materials, such as a single layer of black phosphorus or graphene with deposited Co ad-atoms.

Résumé long en français

Ma thèse est consacrée au développement d'un programme qui permettrait de simuler (de manière générale) le transport quantique de différents canaux (électrons et phonons dans un premier temps) dans des jonctions à l'échelle atomique. Une étape importante consistait à prendre en considération diverses interactions locales des porteurs se propageant au sein de la jonction. Par exemple, nous avons mis en oeuvre des interactions électron-phonon ou spin-spin entre des électrons de conduction et des adsorbats locaux tels que des molécules (éventuellement magnétiques). Cela peut ouvrir la possibilité d'aborder de nombreux phénomènes intéressants, par exemple, la conversion d'énergie et divers effets thermoélectriques, ou des problèmes plus fondamentaux tels que l'intrication quantique entre différents degrés de liberté ou les effets de décohérence quantique dus aux interactions inélastiques. Afin de l'appliquer à des matériaux réels, nous utilisons le cadre de la mécanique quantique basé sur des calculs de structure électronique réalistes effectués avec le package Quantum ESPRESSO basé sur la théorie fonctionnelle de la densité. En conjonction avec les codes WANNIER90 et EPW, il est utilisé pour construire l'hamiltonien sur la base des fonctions de Wannier pour les électrons ou calculer la matrice dynamique pour les phonons ainsi que pour évaluer les constantes de couplage électron-phonon dans la représentation spatiale (en utilisant le code EPW). Ce sont les paramètres clés pour les simulations ultérieures de transport à grande échelle avec notre code. L'hamiltonien total conduisant à la propagation d'une porteuse est composé de deux premiers termes décrivant le réseau parfait et une partie centrale, tandis que le troisième terme contient le couplage entre eux. Les deux premiers termes sont calculés séparément. La supercellule décrivant la partie centrale doit contenir suffisamment de cellules primitives du même type que les électrodes dans toutes les directions pour s'approcher du potentiel de type des électrodes sur les frontières. Ainsi, le code des transports considère les bords de la région centrale comme un prolongement des électrodes.

Sans interactions, c'est-à-dire dans le régime élastique, le transport quantique est régi par l'approche de Landauer-Büttiker reliant de nombreuses propriétés intéressantes (telles que le courant de charge, le courant d'énergie ou encore la chauffe des électrodes) à la fonction de transmission de mécanique quantique pour qu'un porteur se propage à travers la jonction. Le calcul de la transmission élastique peut se faire par plusieurs méthodes. Nous avons travaillé selon deux approches : i) La propagation en temps réel des paquets d'ondes en utilisant soit l'équation de Schrödinger, soit via les polynômes de Chebyshev pour développer l'opérateur d'évolution, et ii) Les fonctions de Green hors équilibre (NEGF), permettant de calculer les spectres de transmission ainsi que les courants de charge, d'énergie ou de chaleur.

La méthode de paquet d'onde consiste en une simulation en temps réel de la propagation d'une porteuse à travers la jonction. Il peut évoluer dans le temps en utilisant différentes techniques de calcul telles que la solution directe de l'équation de Schrödinger ou par l'opérateur d'évolution développé par les polynômes de Chebyschev. La probabilité de transmission est donnée par la portion de la partie transmise du paquet d'ondes. Là où l'équation de Schrödinger permet de prendre en compte un système avec un potentiel dépendant du temps, la méthode des polynômes de Chebyschev permet d'économiser du temps de calcul. Différents types d'interactions peuvent être pris en compte de manière relativement simple et établie dans les deux approches, avec certains avantages et inconvénients. Le couplage électron-phonon est considéré uniquement dans la jonction et est exprimé en base d'orbitale de Wannier pour la partie électronique et en base de mode d'oscillation pour les phonons. Nous utilisons un modèle de dynamique de paquet d'onde similaire à celui proposé par Mingo et al [START_REF] Mingo | Calculation of the inelastic scanning tunneling image of acetylene on cu (100)[END_REF], qui consiste à décrire la fonction d'onde totale du système sur différents canaux phononiques. L'interaction spin-spin intervient lorsqu'une impureté magnétique est introduite dans un système métallique non magnétique, dans certaines conditions et à basse température. Cette impureté est filtrée en raison des interactions avec les électrons de conduction formant un état fortement corrélé. Le filtrage provient d'événements de spin-flip, inversant le spin de l'impureté et de l'électron tout en préservant leur spin total.

La méthode de paquet d'onde a été appliquée aux paquets d'ondes électroniques, alors que la propagation des vibrations atomiques était principalement traitée par des paquets d'ondes classiques. Nous proposons une approche permettant de traiter de manière quantique la propagation de phonon dans une base spatial (comme pour les électrons avec la base d'orbitale de Wannier) représentant le déplacement atomique de chaque atome suivant dans les trois directions de l'espace. Pour Obtenir l'Hamiltonien phononique d'un système, on part de la matrice dynamique calculée à l'aide du package Quantum ESPRESSO en base de déplacement atomique. Afin de pouvoir récupérer l'Hamiltonien, il faut changer de repère et passer en base des modes de vibration et conserver la matrice de passage entre les deux bases.

La pleine puissance du formalisme NEGF apparaît lorsqu'on va au-delà du régime de transport cohérent des transporteurs n'interagissant pas. Diverses interactions à plusieurs corps (telles qu'électron-électron, phonon-phonon ou électronphonon), apparaissant dans la région de diffusion, peuvent être commodément incluses (de manière perturbatrice) en ajoutant les auto-énergies d'interaction correspondantes à celles de contact. Cela entraînera une correction non élastique supplémentaire des courants électriques ou phononiques totaux. À l'inverse, nous pouvons prendre en compte l'interaction électron-phonon pour le calcul des phonons. Les corrections anharmoniques dues aux interactions phonon-phonon dans la jonction peuvent également être incluses de manière similaire en construisant les auto-énergies d'interaction. La plupart des calculs ont été effectués à l'aide de diagrammes de Feynman d'ordre le plus bas impliquant des processus à trois phonons (corrections du second ordre). Il a été constaté que ces effets à plusieurs corps jouent un rôle important dans la réduction significative de la conductance thermique des nanojonctions à l'échelle atomique à haute température.

L'effet des processus de spin-flip sur le transport d'électrons a également été étudié avec l'approche NEGF pour les points quantiques ou pour les chaînes atomiques simulées par un modèle magnétique s-d. Étant donné que le formalisme NEGF devient plutôt compliqué et prend plus de temps, la plupart des calculs jusqu'à présent se limitaient aux corrections d'ordre le plus bas en ce qui concerne les constantes de couplage électron-phonon. Il est conçu pour traiter les systèmes à plusieurs électrons fournis par deux mers de Fermi dans les électrodes gauche et droite. On a implémenté un modèle prenant en compte l'intéraction électronphonon où, la jonction moléculaire est couplée à son bain thermique (modes de vibration moléculaire thermalisants) mais est découplée phononiquement des deux électrodes, seules les fonctions de Green phononiques de la jonction sont nécessaires. Il est donc pratique de travailler sur la base de ses modes de vibration comme dans le cas du paquet d'ondes. Par ailleurs, un second modèle, que nous voulions implémenter, consiste à relier phononiquement la jonction et les électrodes. Dans ce cas, ce problème serait mieux traité en base déplacement avec la matrice dynamique de chaque pièce et leurs couplages.

Ces méthodes ont été testées une première fois sur deux modèles 1D de deux chaînes Ag reliées par une molécule de benzène ou de vanadium. Nous avons observé l'influence de la géométrie des molécules sur le transport des électrons avec ou sans interaction électron-phonon et identifié les modes de vibration les plus importants pour la propagation des électrons et des phonons. Les états moléculaires du benzène participants au transport des électrons ont pu être révélés avec la propagation de paquet d'onde. Nous avons constaté que l'orientation de la molécule permet ou non l'accès des électrons par certains états moléculaires, induisant ainsi une différence de transmission en fonction de l'orientation. Dans un second temps, nous avons testé le couplage électron-phonon et la transmission photonique. Des modes de vibrations du benzène, qui ont été activés par la connexion entre les électrodes et la molécule, participent grandement à l'interaction électron-phonon. De plus, ces modes permettent aussi la propagation de paquet d'onde phononique à travers la jonction pour certaines fréquences d'oscillations. Avec le formalisme des fonctions de Green, nous avons aussi calculé différents types de courant. Cela a permis en autre d'étudier des effets thermoélectriques tel que l'effet Seebeck, Peltier ou encore l'effet Joule. En considérant un couplage électron-phonon, des modes de vibrations induisent des sauts de conductance aux voltages correspondants aux fréquences de ces modes.

Dans le deuxième exemple, nous avons considéré les jonctions Ag/vanadocène afin de tester les interactions spin-spin. Sans tenir compte des interactions entre spins, ce système privilégie la propagation d'électron avec une certaine orientation de spin et permet donc un filtrage de spin. En ajoutant l'interaction spin-spin, cela a induit une probabilité significative de spin-flip modifiant le taux de transmission total.

Enfin, nous avons appliqué les deux méthodes aux matériaux 2D tels qu'une mono-couche de phosphore noir ou de graphène avec des atomes de Co déposés sur ces matériaux. En étudiant le phosphore noir, nous avons tout d'abord écrit un article [START_REF] Lahrache | Effect of bonding and antibonding character of electronic states on their tunneling spectra[END_REF]. Dans lequel, nous expliquons l'implication des orbitales moléculaires non liantes dans les mesures de microscope électronique à balayage (STM) sur un solide de phosphore noir. Les orbitales non liantes, étant plus étendues en dehors du matériau, sont plus sujets au recouvrement avec les orbitales de la pointe d'un STM. Par la suite, nous nous sommes intéressés à modéliser le transport sur une mono-couche. L'idée était de tirer profit de la structure "armchair" du phosphore suivant une direction privilégié et de déposer une ligne d'atome de Co le long de cette direction. Nous voulions étudier le transport suivant la direction orthogonale à cette ligne de Co.

L'avantage de ce genre de système 2D, où la jonction est periodique suivant la direction perpendiculaire au transport, est que nous pouvons appliquer les méthodes NEGF et paquet d'onde dessus. Ce pendant, nous nous sommes aussi intéressés à un modèle dans lequel la jonction contient un défaut ponctuel. Sur une mono-couche de graphène, nous testons la propagation d'un paquet d'onde en deux dimensions. Avec une jonction où l'on dépose un atome de Co sur le graphène, on peut ajouter l'interaction spin-spin dans les calculs.

Dans le futur, en plus du rapprochement d'une molécule en contact avec un bain thermal avec la methode NEGF, nous allons mettre en oeuvre une réalisation avec laquelle la jonction est couplée à des conducteurs non seulement électroniquement, mais aussi par vibrations. Les interactions électron-photon peuvent également être élaborées puisque formellement le formalisme est très similaire au cas déjà mis en oeuvre des interactions électron-phonon. Le calcul électron-phonon avec paquet d'ondes a été introduit avec un exemple d'un canal phonon, ce qui fournit la probabilité d'excitation du phonon. Cependant, nous pouvons approfondir ce modèle. On peut analyser la décohérence des électrons ou des phonons en fonction du nombre de canaux que l'on simule. Nous pouvons étudier les trains de paquets d'ondes via différents canaux et observer comment ils sont distribués entre les canaux. Dans le futur, un projet plus ambitieux consistera à simuler simultanément la propagation d'un phonon et d'un électron pour observer comment ils interagissent dans la jonction ou, à l'inverse, modéliser la manière de dissiper un phonon créé par une interaction avec un paquet d'ondes électroniques. Le calcul électron-phonon avec paquet d'ondes a été introduit avec un exemple d'un canal phonon, ce qui fournit la probabilité d'excitation du phonon. Cependant, nous pouvons approfondir ce modèle. On peut analyser la décohérence des électrons ou des phonons en fonction du nombre de canaux que l'on simule. Nous pouvons étudier les trains de paquets d'ondes via différents canaux et observer comment ils sont distribués entre les canaux. Dans le futur, un projet plus ambitieux consistera à simuler simultanément la propagation d'un phonon et d'un électron pour observer comment ils interagissent dans la jonction ou, à l'inverse, modéliser la manière de dissiper un phonon créé par une interaction avec un paquet d'ondes électroniques.

Chapter 1 Introduction

Motivation and state-of-the-art

The quantum transport field studies propagation and dynamics of various carriers, such as electrons, phonons, magnons, etc. in different materials where quantummechanical effects play an important role. They are important to be taken into account in understanding the behaviour and properties of mesoscopic systems and nanoscale devices with characteristic size smaller than the phase coherence length so that quantum-mechanical nature of carriers is important. It has led to several significant advancements in various fields such as field-effect transistors [3] or molecular spintronics [4]. Modern transistors and integrated circuits were created by the development of quantum transport models. Quantum transport domain is at the origin of studies of Hall effect, magnetic tunnel junctions, quantum interference in atomic and molecular junctions [START_REF] Carlos | Molecular Electronics : An Introduction to Theory and Experiment[END_REF], etc. Quantum transport have a crucial role in understanding and manipulating the behaviour of electrons in quantum dots. The quantum dots have a unique electronic and optical properties which makes them useful in applications such as displays, solar cells [START_REF] Barnham | A new approach to high-efficiency multiband-gap solar cells[END_REF][START_REF] Aeberhard | Microscopic nonequilibrium theory of quantum well solar cells[END_REF], medical imaging, etc. One can mention many other domains such as spintronics, quantum computing or quantum optics, for example, where quantum transport and quantum effects have a determinant role in the advancement.

My PhD thesis is devoted to developing a general playground and computational platform which would allow simulating the quantum transport of different channels (electrons and phonons at the first time) in atomic-scale junctions, as illustrated schematically in Figure 1.1. An important ingredient was to take into consideration various local interactions of propagating carriers within the junction. For example, we have implemented electron-phonon or spin-spin interactions between conduction electrons and local adsorbates such as (possibly magnetic) molecules. This can open the possibility to address many interesting phenomena 11 A gradient of temperature or a voltage can be applied between the leads to drive the motion of carriers. such as, for example, energy conversion and various thermoelectric effects, or more fundamental problems such as quantum entanglement between different degrees of freedom or quantum decoherence effects due to inelastic interactions. In order to be applied to real materials, we use the quantum mechanical framework based on realistic electronic structure calculations performed with Quantum-ESPRESSO (QE) package based on the Density Functional Theory (DFT). In conjunction with Wannier90 and EPW codes, it is used to construct the Hamiltonian in the basis of Wannier functions for electrons or calculate the Dynamical matrix for phonons as well as to evaluate electron-phonon coupling constants in the real space representation (using EPW code). These are the key parameters for subsequent large-scale transport simulations with our code.

Without interactions, i.e. in the elastic (or coherent) regime, quantum transport is ruled by Landauer-Büttiker approach relating many transport properties of interest to the quantum-mechanical transmission function for a carrier to propagate across the junction. For electrons, it allows calculating the particle current as well as the energy or the heat currents (from the Left or Right electrode to the

I L/R = e h dE T (E)[f L/R (E)) -f R/L (E)] I E L/R = 1 h dE E T (E)[f L/R (E)) -f R/L (E)] I Q L/R = -Q L/R = 1 h dE (E -µ L/R ) T (E)[f L/R (E)) -f R/L (E)], (1.1)
where T (E) is the total energy-dependent electron transmission function (including two spin polarizations) and f L/R are Fermi-Dirac electron distribution functions of the Left/Right lead kept at the chemical potential µ L/R . The heat power released in electrodes, Q L/R , is opposite to the heat current going from the corresponding electrode into the scattering region.

For phonons, the transmission function, τ (E), allows calculating thermal (energy) currents by:

I E L/R = 1 h dE E τ (E)[n L/R (E)) -n R/L (E)],
where n L/R are now Bose-Einstein distribution functions of two electrodes. The calculation of elastic transmission can be done by several methods [START_REF] Fan | Linear scaling quantum transport methodologies[END_REF]. Figure 1.2 demonstrates some of them.

Scattering approach: It consists of calculating all transmission and reflection amplitudes, t ij and r ij , for all Bloch states propagating from the Left electrode to the junction. These scattering states can be calculated using wave function matching technique by dividing the whole system into the slabs and constructing the global smooth solution.

Such a strategy is implemented in PWcond code [START_REF] Smogunov | Ballistic conductance of magnetic co and ni nanowires with ultrasoft pseudopotentials[END_REF] (a part of QE) for calculating the electron transport directly after the SCF calculation using a plane wave basis.

Non-equilibrium Green's functions (NEGF):

It is a widely-used method [START_REF] Carlos | Molecular Electronics : An Introduction to Theory and Experiment[END_REF][START_REF] Haug | Quantum kinetics in transport and optics of semiconductors[END_REF][START_REF] Pourfath | The non-equilibrium Green's function method for nanoscale device simulation[END_REF] where the elastic transmission is calculated using the Green's function of the central region and contact self-energies, Σ L/R , describing its coupling to Left/Right electrodes. The method was extensively applied to both electron and phonon transport.

Wave packet method: It consists of a real-time simulation of a carrier propagation across the junction [START_REF] Bokes | Stroboscopic wave-packet description of nonequilibrium many-electron problems[END_REF][START_REF] Konôpka | Wave-packet representation of leads for efficient simulations of time-dependent electronic transport[END_REF]. It can be evolved in time using different calculational techniques such as direct solution of the Schrödinger equation, by evolution operator [START_REF] Tae | Unitary quantum time evolution by iterative lanczos reduction[END_REF][START_REF] Jad | Chebyshev matrix product state approach for time evolution[END_REF] expanded by Chebyshev polynomials, by Suzuki-Trotter time evolution, which decomposes states in Trotter subspace, or by Krylov time evolution [START_REF] Hochbruck | On krylov subspace approximations to the matrix exponential operator[END_REF][START_REF] John Daley | Time-dependent density-matrix renormalization-group using adaptive effective hilbert spaces[END_REF]. The transmission probability is given by the portion of the transmitted part of the wave packet. The approach was applied to electronic wave packets, while propagation of atomic vibrations was mainly treated by classical wave packets.

Different types of interactions can be accounted for in a relatively straightforward and established way into the last two approaches, with certain advantages and drawbacks:

NEGF, interactions:

The full power of the NEGF formalism appears when going, however, beyond the coherent transport regime of noninteracting carriers. Various many-body interactions (such as electron-electron, phonon-phonon or electron-phonon), appearing in the scattering region, can be conveniently included (in a perturbative way) by adding the corresponding interaction self-energies to the contact ones [START_REF] Sergueev | Efficiency of thermoelectric energy conversion in biphenyl-dithiol junctions: Effect of electron-phonon interactions[END_REF][START_REF] Frederiksen | Inelastic transport theory from first principles: Methodology and application to nanoscale devices[END_REF][START_REF] Kim | First-principles quantum transport with electron-vibration interactions: A maximally localized wannier functions approach[END_REF]. That will result in an additional non-elastic corrections to the total charge currents. NEGF can also be used for phonon transport [START_REF] Galperin | Heat conduction in molecular transport junctions[END_REF] with electron-phonon interactions. Anharmonic corrections due to phonon-phonon interactions in the junction can be also included in a similar way by constructing the interaction self-energies. Here, most calculations were performed using the lowest order Feynman diagrams involving three-phonon processes (second-order corrections) [START_REF] Mingo | Anharmonic phonon flow through molecular-sized junctions[END_REF][START_REF] Wang | Nonequilibrium green's function approach to mesoscopic thermal transport[END_REF]. These manybody effects were found to play an important role in reducing significantly the thermal conductance of atomic-scale nanojunctions at high temperatures. Effects of spin-flip processes on electron transport were also investigated with NEGF approach for quantum dots [START_REF] Chen | Spin-flip effects on andreev reflection process[END_REF][START_REF] Fabricio | Spin-polarized current and shot noise in the presence of spin flip in a quantum dot via nonequilibrium green's functions[END_REF] or for atomic chains [START_REF] Hurley | Spin-flip inelastic electron tunneling spectroscopy in atomic chains[END_REF] simulated by a magnetic s-d model. Since the NEGF formalism becomes rather complicated and timeconsuming, most calculations so far were restricted to the lowest order corrections with respect to electron-phonon coupling constants. Advantages: It is designed to deal with many-electron systems provided by two Fermi seas in the Left and Right electrodes.

Drawbacks: Interactions are treated perturbatively so that calculations are normally restricted to the lowest order corrections.

Wave packets, interactions: A time-dependent approach based on wave packets with electron-phonon interactions was presented by Monturet and Lorente [START_REF] Monturet | Inelastic effects in electron transport studied with wave packet propagation[END_REF] to model tight-binding chains with few sites connected to vibrations. The spin-transfer phenomena were also simulated based on wave packets and density matrix approach [START_REF] Kim | How many electrons are needed to flip a local spin[END_REF][START_REF] Kim | Quantum mechanics of spin transfer in coupled electron-spin chains[END_REF] in order to validate a semi-classical approximation. Advantages: Non-perturbative solution; relatively easy to include time-dependent perturbations and some other ingredients, such as magnetic field, for example. Drawbacks: Difficult to simulate more than few carriers at the same time because of entanglement of the full wave function.

Our code aims at implementing both NEGF and wave packet's methods and is going to compliment plane wave transport code PWcond, already integrated in Quantum-ESPRESSO. In the NEGF, we include electron-phonon interactions in the junction, similarly to [START_REF] Sergueev | Efficiency of thermoelectric energy conversion in biphenyl-dithiol junctions: Effect of electron-phonon interactions[END_REF][START_REF] Frederiksen | Inelastic transport theory from first principles: Methodology and application to nanoscale devices[END_REF][START_REF] Kim | First-principles quantum transport with electron-vibration interactions: A maximally localized wannier functions approach[END_REF], with the goal to study their effects on both electron and phonon transport at the same time in a self-consistent way. Within the wave packet approach, we implement a real-time quantum dynamics and propagation of electron and phonon wave packets for 1D (molecular junctions) and also for 2D (Graphene and others) systems. Phonon wave packets, in full analogy with electrons, are constructed as a specific quantum superposition of single-phonon states corresponding to different normal modes. We want to include furthermore the effect of various interactions on the wave packet's propagation. This has been realized, firstly, for an electron interacting with local vibrations on the molecule or with a local spin which follows closely to [START_REF] Monturet | Inelastic effects in electron transport studied with wave packet propagation[END_REF][START_REF] Kim | How many electrons are needed to flip a local spin[END_REF]. Similar interactions can also be considered for a phonon wave packet in the future. Note that in perspective also other kinds of interactions, like electron(or phonon)/photon, may be also introduced in our code.

There are many DFT-based codes with different implementations of quantum transport. Here some of them are listed:

The Quantum-ESPRESSO package [START_REF] Giannozzi | Quantum espresso: a modular and open-source software project for quantum simulations of materials[END_REF][START_REF] Giannozzi | Advanced capabilities for materials modelling with quantum espresso[END_REF] proposes transport calculations with the PWcond code. It is based on plane waves and solves the 3D scattering problem, calculating all transmission and reflection coefficients.

The Dacapo code [START_REF] Sune | An object-oriented scripting interface to a legacy electronic structure code[END_REF] is based on the plane-wave pseudopotential approach in DFT calculations. It allows calculating various transport quantities with NEGF methods such as conductance, current, transmission spectra. From this, Ref. [START_REF] Hjorth Larsen | The atomic simulation environment-a python library for working with atoms[END_REF] develops ASE -the oriented object code in python implementing the DFT.

The code SIESTA [START_REF] García | Siesta: Recent developments and applications[END_REF] is a famous code which performs electronic structure calculations based on DFT. It uses a localized basis set to describe the electronic structure of the system. Transport is calculated with TranSIESTA [START_REF] Brandbyge | Density-functional method for nonequilibrium electron transport[END_REF][START_REF] Papior | Improvements on non-equilibrium and transport green function techniques: The next-generation transiesta[END_REF]) which uses the NEGF method.

Quantum ATK (Atomistix ToolKit) is a commercial software package developed by Synopsis [START_REF] Schneider | Atkforcefield: a new generation molecular dynamics software package[END_REF]. It does DFT calculations employing a localized basis set and more precisely the numerical atomic orbital (NAO) basis set. It performs NEGF calculations [START_REF] Ghorbani-Asl | A computational study of the quantum transport properties of a cu-cnt composite[END_REF]. The electron and electron-phonon calculations are done with DFT, phonon transport calculations are done with Boltzmann equation [START_REF] Gunst | Firstprinciples method for electron-phonon coupling and electron mobility: Applications to two-dimensional materials[END_REF].

Kwant [START_REF] Groth | Kwant: a software package for quantum transport[END_REF][START_REF] Kloss | Tkwant: a software package for time-dependent quantum transport[END_REF] is a Python package for numerical calculations of quantum transport properties using a tight-binding model. It is based on the Landauer-Buttiker formalism and wave packet propagation in multi-terminal setup. It is important to note that Kwant focuses specifically on the numerical simulation of quantum transport phenomena and does not provide electronic structure calculation capabilities like DFT codes.

OpenMX is a DFT code working in plane wave basis set. It can perform several calculations on electronic structure, including current calculations with NEGF method.

KITE is an open-source software coded in python. It uses machine learning to use DFT in tight-binding model without solving Konh-Sham equations [START_REF] Simão | Kite: high-performance accurate modelling of electronic structure and response functions of large molecules, disordered crystals and heterostructures[END_REF]. It calculates wave packet propagation using Chebyshev polynomials.

KNIT [START_REF] Kazymyrenko | Knitting algorithm for calculating green functions in quantum systems[END_REF] is a python library which calculates the transport for multi-electrode systems with NEGF method using tight-binding models.

GOLLUM and SMEAGOL [START_REF] Ferrer | Gollum: a next-generation simulation tool for electron, thermal and spin transport[END_REF][START_REF] Sadeghi | Oligoyne molecular junctions for efficient room temperature thermoelectric power generation[END_REF][START_REF] Sadeghi | Theory of electron, phonon and spin transport in nanoscale quantum devices[END_REF] are transport codes based on DFT packages such as Castep, VASP, ABINIT and Quantum-ESPRESSO interfaced with Wannier90 for the first one, and with SIESTA for the second one. They are both using Hamiltonians in user-defined tight-binding model or in localized basis provided by DFT calculations. Historically, SMEAGOL was created first to perform electron transport with the NEGF method. GOLLUM can calculate the electron and phonon transport in multi-terminal geometries.

TB_Sim is a tight-binding code to compute the structural, electronic, optical and transport properties of various nanostructures using Kubo or NEGF formalisms for transport.

OMEN is a transport code based on semi-empirical tight-binding models. Carrier and current densities are calculated by injecting into the system electrons and holes at different energies and solving their evolutions with the wave function or the NEGF formalisms.

Density Functional Theory

Hohenberg-Kohn theorems and Kohn-Sham equations

To describe the electronic structure of the system, one can begin with the Schrödinger equation:

Ĥsys Ψ sys = E sys Ψ sys .

(1.

2)

The total many-body Hamiltonian for N electrons and M nuclei is the sum of kinetic energy terms and electron-nuclei, electron-electron, and nuclei-nuclei interactions:

Ĥsys = - N i=1 ℏ 2 ∇ 2 i 2m - M A=1 ℏ 2 ∇ 2 A 2m A - N i=1 M A=1 e 2 Z A r iA + 1 2 N i̸ =j e 2 r ij + 1 2 M A̸ =B e 2 Z A Z B R AB (1.3)
As nuclei are much heavier (and so slower) than electrons, one usually applies the adiabatic (or Born-Oppenheimer) approximation, which decouple the motion of electrons and nuclei. One then arrives at the electronic Hamiltonian where nuclei positions enter as external parameters:

Ĥ = - N i=1 ℏ 2 ∇ 2 i 2m - N i=1 M A=1 e 2 Z A r iA + 1 2 N i̸ =j e 2 r ij = T + V + Ŵ (1.4)
containing kinetic T , external single-body V (due to electron-nuclei interactions) and electron-electron Ŵ contributions and the electronic Schrödinger equation:

ĤΨ = EΨ. (1.5)
The solution of the two last equations gives the electronic wave function Ψ and the electronic energy E. The total energy is then the sum of E and the constant energy from nuclei-nuclei potential, since the latters are treated classically. The Density Functional Theory (DFT) is a computational quantum mechanical method used to calculate the electronic structure of many-electron systems and provides a practical approach to study the properties and behaviour of atoms, molecules and solids.

The main variable in DFT is the electron density which can be defined for any many-electron wave function Ψ as:

n(r) = ⟨Ψ| n(r) |Ψ⟩ , n(r) = ψ † (r)ψ(r) (1.6)
using second quantization representation and field operators, ψ † (r) and ψ(r), creating or annihilating an electron at r, respectively.

The theoretical footing of DFT is based on two theorems provided by Hohenberg and Kohn in 1964 [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF][START_REF] Engel | Density functional theory[END_REF]. The first Hohenberg-Kohn theorem states that there are one-to-one mapping between external potential V , the ground state electron wave function Ψ, and, further, the ground state density n(r):

V (r) ⇔ Ψ(r 1 , r 2 , ..., r N ) ⇔ n(r).
(1.7)

Therefore, all observables in the ground state are unique functional of the electron density n(r), for example the wave function Ψ[n] itself, which is much simpler since it depends on only one variable r.

The second Hohenberg-Kohn theorem states that the exact ground state density n 0 (r) corresponding to the external potential V (r) minimizes the energy functional:

∂E V [n(r)] ∂[n(r)] = 0, (1.8)
with

E V [n] = ⟨Ψ[n]| T + Ŵ |Ψ[n]⟩ + V (r)n(r)dr (1.9)
so that the ground state total energy is given by

E 0 = E V [n 0 ].
In 1965, Kohn and Sham [START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF] proposed an efficient way to solve for n(r) mapping a "real" interacting electrons problem to the one of non-interacting electrons. They suggested to rewrite the variational functional E V [n] in the following way:

E V [n] = T s [n] + V (r)n(r)dr + E H + E xc [n].
(1.10)

Here T s [n] is the kinetic energy of non-interacting electrons,

T s [n] = ⟨Ψ s [n]| T |Ψ s [n]⟩, where Ψ s [n]
is the inverse of the last map in Eq. 1.7 for non-interacting electrons and the Hartree energy correspond at

E H = e 2 2 n(r)n(r ′ ) |r -r ′ | drdr ′ .
All remaining corrections are collected in the exchange-correlation functional,

E xc [n],
which is unknown and needs to be approximated, in general.

The variational problem for so written functional is then equivalent to solution of single-particle effective equations:

- ℏ 2 2m ∇ 2 + V eff (r) ϕ i (r) = ϵ i ϕ i (r), V ef f (r) = V (r) + e 2 n(r ′ ) |r -r ′ | dr ′ + V xc (r), (1.11)
where ϕ i are the so-called Kohn-Sham (KS) orbitals of the non-interacting system and V xc (r) = δE xc [n(r)]/δn(r) is the exchange-correlation potential. The ground state electron density is then given by:

n 0 (r) = occ. i |ϕ i (r)| 2 (1.12)
so that the KS equations have to be solved self-consistently. It is very important to realize that the KS scheme would produce the exact ground state density n(r) and, as a consequence, the total energy if one knew the exact exchange-correlation functional.

Comparing two different forms Eq. 1.9 and Eq. 1.10 of the same functional, one can see that the exchange-correlation energy corrects for both kinetic and electron-electron terms:

E xc [n] = T [n] -T s [n] + W [n] - e 2 2 n(r)n(r ′ ) |r -r ′ | drdr ′ . (1.13)
In DFT, the exchange-correlation term is always approximated by different methods. We will discuss later in this chapter several methods of approximation.

Spin polarized systems

The DFT can also be extended to include external magnetic field B(r). The first Hohenberg-Kohn theorem is generalized and establishes one-to-one mappings:

{V (r), B(r)} ⇔ Ψ(r 1 α 1 , r 2 α 2 , ..., r N α N ) ⇔ {n(r), m(r)} (1.14)
where spin indexes α i for the electron i were also introduced. The additional variable, the magnetization density m(r), is calculated as:

m(r) = ⟨Ψ| m(r) |Ψ⟩ , m(r) = α,β=↓,↑ ψ † (r, α)σ α,β ψ(r, β) (1.15)
with Pauli matrices:

σ x = 0 1 1 0 , σ y = 0 -i i 0 , σ z = 1 0 0 -1 . (1.16)
The Kohn-Sham scheme proceeds by introducing the total energy variational functional E V,B [n, m] as follows:

E V,B [n, m] = T s [n, m] + V (r)n(r)dr -B(r)m(r)dr+ e 2 2 n(r)n(r ′ ) |r -r ′ | drdr ′ + E xc [n, m], (1.17)
where exchange-correlation energy functional E xc depends now on both electron and magnetization densities. Konh-Sham equations are written for two-component spinor wave functions and state:

β - ℏ 2 2m ∇ 2 + V eff (r) δ αβ -µ B σ αβ B ef f (r) ϕ i (r, β) = ϵ i ϕ i (r, α) V ef f (r) =V (r) + e 2 n(r ′ ) |r -r ′ | dr ′ + V xc (r) B ef f (r) =B(r) -B xc (r) (1.18)
and exchange-correlation potential and magnetic field are given by V xc (r) = δE xc /δn(r) and B xc (r) = δE xc /δm(r), respectively. The ground state electron and magnetization densities are calculated as:

n 0 (r) = occ. i α |ϕ i (r, α)| 2 , m 0 (r) = occ. i ϕ * i (r, α)σ αβ ϕ i (r, β) (1.19)
where magnetic moment is expressed in units of magneton Bohr, µ B = eℏ/2mc. Note finally that instead of using n(r) and m(r) one can work with another set of four variables, ñαβ (r), which compose the so-called density matrix defined as:

ñαβ (r) = ⟨Ψ| ψ † (r, α)ψ(r, β) |Ψ⟩ (1.20)
and can be calculated in terms of Kohn-Sham orbitals as:

ñαβ (r) = occ. i ϕ * i (r, α)ϕ i (r, β). (1.21)
The relation to the charge and spin magnetizations are provided by:

ñ(r) = Tr[ñ(r)], m(r) = Tr[ñ(r)σ] (1.22)
and vice versa:

ñαβ (r) = 1 2 [n(r)δ αβ + σ αβ • m(r)] (1.23)

Exchange-correlation Energy Functional

Konh-Sham equations describe an auxiliary non-interacting system for a manyelectron problem. It is supposed to describe properly the ground state energy and the charge density if the exchange-correlation (xc) energy is exactly known, which is not true in practice. We will discuss in this section main approximations which are used to calculate it. The first one, the most simple, is the Local Density Approximation (LDA), which assumes that the total xc energy can be expressed as an integral of independent local contributions, ϵ xc , in the form:

E LDA xc = dr n(r)ϵ xc [n(r)], (1.24)
where the local xc energy, ϵ xc (n), can be taken as the one for the infinite homogeneous electron gas with the same electron density n everywhere in space. Since the LDA is the local approximation, it may not be very appropriate for describing systems with rapidly changing electron densities or with important long range nonlocal correlations. Despite this limitation, LDA remains a valuable tool in DFT, and it is often used as a starting point for more complex exchange-correlation functional which take into account other effects beyond the local density. The Local Spin Density Approximation (LSDA) is similar to LDA and is used in spin-polarized calculations. It takes into account both electron density and the local magnetization:

E LSDA xc = dr n(r) ϵ xc [n(r), |m|]. (1.25)
One of the most popular parameterizations was proposed by Perdew and Zunger [START_REF] John | Self-interaction correction to densityfunctional approximations for many-electron systems[END_REF], it consists in interpolating the accurate intermediate values obtained from the quantum Monte Carlo data of Ceperley and Alder [START_REF] David | Ground state of the electron gas by a stochastic method[END_REF]. Generalized Gradient Approximation (GGA) includes the gradient of the electron density in addition to the electron density itself in the integral:

E GGA xc = dr n(r) ϵ xc [n(r), |m|, ∇n(r), ∇m z (r)]. (1.26)
This approximation contains more information than local information due to its gradients ∇n and ∇m. PW91 [START_REF] John | Accurate and simple analytic representation of the electron-gas correlation energy[END_REF] and PBE [START_REF] John P Perdew | Generalized gradient approximation made simple[END_REF] are the most widely used parametrizations of GGA which can be reliably used over a very wide range of materials. Typically (but not always), the GGA is more accurate than LDA. It is believed to describe better, in particular, magnetic systems containing transition metal elements like Co, Ni, or Fe. On the other side, the structural and electronic properties of heavier elements, such as Au, Pt, or Ir are known to be treated better within LDA approximation. We note finally that both (semi)local LDA and GGA approximations are poorly applied to systems where non-local effects, strong correlations or self-interaction corrections may play an important role, such as different types of oxides, molecular systems containing magnetic centers and many others. For such systems more sophisticated (and time-consuming) functionals, such as META-GGA, DFT+U, or different types of orbital-dependent and hybrid functionals have been proposed in a variety of flavours [START_REF] Engel | Density functional theory[END_REF].

Quantum ESPRESSO

Quantum-ESPRESSO (QE) is an open-source software package for electronic structure calculations and materials modelling based on DFT [START_REF] Giannozzi | Quantum espresso: a modular and open-source software project for quantum simulations of materials[END_REF][START_REF] Giannozzi | Advanced capabilities for materials modelling with quantum espresso[END_REF]. It is widely used by the scientific community for predicting the electronic, structural, and dynamical properties of different materials. QE provides a suite of tools for performing various calculations, including ground-state calculations, excited-state calculations, molecular dynamics simulations, and more. It can handle a variety of materials, including solids, surfaces, interfaces, and molecules, and it supports a range of exchange-correlation functionals such LDA, GGA, meta-GGA, and various hybrid functionals. The package includes several codes, such as PWscf (Plane Wave Self-Consistent Field), which is the main program for performing DFT calculations. QE also provides a user-friendly interface, called PWgui, for setting up and running calculations and many post-processing codes such as pp.x (to plot different physical quantities such as the charge density or wave functions) or projwfc.x (to calculate total, projected, and vacuum density of states). It also implements different levels of parallelization on both shared-memory and distributed-memory architectures.

Package PWscf (pw.x)

PWscf uses a plane-wave basis set to expand wave functions of electrons in periodic systems, and uses pseudopotentials (PPs) to describe electron-ion interactions. It employs norm-conserving (NC) or ultrasoft (US) PPs enabling to remove core electrons from calculations and to replace the real valence wave functions by so-called pseudo wave functions. The latter are much smoother and could be expanded using smaller number of plane waves, which improves the computational efficiency of calculations. The PPs can be split into a local part, V loc , and a non-local part, V N L , which can be written by introducing so-called projector functions of well-defined angular momentum l centered at an atomic site i:

V N L = i l,m D i l,m β i l Y i l,m β i l Y i l,m , (1.27)
Figure 1.3: Self-consistent scheme of DFT calculations.

where β i l , Y i l,m , D i l,m are radial components of projector functions, spherical harmonics and pseudopotential coefficients, respectively.

In PWscf, the electronic structure calculations are done using a self-consistent field algorithm, where electronic wave functions and electron density are iteratively updated until a convergence threshold is reached, as illustrated in Figure 1.3. First, we give an initial guess for electron charge density, which allows to construct the effective potential. From this, it solves the Kohn-Sham equations to obtain eigenvalues and eigenfunctions. A new electron density is recalculated from these eigenfunctions using Eq. 1.12. At each iteration, it is verified if the difference between the previous and the new electron density is less than the convergence threshold fixed by the user. If convergence is not achieved, it recalculates again the effective potential with the new electron density, otherwise the system is considered being well described with the calculated charge density and electron levels. The convergence criteria for the self-consistent field iterations can be adjusted by the user.

PWscf can perform a variety of calculations, including total energy calculations, geometry optimizations, electronic spectra, and more. It can also calculate the electronic and magnetic properties of materials (with collinear or noncollinear order), such as band structures, density of states, magnetic moments, and magnetic anisotropy energies related to spin-orbit interactions.

Package PHonon (ph.x)

PHonon is a code included in Quantum-ESPRESSO that calculates the vibration properties of a variety of structures, including crystalline systems as well as isolated molecules. The theory behind PHonon is based on density functional perturbation theory (DFPT) [START_REF] Gonze | Dynamical matrices, born effective charges, dielectric permittivity tensors, and interatomic force constants from densityfunctional perturbation theory[END_REF][START_REF] Baroni | Phonons and related crystal properties from density-functional perturbation theory[END_REF][START_REF] Dal | Density-functional perturbation theory with ultrasoft pseudopotentials[END_REF] which can compute vibrational modes for any wave vector using only self-consistent calculations in a primitive unit cell (contrary to the frozen phonon method employing a supercell geometry). It uses the harmonic approximation, where the atoms oscillate around their equilibrium positions with small amplitudes. Thus, the total energy of the ion-electron system in the adiabatic approximation (Eq. 1.4), known as the Born-Oppenheimer energy surface, is expanded in a Taylor series around equilibrium positions with atomic displacements, x iαR : .28) where i, α, and R stand for the atom in the unit cell, Cartesian direction, and the Bravais lattice vector, respectively. The first derivatives of the total energy correspond to the linear response of the system to external perturbation. They define atomic forces, which can be calculated using the Hellman-Feynman theorem as:

E = E(0) + iαR ∂E x iαR x iαR + 1 2 iαR i ′ α ′ R ′ ∂ 2 E ∂x iαR ∂x i ′ α ′ R ′ x iαR x i ′ α ′ R ′ + ... ( 1
F iαR = - ∂E ∂x iαR = -⟨Ψ| ∂H KS ∂x iαR |Ψ⟩ , (1.29)
where H KS is the Kohn-Sham Hamiltonian and Ψ is the ground state wave function. In the absence of external perturbations, the forces are zero at equilibrium. Second-order terms define the Dynamical matrix given by:

D iαR,i ′ α ′ R ′ = 1 √ m i m i ′ ∂ 2 E ∂x iαR ∂x i ′ α ′ R ′ , (1.30)
where m i are atomic masses. The eigenvalues and eigenvectors of the Dynamical matrix provide phonon frequencies and polarization vectors of the corresponding mode, respectively. In order to calculate the dynamical matrix, ph.x performs SCF calculations for a set of all inequivalent atomic displacements in the primitive cell which is much more time-consuming than the usual SCF calculation (roughly by a factor of 3N, where N is the number of atoms in the unit cell). The Dynamical matrix then can be calculated in the real basis using q2r.x code, and the phonon bands can be plotted on a specified k-path using matdyn.x routine.

Advanced features of the PHonon package include the calculation of third-order energy derivatives, Raman parameters, phonon-phonon interaction coefficients etc. The PHonon package also allows to calculate electron-phonon matrix elements in Konh-Sham basis for electrons and mode basis for phonons. They can be used for evaluating electron-phonon (Eliashberg) spectral function for estimation of superconducting critical temperature. PHonon contains also some capabilities to calculate electron-phonon coefficients in the Wannier basis but we have decided to use for that purpose the EPW code (described below) implementing a calculation of a variety of properties related to electron-phonon interactions.

Package Wannier90

Wannier90 is a software package for calculating maximally localized Wannier functions (MLWFs) [START_REF] Arash | wannier90: A tool for obtaining maximallylocalised wannier functions[END_REF][START_REF] Marzari | Maximally localized wannier functions: Theory and applications[END_REF][START_REF] Pizzi | Wannier90 as a community code: new features and applications[END_REF]] using band structure calculations from plane wave codes such as, for example, Quantum-ESPRESSO. MLWFs provide a real-space representation of the electronic structure, and they are often used to study and analyse the properties of materials, such as the electronic band structure, topological properties, response functions, etc.

Wannier functions have been introduced by Gregory Wannier in 1937 as a way to describe the electronic structure of a periodic crystal in a real space. They are defined by a unitary transformation of the Bloch wave functions:

ψ nk (r) = u nk (r)e ik•r (1.31)
which are eigenstates (associated with band n and wave vector k) of the periodic crystal Hamiltonian. Functions u nk are translationally periodic:

u nk (r + R) = u nk (r)
with R being Bravais lattice vectors. The idea behind Wannier functions is to provide a more intuitive and physical representation of the electronic structure than Bloch wave functions with the inverse transformation:

W nR (r) = V (2π) 3 BZ ψ nk (r)e -ik•R dk, (1.32)
where V is the main volume for normalization of wave functions. It returns wave functions localized at the lattice vectors R. This procedure works relatively well for the case of isolated bands producing one-to-one correspondence. For more complex cases, we have to use a more general transformation matrix U to pass from bands m to orbitals n basis:

W nR (r) = V (2π) 3 BZ m U nm (k)ψ mk (r) e -ik•R dk (1.33)
The transformation is designed to produce a set of functions localized as much as possible in real space, meaning with a well-defined position and sharply peaked around a particular atomic site. For this, they are searched in such a way to reduce the spread Ω defined as the sum of the second moments of all the Wannier functions in a reference cell [START_REF] Marzari | Maximally localized wannier functions: Theory and applications[END_REF]:

Ω = n [⟨0n| r 2 |0n⟩ -⟨0n| r |0n⟩ 2 ], (1.34)
where |0n⟩ are Wannier functions in the cell "0". Wannier functions have a wide range of applications in the study of materials, including the calculation of transport properties, the analysis of topological insulators, and the study of polarization and other response functions. They are also used in the construction of tight-binding models, which provide a simplified description of the electronic structure of a material in terms of a few parameters.

The Wannier90 code works by performing a unitary transformation of the Bloch eigenstates to obtain a set of MLWFs that are maximally localized in real space. This transformation is done using the following steps:

-first, generate a set of initial Wannier functions. These can be chosen arbitrarily, but it is often convenient to use the atomic orbitals of constituent atoms.

-secondly, project the Bloch eigenstates onto the initial Wannier functions to obtain a set of projections.

-thirdly, use the projections to construct a unitary transformation matrix that minimizes the spread of the MLWFs. The spread is a measure of the extent of the Wannier functions in real space, and the minimization procedure ensures that the MLWFs are the most localized.

-finally, apply the unitary transformation matrix to the Bloch eigenstates to obtain the MLWFs.

To use the Wannier90 code, you typically need to provide the code with a set of Bloch eigenstates, calculated with an NSCF run and with usually more bands than in SCF calculation, and a set of crystal coordinates that define the periodicity of the system. The code also requires various input parameters that control the calculation, such as the number of MLWFs to be generated and the convergence criteria for the minimization procedure. The resulting Wannier functions can be used to calculate a wide range of properties, such as charge density, band structures, density of states and many others. The Wannier90 code also includes tools for visualizing the MLWFs (with the software xcrysden in the format "xsf" for example) and calculating various transport properties.

Package EPW (epw.x)

Electron-phonon coupling appears when the self-consistent DFT potential is expanded around equilibrium atomic positions:

H KS (X) = H 0 KS + iαR ∂V SCF ∂X iαR x=0 X iαR (1.35)
where the first term, H 0 KS = -(ℏ 2 /2m)∇ 2 + V SCF (X = 0), is the Kohn-Sham Hamiltonian at the equilibrium and X iαR are (reduced, see Eqs. 2.47) atomic displacements of the atom i at the lattice vector R in the Cartesian direction α. In the second quantization, the displacements become operators (as we will discuss later) given by Eq. 2.73:

XiαR = 1 √ N qν ℏ 2ω qν V qν iα e iqR bqν + b † -qν , (1.36)
where bqν and b † -qν are phonon's annihilation and creation operators of the qν vibrational mode described by its polarization vector V qν iα . Inserting this expression, the electron-phonon interaction term in the second quantized form can be cast into the form:

Ĥe-ph = kmn,qν

g mn,ν (k, q)c † k+q,m c k,n bqν + b † -qν , (1.37)
where c † k,n and c k,n are creation and annihilation operators for an electron of the wave vector k and the band n.

Electron-phonon coupling constants, g mn,ν (k, q), can be written using derivatives of the KS potential along a given vibrational mode, ∂ qν V SCF , as follows:

g mn,ν (k, q) = ℏ 2ω qν ⟨ψ m,k+q | ∂ qν V SCF |ψ n,k ⟩ . (1.38)
with ∂ qν V SCF related to derivatives in atomic basis by:

∂ qν V SCF = 1 √ N iαR V qν iα e iqR ∂ X iαR V SCF (1.39)
which follows from Eq. 1.36.

In order to evaluate accurately various material properties related to electronphonon interactions such as superconducting temperatures, scattering rates, electron or phonon lifetimes, and so on, a very fine mesh in reciprocal space is often required.

The EPW (Electron-Phonon Interactions using Wannier functions) [START_REF] Noffsinger | Epw: A program for calculating the electron-phonon coupling using maximally localized wannier functions[END_REF][START_REF] Poncé | Epw: Electron-phonon coupling, transport and superconducting properties using maximally localized wannier functions[END_REF][START_REF] Lee | Electron-phonon physics from first principles using the epw code[END_REF] is designed to accomplish this task. It makes use of Quantum-ESPRESSO package to calculate the electronic structure and phonons (with ph.x code) and Wannier90 to construct maximally localized Wannier functions (MLWFs). The strategy of EPW is then to construct electron-phonon coupling elements in a real space representation for both electrons (in the basis of Wannier functions) and phonons (in the basis of atomic displacements):

g mn,iα (R ′ , R) = ⟨m0| ∂ X iαR V SCF |nR ′ ⟩ (1.40)
where |nR ′ ⟩ are electronic Wannier function n in the cell R ′ . First, coupling matrices in reciprocal space, g(k, q) (Eq. 1.38), are calculated on a coarse grid of k and q from the data provided by pw.x and ph.x (providing also the derivatives of the SCF potential) runs. Then real space coupling matrices, g(R ′ , R) (Eq. 1.40), are constructed from reciprocal space ones by Fourier transformation [START_REF] Giustino | Electron-phonon interaction using wannier functions[END_REF][START_REF] Poncé | First-principles predictions of hall and drift mobilities in semiconductors[END_REF]. In the second step, the inverse Fourier transformation is applied back to real space matrices in order to interpolate g(k, q) on a much denser mesh of k and q vectors.

For our transport code, we need to save real space couplings g mn,iα (R ′ , R) calculated by EPW at the intermediate stage. These parameters, written in appropriate format, will be next used in our transport calculations.

INTRODUCTION TO OUR TRANSPORT CODE AND MY CONTRIBUTION29

Introduction to our transport code and my contribution

The choice of computational path through different codes presented in previous sections depends on the type of system (magnetic or not) and the carrier transport we want to study. As it is shown in Fig. 1.4, it starts with SCF calculation followed by either NSCF (for electrons) or PHonon (for phonons) runs. Our model is based on a system composed of a periodic lattice (lead), large enough to propagate a carrier, and a junction (a molecule, for example) introduced in the centre of the system. The total Hamiltonian driving the propagation of a carrier is then given by:

Ĥ = Ĥlead + Ĥcenter + Ĥlead-center (1.41)
with two first terms describing the perfect lattice and the central region, while the third term contains the coupling between them. Two first terms can be calculated separately. The supercell describing the central part should contain enough of the lead-like units cells in all directions to approach the lead-like potential on the boundaries. Thus, the transport code considers the central region edges as a prolongation of the lead.

The calculations with electron's Hamiltonian and Dynamical matrix (governing the phonon dynamics) are already implemented. The electron-phonon couplings can be extracted from EPW calculations, and exchange spin-spin interactions in the case of magnetic junctions can be calculated from spin-dependent DFT Hamiltonians. It is possible in the perspective to add some other types of carriers and other interactions. In next chapters, we will discuss in more detail different types of transport calculations implemented in our code and present some test examples.

I discuss finally my contributions to the code. On the time of my arrival, the transmission calculation via the NEGF method was already implemented for electrons and phonons. Moreover, it was already possible to simulate the dynamics of wave packets with the Schrödinger method. I have started by implementing the wave packet propagation with Chebyschev method. This algorithm works much faster than the Schrödinger method, reducing the computation time by the factor of 100. I have also suggested to use EPW to extract electron-phonon couplings in order to work in the same electron basis as provided by Wannier90 code. I contributed to the calculation of the phonon Hamiltonian too. I have also implemented the spin-spin interaction for magnetic junctions and developed the formalism of the 2D wave packets.

For NEGF calculations, our initial goal was to calculate the electron and phonon currents in a self-consistent way taking into consideration electron-phonon interactions in the junctions affecting both electron and phonon channels -a "full approach" implying a self-consistent loop over electron and phonon Green's functions. I have started by working on the elastic currents implementation from the transmission functions already implemented in the code, for both electrons and phonons. I contributed next to the development of the current calculation with the electron-phonon interactions for electrons in the model of a junction (molecule) coupled to its own thermal bath and described therefore by a bare phonon Green's function. I have further started to implement the phonon's current calculation in inelastic regime and also to model the electron current by renormalising electronic Green's functions with a dressed phonon Green's functions due to phonon contacts. We still need to test carefully these preliminary implementations, and few further steps will be needed to achieve our final objective.

I would emphasise again that our code implements a quantum description of electron and phonon carriers based on realistic parameters extracted from ab initio

INTRODUCTION TO OUR TRANSPORT CODE AND MY CONTRIBUTION31

DFT calculations and should therefore apply to a variety of real materials such as 1D molecular junctions, 2D systems or material interfaces, to study different kinds of electron and phonon related phenomena.

Chapter 2

Wave packet

Wave packet (WP) methods are well studied and widely used approaches in transport field [START_REF] Bokes | Stroboscopic wave-packet description of nonequilibrium many-electron problems[END_REF]. WP dynamics can be simulated by different methods. The simplest one consists in applying directly the Schrödinger equation [START_REF] Konôpka | Stroboscopic wave packet description of time-dependent currents through ring-shaped nanostructures[END_REF][START_REF] Konôpka | Wave-packet representation of leads for efficient simulations of time-dependent electronic transport[END_REF] to evolve in time the wave function. It can also be resolved by using Green's functions [START_REF] Gaury | Numerical simulations of time-resolved quantum electronics[END_REF]. Another way consists in calculation of an evolution operator using different methods. Halimed et al. [START_REF] Jad | Chebyshev matrix product state approach for time evolution[END_REF] compare the Chebyshev method with the Krylov's one, and the Suzuki-Trotter time evolution and conclude that Chebyshev method may be more efficient and fast. Some codes already implement this scheme with and without including interactions. So far, the wave packet method was developed to study phonon transport either with Boltzmann equation [START_REF] Triozon | Diffusion quantique et conductivité dans les systèmes apériodiques[END_REF][START_REF] Yu | Wave packet dynamics in periodically kicked nonlinear systems[END_REF] or within classical or semi-classical approaches. Moreover, electron-phonon interactions were included mainly for model systems as in Ref. [START_REF] Monturet | Inelastic effects in electron transport studied with wave packet propagation[END_REF]. Our contribution consists in proposing a quantum-mechanical treatment of a phonon dynamics using a single-phonon wave packets constructed in the real space. We also developed a framework to perform calculations on the same footing for electrons and phonons, including electronphonon or spin-spin interactions, with realistic parameters describing real 1D and 2D systems.

Wave packet model

In a periodic lattice, quasi-particles as electrons or phonons are delocalized and can be treated as Bloch waves characterized by their wave vector k. The localized state of the particle can be described by a function describing the probability's amplitude to find a carrier in a specific point in space. The WP will typically have a finite width, which reflects the uncertainty in the position of the particle. The wave function describing our propagating particle can be defined as:

ψ(r i ) = F (r i )ϕ k (r i ) (2.1)
where r i runs over all the orbitals of the system. This function is the product of the envelope function F (r i ) defining the shape and localization of WP and the Bloch wave function ϕ k (r i ) encoding all the information on the direction of propagation, the velocity or the phase of the wave function. The total number of orbitals N tot includes the orbitals of the central part N C and of the lead N L . As the lead is a periodic lattice, it can be defined by the primitive cell, which reduces the memory allocation for calculations. In this way, the total number of orbitals in the lead N L is expressed with the number of orbitals per cell, N orb , and the number of cells in the lead, N cell , as N L = N orb ×N cell . With WP, we work with large scale systems, represented in Figure 2.1. The probability to find a carrier in a cell is represented by |ψ| 2 , which has a Gaussian shape due to the envelope function F (r). The white boxes represent infinite lead's cells. The cells overlapping with molecular junction supercell are represented in grey. In green, we represent the cell R mol related to the molecule and few possible orbitals of the lead perturbed my the molecule. We can notice that the model is valid for 2D systems as well as for 1D system if the number of cells in y direction is set to one. The total number of orbitals can be expressed as:

N tot = N orb (N cell -1) + N C .
Then, all carriers coming from the lead evolve according to the lead Hamiltonian ĤL . The first property we can extract from this Hamiltonian is the band structure, which gives the energy, E(k), as a function of wave vector k. To do this, one makes a Fourier transformation of the lead Hamiltonian to reciprocal space:

H L mn (k) = R ĤL 0m,Rn e ikR , (2.2)
where the first index was fixed at the cell 0 and the sum runs over all cells, R, of the periodic lead. It can be noticed that the matrix H L (k) has the dimension of N orb × N orb and is Hermitian. It can be diagonalized as:

n H L n ′ n (k)U nm (k) = E m (k)U n ′ m (k)
where its eigenvalues, E m (k), give energy bands m and eigenvectors, U nm (k), define wave functions on the orbital n. The wave function, corresponding to the state km, can be then written as the Bloch function:

ϕ km (R, n) = U nm (k)e ikR (2.3)
at the orbital n in the cell R. The band structure (or phononic dispersion for phonons) is used to select the Bloch state at the energy of interest, determining possible k and m as transport channels. The initial wave function of Eq. (2.1) can be then constructed as:

ψ(R, n) = 1 c F (R)ϕ km (R, n), (2.4)
where c is the normalizing constant, and should be evolved in time with the full Hamiltonian. This wave function is associated with a carrier in the lead with the selected energy we want to propagate. In the majority of cases, the envelope function F to be used has a simple Gaussian shape:

F (R) = e -|R-R 0 | 2 σ 2 (2.5)
and is defined by two parameters, the initial cell position, R 0 , and the width of the packet controlled by σ.

After defining the wave function, we focus now on the propagation properties. The velocity and the direction of an electron in the lead is required to know where to send the electrons. In theory, the group velocity is determined by the lead band dispersion, the derivative of the band energy with respect to the wave vector, V k = ∇ k E(k)/ℏ. However, computing this expression involves working on a dense k-grid and is not numerically very accurate. Instead, one can use the definition of the velocity operator as the time derivative of the position operator expressed in the Heisenberg representation r(t), which leads to:

V = dr(t) dt = 1 iℏ [r, ĤL ] (2.6)
Consider now the Bloch function in Eq. 2.3. Its velocity in the x direction will be given by:

V x = ⟨ϕ| Vx |ϕ⟩ = - i ℏ ⟨ϕ| x ĤL -ĤL x |ϕ⟩ (2.7)
The expression for the velocity can be written in terms of orbital positions and Hamiltonian elements:

V x = - i ℏ ni ′ ϕ(n) * x n H L ni ′ -H L ni ′ x i ′ ϕ(i ′ ), (2.8)
where n is restricted to the unit cell 0 while i ′ runs over all the orbitals of the lead. Using now explicit Bloch form in Eq. 2.3 one can rewrite the expression in the final form:

V x = - i ℏ nn ′ ,R ϕ(0, n) * x 0,n H L 0n,Rn ′ -H L 0n,Rn ′ (R x + x 0,n ′ ) ϕ(0, n ′ )e ikR , (2.9)
where now n, n ′ run over the unit cell 0.

In this section, we constructed our wave packet in the lead part as a propagating Bloch function enveloped with a Gaussian function. Moreover, we defined some tools for Wave packet propagation as the band structure and the velocity which allow us to select the state at the specific energy E m (k). In next, we will focus on how to propagate the wave function we constructed. The most simple way to propagate a Wave packet is to compute on several time steps the wave function. As shown in Fig. 2.2, we discretize the time t n to which we associate a wave function at this time ψ(t n ). We pass from a time to the next one by the time step t n+1 = t n + dt n (in our case, we consider a constant step in time dt n = dt). We use the same process for wave functions:

Carrier propagation

Evolution by Schrödinger equation method

ψ(t n + 1) = ψ(t n ) + dψ(t n ).
(2.10)

To calculate the partial derivative of the wave function over the time, we use the time dependent Schrödinger equation:

iℏ dψ(t n ) dt = Ĥψ(t n ) ⇐⇒ dψ(t n ) = - i ℏ Ĥψ(t n )dt.
(2.11)

Starting from the wave packet we initialized in the previous section, we calculate iteratively the wave function at each step from the previous time with the expression:

ψ(t n + 1) = ψ(t n ) - i ℏ Ĥψ(t n )dt . (2.12)
This method is straightforward. In case if we want to apply a Hamiltonian with a time dependence, we can modify at each iteration the Hamiltonian. However, it is time-consuming to calculate the product Ĥψ(t n ) for huge systems and at the fine grid of time steps. Moreover, the dependency on the previous step makes it difficult to implement parallel calculations.

Evolution with Chebyshev polynomials

There are many publications on time evolution of quantum wave packets. Various methods, improving on direct use of the Schrödinger equation, have been proposed implementing the computation of the evolution operator:

Ŝ(t 2 , t 1 ) = e -i ℏ Ĥ(t 2 -t 1 ) . (2.13)
Chebyshev polynomials have various applications in mathematics and physics, including the computation of the evolution operator. In particular, Chebyshev method was used in many studies on spin polarized systems [START_REF] Bruognolo | Dynamic structure factor of the spin-1 2 xxz chain in a transverse field[END_REF] and with electronphonon interactions too [START_REF] Ishii | Quantum transport properties of carbon nanotube field-effect transistors with electron-phonon coupling[END_REF]. This method allows to calculate the diffusion and the conductance too [START_REF] Ishii | Order-n electron transport calculations from ballistic to diffusive regimes by a time-dependent wave-packet diffusion method: Application to transport properties of carbon nanotubes[END_REF][START_REF] Tae | Unitary quantum time evolution by iterative lanczos reduction[END_REF]. With this method, one is able to calculate the wave function |ψ(t n )⟩ directly from the initial wave function |ψ(t 0 )⟩ at any time of interest:

|ψ(t n )⟩ = Ŝ(t n , t 0 ) |ψ(t 0 )⟩ . (2.14)
The evolution operator is often used to solve differential equations numerically. Thereafter, we will consider the operator of interest as a complex function of the Hamiltonian Ŝ = f ( Ĥ) (the evolution operator being the exponential function).

To calculate it we approximate the operator as the polynomial sum up to the order N as f ( Ĥ) ≈ R N ( Ĥ). First, we define orthogonal polynomials as a basis of polynomials P n with the hermitian product between them defined as:

(P m |P n ) = P m (E)P n (E)n(E)dE = δ mn (2.15)
with some specific density function, n(E). These polynomials satisfy recurrence relations:

E P n (E) = a n P n (E) + b n-1 P n-1 (E) + b n P n+1 (E) , (2.16)
with some a n and b n coefficients being real and positive. The goal is therefore to find the best polynomial approximation of an operator f ( Ĥ), i.e. the one that gives the best accuracy at a given numerical cost. For an order N , the polynomials which minimize the difference between R N ( Ĥ) and f ( Ĥ) have the form [START_REF] Triozon | Diffusion quantique et conductivité dans les systèmes apériodiques[END_REF]:

R N ( Ĥ) = N n=0 (P n |f )P n ( Ĥ). (2.17)
In the orthogonal polynomial's theory, the polynomial R N ( Ĥ) converge to the operator f ( Ĥ) when N → +∞. From this, we can try to solve Eq. 2.14 by developing the product of f ( Ĥ) and a state |ψ⟩:

f ( Ĥ) |ψ⟩ = lim N →+∞ R N ( Ĥ) |ψ⟩ = lim N →+∞ N n=0 (P n |f )P n ( Ĥ) |ψ⟩ . (2.18)
Let N be the degree of the expansion for which we estimate that the accuracy is sufficient. We then have to calculate numerically the sum:

f ( Ĥ) |ψ⟩ ≈ N n=0 (P n |f )P n ( Ĥ) |ψ⟩ . (2.19)
The coefficients (P n |f ) depend on the kind of polynomials we use. This part of the calculation of Eq. 2.19 is not computationally expensive. On the contrary, the calculation of |ψ n ⟩ = P n ( Ĥ) |ψ⟩ will run over N tot × N tot elements for each N vector. The computer memory will be able to store few of them in the same time. Thus, we have to save the sum of them at each iteration. |ψ n ⟩ can be expressed from Eq. 2.16 which gives:

b n |ψ n+1 ⟩ = Ĥ |ψ n ⟩ - N m C P |ψ m ⟩ , (2.20)
where C P is a tri-diagonal matrix composed of a and b coefficients:

C P =       a 0 b 0 0 ... ... 0 b 0 a 1 b 1 ... ... ... ... \ \ \ ... ... 0 ... b n-1 a n b n 0 ... ... ... \ \ \       . (2.21)
We initialize the 0 and 1 st orders, and we can write the general expression for the term n + 1 as a function of n and n -1 terms:

|ψ 0 ⟩ = P 0 ( Ĥ) |ψ⟩ = |ψ⟩ |ψ 1 ⟩ = Ĥ b 0 |ψ 0 ⟩ - a 0 b 0 |ψ 0 ⟩ |ψ n+1 ⟩ = Ĥ b n |ψ n ⟩ - a n b n |ψ n ⟩ - b n-1 b n |ψ n-1 ⟩ . (2.22)
In this way, we need to store only 3 vectors at the same time for each n+1 iteration which allows to reduce a memory consummation.

In the case of Chebyshev polynomials, the density function is of the form,

n(E) ∝ 1/ [1 -(E -a) 2 /4b 2 ].
Recurrence coefficients are constants, a n = a and b n = b, and are defined by the energy window W of the system's spectrum, W = max( Ĥ) -min( Ĥ), a corresponds to the middle point of the spectrum and 2b = W/2. To be sure to converge, we take a second energy window:

W = 2 max( Ĥ) -min( Ĥ) .
(2.23)

For Chebyshev polynomials, (P n |f ) for exponential function are expressed in terms of Bessel functions J n (x). These coefficients are initialized at the beginning as:

(P 0 |f ) = e -i ℏ at J 0 (- 2bt ℏ ) , (2.24)
and are given for other orders by:

(P n |f ) = (-i) n √ 2 e -i ℏ at J n (- 2bt ℏ
) .

(2.25)

In our code, we put as a parameter the convergence criteria α defined by

α = e 1 W t 4ℏ(N + 1) , (2.26)
which takes values in the interval 0 ≤ α ≤ 1. These criteria indicate if the N polynomial degrees are adequate to describe the evolution operator for a time difference t and an energy windows W . The lower is α, the better is the precision of calculation, and thus the wave packet do not diverge in time. We extract the power N from Eq. 2.26 as a function of α and a maximum time t:

N = eW t 4ℏα -1. (2.27)
The strength of the Chebyshev method is in its efficiency and the rapidity of calculation. The heaviest part is the calculation of |ψ n ⟩, it has to be done once for different time steps evolution. The time dependent part is in coefficient (P n |f ). We note, however, that within this method the Hamiltonian has to be time-independent which limits its applicability to some extent, but for majority calculations, the Chebyshev method is widely used and very much appreciated. Now, we will discuss different models of propagation with interactions.

Different models of propagation

Electron-phonon interaction

The first model of propagation we implement, similar to Mingo et al. [START_REF] Mingo | Calculation of the inelastic scanning tunneling image of acetylene on cu (100)[END_REF], is based on the Frölich-like Hamiltonian:

Ĥ = Ĥe + Ĥph + Ĥep , (2.28)
where Ĥe is the electronic Hamiltonian, Ĥph represents the phonon system and Ĥep is the electron-phonon coupling term. Electron Hamiltonian can be expressed in the tight-binding form as follows:

Ĥe = Ntot jj ′ t j ′ j ĉ † j ′ ĉj (2.29)
where ĉ † j and ĉj are creation and annihilation operators at the site j and t j ′ j are hopping parameters.

For phonons, Ĥph can be expressed in vibration mode basis. As for electrons, this Hamiltonian is composed of creation ( b † ) and annihilation ( b) operators, for the mode λ:

Ĥph = λ ℏω λ b † λ bλ + 1 2 (2.30)
Electron-phonon couplings, Eq. 1.40, are written in the basis of Wannier orbitals and atomic displacements. As we consider interactions only in the junction, the periodicity brought by the displacement basis is less appropriate than the vibration mode basis. With the expression of the displacement operator Eq. 1.36, the electron-phonon couplings have the following expression in the vibration mode basis [START_REF] Frederiksen | Inelastic transport theory from first principles: Methodology and application to nanoscale devices[END_REF]:

g λ j ′ j = iα ℏ 2ω λ V λ iα g j ′ j,iα (0; 0), (2.31)
where λ should replace {qν} and the translational vectors, R, R ′ , for the supercell of the junction are set to 0 to single out all the displacements and Wannier functions in the unit cell "0". From this, the interaction term Ĥep can be rewritten as: This expression is visualized in Figure 2.3. Blue arrows describe the electron coming from the site j and departing from the site j ′ after the interaction. During the interaction, a phonon with a mode λ, represented in a green wave, can be created or annihilated. This interaction is characterized by the electron-phonon element g λ j ′ j . It can be noticed that the energy of the electron is different between the site j and the site j ′ .

Ĥep = j ′ jλ g λ j ′ j ĉ † j ′ ĉj ( b † λ + bλ ) (2.32) j j ′ j ′ j
The total wave function of the system |Ψ⟩ described by the Hamiltonian (2.28) can be generally written as a superposition:

|Ψ⟩ = jn ϕ n (j) • |j⟩ ⊗ |n⟩ , (2.33)
where |j⟩ and |n⟩ describe an electron on the site j and a phononic state with n phonons in a chosen mode λ, respectively, the products of which make up the basis states for the total system. For simplicity, we restrict here our discussion to only one active mode λ but the approach is straightforwardly generalized to the multi-mode case. This state can not be in general represented as a simple product of the electronic wave function |ψ e ⟩ and the phononic wave function |ψ ph ⟩:

|ψ e ⟩ = j A j |j⟩ |ψ ph ⟩ = n B n |n⟩ , (2.34)
because ϕ n (j) can not be generally factorized:

ϕ n (j) ̸ = A j B n
representing generally an "entangled" state of the composed electron-phonon system.

Note that even if it were the case and:

|Ψ⟩ = |ψ e ⟩ × |ψ ph ⟩ . (2.35)
at the initial moment, the state will get "entangled" at later times due to el-ph interactions and will evolve into the general form of Eq. 2.33.

In our model, represented in Fig. 2.4, the electronic wave packet ϕ 0 (j) represented by (0) is initially on the chain with no phonons. It propagates to the molecular junction and at (1) interacts with a molecule and will be able to create or annihilate phonons with a certain probability, led by interaction term Ĥep . At later times, the probability to find an electron (2) is redistributed over different phononic channels. Each channel corresponds to a different number of phonons, and an electron is described with the wave function ϕ n (j). Note that the total energy, E = E e + E ph , is conserved, therefore electron wave packets propagate slightly slower when the number of phonons increase.

We note finally that the probability to have n phonons in the system, P n , after interactions is given by:

P n = ⟨n| Ψ |n⟩ = j |ϕ n (j)| 2
(2.36) Similarly, the probability of an electron to be found at the site i, P i , is expressed as:

P i = ⟨i| Ψ |i⟩ = n |ϕ n (i)| 2 (2.37)
because of orthogonality of phononic states |n⟩. It means, for example, that the only electron subsystem after all the interactions is not in a pure state but is rather described by a statistical mixture of electronic states |ϕ n ⟩.

Spin-Spin interaction

We now turn to the second model of interactions, namely spin-spin interactions between the conduction electron and the local spin (of a magnetic molecule, for example). In Quantum-ESPRESSO, we have a possibility to calculate systems polarized in spin by inducing a starting magnetization on some atoms and letting the system relax during SCF calculation. The magnetic interaction we implement is inspired by the famous s-d or Kondo model. It includes ingredients needed to study Kondo effect, which was first described by Jun Kondo in 1964. When a magnetic impurity is introduced in a non-magnetic metallic host, at certain conditions and at low temperature it becomes screened due to interactions with conduction electrons forming a collective spin singlet state. The screening comes from spin-flip events which can change the spin projections of the impurity and of the electron at the same time, preserving the projection of the total spin. We mention some previous works such as [START_REF] Kim | How many electrons are needed to flip a local spin[END_REF][START_REF] Kim | Quantum mechanics of spin transfer in coupled electron-spin chains[END_REF][START_REF] Pertsova | Electrical control of spin dynamics in finite one-dimensional systems[END_REF] where wave packet approach was also used to explore spin-spin interactions. Figure 2.5: Scheme of two spin-flip paths contributing to the s-d model. The states of the conduction electron i and j are represented by upper red lines, and an impurity state is represented by the line beneath. In a), the electron arrives with the spin opposite to the spin of the localized state. From this, there are two options: in b) electron first jumps onto the impurity level which gets doubly occupied; in c) an electron first jumps out of the impurity level and then another electron fills it. Both processes contribute to the same final state in d) where the spin exchange between the impurity and the conduction electron is observed.

The relevant total Hamiltonian is written as:

Ĥ = ijσ t ij ĉ † iσ ĉjσ + ij J ij ⃗ s ij ⃗ S, (2.38)
where t ij are hopping elements between orbitals i and j. The interaction of the conduction electron and the local spin is provided by the second term of Eq. 2.38. The conduction electron spin operators are described using Pauli matrices ⃗ σ:

⃗ s ij = αβ ĉ † iα ⃗ σ αβ 2 ĉjβ (2.39)
and local spin is represented by the spin vector operator ⃗ S. The matrix J represents the coupling constants for exchange interactions between sites i and j. We note that conduction electron spin operators, ⃗ s ij , refer to all the Wannier orbitals ij (including those participating in the local spin formation), which is somewhat different from the standard Kondo model formulation.

The commutation rules for spin operators are as follows:

[ Ŝx , Ŝy ] = i Ŝz , [ Ŝy , Ŝz ] = i Ŝx , [ Ŝz , Ŝx ] = i Ŝy .
It is convenient to define ladder operators as:

Ŝ+ = Ŝx + i Ŝy and Ŝ-= Ŝx -i Ŝy .
The Kondo Hamiltonian (2.38) can be rewritten, using these ladder operators:

Ĥ = Ĥ0 + ij J ij ŝz ij Ŝz + 1 2 ŝ+ ij Ŝ-+ ŝ- ij Ŝ+ , (2.40)
where Ĥ0 is the electron Hamiltonian of the system without spin-spin interactions.

The spin of an electron or magnetic impurity is represented by the spin quantum number S determining the magnitude of the spin. For conduction electron it is of course 1/2, but it could take any positive integer or half-integer values for the local spin S = 0, 1/2, 1, 3/2, 2, etc. S goes with the quantum number m s which refers to the spin projection on the z direction and takes values {-S, -S +1, .., S -1, S}. Quantum numbers S and m s define the Zeeman basis, |S, m s ⟩, where spin operators act as follows:

Ŝz |S, m s ⟩ = m s |S, m s ⟩ Ŝ+ |S, m s ⟩ = S(S + 1) -m s (m s + 1) |S, m s + 1⟩ Ŝ-|S, m s ⟩ = S(S + 1) -m s (m s -1) |S, m s -1⟩ . (2.41)
The appearance of spin-spin interactions in the Kondo model can be understood (and derived) from another well-known model -the Anderson impurity model, which describes the conduction electron propagating through a site with a strong Coulomb interaction. Two processes contributing to spin-flip contributions are illustrated in Figure 2.5. The first one goes through the doubly occupied impurity level, while the second process involves only single occupation of the level. The both paths contribute to the spin-flip of the conduction electron between initial and final states. Taking into consideration these two processes, one can set up the effective Kondo model with parameters J related to those of the Anderson impurity model [START_REF] Cyril | The Kondo problem to heavy fermions[END_REF][START_REF] Yosida | Theory of magnetism[END_REF].

Magnetic defect

Pertubed cells

Ideal unperturbed cells

Figure 2.6: General model of the system with spin-spin interactions. Each square represents a cell of the system. In red, the magnetic impurity; in green, the lead cells perturbed by the magnetic impurity; in blue, the lead cells considered too far from the junction to be affected by the magnetization and having therefore J ij = 0.

To construct the junction's Hamiltonian Eq. (2.40), we need to know the Hamiltonian Ĥ0 and parameters J. For this, we need the Hamiltonians for spin up and spin down electrons separately constructed using the Wannier90 code. In DFT calculations, the local spin appears to have its maximum value m s = S so the DFT Hamiltonians for spin up or down should be compared with:

Ĥ↑↓ = Ĥ0 + ij J ij s ↑↓ ij S. (2.42)
Knowing the local spin S from DFT calculations, we can get Ĥ0 and couplings J ij from following equations:

Ĥ0 = 1 2 ( Ĥ↑ + Ĥ↓ ) (2.43)
and

J ij = ( Ĥ↑ ij -Ĥ↓ ij )/S, (2.44)
assuming conduction electron's spin s z = ±1/2 for up and down channels.

We can calculate now the spin-flip probability of a propagating electron as shown in Figure 2.6. We consider an electron coming from a non-magnetic lead (shown in blue) and going through the junction (shown in green). The spin of the electron will interact with the local spin (shown in red) during its propagation with exchange parameters J ij . determined as discussed above.

The total wave functions have the form:

|Ψ⟩ = j,ms ψ ms (j) • |j⟩ ⊗ |m s ⟩
and is described by electron wave functions, ψ ms (j), on 2S + 1 chains associated with different z-components of the local spin S (finite in contrast to the phononic case presented before).

Phononic propagation

We will start our discussion of the second quantization of atomic vibrations on the simplest case of a finite system, such as a molecule, for example. This will allow us to introduce the notion of phonons as quanta of vibrational excitations. Then we generalize the formalism to periodic crystals, making a direct analogy of the procedure with the previous case of a finite system.

Finite system

The ensemble of N atoms of a molecule is described in quantum mechanics by the Hamiltonian with the kinetic and potential energy parts:

Ĥ = T + Ê({x iα }) (2.45)
where x iα stand for displacements in Cartesian coordinates (α = x, y, z) of the i-th atom (i = 1..N ) from equilibrium positions. If we develop the potential energy in a Taylor expansion, the first constant term, E(0), can be put to zero, the firstorder terms, ∂E ∂x iα , giving the forces, are also zero at equilibrium. If we keep only second-order terms (harmonic approximation) the Hamiltonian takes the following form:

Ĥ = iα p2 iα 2m i + 1 2 iα,i ′ α ′ ∂ 2 E ∂x iα ∂x i ′ α ′ xiα xi ′ α ′ (2.46)
where conjugate displacement and momentum operators satisfy canonical commutation relations:

[x iα , pi ′ α ′ ] = ihδ ii ′ δ αα ′
It is convenient to re-scale displacement and momentum operators like:

piα √ m i - → Piα ; √ m i xiα - → Xiα . (2.47)
The canonical commutation relations are preserved:

[ Xiα , Pi ′ α ′ ] = ihδ ii ′ δ αα ′
and the Hamiltonian takes the following form:

Ĥ = 1 2 iα P 2 iα + 1 2 iα,i ′ α ′ Xiα D iα,i ′ α ′ Xi ′ α ′ (2.48)
with the dynamical matrix D defined by the second derivative of the potential energy with respect to atomic displacements as follows:

D iα,i ′ α ′ = 1 √ m i m i ′ ∂ 2 E ∂x iα ∂x i ′ α ′ = ∂ 2 E ∂X iα ∂X i ′ α ′ (2.49)
We can notice that the first and second terms are of the form of ⟨P |P ⟩ and ⟨X| D |X⟩, constructed on vector-columns P iα and X iα , respectively. One can therefore proceed by solving the eigenvalue problem for the dynamical matrix:

i ′ α ′ D iα, i ′ α ′ φ ν i ′ α ′ = ω 2 ν φ ν iα (2.50)
which defines normal modes ν with their frequencies ω ν and corresponding polarization vectors φ ν . Since the D matrix is by construction real and symmetric, its eigenvectors φ ν are also real-valued. Now we can insert the identity resolution and the D-matrix representation in the normal mode's basis φ ν which we simply denote as |ν⟩:

Î = ν |ν⟩ ⟨ν| , D = ν ω 2 ν |ν⟩ ⟨ν| (2.51)
into the kinetic and potential energy terms, respectively, of the Hamiltonian which can be thus rewritten as:

Ĥ = 1 2 ν P ν ν P + 1 2 ν ω 2 ν X ν ν X = 1 2 ν P 2 ν + 1 2 ν ω 2 ν X2 ν (2.52)
with new displacement and momentum operators referring to normal modes:

Xν = iα φ ν iα Xiα Pν = iα φ ν iα Piα (2.53)
and the inverse transformation:

Xiα = ν φ ν iα Xν Piα = ν φ ν iα Pν (2.54)
New operators are Hermitian and satisfy again the canonical commutation relations:

[ Xν , Pν ′ ] = ihδ νν ′ (2.55)
The Hamiltonian from Eq. 2.52 describes the set of independent oscillators corresponding to different normal modes. Therefore, it can be conveniently rewritten as:

Ĥ = ν ℏω ν b † ν bν + 1 2 (2.56)
by introducing creation and annihilation operators for each mode ν:

bν = ω ν 2ℏ Xν + i ω ν Pν b † ν = ω ν 2ℏ Xν - i ω ν Pν (2.57)
These creation operators act on the vacuum state and create an excited state in a specific mode, called a phonon. From Eqs. 2.55, 2.57 one sees that creation and annihilation operators satisfy the commutation relations:

[ bν , b † ν ′ ] = δ νν ′ (2.58)
so that phonons satisfy the boson statistics.

The basis states for the system can be characterized by the number of phonons in each mode, |n 1 , .., n N ⟩. Applying the creation or annihilation operator will increase or decrease the number of phonons in a specific mode, like:

b † ν |n 1 , .., n ν , .., n N ⟩ = √ n ν + 1 |n 1 , .., n ν + 1, .., n N ⟩ bν |n 1 , .., n ν , .., n N ⟩ = √ n ν |n 1 , .., n ν -1, .., n N ⟩ (2.59)
Note, finally, that displacement operators can be expanded in terms of creation and annihilation operators as follows:

Xiα = ν ℏ 2ω ν φ ν iα ( bν + b † ν ) . (2.60)
In order to describe the quantum propagation of a phonon, we need to define local creation and annihilator operators in the displacement basis as:

biα = ν φ ν iα bν b † iα = ν φ ν iα b † ν (2.61)
which corresponds to the unitary rotation of the basis. In this way, we can think of creating or annihilating a phonon on a specific site iα by applying the operator b † iα or biα :

b † iα |0⟩ = |ϕ iα ⟩ biα |ϕ iα ⟩ = |0⟩ (2.62)
The general one-phonon state can be then expanded as:

|Ψ⟩ = iα φ iα |ϕ iα ⟩ (2.63)
where coefficients φ iα can be seen as the real-space wave function of the phonon. Its evolution in time will be governed by the Hamiltonian matrix:

H iα,i ′ α ′ = ⟨ϕ iα | Ĥ |ϕ i ′ α ′ ⟩ = ν ℏω ν φ ν iα φ ν i ′ α ′ (2.64)
which simply means that it is proportional to the square root of the Dynamical matrix, H = ℏ √ D. This is analogous to the tight-binding Hamiltonian for electrons. The phononic wave function, φ iα , can be then propagated with either Schrödinger equation or by the Chebyshev method, as in the electronic case.

Note finally that the wave function for the phonon state b † ν |0⟩ in the real space representation is simply φ ν iα .

Periodic system

The Hamiltonian for a periodic system can be generally written as:

Ĥ = 1 2 iαR P 2 iαR + 1 2 iαR, i ′ α ′ R ′ XiαR D iαR, i ′ α ′ R ′ Xi ′ α ′ R ′ (2.65)
where X iαR stand for atomic displacements in Cartesian coordinates (α = x, y, z) of the i-th atom in the cell defined by the translation R.

Again, first and second terms have the form of ⟨P |P ⟩ and ⟨X| D |X⟩, respectively. We diagonalize therefore the dynamical matrix and search for its eigenfunctions. As for electrons, the standard choice for periodic dynamical matrix are the Bloch functions:

i ′ α ′ R ′ D iαR, i ′ α ′ R ′ φ qν i ′ α ′ R ′ = ω 2 qν φ qν iαR φ qν iαR = 1 √ N V qν iα e iqR i ′ α ′ D iα, i ′ α ′ (q)V qν i ′ α ′ = ω 2 qν V qν iα where D iα, i ′ α ′ (q) = i ′ α ′ R ′ D iαR, i ′ α ′ R ′ e iq(R ′ -R)
is the Fourier transforming of the dynamical matrix. Note that for -q one always has:

ω -qν = ω qν ; V -qν iα = V qν * iα due to D iα, i ′ α ′ (-q) = D iα, i ′ α ′ (q) * .
Inserting again the unity resolution and the spectral representation of the dynamical matrix in Eq. 2.65 one obtains:

Ĥ = 1 2 qν Pqν P † qν + 1 2 qν ω 2 qν X † qν Xqν (2.66)
where

Pqν = ⟨P |qν⟩ = iαR PiαR φ qν iαR Xqν = ⟨qν|X⟩ = iαR XiαR φ qν * iαR (2.67)
with inverse transformation:

PiαR = qν Pqν φ qν * iαR XiαR = qν Xqν φ qν iαq (2.68)
They still represent the conjugated variables:

[ Xqν Pq ′ ν ′ ] = iℏδ qq ′ δ νν ′
but are not self-adjoint operators, namely:

P † qν = P-qν X † qν = X-qν (2.69)
so that the Hamiltonian can be also written as:

Ĥ = 1 2 qν Pqν P-qν + 1 2 qν ω 2 qν X-qν Xqν (2.70)
In Appendix 4 we show how Hamiltonian Eq. 2.70 can be written in the similar form as for finite systems:

Ĥ = qν ℏω qν b † qν bqν + 1 2 , (2.71)
where creation and annihilation operators for the mode qν adopt well-known form:

bqν = ω qν 2ℏ Xqν + i ω qν P-qν b † qν = ω qν 2ℏ X-qν - i ω qν Pqν . (2.72)
Finally, real space displacement operators can be also cast into well-known expression: and define phononic basis in the real space:

XiαR = qν ℏ 2ω qν φ qν iαR bqν + b † -qν . ( 2
|ϕ iαR ⟩ = b † iαR |0⟩ . (2.75)
The general one-phonon state can then be expanded as:

|Ψ⟩ = iαR ϕ iαR |ϕ iαR ⟩ , (2.76)
where coefficients ϕ iαR can be considered as the real-space wave function of the phonon. The time evolution is defined by the Hamiltonian matrix:

H iαR,i ′ α ′ R ′ = ⟨ϕ iαR | Ĥ |ϕ i ′ α ′ R ′ ⟩ = qν ℏω qν φ qν iαR φ qν * i ′ α ′ R ′ .
(2.77)

Note again that the wave function for the Bloch phonon state b † qν |0⟩ in the real space representation is simply φ qν iαR . Weighted by the Gaussian function, it can be considered as a starting wave function for a phononic wave packet which can be propagated using the Schrödinger equation or the Chebyshev method.

Chapter 3

Non-equilibrium Green Function Method

Green's Functions

The Green's functions method is widly-used in transport calculation codes. Some codes perform current calculations as Quantum ATK [START_REF] Schneider | Atkforcefield: a new generation molecular dynamics software package[END_REF], openMX, KNIT [START_REF] Kazymyrenko | Knitting algorithm for calculating green functions in quantum systems[END_REF], GOLLUM [START_REF] Ferrer | Gollum: a next-generation simulation tool for electron, thermal and spin transport[END_REF], etc. Within the NEGF formalism one can also consider interactions in the Central region, as in Refs. [START_REF] Sergueev | Efficiency of thermoelectric energy conversion in biphenyl-dithiol junctions: Effect of electron-phonon interactions[END_REF][START_REF] Aeberhard | Microscopic nonequilibrium theory of quantum well solar cells[END_REF][START_REF] Frederiksen | Inelastic transport theory from first principles: Methodology and application to nanoscale devices[END_REF], which include electron-electron, electron-phonon, phonon-phonon, or electron-photon interactions. It is important however that interactions are limited by the Central "active" region while electrodes are described by noninteracting Hamiltonians. With our code, we include the possibility to add the electron-phonon interactions in the current calculations. Our initial goal was to calculate self-consistently the inelastic corrections to both electron and phonon Green's functions due to electron-phonon interactions in the Central region (we have had time however to advance only on the first part of the task -on the effect of interactions on the electron transport). This full treatment of two connected electron-phonon transport channels is not yet implemented to the best of our knowledge in any code.

We start our discussion of Green's functions approach to the transport with an equilibrium case, and then present the formalism for a non-equilibrium situation.

Equilibrium Case

In physics, Green's functions are an important concept used in different domains to solve, for example, differential equations or to study many-body systems. It provides a way to calculate various properties for a given Hamiltonian with possible interactions. There are many references discussing different aspects of the Green's 55 functions formalism and the variety of its applications [START_REF] Eleftherios | Green's functions in quantum physics[END_REF][START_REF] Gerald D Mahan | Many-particle physics[END_REF][START_REF] Haug | Quantum kinetics in transport and optics of semiconductors[END_REF][START_REF] Bruus | Introduction to Many-body quantum theory in condensed matter physics[END_REF][START_REF] Balzer | Nonequilibrium Green's Functions Approach to Inhomogeneous Systems[END_REF][START_REF] Pourfath | The non-equilibrium Green's function method for nanoscale device simulation[END_REF][START_REF] Carlos | Molecular Electronics : An Introduction to Theory and Experiment[END_REF].

One usually starts by discussing three different pictures in quantum mechanics. In the Schrödinger picture, a quantum state is a function of time while operators are time-independent:

|Ψ S (t)⟩ = e -i Ĥt/ℏ |Ψ S (0)⟩ .

(3.1)

ÔS (t) = ÔS (3.2)
In the Heisenberg picture, on the contrary, the operators evolve in time while a state is time-independent:

|Ψ H ⟩ = |Ψ S (0)⟩ . (3.3) ÔH (t) = e i Ĥt/ℏ ÔS e -i Ĥt/ℏ (3.4)
Finally, in the intermediate representation, called the interaction picture, both states and operators depend on time. It is useful when the total Hamiltonian is decomposed as:

Ĥ = Ĥ0 + V , (3.5)
where Ĥ0 is a non-interacting known part and V can be time-dependent and include interactions. The states and operators evolve as:

|Ψ I (t)⟩ = e i Ĥ0 t/ℏ |Ψ S (t)⟩ , (3.6)
ÔI (t) = e i Ĥ0 t/ℏ ÔS e -i Ĥ0 t/ℏ . (3.7)
The time evolution is therefore driven by the interaction part of the Hamiltonian:

iℏ∂ t |Ψ I (t)⟩ = VI (t) |Ψ I (t)⟩ , (3.8)
where VI (t) is the interaction term in the interaction picture. It is important that all the three pictures, dividing differently the time evolution between states and operators, give the same expectation values, ⟨Ψ| O |Ψ⟩. From Eq. 3.6 one can find the relation between |Ψ I (t)⟩ and |Ψ I (t 0 )⟩:

|Ψ I (t)⟩ = e i Ĥ0 t/ℏ e -i Ĥ(t-t 0 )/ℏ |Ψ S (t 0 )⟩ = e i Ĥ0 t/ℏ e -i Ĥ(t-t 0 )/ℏ e -i Ĥ0 t 0 /ℏ |Ψ I (t 0 )⟩ (3.9)
and define the operator Ŝ(t, t 0 ) from:

|Ψ I (t)⟩ = Ŝ(t, t 0 ) |Ψ I (0)⟩ (3.10)
which is called the evolution operator and has the following properties:

• Ŝ † (t, t 0 ) Ŝ(t, t 0 ) = Ŝ(t, t 0 ) Ŝ † (t, t 0 ) = 1 • Ŝ † (t, t 0 ) = Ŝ-1 (t, t 0 ) • Ŝ(t, t 0 ) Ŝ(t 0 , t) = 1 • Ŝ(t 0 , t 1 ) Ŝ(t 1 , t 2 ) = Ŝ(t 0 , t 2 ).
Note that the evolution operator allows to relate operators in the Heisenberg picture in Eq. 3.4 with those in the interaction picture:

ÔH (t) = Ŝ(0, t) ÔI (t) Ŝ(t, 0) (3.11)
Now, the (time-ordered) electronic Green's function (at zero temperature) is defined as:

iℏG(r ′ , t ′ ; r, t) = ⟨Ψ 0 | T t ψH (r ′ , t ′ ) ψ † H (r, t) |Ψ 0 ⟩ ⟨Ψ 0 |Ψ 0 ⟩ , (3.12)
where Ψ 0 is the interacting ground state and the time-ordering operators (for electrons) are defined by:

T t Â(t 1 ) B(t 2 ) = θ(t 1 -t 2 ) Â(t 1 ) B(t 2 ) -θ(t 2 -t 1 ) B(t 2 ) Â(t 1 ). (3.13)
The Green's function G(r ′ , t ′ ; r, t) represents the amplitude of probability to create an additional electron in the ground state and to observe it at the later time. The time-ordering operator allows to describe at the same time both electrons and holes propagation: when t ′ > t it describes an electron propagating from r to r ′ , when t ′ < t -a hole propagating from r ′ to r.

One usually introduces also four other Green's functions called retarded, advanced, inferior and superior by following definitions:

G r (r ′ , t ′ ; r, t) = -iθ(t ′ -t)⟨{ψ H (r ′ , t ′ ), ψ † H (r, t)}⟩ G a (r ′ , t ′ ; r, t) = iθ(t -t ′ )⟨{ψ H (r ′ , t ′ ), ψ † H (r, t)}⟩ G < (r ′ , t ′ ; r, t) = i⟨ψ † H (r, t), ψ H (r ′ , t ′ )⟩ G > (r ′ , t ′ ; r, t) = -i⟨ψ H (r ′ , t ′ ), ψ † H (r, t)⟩ (3.14)
where {} are anticommutators. These functions are used to describe different kinds of correlations and allow to calculate various response functions (retarded functions). These four Green's functions are not independent of each other, obeying the following relationship:

G r -G a = G > -G < . (3.15)
From Eq. 3.8 one finds the equation for the evolution operator:

iℏ ∂ t Ŝ(t, t 0 ) = VI (t) Ŝ(t, t 0 ). (3.16)
Integrating this equation, we obtain a perturbation expansion:

Ŝ(t, t 0 ) = Ŝ(t 0 , t 0 ) - i ℏ t t 0 dt 1 VI (t 1 ) Ŝ(t 1 , t 0 ) = ∞ n=0 - i ℏ n t t 0 dt 1 ... t n-1 t 0 dt n VI (t 1 ) ... VI (t n ).
(3.17)

Using the time-ordering operator, this equation can be conveniently written as the time-ordered exponential:

Ŝ(t, t 0 ) = ∞ n=0 1 n! - i ℏ n t t 0 dt 1 ... t t 0 dt n T t VI (t 1 ) ... VI (t n ) = T t exp - i ℏ t t 0 dt ′ VI (t ′ ) . (3.18)
The Green's function in Eq. 3.12 is defined on the ground state |Ψ 0 ⟩ of an interacting system. According to the Gell-Mann and Low theorem [START_REF] Gell-Mann | Quantum electrodynamics at small distances[END_REF], it can be constructed from the ground state |Φ 0 ⟩ of a non-interacting system with the help of the evolution operator in two ways:

|Ψ 0 ⟩ = Ŝ(0, -∞) |Φ 0 ⟩ , (3.19)
e -iθ |Ψ 0 ⟩ = Ŝ(0, ∞) |Φ 0 ⟩ , (3.20)
assuming adiabatic switching on or off interactions at ±∞, where e iθ is some unknown phase. Inserting these expressions for the interacting ground state, we can rewrite the Green's function as following:

iℏG(r ′ , t ′ ; r, t) = ⟨Φ 0 | Ŝ(∞, 0)T t ψH (r ′ , t ′ ) ψ † H (r, t) Ŝ(0, -∞) |Φ 0 ⟩ ⟨Φ 0 | Ŝ(∞, 0) Ŝ(0, -∞) |Φ 0 ⟩ = ⟨Φ 0 | T t Ŝ(∞, -∞) ψI (r ′ , t ′ ) ψ † I (r, t) |Φ 0 ⟩ ⟨Φ 0 | Ŝ(∞, -∞) |Φ 0 ⟩ , (3.21)
where the field operators in the last line were transformed to the Interaction picture according to Eq. 3.11, or explicitly:

iℏG(r ′ , t ′ ; r, t) = ⟨Φ 0 | T t exp -i ℏ ∞ -∞ dt 1 VI (t 1 ) ψI (r ′ , t ′ ) ψ † I (r, t) |Φ 0 ⟩ ⟨Φ 0 | T t exp -i ℏ ∞ -∞ dt 1 VI (t 1 ) |Φ 0 ⟩ (3.22)
This final expression is used to calculate equilibrium Green's function using perturbation expansion with respect to the potential V . It contains an integral over the time axis which should be modified if the system is driven out of equilibrium, as we now will discuss.

Non-Equilibrium case

In the previous equation the integration ran from -∞ to +∞ because the Gell-Mann and Low theorem allowed to relate ground states of interacting and noninteracting systems assuming adiabatic switching on and off interactions. In nonequilibrium situation (in applied voltage, for example) it is not any more true for the time at +∞, the system is not guaranteed to return to its non-equilibrium noninteracting state after switching off interactions. Therefore, the only well-defined state is the one at -∞ and one is forced to use the contour C with two branches for the evolution operator S as it is demonstrated in Fig. 3.1, one going from -∞ to +∞ (the contour part C + ) and another one -from +∞ to -∞ (the contour part C -) with the contour parameter along the C denoted as τ . The advantage of this contour is that one avoids to refer explicitly to unknown state at +∞, and only the reference state is the one at -∞ which is well known. This formalism, often called Keldysh or Non-equilibrium Green's function (NEGF), was foreshadowed by the work of Julian Schwinger and proposed almost simultaneously by Keldysh [START_REF] Keldysh | Diagram technique for nonequilibrium processes[END_REF] and separately by Kadanoff and Baym [START_REF] Leo | Quantum statistical mechanics[END_REF]. It also allows to include naturally time-dependent perturbation potentials V (t) in the scheme.

One then defines the contour-ordered Green function as:

iℏG(1 ′ ; 1) = ⟨Ψ 0 | T C [ ψH (1 ′ ) ψ † H (1)] |Ψ 0 ⟩ (3.23)
using the contour ordering operator T C and shorthand notations (1) referring to (r 1 , τ 1 ) on the contour. This contour-ordered Green's function plays an analogous role in NEGF as the causal (or time-ordered) Green's function in equilibrium theory. It has for example a similar perturbation expansion:

iℏG(1 ′ ; 1) = ⟨Φ 0 | ŜC T C ψI (1 ′ ) ψ † I (1) |Φ 0 ⟩, (3.24)
Figure 3.1: Keldysh contour composed of two branches going from -∞ to +∞ and back from +∞ to -∞.

where the evolution operator on the contour is given by:

ŜC = T C exp - i ℏ C dτ VI (τ ) . (3.25)
Since time labels lie now on the contour with two branches, one has now four different Green's functions defined in the time domain:

G(1, 2) =        G ++ (1, 2) t 1 , t 2 ∈ C + G --(1, 2) t 1 , t 2 ∈ C - G < (1, 2) t 1 ∈ C + , t 2 ∈ C - G > (1, 2) t 1 ∈ C -, t 2 ∈ C + (3.26)
From this one sees that G < , G > , G ++ , and G --correspond to "lesser", "greater", time-ordered (or causal), and anti time-ordered Green's functions, respectively:

iℏG < (1 ′ , 1) = -⟨Ψ 0 | ψ † H (1) ψH (1 ′ ) |Ψ 0 ⟩ iℏG > (1 ′ , 1) = ⟨Ψ 0 | ψH (1 ′ ) ψ † H (1) |Ψ 0 ⟩ iℏG ++ (1 ′ , 1) = ⟨Ψ 0 | T t [ ψH (1 ′ ) ψ † H (1)] |Ψ 0 ⟩ iℏG --(1 ′ , 1) = ⟨Ψ 0 | T t[ ψH (1 ′ ) ψ † H (1)] |Ψ 0 ⟩ , (3.27)
where T t is the anti time-ordering operator. One can define, moreover, retarded and advanced Green's functions as:

G r =G ++ -G < G a =G ---G > . (3.28)
These definitions can be in fact applied to any function on the contour. From Eqs. 3.27 it follows that four Green's functions are not independent, in particular, one has:

G ++ + G --= G > + G < G r -G a = G > -G < . (3.29)
We also list some other useful relations:

G ++ (1, 2) = θ(t 1 -t 2 )G > (1, 2) + θ(t 2 -t 1 )G < (1, 2) G --(1, 2) = θ(t 2 -t 1 )G > (1, 2) + θ(t 1 -t 2 )G < (1, 2) G r (1, 2) =θ(t 1 -t 2 ) G > (1, 2) -G < (1, 2) G a (1, 2) =θ(t 2 -t 1 ) G < (1, 2) -G > (1, 2) .
(3.30)

We discuss now how to express Green's functions with the perturbation approach. The Wick decomposition allows a perturbation expansion of Green's functions starting from Eqs. 3.24, 3.25. In particular, if the interaction Hamiltonian includes external potential (due to coupling to electrodes, for example) and many-body interactions, H = H 0 + V ext + W int , the contour-ordered Green's will satisfy the Dyson equation:

G(1 ′ , 1) = g(1 ′ , 1) + C dτ 2 g(1 ′ , 2)V ext (2)G(2, 1) + C dτ 2 C dτ 3 g(1 ′ , 2)Σ(2, 3)G(3, 1) (3.31)
where g is the Green's function for non-interacting system, H = H 0 , and manybody interactions contribute to the (irreducible) self-energy functional, Σ. It can be calculated for e -e or e -ph interactions taking into account various Feynman diagrams. Fig. 3.2 presents some key ingredients, like interaction vertices and Green's functions, needed to construct them. Note that there is also symmetric Dyson equation:

G(1 ′ , 1) = g(1 ′ , 1) + C dτ 2 G(1 ′ , 2)V ext (2)g(2, 1) + C dτ 2 C dτ 3 G(1 ′ , 2)Σ(2, 3)g(3, 1) (3.32)
with non-interacting g on the right side of integrals.

In considering the Dyson equation on the contour one encounters different types of contour integrals such as for example:

F (τ, τ ′ ) = C dτ 1 A(τ, τ 1 )B(τ 1 , τ ′ ), (3.33)
where we skip all space variables for simplicity, the integration over them is assumed otherwise.

We briefly discuss now how to calculate the corresponding real-time functions, F </> or F r/a . The Langreth theorem proposes to deform the real timeline into a double contour shape, as it is shown in Fig. 3

.3. Lets calculate the "lesser" function

Non-interaction Green's function F < (t, t ′ ). We consider C 1 as a twisted contour C + of the contour in Fig. 3.1 and C 2 as a twisted contour C -:

F < (t, t ′ ) = C 1 dt 1 A ++ (t, t 1 )B < (t 1 , t ′ ) + C 2 dt 1 A < (t, t 1 )B ++ (t 1 , t ′ ).
(3.34)

In the first part, t 1 and t are on the same contour. That is why we get A ++ . On the contrary, we have B < because t 1 and t ′ are not on the same contour. Similarly, for the second integral, where t 1 is in the same contour as t ′ . The Langreth theorem allows to calculate A ++ in Eq. 3.30 by splitting the contour C 1 into two parts as a casual contour. Then, the first part can be calculated as a function of the retarded terms of A and B:

C 1 dt 1 A ++ (t, t 1 )B < (t 1 , t ′ ) = t -∞ dt 1 A > (t, t 1 )B < (t 1 , t ′ ) + -∞ t dt 1 A < (t, t 1 )B < (t 1 , t ′ ) = +∞ -∞ dt 1 A r (t, t 1 )B < (t 1 , t ′ ) . (3.35)
By prolongation, we can extend the integral window from -∞ to +∞. Repeating the same operation with the second part, we obtain the final expression:

F < (t, t ′ ) = +∞ -∞ dt 1 A r (t, t 1 )B < (t 1 , t ′ ) + A < (t, t 1 )B a (t 1 , t ′ ) . (3.36)
For the greater term F > , we obtain the same result by replacing "<" terms by ">".

To calculate the retarded term, we can start again with Eqs. 3.28:

F r (t, t ′ ) = θ(t ′ -t) F > (t, t ′ ) -F < (t, t ′ ) . (3.37)
Then, we get from Eq. 3.36:

F r (t, t ′ ) = θ(t -t ′ ) +∞ -∞ dt 1 A r (t, t 1 )(B > (t 1 , t ′ ) -B < (t 1 , t ′ )) + (A > (t, t 1 ) -A < (t, t 1 ))B a (t 1 , t ′ ) = θ(t -t ′ ) t -∞ dt 1 (A > (t, t 1 ) -A < (t, t 1 )) (B > (t 1 , t ′ ) -B < (t 1 , t ′ )) + t ′ -∞ dt 1 (A > (t, t 1 ) -A < (t, t 1 )) (B < (t 1 , t ′ ) -B > (t 1 , t ′ )) . (3.38)
Reversing the integral from -∞ to t ′ , we obtain finally:

F r (t, t ′ ) = t t ′ dt 1 A r (t, t 1 ) B r (t 1 , t ′ ) = ∞ -∞ dt 1 A r (t, t 1 ) B r (t 1 , t ′ ). (3.39)
For a stationary case (a constant applied voltage, for example) all the functions depend only on the difference, t -t ′ , so after Fourier transform, one gets a simple relation in the energy (or frequency) domain:

F r (E) = A r (E)B r (E). (3.40)
Using similar rules for multi-variable contour integration and one-body potential one can re-write the Dyson equation 3.31 in the energy domain for different functions as:

G r/a =g r/a + g r/a V ext G r/a + g r/a Σ r/a G r/a G < =g < + g < V ext G a + g r V ext G < + g < Σ a G a + g r Σ < G a (E) + g r Σ r G < (3.41)
Greater Green's function is given by the last equation simply replacing "<" by ">". By using this method, we can describe different types of self-energy and Green's functions. We will see in the next section that sometimes it is easier to use a simpler contour. Our goal in this section is to calculate the current, taking into account (or not) inelastic interactions between electrons and phonons. Note that several studies were done to calculate electron-boson interactions, as Aeberhard et al. [START_REF] Aeberhard | Microscopic nonequilibrium theory of quantum well solar cells[END_REF] who presented calculations in the presence of electron-photon interactions by NEGF methods. Batge et al. [START_REF] Bätge | Nonequilibrium open quantum systems with multiple bosonic and fermionic environments: A hierarchical equations of motion approach[END_REF] used time-ordering to calculate quantum electron and heat currents, taking into account interactions with bosons. Subsequently, we will focus on the inelastic terms in the lowest-order of perturbation and discuss the influence of electron-phonon interactions on the electronic charge current I L/R and the heating Q L/R in Eq. 1.1. In addition, we will present similar formalism to calculate the energy current due to phonons taking into consideration electronphonon interactions.

Current calculation

Electron current

We now derive a general expression for different currents in terms of Green's functions. The system consists of the Central region (C) (which could be a molecule, atomic chain etc.) connected to two electrodes, left L and right R. The Hamiltonian matrix has therefore the block-wise form:

H 0 + V ext =   H LL V LC 0 V CL H CC V CR 0 V RC H RR   , (3.42)
where H 0 has only diagonal blocks and describes decoupled systems and V ext couples the Central region with two leads. Total Hamiltonian can also contain manybody interactions:

Ĥ = Ĥ0 + Vext + Ŵint (3.43)
which we limit to the Central region. Different Green's function for isolated leads at equilibrium can be derived from their definitions, Eq. 3.14. In the basis of eigenstates, denoted by k, they take the form:

g < k (t, t ′ ) = - 1 iℏ f (ε k ) e i ℏ ε k (t ′ -t) g > k (t, t ′ ) = 1 iℏ (1 -f (ε k ))e i ℏ ε k (t-t ′ ) g r k (t, t ′ ) = 1 iℏ θ(t -t ′ ) e i ℏ ε k (t-t ′ ) g a k (t, t ′ ) = - 1 iℏ θ(t ′ -t) e i ℏ ε k (t ′ -t) , (3.44)
where f (E) is the Fermi-Dirac distribution function.

In the energy space it then gives:

g < k (E) = 2πi f (ε k )δ(E -ε k ) g > k (E) = 2πi 1 -f (ε k ) δ(E -ε k ) (3.45) g r/a k (E) = 1 E ± iη -ε k . (3.46)
One can see that the "<" and ">" functions can be expressed in terms of retarded and advanced functions as follows:

g < k (E) = -f (E)[g r k (E) -g a k (E)] g > k (E) = -i(1 -f (E))[g r k (E) -g a k (E)]. (3.47)
Eqs. 3.46, 3.47 are valid of course in any basis such as the basis of atomic orbitals used in Eq. 3.42 so that, in matrix notations, we have:

g r/a (E) = [E ± iη -H LL ] -1 g < (E) = -f (E)[g r (E) -g a (E)] g > (E) = -i(1 -f (E))[g r (E) -g a (E)], (3.48)
where the two last relations are valid for any system in equilibrium (due to so-called fluctuation-dissipation theorem [START_REF] Haug | Quantum kinetics in transport and optics of semiconductors[END_REF]).

The electric current flowing from the Left lead to the Central region is calculated as follows:

I L (t) = -e d NL,H (t) dt = - ie ℏ [H, NL,H (t)] , (3.49)
where the operator (in the Heisenberg representation) of the number of electrons in the Left lead is:

NL = k c † k c k = i c † i c i (3.50)
and is defined using either eigenstates k or local orbitals i of the Left lead. Non-zero contributions in the commutators come from terms coupling the Left lead and the Central region:

I L = ie ℏ j∈L,n∈C V jn c † j (t)c n (t) -V nj c † n (t)c j (t) (3.51)
and using "lesser" Green's functions:

I L = e j∈L,n∈C [G < nj (t, t)V jn -V nj G < jn (t, t)] = e • Tr[G < CL (t, t)V LC -V CL G < LC (t, t)], (3.52)
where in the last line matrix multiplication notations are assumed for shortness and the trace is taken over Central region indices.

Passing to the energy space by Fourier transformation, we have:

I L = e h dE Tr[G < CL (E)V LC -V CL G < LC (E)]. (3.53)
The G < CL and G < LC terms can be developed by Dyson's equations (3.32) and (3.31), respectively, to produce:

I L = e h dET r[(G < CC V CL g a LL + G r CC V CL g < LL )V LC -V CL (g < LL V LC G a CC + g r LL V LC G < CC )] = e h dE Tr[(G r CC -G a CC )V CL g < LL V LC -V CL (g r LL -g a LL )V LC G < CC ]. (3.54)
Using relations (3.28) between different Green's functions, one can rewrite the expression as follows:

I L = e h dE Tr[V CL g < LL V LC (G > CC -G < CC ) -V CL (g > LL -g < LL )V LC G < CC ]. (3.55)
Now, one can introduce various contact self-energies of the Central part due to its coupling to the Left electrode:

Σ r/a L (E) = V CL g r/a LL (E)V LC Σ < L (E) = V CL g < LL (E)V LC = -f L (E)V CL (g r LL (E) -g a LL (E))V LC = if L (E)Γ L (E) Σ > L (E) = V CL g > LL (E)V LC = -i(1 -f L (E))Γ L (E) (3.56)
with the coupling matrix defined as:

Γ L (E) = i[Σ r L -Σ a L ]. (3.57)
Using these "</>" self energies, the current can be expressed finally as:

I L = e h dE Tr[Σ < L (G > CC -G < CC ) -(Σ > L -Σ < L )G < CC ] = e h dE Tr[Σ < L (E)G > CC (E) -Σ > L (E)G < CC (E)] . (3.58)
Using instead of the particle operator, NL , the energy, ĤL = k ε k c † k c k or ĤL -µ L NL , one can arrive at the energy and heat currents flowing from the Left electrode:

I E L = 1 h dE E Tr[Σ < L (E)G > CC (E) -Σ > L (E)G < CC (E)] I Q L = 1 h dE (E -µ L ) Tr[Σ < L (E)G > CC (E) -Σ > L (E)G < CC (E)]. (3.59)
Lesser and greater Green's functions of the Central region, entering Eq. 3.58, can be derived from the Dyson equation (see Refs. [START_REF] Carlos | Molecular Electronics : An Introduction to Theory and Experiment[END_REF][START_REF] Haug | Quantum kinetics in transport and optics of semiconductors[END_REF]):

G </> CC = G r CC (Σ </> L + Σ </> R + Σ </> int )G a CC (3.60) with G r/a CC = 1 (E -H CC -Σ r/a L -Σ r/a R -Σ r/a int ) . (3.61)
Here, interaction self-energies, Σ int , resulting from interactions in the Central region (el-el or el-ph) are added to corresponding contact ones. Inserting Eqs. 3.60 in Eq. 3.58, the total current can be readily written as the sum of elastic and inelastic contributions:

I L = I L,elastic + I L,inelastic (3.62)
I L,elastic = e h dE(f L -f R ) Tr Γ L G r CC Γ R G a CC (3.63) I L,inelastic = e h dE Tr i(1 -f L )Γ L G r CC Σ < G a CC + if L Γ L G r CC Σ > G a CC . (3.64)
If no interactions take place in the Central region, the inelastic contribution disappears and one is left with the well known Landauer expression for the current:

I L = e h (f L -f R )T (E)dE, (3.65)
where the transmission function is calculated as:

T (E) = Tr Γ L (E)G r CC (E)Γ R (E)G a CC (E) . (3.66)
With the general Eq. 3.58, we can calculate the current, including any type of interactions occurring in the Central region.

For phonons, we define Green's functions differently, but we obtain an equivalent general expression for the current.

Phonon current calculation

For atomic vibrations the Green's functions are constructed using displacement operators (in the Heisenberg picture) [START_REF] Mingo | Anharmonic phonon flow through molecular-sized junctions[END_REF][START_REF] Gerald D Mahan | Many-particle physics[END_REF][START_REF] Galperin | Heat conduction in molecular transport junctions[END_REF]:

iℏD i ′ α ′ R ′ ,iαR (τ ′ , τ ) = ⟨T C { Xi ′ α ′ R ′ (τ ′ ) XiαR (τ )}⟩, (3.67)
where the contour-ordering of two operators is now defines as:

T C X(τ 2 ) Ŷ (t 1 ) = θ(τ 2 -τ 1 ) X(τ 2 ) Ŷ (τ 1 ) + θ(τ 1 -τ 2 ) Ŷ (τ 1 ) X(τ 2 ). (3.68)
As in the case of electrons, the Dynamical matrix, playing the role of the Hamiltonian in the NEGF formalism, is divided into blocks, where W couples the Central region with two leads as:

D 0 + D ext =   D LL W LC 0 W CL D CC W CR 0 W RC D RR   (3.69)
and possible interactions can be again considered in the Central region.

For isolated leads, governed by D LL or D RR , we can expand displacement operators in terms of phonon operators as in Eq. 2.73 which leads us to:

d i ′ α ′ R ′ ,iαR (τ ′ , τ ) = qν φ qν i ′ α ′ R ′ φ qν * iαR 2ω qν dqν (τ ′ , τ ), (3.70)
where phononic Green's functions are defined as follows:

dqν (τ ′ , τ ) = 1 i ⟨T C { Âqν (τ ′ ) Â † qν (τ )}⟩ (3.71)
and are constructed on operators:

Âqν = bqν + b † -qν . (3.72)
Non-interacting phononic Green's functions are given by:

d< qν (t ′ , t) = -i⟨ Â † qν (t) Âqν (t ′ )⟩ d> qν (t ′ , t) = -i⟨ Âqν (t ′ ) Â † qν (t)⟩ dr qν (t ′ , t) = -iθ(t ′ -t) ⟨[ Âqν (t ′ ) Â † qν (t)]⟩ da qν (t ′ , t) = iθ(t -t ′ ) ⟨[ Âqν (t ′ ) Â † qν (t)]⟩ (3.73)
and in the frequency domain:

d< qν (ω) = -2πi [(1 + n qν )δ(ω + ω qν ) + n qν δ(ω -ω qν )] d> qν (ω) = -2πi [(1 + n qν )δ(ω -ω qν ) + n qν δ(ω + ω qν )] dr,a qν (ω) = 1 ω ± iη -ω qν - 1 ω ± iη + ω qν = 2ω qν (ω ± iη) 2 -ω 2 qν , (3.74)
where n qν are occupation numbers of the mode qν given, in equilibrium, by the Bose-Einstein distribution.

One can again express "<" and ">" functions in terms of retarded and advanced functions as follows [START_REF] Galperin | Heat conduction in molecular transport junctions[END_REF]:

d< qν (ω) = F (ω)[ dr qν (ω) -da qν (ω)] d> qν (ω) = F (-ω)[ dr qν (ω) -da qν (ω)], (3.75)
where

F (ω) = n(ω), ω > 0 1 + n(|ω|), ω < 0 (3.76)
Going back to the real space, one arrives (in matrix notations) at:

d r,a (ω) = [(ω ± iη) 2 -D LL ] -1 d < (ω) = F (ω)[d r (ω) -d a (ω)] d > (ω) = F (-ω)[d r (ω) -d a (ω)] (3.77)
which are analogous to Eqs. 3.48 for electrons. We can introduce contact self-energies of the Central part due to its coupling to the Left electrode also for phonons:

Π r/a L (ω) = W CL d r/a LL (ω)W LC Π <,> L (ω) = W CL g <,> LL (ω)W LC .
(3.78)

For phonons, the energy current flowing from the Left lead to the Central region is obtained from the time derivative of the energy of the Left lead, qν ℏω qν b † qν b qν . Passing through similar steps as for electrons one can derive the following expression for phononic energy (or thermal) current:

J L = ℏω Tr[Π < L (ω)D > CC (ω) -Π > L (ω)D < CC (ω)]dω. (3.79)
The needed Green's functions are calculated from similar equations:

D </> CC = D r CC (Π </> L + Π </> R + Π </> int )D a CC (3.80) with D r/a CC = 1 ((ω ± iη) 2 -D CC -Π r/a L -Π r/a R -Π r/a int ) . (3.81)
Again, interaction self-energies, Π int , resulting from interactions in the Central region (el-ph or ph-ph) can be added to corresponding contact ones.

Electron-phonon coupling

Electrons: Hartree self-energy

In Section 3.1.2 we discussed irreducible self-energies Σ. For electron-phonon coupling two terms, Hartree and Fock (exchange), appear in the lowest order (which are second-order terms, i.e. include two electron-phonon interactions). In the case when the molecular junction is coupled to its thermal bath (thermalizing molecular vibrations) but is decoupled vibrationally from the two electrodes, only phononic Green's functions of the junction are needed. Therefore, it is convenient to work in the basis of its vibrations as in the wave packet case, Eq. 2.31. We first discuss the Hartree diagram shown in Fig. 3.4 for interaction with a specific mode λ in the junction.

Σ H (τ 1' ,τ 1 ) G (τ',τ 1' ) G (τ 1 ,τ) D λ (τ 2 ,τ 1 ) G (τ 2 ,τ 2 ) τ 2 τ 1 ,τ 1' g λ g λ
The loop in the diagram corresponds to the electronic charge, n = -iℏG < (0) so that:

Σ H (τ ′ 1 , τ 1 ) iℏ = δ(τ ′ 1 , τ 1 ) g λ iℏ C dτ 2 iℏD λ (τ 2 , τ 1 ) [-iℏG < (0)] g λ iℏ , (3.82)
where it is assumed that g λ and G < are matrices in the Wannier basis while phononic Green's function D λ is a scalar. The self-energies Σ >,< H = 0 because two times on the contour should be the same. On the other side, for Σ ++ H one has :

Σ ++ H (t ′ 1 , t 1 ) = -iℏδ(t ′ 1 , t 1 )g λ +∞ -∞ D ++ λ (t 2 , t 1 )dt 2 + -∞ +∞ D +- λ (t 2 , t 1 )dt 2 G < (0)g λ . (3.83)
Combining two integrals (coming from two contour branches) together, one gets:

Σ ++ H (t ′ 1 , t 1 ) = -iℏδ(t ′ 1 , t 1 )g λ +∞ -∞ D r λ (t 2 , t 1 )dt 2 G < (0)g λ , (3.84)
where

D r λ = D ++ λ -D +- λ .
Going to the energy domain for all the quantities and using Σ r H = Σ ++ H -Σ +- H = Σ ++ H one finally gets:

Σ r H (E) = -iℏ g λ D r λ (0) +∞ -∞ dE ′ 2π G < (E ′ )g λ . (3.85)
Summing up contributions from all the modes λ, we get the total Hartree self energy:

Σ r H (E) = -iℏ λ g λ D r λ (0) +∞ -∞ dE ′ 2π G < (E ′ )g λ . (3.86)
One can notice that the Hartree self-energy does not depend on energy and provides just the energy renormalization. It is therefore usually neglected [START_REF] Carlos | Molecular Electronics : An Introduction to Theory and Experiment[END_REF] and we will also do that. Next, we discuss the calculation of the Fock self-energy.

Electrons: Fock (exchange) self-energy

Contrary to the Hartree self-energy, the Fock self-energy is usually used in the inelastic current calculations. In this subsection, we derive the inelastic self-energies Σ we use in our code. Globally, in this subsection we will do the same step as in the previous one. The Fock diagram, shown in Figure 3.5(a), is given by: This equation comes from the Dyson equation (3.32). Like before, we have two electron-phonon interactions, but the electronic Green's function connects now different vertices. As previously, it is assumed that g λ and G are the matrices in the Wannier basis while phononic Green's function D λ is a scalar. We do not need to use the Langreth rules, but we can deduce them with a simple contour. The time on the contour corresponds to the self-energy we derive. The relations (3.26), illustrated in Figure 3.5(b), show how τ 1 and τ 2 should be placed on the contour. The two times in green, τ 1 and τ 2 , correspond to the " + +" term, and should be both on the C + part of the contour. In red, τ 1 on C -and τ 2 on C + correspond to the " < " term. In blue, τ 1 on C + and τ 2 on C -correspond to the " > " term. For the " --" case, both times are on the C -contour in purple. With this, we can determine which Green's function we have to use.

The simplest real-time self-energies are Σ < and Σ > which are just the products so that one has:

Σ < (t 2 , t 1 ) = iℏg λ G < (t 2 , t 1 ) D < λ (t 2 , t 1 )g λ Σ > (t 2 , t 1 ) = iℏg λ G > (t 2 , t 1 ) D > λ (t 2 , t 1 )g λ . (3.88)
For the retarded self-energy, we begin by using Eq. 3.28 which allows us to express it in terms of "++" and "<" terms:

Σ r (t 2 , t 1 ) = Σ ++ (t 2 , t 1 ) -Σ < (t 2 , t 1 ). (3.89)
It is easy to develop this term with the help of Figure 3.5(b):

Σ r (t 2 , t 1 ) = iℏg λ G ++ (t 2 , t 1 ) D ++ λ (t 2 , t 1 ) -G < (t 2 , t 1 ) D < λ (t 2 , t 1 ) g λ .
(3.90) One can then use again Eq. 3.89 for "++" term and arrive at:

Σ r (t 2 , t 1 ) = iℏg λ (G r (t 2 , t 1 ) + G < (t 2 , t 1 ))(D r λ (t 2 , t 1 ) + D < λ (t 2 , t 1 )) -G < (t 2 , t 1 )D < λ (t 2 , t 1 ) g λ = iℏg λ G r (t 2 , t 1 )D r λ (t 2 , t 1 ) +G r (t 2 , t 1 )D < λ (t 2 , t 1 ) + G < (t 2 , t 1 )D r λ (t 2 , t 1 ) g λ .
(3.91)

In the stationary case (dependence only on (t 2 -t 1 )), we perform the Fourier transformation to pass to the energy domain:

Σ r (E) = iℏ λ g λ +∞ -∞ dΩ 2π G r (E -Ω) D r λ (Ω) +G r (E -Ω) D < λ (Ω) + G < (E -Ω) D r λ (Ω) g λ , (3.92)
where we have also reintroduced the summation over all modes λ. For " < " and " > " self-energies, we get similar equations in the energy domain:

Σ < (E) = iℏ λ g λ +∞ -∞ dΩ 2π G < (E -Ω) D < λ (Ω) g λ Σ > (E) = iℏ λ g λ +∞ -∞ dΩ 2π G > (E -Ω) D > λ (Ω) g λ .
(3.93)

These expressions can be found in many studies of electron-phonon interactions, as in Refs. [START_REF] Haug | Quantum kinetics in transport and optics of semiconductors[END_REF][START_REF] Pourfath | The non-equilibrium Green's function method for nanoscale device simulation[END_REF][START_REF] Frederiksen | Inelastic transport theory from first principles: Methodology and application to nanoscale devices[END_REF]. Moreover, the NEGF approach can be also applied to other electron-boson interactions such as electron-photon [START_REF] Aeberhard | Microscopic nonequilibrium theory of quantum well solar cells[END_REF] in a similar way.

In our code, the calculations are done in orbital basis. To understand the element matrix calculations, we consider for the event at τ 1 the interaction between sites i and j and for the event at τ 2 the interaction between sites i ′ and j ′ . We can consider two kinds of setup. The first one, as discussed by, for example, Sergueev et al. [START_REF] Sergueev | Efficiency of thermoelectric energy conversion in biphenyl-dithiol junctions: Effect of electron-phonon interactions[END_REF], considers an isolated molecule kept in contact with its own thermal bath and can be therefore described conveniently in the phonon mode basis λ. The second situation consists of connecting the junction and leads phononically. This problem would be better treated in displacement basis with the dynamical matrix of each part and their couplings as shown in Eq. 3.69.

We can express the electron-phonon self-energy in Eqs. 3.92, 3.93 in Wannier basis for electrons and displacement basis for phonons b (including atomic site and three Cartesian directions) by using relation (3.70) between phonon Green's functions and electron-phonon couplings, see Eq. 1.40, in displacement basis. For example, the lesser term can be expressed as:

Σ < j ′ j (E) = iℏ b ′ b g b ′ j ′ i ′ +∞ -∞ dΩ 2π G < i ′ i (E -Ω) D < b ′ b (Ω) g b ij (3.94)
and similar for other self-energies.

Phonons: lowest order self-energy

In this section, we were inspired by Galperin et al. [START_REF] Galperin | Heat conduction in molecular transport junctions[END_REF] in which they calculated thermal transport with phonon Green's functions. Here, we describe how electron-phonon interactions renormalise the phonon's Green's function. Whereas for electrons there are two self-energies at the lowest order of perturbation, for phonons there is only one possibility to combine two interactions. At each interaction, phonons are created or annihilated.

Π b ′ b (τ 2 , τ 1 ) =iℏ iji ′ j ′ g b ′ j ′ i ′ G i ′ i (τ 2 , τ 1 )g b ij G jj ′ (τ 1 , τ 2 ) =iℏTr g b ′ G(τ 2 , τ 1 )g b G(τ 1 , τ 2 ) . (3.95)
We omit in the following the trace over Wannier functions to simplify equations. The first thing we can notice in this equation is that the phonon self-energy is completely defined by electron Green's functions without phonon propagators. The second thing we notice is the inversion of times of the second Green's function. Similar to electronic case, we use the contour in Figure 3.6(b) to determine which Green's function is used. We can get quite straightforwardly the expressions for "</>" self-energies:

Π < b ′ b (t 2 , t 1 ) = iℏ g b ′ G < (t 2 , t 1 )g b G > (t 1 , t 2 ) Π > b ′ b (t 2 , t 1 ) = iℏ g b ′ G > (t 2 , t 1 )g b G < (t 1 , t 2 ).
(3.96)

For the retarded term it is now convenient to use Eq. 3.30:

Π r b ′ b (t 1 , t 2 ) = θ(t 1 -t 2 ) ]Π > (t 1 , t 2 ) -Π < b ′ b (t 1 , t 2 ) . (3.97)
By developing the two terms with Eq. 3.30 we obtain:

Π r b ′ b (t 2 , t 1 ) = iℏ θ(t 2 -t 1 ) g b ′ G > (t 2 , t 1 )g b G < (t 1 , t 2 ) -g b ′ G < (t 2 , t 1 )g b G > (t 1 , t 2 ) = iℏθ(t 2 -t 1 ) g b ′ G r (t 2 , t 1 )g b G < (t 1 , t 2 ) + g b ′ G < (t 2 , t 1 )g b G a (t 1 , t 2 ) -g b ′ G a (t 2 , t 1 )g b G < (t 1 , t 2 ) -g b ′ G < (t 2 , t 1 )g b G r (t 1 , t 2 ) , (3.98)
where in the last step we used Eq. 3.29 to replace "greater" functions,

G > = G r -G a + G < .
When we develop the two last terms and, more specifically, G a (t 2 , t 1 ) and G r (t 1 , t 2 ) as:

G a (t 2 , t 1 ) = θ(t 1 -t 2 ) G < (t 2 , t 1 ) -G > (t 2 , t 1 ) G r (t 1 , t 2 ) = θ(t 1 -t 2 ) G > (t 1 , t 2 ) -G < (t 1 , t 2 ) (3.99)
we see that the two θ functions will cancel out them. The final expression of the retarded self-energy will therefore take the form:

Π r b ′ b (t 2 , t 1 ) = iℏ g b ′ G r (t 2 , t 1 )g b G < (t 1 , t 2 ) + g b ′ G < (t 2 , t 1 )g b G a (t 1 , t 2 ) . (3.100)
As before, in the stationary case, we can go to the energy space by a Fourier transform and putting back the trace notation we obtain the final expression:

Π r b ′ b (Ω) = iℏTr +∞ -∞ dE 2π g b ′ G r (E + Ω)g b G < (E) + g b ′ G < (E + Ω)g b G a (E) .
(3.101) Applying the same Fourier transformations to " < " and " > " terms, we obtain:

Π < b ′ b (Ω) = iℏTr +∞ -∞ dE 2π g b ′ G < (E + Ω)g b G > (E) Π > b ′ b (Ω) = iℏTr +∞ -∞ dE 2π g b ′ G > (E + Ω)g b G < (E) .
(3.102)

In this section, we have derived inelastic self-energies needed to calculate the currents for electrons and phonons. In principle, all the Green's functions, electronic and phononic, entering the diagrams could be calculated self-consistently in the real space basis. That would make up a complete scheme for calculating the properties of coupled electron-phonon systems (second model). We have however implemented only the first step at the moment, where phononic Green's functions in the junction are those of an isolated molecule kept in contact with its thermal bath (the first model).

Chapter 4 Test on different systems

We test various implementations in our transport code on several simple models. First, we consider 1D junctions consisting of two Ag atomic chains connected by Benzene or Vanadocene molecules to model electron-phonon or spin-spin interactions, correspondingly. Then several 2D systems, such as Black Phosphorus or Graphene with Co ad-atoms in different geometry will be discussed.

1D system: Ag wire

Ag perfect wire

We consider first a linear chain (a wire) of silver atoms, which will be used as a lead in further transport calculations. It is interesting due to the presence of only one s-like band at the Fermi energy (E F ), as we will see, and has also a gap close to E F which could be also of potential interest.

For electronic structure calculations we used the first principles DFT software Quantum-ESPRESSO (QE), discussed in Section 1.3, which is based on plane wave expansion of electronic wave functions. We build a primitive cell, composed of an Ag atom placed at the origin (0, 0, 0), which is large enough, of 10 Å, in the x and y directions, to avoid spurious chain-chain interactions. We have started by minimizing the total energy using the pw.x code with respect to the cell parameters a in the z direction and have found the equilibrium distance of about 2.6 Å.

SCF calculations were done using Perdew-Burke-Ernzerhof (PBE) functional for exchange-correlation potential, which is one of the most popular GGA parametrizations. Ultrasoft pseudopotentials were employed to describe electron-ion interactions with energy cutoffs of 30 and 300 Ry for wave functions and charge density, respectively. The mesh of 1 × 1 × 40 k-points were adopted and the smearing of 0.01 Ry was used for integrating over 1D Brillouin zone. We can determine the orbital contributions to different bands using the information from the DOS shown in Figure 4.2 calculated with a NSCF run using 1 × 1 × 80 k-points in order to obtain smooth-shape curves. With projwfc.x we can calculate the total DOS as well as the DOS projected on different Ag orbitals For Wannier calculations with the code Wannier90, we select an s-orbital in the middle of the Ag-Ag bond as well as five d-orbitals centered on Ag atoms, as illustrated in Figure 4.3. The frozen energy window spreads from -∞ to 3 eV which allows to reproduce well all the bands discussed above from the calculated Wannier Hamiltonian. The bands are calculated taking into account orbital-orbital hoppings up to 18.5 Å, corresponding to the 7th neighbour, which is especially important for correct description of the s band crossing the Fermi energy. The phonon carriers have their own energy dispersion, represented by their frequencies ω as a function of a wave vector q. 4.4) indicating at the transverse instability of a linear chain, presumably towards the formation of a zigzag configuration. In our transport code, we have also the possibility to compute the dispersion relations by diagonalizing the phonon's Hamiltonian calculated as a square root of the Dynamical matrix. We took here a cutoff for interaction distance smaller than for electrons, of 12 Å. From Figure 4.4 we can notice that we reproduce correctly the longitudinal band with the Hamiltonian while all imaginary frequencies have been set to zero in order to avoid their impact on subsequent phonon transport calculations.

After the detailed analysis of the Ag chain, we present now different types of molecular junctions where Ag chains play a role of electrodes.

Ag wire Benzene junctions

We begin with a simple benchmark case of a benzene molecule, which represents a basic molecule widely studied in literature. We mention few publications: in Ref. [START_REF] Mark | Electrode-molecule energy level offsets in a gold-benzene diamine-gold single molecule tunnel junction[END_REF] and [START_REF] Aranda | An approach to the electronic structure of molecular junctions with metal clusters of atomic thickness[END_REF] they simulated the adsorption of a benzene molecule on a Ag wire and studied the impact of charge transfer on molecular states. Ref. [START_REF] Kim | First-principles quantum transport with electron-vibration interactions: A maximally localized wannier functions approach[END_REF] reported on NEGF calculations based on DFT for Carbon atomic chain or nanotube molecular junctions with a benzene molecule. They looked at electron-phonon interactions too, which will be also the subject of our study.

Elastic transport

Inspired by the recent publication [START_REF] Li | Symmetry aspects of spin filtering in molecular junctions: Hybridization and quantum interference effects[END_REF] discussing geometry-dependent interference effects for spin-filtering in Benzene-based junctions, we compare two possibles ori-entations of Benzene connected to Ag wires as illustrated in Figure 4.5. The "flat" configuration puts the Benzene plane perpendicularly to the wire's axis while the "tilted" configuration has only one mirror symmetry with respect to the xz plane -the plane of the page in figure. These two configurations can be realized depending on the distance between Ag wires. The supercell in the z direction is made of 10 Ag atoms on the left of Benzene and 9 on the right, which is enough to match the SCF potential to the one of the perfect Ag wire on both sides. During the relaxation, we fixed two Ag atoms on both sides (contacting the leads) and relax the other Ag atoms only in the z direction. We observed that the tilting angle of Benzene increases with the distance between the Ag apex atoms. The flat and tilted configurations were obtained for apex-apex distance of 4.68 Å and 5.66 Å, respectively. The SCF calculations were done using the 1x1x4 k-mesh. To select the energy window for Wannier90 calculations, the DOS calculations are useful. The total DOS for both configurations is shown in Figure 4.6 and was calculated with 1x1x20 k points. Figure 4.2(a) and (b) present the DOS for the flat and tilted configurations, respectively. We can recognize, in both of them, Ag-related structures with the contribution of Ag d orbitals in the energy window from about -4 to -2 eV and the s-like structure around the Fermi energy. In contrast with the Ag wire, many peaks appear at lower energy due to Benzene molecular states. Figure 4.6(c),(d) present the contribution of the molecule to the total DOS, projected DOS (PDOS). We can identify peaks corresponding to molecular levels of the benzene. Molecular PDOS in both cases has a non-zero contribution around the Fermi energy, which indicated that it will participate in electron transport, which should be mediated by frontier orbitals, mostly by LUMO (lowest unoccupied molecular orbital) state located at about 1.6 eV above the Fermi level. To describe correctly LUMO states with Wannier90, the energy window [-25; 3] eV was used. To achieve a good localization of Wannier orbitals, we add two s orbitals on Ag apex atoms while for Benzene, we used 18 orbitals: π-orbitals on C, s orbitals on H, and s orbitals in the middle of C-C bonds. In both configurations, the benzene frontier orbitals hybridize with Ag s-states around the Fermi energy producing a finite PDOS at the Fermi energy, which is noticeably higher for the "tilted" configuration, indicating that probably more molecular states can hybridize in this case with Ag states.

The representation of molecular states of Benzene, shown in Figure 4.7, allows understanding of the composition of electronic wave functions on each site. We can notice that both HOMO (highest occupied molecular orbital) and LUMO states are two-fold degenerate with an energy of ϵ 0 ∓ t, respectively, where ϵ 0 is the onsite energy and -t is the nearest neighbour hopping element. In the case of the junction, we have to take into consideration the overlap between the s-orbital of Ag atoms and molecular states of the Benzene. For the "flat" configuration, only HOMO-1 state can have a non-zero overlap because of opposite alternating phases of electronic wave functions for other states on the Benzene plane (cancelling out). For the tilted configuration, the axial symmetry with respect to the z axis is broken. It allows therefore to have finite overlaps for one of HOMO (the one on the right) and one of LUMO (the one on the left) states. Figure 4.8(b) present transmission functions for two configurations of Benzene calculated by Green's function method discussed in Section 3.2 with the Hamiltonian constructed on Wannier's functions. In black, the ideal transmission of Ag infinite atomic chain is shown given by the number of bands at specific energy, it is integer and is always higher than the transmission of a molecular junction. For the "flat" configuration, as was argued above, only one symmetric Benzene orbital (denoted as "1") can mediate the transport from the Ag s-band, resulting in a finite tunnelling transmission at the Fermi energy since this orbital lies well below the Fermi energy (at around -5 eV). For the "tilted" configuration, on the contrary, two more molecular orbitals, denoted as "2" and "3" and coming from HOMO and LUMO doublets, get opened producing a destructive interference with the orbital "1". This is translated to the very small transmission in a wide range of energies around the Fermi level. Near 1.6 eV, the transmission goes to 1 due to the resonant tunnelling through the LUMO orbital "3".

In order to confirm this interpretation of NEGF results, we perform a realtime dynamics of electronic wave packets, shown in Figure 4.9. We propagate the wave packet constructed from the Bloch function at the Fermi energy. The corresponding k-point was k z = 0.25 (in units of 2π a ). The group velocity of an electron was found to be about V z ≈ 640 km/s. At t 0 , the wave packet is coming from the left lead. At t 1 , an electron arrives at the junction. The noisy peaks at this point are due to interference between the electron wave function coming from the left and the reflected from the junction part. The contributions from all Benzene molecular orbitals (18 altogether), shown by columns in the middle of panels, allow to confirm previous analysis. In the "flat" case, an electron passes only through the HOMO-1 state, while in the "tilted" case, the three orbitals are seen to be involved. These three pathways produce destructive inference, discussed above, resulting in complete reflection of an electron as in Ref. [START_REF] Sengul | Electrode effects on the observability of destructive quantum interference in single-molecule junctions[END_REF]. The wave packet can be also studied at the LUMO energy (≈ 1.6 eV) to reveal the origins of the transmission peak in the tilted case. We have seen that an electron has a negligible probability to pass through HOMO-1 or HOMO states at this energy, and the main contribution comes from the LUMO state.

At t 2 , an electron has passed the junction and the transmitted and reflected portions of the wave packet should be related to transmission (T ) and reflection (R) probabilities, respectively. They can be estimated as follows:

T = i∈R |ϕ(i)| 2 ; R = i∈L |ϕ(i)| 2 , (4.1)
where ϕ(i) is the electron wave function at the site i. consistency between the two methods: from NEGF calculations, one had at the Fermi energy T ≈ 0.62 and R ≈ 0.38 for the "flat" case and T ≈ 1.00 and R ≈ 0.00 for the "tilted" one. As we have a unique s-band crossing the Fermi energy, the wave packet propagated at the corresponding k-point gives the total transmission as provided by NEGF in Figure 4.8.

On this simple example, we can see a utility of wave packet method. The fact that we can observe the propagation of electrons across the junction in real time and analyze contributions from different molecular orbitals is complementary to the NEGF approach, which allows to obtain the total transmission as a function of energy.

Electron-phonon interactions

We now add electron-phonon interactions. We consider first the isolated Benzene molecule. From EPW calculations we can extract the Dynamical matrix, the electronic Hamiltonian in Wannier basis and electron-phonon coupling matrix in real space. We then diagonalize both the Dynamical matrix and the Hamiltonian and rotate the coupling matrix in order to get electron-phonon elements in the basis of molecular vibrations and molecular orbitals.

In order to identify the most important for transport vibrational mode, we plot in Figure 4.10 for each mode its coupling strength to the LUMO orbital (any of two of which were found to be equal) which is expected to play the most important role in electron transport around the Fermi energy. The mode which interacts most configuration Reflection in junction transmission total flat 0.61781 0.0000 0.38219 1.0000 tilted 0.99733 0.0000 0.00267 1.0000 Table 4.1: Contributions of the final wave packet in the two leads for flat and tilted configurations, which can be related to the transmission and reflection coefficients.

Flat configuration

Tilted configuration X 5000 X 5000 X 5000 X 5000

Figure 4.11: Evolution of electronic wave packet with inclusion of electron-phonon interaction term with the mode "30" of the Benzene. Again, the lines on the left and right sides as well as the columns in the middle describe the contributions to the left and right chains and on molecular levels, correspondingly.

strongly with the LUMO (marked by the red circle) is the mode "30" characterized by its energy of 0.197 eV. For this mode, which we can denote as the mode "0", we add to the electron Hamiltonian the corresponding electron-phonon term as in Eq. (2.32):

Ĥep ij,0 = g λ=30 ij ( b † 0 + b0 )ĉ † i ĉj , (4.2)
where ij run over (all) molecular orbitals. Figure 4.11 shows the wave packet propagation in the presence of electronphonon interactions with the chosen mode considering the simple case of one phonon channel in the model, Figure 2.4. The probability function in the upper channel (with one excited phonon) is multiplied by 5000. At the time t 0 , the wave packet comes from the left lead and no phonon is present in the molecule initially. At the time t 1 , an electron arrives at the molecule and starts to interact with the mode.

From the later time t 2 , one can get the overall probability of phonon excitation. The total wave function is of the form:

|Ψ⟩ = i ϕ 0 (i) |0⟩ + i ϕ 1 (i) |1⟩ , (4.3)
where |0⟩ and |1⟩ are phononic states without phonon and with one phonon, respectively. The probability to excite a phonon is then given by the norm of the wave function in the excitation channel:

P 1 = i |ϕ 1 (i)| 2 . (4.4)
One can see that the "flat" configuration has a negligible probability to excite a phonon. On the contrary, this probability is much higher, of the order of ≈ 10 -4 , for the "tilted" configuration. It is consistent with our previous finding: since an electron does not pass across the LUMO orbital in the "flat" configuration due to symmetry considerations, it can not neither excite a phonon. In the "tilted" case, the LUMO is active, and the phonon can be now excited. We performed similar calculations for all the modes, Figure 4.12, in order to see if the probability to excite a phonon correlates with electron-phonon coupling calculated for the LUMO orbital, and to verify if the mode "30" gives indeed the highest excitation rate. For the "tilted" case, some of the modes get a non-zero excitation rate, but the Highest rate is still related to the mode "30". In the "flat" case, as the LUMO state is not acting in the electron transport, relevant, and we find indeed that the only mode which can be excited is the mode "18". This mode is rotationally symmetric and therefore should be activated by an electron passing through the HOMO-1 orbital, which is the only one participating in the transport. Until now, the electron-phonon couplings were extracted from the isolated Benzene molecule. We compare now these calculations with the electron-phonon calculations of the full Ag/Benzene junction. Figure 4.13 presents in black the frequencies of all 93 modes of the full system (containing 19 Ag atoms and 12 Benzene atoms) for each geometry as well as those obtained from the sub-matrix of the total Dynamical matrix, considering only displacements of Benzene atoms (in red). The highest in frequency modes are matching perfectly for two types of calculations. These modes are intrinsic Benzene modes involving relative displacements of Benzene atoms. The 6 highest modes corresponds to Hydrogen oscillations. On the contrary, Ag atoms are much heavier than Carbon or Hydrogen. Therefore, the modes corresponding to Ag wires are situated at lower frequencies, contributing to a hump going from 0 to 30 meV. Between these extreme cases, the modes are composed of displacements of silver and benzene atoms. In the case of isolated benzene, the 6 lowest-frequency modes, corresponding to rigid displacements and rotations of the Benzene as a whole, should in fact have zero frequencies since do not require any energy. When the Benzene is connected to Ag wires, some of these modes get non-zero frequencies due to interaction with Ag contacts. In the "flat" case, the contacts activate a "breathing" mode with a rigid oscillation of the molecule along the z direction. For the "tilted" one -a "swing" oscillation in the zy plane as shown in Figure 4.13. We can notice that these oscillations fall into the frequency window of the Ag wire vibrations, of [0, 0.027] eV, so they should be active for phonon transport across the molecule, which as we will see is indeed the case.

Flat configuration

Tilted configuration 4.14 displays the excitation rates for those phonons calculated from the molecule sub-matrix. For the "tilted "case, the mode "30" (marked as "5" in figure) still has significant excitation rate, close in fact to the one found previously for the isolated molecule (see Figure 4.12). However, two other modes, the mode "11" and especially the "swinging" mode (marked as "3") strongly benefit from the coupling to Ag wires and may get highly excited. For the "flat" configuration, one finds two modes with relatively high excitation rates -the "breathing" mode (marked as "1") and another one ("2"). One can conclude, therefore, that for both configurations the most interesting modes, easy to be excited, are the modes activated by interaction with Ag chains -the "breathing" and the "swinging" modes. The electron transport for Ag/Benzene junctions calculated with our code has brought therefore consistent results. We test now the phonon transport.

Phonon transport

We calculate first the phonon transmission function by NEGF method. Figure 4.15 presents the phonon's transmission for both geometries, together with the phononic DOS projected on Benzene and different displacements of the apex Ag atom. The "flat" configuration has a wide transmission feature, centred at 7 meV, with transmission close to 1. The corresponding PDOS indicates that this feature should be attributed to the "breathing" mode. Similarly, for the "tilted" configuration, the peak at 12 meV should be associated to the "swinging" mode. To confirm the role of these low-frequency vibrations, we propagate the phononic wave packets at specific energies chosen close to transmission function maxima. The results are shown in Figure 4.16. In the "flat" case, we choose the wave vector q z = 0.076 ( 2π a ) corresponding to the energy ℏω = 7 meV, which gives the group velocity of V z ≈ 1102 m/s. Observing the displacement of the Benzene atoms, we can confirm that the main contribution comes from the "breathing" mode. In the "tilted" case, we work with the wave vector of q z = 0.126 ( 2π a ) at the energy ℏω = 12 meV, which correspond to the velocity V z ≈ 6622 m/s. We see again that the main vibration mediating the wave packet is the "swinging" mode.

NEGF calculations

We now present the calculations of various current-voltage characteristics. They are most relevant since they can be directly compared with experimental data. First, we will consider the elastic regime without electron-phonon interactions. As the "tilted" configuration has a very low transmission around the Fermi energy, we will focus most on the "flat" geometry. Figure 4.17 symmetrically too, +V /2 and -V /2 on the left and right leads, whereas molecular levels energies were kept at their equilibrium values.

The charge current (in green) evolves linearly as a function of voltage. As it is an integral of the transmission between ±V /2 and because the transmission for the "flat" case can be considered as linear around the Fermi energy (Figure 4.8), the current is linear until the voltage reaches the edge of the conduction band at -1.2 eV below the Fermi energy. After this, the charge current saturates due to the gap in the band structure of Ag wires. The associated heating of the left lead (in red) and of the right one (in black) are not symmetric too, describing how much each lead receives the thermal energy due to electron thermalization to the corresponding chemical potential.

One can also model some interesting thermoelectric effects. The Joule effect refers to the heat dissipation that appears when an electric current passes through the resistive junction. Classically, it is due to collisions of electrons with atoms or molecules. These collisions cause the electrons to lose their kinetic energy, resulting in the generation of heat. The greater the resistance of the material, the more collisions occur, and thus, more heat is produced. The amount of the total heat Q generated in the system is given by the Joule's law:

Q = I × V = V 2 /R. (4.5)
Without inelastic effects, the Joule heat should be equal to the sum of two contributions of the left and right heating (from Eq. (3.59)).

Q = Q L + Q R = - 1 h dE (E -eV /2)(f L -f R )T (E, V ) + dE (E + eV /2)(f R -f L )T (E, V ) = IV. (4.6)
Adding two contributions in the Left and Right leads, we indeed recover the Joule's law as shown in Figure 4.17(b).

The most important thermoelectric effect is the Seebeck effect. It consists in the voltage generation due to temperature gradient between two electrodes. In general, electrons tend to diffuse from the hotter lead to the colder one. This movement of electrons is translated into induced voltage between the two electrodes needed to stop electron's flux. The magnitude and the sign of the generated voltage depends on the properties of the material. The Seebeck coefficient represents the voltage induced by the gradient of temperature between the two regions. In practical applications, the Seebeck effect is commonly used in thermoelectric devices, such as thermocouples and thermoelectric generators. A thermocouple consists of two different metals or semiconductors linked together by a junction. When a junction connects a higher temperature reservoir with a lower temperature one, a voltage is generated across the thermocouple. Thermoelectric generators are devices that convert heat energy directly into electrical energy. They typically consist of an array of thermocouples connected in series. The Seebeck effect has applications in various fields, including power generation or waste heat recovery. It can be used to power electronic devices or in charge batteries.

Figure 4.18(a) helps to understand the phenomenon. The gradient of temperature induces a difference in electron occupations in left and right leads. If the temperature of the left lead T L is higher than the temperature of the right one T R (in figure, T R = 0 K), electrons above the Fermi level (equal for two leads) will propagate from the left to the right. On the contrary, electrons below the Fermi level will propagate in the opposite direction -from the right to the left. The two currents will cancel out if the transmission function were constant around the Fermi energy. If not, the net non-zero current will establish in the system, of the sign depending on the slope of the transmission function. shows the calculated current at zero voltage as a function of temperature difference between electrodes. As we observe a slow variation of the transmission in the "flat" configuration, the effect is very small too. Notice that the induced current is negative, which means that it is mediated by holes. It is due to the fact that the transmission function has a negative (though quite small) slope around the Fermi energy as seen in Figure 4.8(b). Another important effect is the Peltier effect, which is reciprocal to the Seebeck one. This phenomenon is related to the cooling of one electrode by passing an electric current while keeping both electrodes at the same temperature. The primary application of the Peltier effect is in thermoelectric cooling and heating systems. By using the Peltier effect, these systems can provide cooling or heating abilities without the need for refrigerants or mechanical compressors. They are commonly used for cooling electronic components, such as computer chips, laser diodes or sensors. Figure 4.18(b) illustrates the effect. Two electrodes are kept at the same temperature while an applied voltage shifts their chemical potentials µ, in figure µ L > µ R . The heat currents into the electrodes are defined by Eq. (1.1). Since f L -f R > 0 electrons will propagate from the left to the right at any energy. For the left electrode, only electrons leaving it at E > µ L will produce the cooling (shown in blue) since the remaining holes should be filled from the chemical potential µ L . At E < µ L , the thermalization will cause the heating of the electrode (red region). For the right electrode, on the contrary, only electrons arriving at E < µ R will produce the cooling (shown in blue) since they should thermalize to higher lying µ R . At E > µ R , the thermalization will cause the heating of the electrode (red region).

As seen from Figure 4.18(b) only a little fraction of electrons will contribute to the Peltier cooling in electrodes. That is why the transmission should have a huge variation around the Fermi energy in order to bias enough the cooling contribution. In our case of rather constant in energy transmission, as seen in Figure 4.18(b), we can not observe Peltier effect -both leads are heating positively at any voltage. To enhance the Peltier effect, a possible solution would be to apply a gate on the benzene molecule to bring a specific molecular orbital close to the Fermi energy. We model this situation for the "tilted" configuration, where the important LUMO orbital is positioned at 1.6 eV above the Fermi energy. Figure 4.19(a) presents the transmission function at applied gate of -1.4 Volt. This will shift all molecular levels. We can observe in particular the peak of transmission from LUMO shifted towards the Fermi energy with a strong variation in energy. Figure 4. 19(b) gives the heating powers calculated for this situation. One can see indeed that at small positive voltage the left lead is getting cooled while the right one is heated. For negative voltages, we have an inverse situation. The strong asymmetry in transmission, due to the LUMO placed closely to the Fermi level, allows therefore to amplify significantly the Peltier effect. We can also calculate the thermal flux from elastic current calculations. The thermal flux is the heat transport induced by a gradient of temperatures at zero applied voltage, similar to Seebeck effect. The different types of carriers can transfer the thermal energy. Figure 4.20 compares the thermal flux carried by electrons and by phonons. For electrons, the thermal flux is represented by the heating of one electrode and the cooling of another one, the sum of the two being zero since IV = 0. For phonons, it is characterized by the energy current. We see that the electronic contribution is in the range of nW atts while it is in pW atts for phonons. The factor of 1000 between the two thermal currents indicates that electrons conduct the heat much more efficiently than vibrations in the considered system.

Finally, we discuss the effect of electron-phonon interactions. The electron- phonon coupling in NEGF formalism can be treated with two approximations. First, we can consider a junction connected to two left and right phonon reservoirs (leads). The second model consist of connecting a junction to another thermal bath [START_REF] Sergueev | Efficiency of thermoelectric energy conversion in biphenyl-dithiol junctions: Effect of electron-phonon interactions[END_REF] at a certain temperature in order to thermalize the junction vibrations. The first model is not yet implemented, so we will discuss here the second (more simple) approach. Figure 4.21(a) and (b) present the differential conductance, both elastic and the total one including inelastic corrections, for the "flat" and "tilted" cases. The elastic conductance of the "flat" configuration is higher than for the "tilted" case because of the strong difference in transmission around the Fermi energy. In the "flat" case, the inelastic conductance is slightly shifted with respect to the elastic one. This shift comes mainly from the elastic correction [START_REF] Carlos | Molecular Electronics : An Introduction to Theory and Experiment[END_REF],

δI elastic , due to possibility for an electron to take energy from phonons and to give it back in the junction without changing its final energy. For the tilted case, the conductance shows steps at specific voltages corresponding to frequencies of vibration modes presented previously. Such steps in conductance are often obtained in simulations (for example, in Refs. [START_REF] Frederiksen | Inelastic transport theory from first principles: Methodology and application to nanoscale devices[END_REF][START_REF] Sergueev | Electron-phonon interactions in molecular electronic devices[END_REF][START_REF] Woo Youn | Application of quantum chemistry to nanotechnology: electron and spin transport in molecular devices[END_REF][START_REF] Galperin | Resonant inelastic tunneling in molecular junctions[END_REF]) and experimentally [START_REF] Troisi | Molecular signatures in the transport properties of molecular wire junctions: what makes a junction "molecular[END_REF][START_REF] Troisi | Molecular signatures in the transport properties of molecular wire junctions: what makes a junction "molecular[END_REF]. The steps are attributed to the opening of another transport channel with an excited phonon. We observe that the first step is the highest one and corresponds to the "swinging" vibrational mode. The height of this step confirms the importance of the mode in inelastic transport, as it has been also seen in the simulations with wave packets. The step at 0.2 eV is associated with the mode 30 of the benzene.

Finally, we can also calculate the amount of energy transferred to the molecule. In the elastic case, the total energy current to the molecule is always zero, I E L +I E R = 0 because of the energy conservation. In the inelastic case, it is not the case producing the heating of the molecule, Q = I E L + I E R . Figure 4.21(c) demonstrates that the heating of the benzene is more important in the "tiled" configuration, which is consistent with our wave packet studies showing more important phonon excitation rates for this case.

This section allowed to demonstrate a majority of functionalities of our transport code on a simple one dimensional system. We simulated the transmission and the current/voltage characteristics with the NEGF method. The wave packet method confirmed NEGF results and allowed to get insight and more details on molecular state resolution and time dynamics. We include the electron-phonon interactions in the current and in the wave packet dynamics. We will discuss now another kind of interaction which has been implemented in the code, namely magnetic exchange interactions between the local spin and the spin of the conduction electron.

Ag wire with a Vanadocene

Metallocenes are a class of organometallic molecules which consist of a transition metal atom trapped between two cyclopentadienyls (C5H5). The vanadocene belongs to metallocene group, it contains the vanadium as the central atom. The vanadocene exhibits interesting electronic and magnetic properties due to the presence of the vanadium atom and has the total spin S = 3/2.

The Ag/Vanadocene junctions have been recently studied in our group in collaboration with experimentalists in view of single-molecule spin-filters based on quantum interference effects [START_REF] Nath Pal | Nonmagnetic single-molecule spin-filter based on quantum interference[END_REF]. We took therefore a model junction consisting of Ag wires connected by a vanadocene molecule, shown in Figure 4.22, as a first simple system to test spin-spin interactions implemented in our transport code. We have started by performing the structural relaxation using spin-polarized cal- culations, where spin moment was initialized at the beginning only on the Vanadium atom. In output, the total magnetization of the system was found to be about 2.52µ B , which means that the system roughly preserves the total spin of the Vanadocene molecule, which is S = 3/2. The small difference comes from (negative) co-polarization of Ag chains.

The SCF calculations were done with PBE exchange-correlation potential (GGA approximation) using 1x1x6 k-points. We used 1x1x2 k-mesh in the NSCF calculation needed to construct Wannier functions. Another NSCF calculation was done with a finer mesh of 1x1x40 in order to obtain the smooth DOS shown in Figure 4.23. Figure 4.23(a) presents the DOS calculated with projwfc.x and allows to compare the total DOS (in black) with the PDOS projected on the Vanadocene molecule (in red). We observe that only Vanadocene orbitals contribute to the differences between spin up and spin down channels. The cyan dashed curve correspond to the DOS calculated with Wannier90. This DOS matches perfectly to the one calculated with projwfc.x in the energy window [-6; 3] eV indicating that we describe well all electronic states around the Fermi energy with the Wannier basis. We select the orbitals for Wannier calculations by analyzing As previously, we first calculate the spin-resolved transmission function, which is shown in Figure 4.24. One can observe that spin up electrons do not transmit at the Fermi level while spin down transmission reaches almost the value of 1 which corresponds to spin-dependent PDOS of d z 2 orbital of the Vanadium which is expected to dominate the transport due to its orientation and strong coupling to leads. One can see that main contributions come from d-orbitals, while contributions of s-and p-orbitals are relatively small. For the Wannier basis we choose therefore five d-orbitals on Vanadium, and we add, as for the Benzene before, π-orbitals on C atoms and s-orbitals on H and on each C-C bond. Altogether, it makes therefore 35 orbitals for Vanadocene molecule. We have found that 3 eV above the Fermi level is enough for the frozen energy window. The perfect matching of both DOS curves in the selected energy window, discussed above, confirms our choice of Wannier functions. Figure 4.23(c) presents the PDOS on different d-orbitals. For spin up, d z 2 , d xz and d yz states are occupied while for spin down only a small portion of the d z 2 is below the Fermi energy. Thus, we count 3 states occupied, which results in the total spin close to S = 3/2. Moreover, the d z 2 orbital (relatively strongly hybridized with Ag chains) is located right above the Fermi energy for spin down, while it is down shifted for the spin up case. This energy shift should play a role in the electron transport around the Fermi energy, as we will see now.

We can study an electron propagation through the Vanadocene with the wave packet method. Figure 4.25 shows wave packet propagation at the Fermi energy for both spins independently. The group velocity was found to be about 650 km/s. At the time t 0 , the wave packet is coming from the left lead and reaches the junction at t 1 . We see that in both spin cases, the d z 2 Vanadium orbital is strongly contributing to the electron transport. In the spin down case, the electron can pass through with a high probability. For spin up, the wave packet is totally reflected. spin states Reflection in junction transmission total ↑, The transmission and reflection probabilities extracted from wave packet simulations are presented in Table 4.2 and matches well to NEGF transmissions. The transmission in spin up case is very close to 0, while the spin down one is about 0.87 where with the NEGF method we obtained 0.86. These results agree pretty well with NEGF transmission curves presented in Figure 4.24.

Once the spin up and spin down Hamiltonians are calculated with DFT and transformed to Wannier basis, we can construct H 0 and J matrices from Eqs. (2.43) and (2.44) knowing the local spin S = 3/2. We can verify if exchange parameters J are localized on Vanadocene and what are their spread into the junction. Fig 4 .26 shows diagonal elements of J. In the lead part, one can observe some small values which are negligible compared to those on the molecule and specially on Vanadium d orbitals. We can thus deduce that the spin exchange events will occur mainly near the Vanadium, which is quite reasonable.

In order to verify if our calculations with spin exchange are consistent, we can propagate wave packets with H 0 and J but deactivating the spin-flip terms in the total Hamiltonian (2.40). We observe exactly the same evolution as in spin-polarized DFT simulations presented in Figure (4.25).

There are no spin-flip events for |↑, 3/2⟩ and |↓, -3/2⟩ initial states where we should recover the DFT result. On the contrary, the state |↓, 3/2⟩ can mix with the state |↑, 1/2⟩ by rising and lowering operators in Eq. ((2.40)).

Figure 4.27 presents the evolution of wave packets including this interaction. We initialize an electron coming from the left lead with a spin state |↓, 3/2⟩ at t 0 . At t 1 , the electron interact with the local spin Of the Vanadocene.

At the time t 2 , an electron is reflected or transmitted with a spin up or down. We could also take another local spin state rather than S = 3/2 for the initial wave packet. Table 4.4 represents spin-flip probabilities for different combinations of conduction electron and local spins. We can see that the probabilities are quite close and are moreover symmetric with respect to the change of sign of two spins, as it should be. The most probable spin-flip is found for the pair of states |↑, -1/2⟩ and |↓, 1/2⟩.

We have tested electron-phonon and spin-spin interactions on model systems based on Ag wire leads. We now discuss our first applications to several 2D materials. 4.2 2D system: Black phosphorus

Effect of bonding and antibonding character on STS spectra

We introduce the 2D materials with a very representative example of Black Phosphorus (BP). One of the 2D layer-structured materials which has attracted recently a great attention [START_REF] Jv Riffle | Impact of vacancies on electronic properties of black phosphorus probed by stm[END_REF][START_REF] Kiraly | Probing single vacancies in black phosphorus at the atomic level[END_REF][START_REF] Zhang | Surface structures of black phosphorus investigated with scanning tunneling microscopy[END_REF] due to possible superconductivity at low temperature, high carrier mobility at room temperature or his semiconductor properties with a tunable band gap [START_REF] Ling | The renaissance of black phosphorus[END_REF][START_REF] Guo | Linear scanning tunneling spectroscopy over a large energy range in black phosphorus[END_REF]. The single layer of BP is named as Phosphorene. We construct its primitive cell with the lattice parameters of a = 3.348 Å and b = 4.587 Å along the x and y direction, respectively. The vacuum region of 20 Å was used to separate the slabs in the z direction and to avoid their artificial interaction. This cell contains 4 Phosphorus atoms (two atoms per sublayer). The DFT calculations were performed within the PBE parametrization for the exchange-correlation functional and using cutoffs of 30 and 300 Ry for wave functions and charge density, respectively.

The band structure of Phosphorene is shown in To investigate the effect of B/AB character of states on their extension out of BP surface and therefore on observed STS spectra, we have performed the calculation of a slab of ten Phosphorene layers stacked in an ABA ... sequence along the z direction. The B layer is shifted by half of a unit-cell period along the X axis with respect to the A layer, which corresponds to the minimum energy configuration. Figure 4.29(a) shows the total DOS integrated over the full 2D Brillouin zone for a 10 layer slab. One sees that the p z -DOS displays a rather symmetric shape to the right and to the left of the gap. On the contrary, the LDOS in the vacuum, Figure 4.29(b), is clearly asymmetric, with a large increase at positive energies where AB states are present compared to negative energies dominated by B states. This asymmetry is a robust result and becomes even more pronounced with increasing the distance from the surface, from 4 to 6 Å.

Experimentally, this asymmetry is observed by STM, as shown in Fig. 4.29(d). The dI/dV spectrum of bulk BP around the Fermi level exhibits a 0.3 V gap corresponding to the gap of bulk BP [START_REF] Kiraly | Probing single vacancies in black phosphorus at the atomic level[END_REF][START_REF] Qiu | Resolving the spatial structures of bound hole states in black phosphorus[END_REF][START_REF] Bouatou | Black phosphorus for directed molecular assembly with weak electronic coupling[END_REF]. Outside the gap, the dI/dV signal exhibits a stronger increase in the conduction band than in the valence band. At larger sample bias this asymmetry still holds with a much larger signal at positive voltage than at negative voltage. This experimental spectrum compares very well with the calculated LDOS, which therefore confirms our theoretical predictions. The discussed here results of this joint work have been published recently in Ref. [START_REF] Lahrache | Effect of bonding and antibonding character of electronic states on their tunneling spectra[END_REF]. I note that Phosphorene has been considered as one of possible 2D systems (along with of course Graphene) to test our transport code for 2D wave packets dynamics. We proceed now with a junction built on 2D Phosphorene. The idea was to take advantage of the armchair structure of Phosphorene in the z direction and to deposit a line of Co atoms along z as shown in Figure 4.30. We want then to study the transport in the y direction and to see, in particular, if it displays any spin-polarization. The relaxation of the junction was done by fixing all P atoms and using a mesh of 1 × 3 × 12 k-points. We have found a total magnetization of 1.02 µ B . Looking at the DOS, shows the PDOS on different d-orbitals of Co. Among them, we observe that only the d z 2 -orbital (for spin down) is located above the Fermi level. For spin up states, all the d-orbitals are placed at the edge of the Phosphorene valence band. This will result in spin unbalance corresponding to the total magnetic moment of about 1 µ B which is found in calculations. Moreover, the d-orbitals for spin down are located in the gap region of Phosphorene, [-0.6; 0.3] eV. In the spin up case, on the contrary, no Co orbital contributes to the DOS at the Fermi energy. This indicates at the full spin-polarization of Co states at the Fermi level, which could be important for spin-filter properties. We have therefore calculated the spin-polarized transmission in the y direction perpendicular to the Co chain, which is shown in not find any strong spin-polarization, there is a partial effect at the valence band edge, at energies around -0.9 eV, but it is not so pronounced. We can argue, however, that due to fully spin-polarized Co states at the Fermi energy one can expect good spin-filter transport properties, but in the Co chain direction z. Since the Fermi energy falls into the gap of Phosphorene these Co states are inactive for the transport in the y direction where the transmission is always zero in the gap. To check this point, we present in Figure 4.33 the band structure of the junction in the z direction and compare it with the one for freestanding Co chain, calculated with the same interatomic Co-Co distance along the z direction. We can see that for the Co wire, the spin up has only the s band crossing the Fermi level while for spin down -three bands. For deposited Co wire, we find that only one spin down band intersects the Fermi energy which should provide therefore fully spin-polarized transport channel, as was argued above.

Phosphorene with a Co linear chain

In this section, we showed an example of transport calculations for 2D material where however translation symmetry is preserved along one direction, along the z in the case of Co chains on Phosphorene. We will consider now another 2D case where translation symmetry is missing in both in-plane directions.

Graphene with a Cobalt ad-atom

Graphene is a well studied 2D material constituted of a single layer of carbon atoms arranged in a two-dimensional honeycomb lattice. It is considered as a very promising system due to its exceptional properties such as an excellent electrical conductivity because of the 0-band gap, chemical stability and resistance to chemical reactions, transparency, good thermal conductivity, etc. Therefore, we have started with the Graphene as a basic material to simulate a propagation of wave packets in two dimensions. First, we have calculated the lattice parameter by minimizing the total energy, and we have found a value of a = 2.44 Å. behaviour. The p x (in purple), p y (in cyan dot) and s (in green) orbitals are not contributing to electronic states around the Fermi energy, where only p z -orbitals (in red) are important. Note that since the unit cell is composed of two carbons, the PDOS should be multiplied by 2 to produce a perfect matching to the total DOS.

The choice of Wannier orbitals is multiple. One can use only p z orbitals on carbons which will allow reproducing well the band structure in a small range of energies around the Fermi energy, [-3; 3] eV. However, in the perspective to implement a magnetic scattering centers, it is better to include more orbitals in order to describe well electronic states at lower energies too. The sp 2 orbitals converged very well in Wannier calculations, but the band structure wasn't well reproduced, we did not get, in particular, the zero band gap. We choose therefore to use s-orbitals in the middle of each C-C bond in addition to Carbon p z -orbitals. This choice allowed to reproduce correctly the QE bands in a large energy window, [-25; 3] eV, as shown by the red curve in Figure 4.34(a), with reduced number of orbitals per unit cell. The calculations were done with an orbital-orbital cutoff of 6.5 Å, which corresponds to 2nd cell neighbours on the Graphene lattice. With our transport code, we can calculate first the 2D band structure along a specific k-path or on some 2D mesh of k-points and then to choose a specific Bloch state to be propagated at the energy of interest. Figure 4.35(a) is obtained with a grid of k-points centred at K-point. Then we filter out the k-points around the energy E = 0.2 eV with a range of ±0.01 eV which results in Figure 4. 35(b). For this fixed energy cut at E = 0.2 eV, the Dirac cone produces a circle, where each k-point has its own group velocity, which can be calculated using equations like Eq. (2.9) and should also correspond to V g ∼ ∇ k E(k). In the case of Dirac cones, the velocity direction should follow the radial direction going outside the circle. Our result confirms this and allows selecting a propagation direction for our wave packet.

From the lead, we built our junction supercell with 7 × 6 unit cells of graphene and put a Co atom at the hollow position above the carbon hexagon as illustrated in Figure 4.36. This junction will be inserted in a larger plane of graphene, considered as a lead far away from the Co in transport calculations. As the interactions involve two neighbours in the lead (defined by our cutoff parameter) we consider two unit cells on the perimeter of the junction as being unperturbed and therefore matching to ideal graphene. The Co has an electronic configuration [Ar]3d 7 4s 2 , which corresponds to the spin S = 3/2. We relaxed the junction by fixing all the carbon atoms except of those closest to the Co ad-atom, and observed that the Co approaches the graphene layer before taking up a final position above the middle of the hexagon.

From SCF calculations, done with 2 × 2 k-points, we obtain the total magneti-

DOS (arb. unit)

.. As in the case of BP/Co discussed above, in order to understand which orbitals of Co we have to take into account and which energy window is necessary in Wannier calculations, we plot the DOS in a larger energy window as shown in Figure 4.38. Lowest levels are composed of 3s and 3p orbitals of Co. On the contrary, the 4s and 3d Co orbitals are contributing around the Fermi energy. One can guess that carbon orbitals start to contribute around -20 eV. Thus, we take all electronic states in an energy window from -20 eV to 15 eV above the Fermi energy to be used for constructing Wannier orbitals. The Wannier calculations were done, as before, with a frozen window up to ±3 eV above the Fermi energy, and we choose as Wannier functions p z Carbon orbitals and s-orbitals at C-C bonds, as was discussed before for ideal Graphene, as well as five Co d-orbitals.

s QE DOS Wannier90

We show in Figure 4.39 the propagation of a 2D wave packet for spin up and down electrons (firstly, without spin-flip terms). This is quite heavy calculation because of huge number of unit cells. In 1D models, presented before, we worked with 800 unit cells of Ag wire along the z direction with 6 orbitals per cell. In Figure 4.39, a simulation grid was 300 × 300 unit cells with 5 orbitals per unit cell of Graphene. This model allows to see difficulties of wave packet propagation, specifically the dispersion of the wave packet in time if it is not wide enough (due to slight nonlinearity of the energy dispersion). We note however that the current implementation of the code may be still significantly optimized and moreover parallelized (it is a serial implementation at the moment) in the future. For spin up and spin down cases, wave packets start to propagate at the same position at t 0 . They propagate towards the Co ad-atom and interact with the perturbed cell at t 1 . After the interaction, the wave packet tends to diffuse in time. Very small difference is seen between the two spin cases with a slightly larger portion localized on the Co cell in spin up channel. We notice that for more massive perturbations, like magnetic molecules for example, the difference may be more pronounced.

We discuss now the effect of spin-spin interactions. As for the Ag/Vanadocene case, we calculate the H 0 and J matrix from spin up and spin down DFT Hamiltonians. We can verify the localization of exchange parameters plotting diagonal elements as it is shown in Carbon atoms are collected on the right. We observe again large values of J on Cobalt d-orbitals and much smaller values on closest Carbon atoms. We can thus see again that spin-spin exchange interactions have a rather local character and are mainly concentrated on the atoms around and especially on the Co ad-atom.

We propagate now the wave packet as before and visualize it on a map having 300 × 300 unit cells of Graphene. The wave packet is initially in the spin down channel and at the same position as before at time t 0 . It propagates in the direction of the central cell with Co atom and starts to interact with it at the time t 1 , opening the second spin up channel due to spin-flip events. At time t 2 , the wave packet in the spin down channel tends to diffuse but conserves its preferable transport direction. For spin up channel, on the contrary, the electronic wave function seems to spread in a more isotropic way around the Co ad-atom, showing a more spherical ring. We have found the probability of the spin-flip to be about ≈ 0.05. We note that it is much smaller than for Ag/Vanadocene case due to presumably 2D character of the system and therefore weaker global impact of local spin-flip interactions on a 2D wave packet.

Conclusions and perspectives

During My PhD, I was developing a methodology and computational tools for simulating quantum transport based on the Quantum-ESPRESSO package and the Wannier90 and EPW codes. We implemented transport calculations in a realistic tight-binding approach with the NEGF and Wave packet methods. The code was tested on different 1D models such as molecular junctions made of Ag wires and Benzene or Vanadocene molecules, but also on 2D systems like a Phosphorene or a Graphene layer. In these systems, we tested the wave packet propagation and the NEGF transmission calculations for both electron and phonon transport, including electron-phonon or spin-exchange interactions.

The complementarity of the NEGF and wave packet methods is exemplified in different calculations of 1D systems. Whereas the NEGF method gives a global transport characteristics of a system, the wave packet method gives a detailed vision on the real dynamics of carriers. Section 4.1.2 illustrates this by identifying molecular states acting in the destructive quantum interference or the phonon modes mediating the thermal transport. In Section 4.1.3, a wave packet analysis allowed to understand the role of d z 2 orbital of the Vanadium in the spin-filter and spin-flip properties.

The calculations at finite voltage are directly linked to different thermoelectric phenomena such as Seebeeck, Joule or Peltier effects. Our code verifies these properties in the elastic regime and also provides consistent results for inelastic conductance. In the future, in addition to the approximation of a molecule coupled phononically to its thermal bath, we are going to implement a realization where the junction is coupled to leads not only electronically but also by vibrations and therefore electron and phonon Green's functions should be calculated self-consistently in the presence of electron-phonon interactions on the molecule. Moreover, other kinds of interactions, such as electron-photon for example, can be also worked out since formally the formalism is very similar to already implemented case of electron-phonon interactions.

The electron-phonon calculation with wave packet was introduced with an example of one or several phonon channels, which provides the calculation of the phonon excitation probability. However, we can explore deeper this model. We can analyse, for example. the decoherence of electrons or phonons after their interaction at the junction. We can study, moreover, a sequence of wave packets -wave packet train -and explore the energy transfer to the molecule vibrations. In the future, a more ambitions project would be to simulate a phonon and an electron propagation simultaneously to observe how they interact in the junction or, on the contrary, to model how a phonon created by an interaction with an electron will dissipate due to its propagation deep into the leads.

We have tested our code on 2D systems such as a single-layer of Black Phosphorus, a Phosphorene, in Section 4.2 and a Graphene in Section 4.3. We are going to put more efforts in this direction in order to optimize the code and to perform some more advanced simulations in 2D materials, including different types of interfaces or molecules.

We also collect and catalogue input and output data for different systems, needed to construct an efficient Wannier basis and to perform subsequent transport simulations. Finally, our transport code is expected to be integrated into the planewave transport code PWcond which is already included in Quantum-ESPRESSO package.

Titre:

  Simulation ab initio du transport et de la dynamique quantique de differents porteurs avec interactions Mots clés: Quantum-ESPRESSO, transport quantique, spin-flip, NEGF, electronphonon, paquet d'onde Résumé:

Figure 1 . 1 :

 11 Figure 1.1: Physical setup for our transport code: it represents a nanoscale junction connected to two leads. Various carriers, like electrons or phonons, propagate through the junction and can moreover interact inside it. A gradient of temperature or a voltage can be applied between the leads to drive the motion of carriers.

Figure 1 . 2 :

 12 Figure 1.2: Three main methods to compute the transmission and reflection functions of different carriers.

Figure 1 . 4 :

 14 Figure 1.4: Calculational scheme of our transport code.

Figure 2 . 1 :

 21 Figure 2.1: Schematic representation of a simulated system with a propagating carrier scattered on a perturbed central region.

Figure 2 . 2 :

 22 Figure 2.2: Algorithm of wave function propagation by the Schrödinger equation.The time, represented under the evolution axis, is discretized in N + 1 steps. Each time corresponds to a wave function ψ above the evolution axis. The time step t n+1 -t n is dt n represented in red. The differential of the wave function between time t n and time t n+1 is dψ(t n ) (represented in red).

Figure 2 . 3 :

 23 Figure 2.3: Vertex describing the scattering of an electron from the state j to the state j ′ and excitation/absorption of a phonon in the mode λ.

Figure 2 . 4 :

 24 Figure 2.4: Schematic illustration of the model for the electron wave packet propagation, including the electron-phonon interactions at the junction. j represents electronic sites of the system, j = [1; N ], while n is the number of phonons in the mode of interest.

. 73 )

 73 Similar to the finite system, one can introduce local phonon basis by means of operators:

Figure 3 . 2 :

 32 Figure 3.2: Source: [11]. Different key quantities entering Feynman diagrams with electron-electron or electron-phonon interactions.

Figure 3 . 3 :

 33 Figure 3.3: Keldysh contour is twisted and divided into two sub-contours C 1 and C 2 . The right-side edge of each contour corresponds to the two arguments of the object F defined by Eq. 3.33. Both contours are linked at -∞.

Figure 3 . 4 :

 34 Figure 3.4: Feynman diagram of the Hartree term.

Figure 3 . 5 :

 35 Figure 3.5: a) Feynman diagram of the Fock term. b) Scheme of the contour corresponding to each type of Fock's self-energy Σ. Each of them has a fixed time τ 1 and τ 2 on the contour.

Figure 3

 3

Figure 3 . 6 :

 36 Figure 3.6: a) Feynman diagram of the phonon self-energy term. b) Scheme of the contour corresponding to each type of self-energy Π. Each of them has a fixed time τ 1 and τ 2 on the contour.

Figure 4 . 1 :

 41 Figure 4.1: Band structure of an infinite Ag linear atomic chain calculated with QE (black) and using Wannier functions (red). All the bands are labelled by their main atomic orbital characters.

Figure 4 . 2 :

 42 Figure 4.2: a) Total DOS of an Ag infinite atomic chain compared to the PDOS on d and s orbitals. b) PDOS resolved over five d orbitals.

Figure 4 . 3 :Figure 4 . 4 :

 4344 Figure 4.3: Real-space images of different Wannier orbitals sampling the electronic structure of the Ag wire and used further in the transport code. Red and blue colours correspond to positive and negative sign, respectively. a) shows s-orbital centred at the Ag-Ag bond; b) represents different d-orbitals.

Figure 4 . 5 :

 45 Figure 4.5: Two studied configurations of Benzene molecule joining two Ag wires. They can be realized at different Ag-Ag wire separations.

Figure 4 .

 4 4 presents the phononic bands of the Ag wire along the z direction. The three bands correspond to three directions of oscillations of Ag atoms: longitudinal one (oscillations along the wire axis z) shows positive energies going up to 30 meV while two transverse bands (oscillations in the xy plane) are degenerate by symmetry and have imaginary frequencies (shown by negative values in Figure

  DOS (arb.unit) PDOS (arb.unit)

Figure 4 .

 4 Figure 4.6: a) and b) represent the total DOS of tilted and flat configurations, respectively, as a function of energy. c) and d) show the projected DOS onto Benzene states (molecular PDOS) for two configurations.

Figure 4 . 7 :

 47 Figure 4.7: Source of the image [5]: (a) Energy diagram of a tight-binding model of Benzene molecule with a hopping elements -t between the nearest Carbon atoms. On-site energy is denoted as ϵ 0 . (b) representation of six molecular orbitals and their energy levels. Green and blue colors represent the positive and negative signs of electronic wave functions, respectively.

Figure 4

 4 Figure 4.8: a) Propagation pathways for electron transport around the Fermi energy for Ag/Benzene junctions in two configurations. b) Transmission as a function of energy for two configurations and for the perfect Ag wire calculated by NEGF method.

Figure 4 . 9 :

 49 Figure 4.9: Propagation of the wave packet across Ag/Benzene junction in flat and tilted configurations. On each panel the lines on the left and right sides show the probability function on the left and right electrodes, respectively, while the columns in the middle describe the contributions from the Benzene molecular orbitals (which are 18 in total).

Figure 4 .

 4 Figure 4.10: a) Electron-phonon coupling elements with the LUMO state for different vibration modes of an isolated Benzene molecule. b) Representation of the important mode "30" mostly coupled to the LUMO. The red arrows represent the direction and the amplitude of atomic oscillations.

Figure 4 . 12 :

 412 Figure 4.12: Excitation rates for each vibration mode for flat and tilted configurations. The mode with the highest rate in each case is represented on inset.

Figure 4 Figure 4 . 13 :

 4413 Figure 4.13: Vibrational frequencies as a function of the mode number for the case of Ag/Benzene junction in two configurations. The frequencies for the full junction and for the Benzene only (fixing all Ag atoms) are both shown by black and red lines, respectively. For both configurations, two lowest frequency modes which are activated due to coupling to Ag chains are represented on insets.

Figure 4 .

 4 Figure 4.14: a) Excitation rates for each vibration mode in two configurations, considering the Benzene sub-matrix blocks in the Dynamical matrix. b) Important modes with the highest rates for each case.

Figure

  Figure 4.14 displays the excitation rates for those phonons calculated from the molecule sub-matrix. For the "tilted "case, the mode "30" (marked as "5" in figure) still has significant excitation rate, close in fact to the one found previously for the isolated molecule (see Figure4.12). However, two other modes, the mode

Figure 4 .

 4 Figure 4.15: a) Phonon transmission as a function of energy for the perfect Ag wire and for the two configurations of Ag/Benzene junction. b) Phonon's PDOS on the Benzene and on the apex Ag atom.

Figure 4 . 16 :

 416 Figure 4.16: Phononic wave packet at the moment of its passage across the molecule for two configurations at specific energies.

  (a) represents the charge current as well as the heating power dissipated in two electrodes as a function of voltage at zero temperature. In our symmetric geometry with equal coupling of molecular orbitals to the left and right electrodes the voltage was assumed to apply Different type of current (arb. units) QR Heating in right lead QL Heating in left lead IL.V Power dissipated QR +QL Total heating Voltage (V)

Figure 4 .

 4 Figure 4.17: a) Charge current and heating of two electrodes as a function of voltage for Ag/Benzene junction in the "flat" configuration. b) Joule effect: the total heating of electrodes matches perfectly to the total dissipated power given by IV .

Figure 4 Figure 4 . 18 :

 4418 Figure 4.18: Thermoelectric effects for the Ag/Benzene junction in the "flat" configuration. Illustrations of the Seebeck effect (a) and of the Peltier effect (b). c) Induced charge current as a function of temperature of the left electrode at zero voltage. d) Heating of the two electrodes as a function of applied voltage at constant (room) temperature.

Figure 4 .

 4 Figure 4.19: a) Transmission function in the "tilted" configuration with a gate of -1.4 V applied locally on the Benzene molecule. b) The corresponding heating as a function of voltage.

Figure 4 . 20 :

 420 Figure 4.20: Thermal flux as a function of temperature gradient in the "flat" configuration. a) Electron's heat current at zero voltage. b) Phonon's energy current.

Figure 4 . 21 :

 421 Figure 4.21: Differential conductance as a function of voltage in elastic and inelastic regimes for the "flat" (a) and "tilted" (b) configurations. c) Heating power to the molecule as a function of voltage for both cases.

Figure 4 . 22 :

 422 Figure 4.22: Model of a magnetic junction: Ag wires connected by a Vanadocene molecule.

Figure 4 .

 4 Figure 4.23: a) DOS as a function of energy for the Vanadocene connected to Ag wires. Total DOS calculated by different methots as well as the molecular PDOS are shown. b) PDOS on Vanadium d-, p-and s-orbitals. c) PDOS on different Vanadium d-orbitals.

Figure 4 .

 4 23(b) where PDOS on different Vanadium orbitals are plotted. Transmission Perfect Ag wire Spin up Spin down

Figure 4 . 24 :

 424 Figure 4.24: Transmission as a function of energy for spin up and spin down channels through a Vanadocene molecular junction connecting Ag wires.

Figure 4 . 25 :

 425 Figure 4.25: Evolution of a wave packet through the Ag/Vanadocene junction for spin up and down channels. The columns in the middle of each panel display now the contributions from Wannier functions of the Vanadocene molecule (35 in total).

Figure 4 .

 4 Figure 4.26: On-site exchange interaction elements for all Wannier orbitals in the Vanadocene connected ro Ag wires. Ag wire's orbitals far from Vanadocene are all represented on the left, while the molecule's orbitals -on the right. spin states Reflection in junction transmission total |↑, 1/2⟩ 0.1345 0.0000 0.1155 0.2500 |↓, 3/2⟩ 0.3688 0.0000 0.3812 0.7500 Total 0.5033 0.0000 0.4967 1.0000

Figure 4 . 27 : 2 m 4 :

 42724 Figure 4.27: Time propagation of |↓, 3/2⟩ wave packet through the Vanadocene molecule when spin-flip interactions are switched on. electron spin m S = 3 2 m S = 1 2 m S = -1 2 m S = -3 2

Figure 4 .

 4 28(a) and it was calculated with (20 × 20 × 1) k-points in the SCF calculation. It shows a gap

Figure 4 .

 4 Figure 4.28: a) Band structure of Phosphorene (a monolayer of BP) along the k-path ΓXW Y Γ. b) Images of the bonding and antibonding states at the Γ point. c) bonding and antibonding distribution probabilities along the z direction while integrated in the xy plane.

Figure 4 .

 4 28(c) presents the distribution probabilities of these states along the z direction perpendicular to the Phosphorene plane. One can see that the two states have rather different behavior: while the B state (corresponding to the valence band) shows a strong localization in between the two sites (making a bond), the AB state (corresponding to the conduction band) presents a node there and spreads much more on both sides in the vacuum region. These features can be readily reproduced by a simple model with two Gaussians simulating the two atomic wave functions with nonzero overlap S.In collaboration with the MPQ group (Paris Diderot University), we have worked on simulation of STM (Scanning tunnelling microscope) images, in par-

Figure 4 .

 4 Figure 4.29: a) Total DOS and PDOS on different p-orbitals of 10 layer BP slab integrated over full 2D Brillouin Zone. b) LDOS in the vacuum plates at different heights D of 10 layer slab integrated over full 2D Brillouin Zone. c) STM image with atomic resolution of the area where the spectrum was measured. d) dI/dV spectrum measured over a larger sample bias range.

Figure 4 . 30 :

 430 Figure 4.30: Model of a junction consisting of a periodic line of Co atoms adsorbed on the Phosphorene.

Figure 4 .Figure 4 .

 44 Figure 4.31: a) DOS of the Co line ad-atoms on Phosphorene, total DOS as well as different Co contributions are shown. b) DOS decomposed over s and d Co orbitals.

Figure 4 .

 4 Figure 4.31(b) shows the PDOS on different d-orbitals of Co. Among them, we observe that only the d z 2 -orbital (for spin down) is located above the Fermi level. For spin up states, all the d-orbitals are placed at the edge of the Phosphorene valence band. This will result in spin unbalance corresponding to the total magnetic moment of about 1 µ B which is found in calculations. Moreover, the d-orbitals for spin down are located in the gap region of Phosphorene, [-0.6; 0.3] eV. In the spin up case, on the contrary, no Co orbital contributes to the DOS at the Fermi energy. This indicates at the full spin-polarization of Co states at the Fermi level, which could be important for spin-filter properties. We have therefore calculated the spin-polarized transmission in the y direction perpendicular to the Co chain, which is shown in Figure4.32. The total transmission is obtained by integrating the k-resolved transmission over 200 k z -points. Unfortunately, we do

Figure 4 . 32 .Figure 4 . 32 :

 432432 Figure 4.31(b) shows the PDOS on different d-orbitals of Co. Among them, we observe that only the d z 2 -orbital (for spin down) is located above the Fermi level. For spin up states, all the d-orbitals are placed at the edge of the Phosphorene valence band. This will result in spin unbalance corresponding to the total magnetic moment of about 1 µ B which is found in calculations. Moreover, the d-orbitals for spin down are located in the gap region of Phosphorene, [-0.6; 0.3] eV. In the spin up case, on the contrary, no Co orbital contributes to the DOS at the Fermi energy. This indicates at the full spin-polarization of Co states at the Fermi level, which could be important for spin-filter properties. We have therefore calculated the spin-polarized transmission in the y direction perpendicular to the Co chain, which is shown in Figure4.32. The total transmission is obtained by integrating the k-resolved transmission over 200 k z -points. Unfortunately, we do

Figure 4 . 33 :

 433 Figure 4.33: Band structure along the k z direction of an isolated Co wire (left) and of Phosphorene/Co junction (right). For isolated Co wire the same lattice parameter (in the z direction) as the one in the Phosphorene/Co junction was used.

Figure 4 .

 4 34(a) presents the band structure of the Graphene along the usual k-point path. The black curve shows the QE bands after performing an SCF calculation with 12 × 12 × 1 k-points. The Graphene is known for its famous Dirac cones, at K = (2/3, 2/3) (in crystal units), around which the bands show quasilinear dispersion as seen in Figure4.34(a). The DOS of the Graphene around the Fermi level is shown in Figure4.34(b). We observe the zero band gap DOS (arb. unit)

Figure 4 . 34 :

 434 Figure 4.34: Graphene: a) band structure along k-path ΓM KΓ calculated by QE as well as with Wannier functions; b) total DOS and PDOS resolved on different s and p orbitals.

Figure 4 . 35 :

 435 Figure 4.35: Graphene states for transport: a) 2D band structure calculated with our code with a grid of k-points centred at the K-point; b) velocity vectors of Graphene on the Dirac cone. The purple arrows represent the velocity in the zy plane, the length of the arrows is proportional to the velocity magnitude.

Figure 4 . 36 :

 436 Figure 4.36: Model of a Co ad-atom deposited on Graphene layer with a supercell periodicity of 6 × 7 unit cells of Graphene in the xy plane.

Figure 4 . 37 :

 437 Figure 4.37: Co ad-atom/Graphene system: spin-resolved total DOS and the PDOS on s and d Co orbitals. Also, the total DOS calculated with Wannier90 is shown (in cyan).

Figure 4 . 38 :

 438 Figure 4.38: Co ad-atom/Graphene system: total DOS and PDOS on different Co orbitals in a large energy window used to chose Wannier90 parameters.

Figure 4 . 39 :

 439 Figure 4.39: Simulation of a 2D wave packet propagation in a Co adatom/Graphene system for spin up and spin down electrons (without spin-flip).

Figure 4 . 40 .

 440 Graphene orbitals which are far from the Co atom are all represented on the left part of the figure while the perturbed Exchange interaction strength (eV)

Figure 4 .

 4 Figure 4.40: Spin-exchange interaction elements on Wannier orbitals for a Co adatom/Graphene.

Figure 4 . 41 :

 441 Figure 4.41: Simulation of a 2D wave packet propagation with spin exchange interactions on a Co ad-atom/Graphene. The wave packet is initialized at the time t 0 in the spin down channel.

Table 4

 4

	.1 demonstrates the

Table 4 .

 4 2: Different contributions in final wave packet for spin up and spin down electrons.

	3 2 ↓, 3 2	0.9964 0.1244	0.0000 0.0000	0.0036 0.8756	1.0000 1.0000

Table 4 .

 4 3 gives a probability of having spin-flip, ≈ 0.22, which is quite high (much

	exchange interaction strength (eV)
	Wannier orbitals

Table 4 . 3

 43

		0.1345	0.0000	0.1155	0.2500
	|↓, 3/2⟩	0.3688	0.0000	0.3812	0.7500
	Total	0.5033	0.0000	0.4967	1.0000

: Different contributions in the final wave packet, taking into account spin-flip events between |↓, 3/2⟩ and |↑, 1/2⟩ states. larger than the probability to excite a phonon in the previous section). Moreover, comparing this result with Table

4

.2, we see that the spin interaction decreases the reflection probability by a δT ≈ 0.512 what is more than 50%.

Acknowledgements / Remerciements

Appendix

In order to proceed in analogy with the finite system's analysis we need realvalue eigenfunctions φ, not complex, which could be always constructed since the dynamical matrix is always real and symmetric by construction, even for periodic systems. To do that we pass to notation {q, -q, ν} → {Q, p, ν} where Q number all the {q, -q} pairs and p = 1, 2 indicates the sign, +/-. In these new notations:

iαR so these pairs with p = 1, 2 correspond to the same eigenvalue of the dynamical matrix and can therefore be combined to produce real-value eigenfunctions:

with the unitary matrix:

This is similar to the transformation from two conjugated exponentials (complex) to cos and sin functions (real). Inverse transformation relates the real function basis to the standard one:

Using this real function basis in the unity resolution and in the spectral representation of the dynamical matrix yields:

where again

123 with inverse transformation:

These operators refer again to the conjugate variables:

but are now Hermitian:

Hamiltonian can be now written as follows:

in terms of standard creation and annihilation operators:

From Eq. ( 10) one finds for displacement operators:

Rotation of creation and annihilation operators to the standard basis:

and inverse transformations:

The Hamiltonian in the standard basis keeps the same form:

The expression for creation operators in terms of displacement operators:

The expression of displacement operators:

We can go now to original notations, {Q, p, ν} → {q, -q, ν}.