Groningin Molecular

Khon-Sham Ks

Srff

Quantum Interference EDL Electronic Double Layer vdW van der Waals TMDCs Transition Metal Dichalcogenides HICD High Interfacial Charge Density FETs Field-Effet Transistors FFs Force Fields ILs Ionic Liquids CPU Central Prossessing Unit CPMD

List of Figures

interface produced by classical MD with a surface charge of 73 µC/cm 2 . The electron density is shown at an isosurface value of 0.015 e/ a.u. 3 . Color code is cyan for Mo, yellow for S, organge fo BF 4 ions and blue for EMIM ions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.9 EMIM-BF 4 ions density and charge density profiles with respect to the distance from each electrodes with surface defects of 10 and 20 S vacancies. des liaisons hydrogènes entre l'atome H2 du cation du IL et les atomes F de l'anion BF - 4 en fonction de la longueur de la chaîne alkyle du cation (DMIM + EMIM + , BMIM + et HMIM + ) obtenue par MD classique à température ambiante. Line inférieure gauche: fonction de distribution combinée montrant la liaison hydrogène obtenue par FPMD. À droite : étiquetage des principaux atomes cités dans le texte principal et des vecteurs définissant l'angle pour lequel les diagrammes de densité de distribution sont affichés . 6.8 À gauche : étapes initiales par MD classique du processus de formation de l'interface constituée de 500 paires de EMIM-BF 4 en contact avec une bicouche de MoS 2 . A droite : modèle structural final typique de l'interface EMIM-BF 4 @MoS 2 équilibrée à 300 K. . . . . . . . . . . . . . . . . . . . . . 6.9 Distributions de la densité électronique des interfaces EMIM-BF 4 @MoS 2 réalisées avec la surface sans défauts de MoS Note that the charge values of some equivalent atoms have been averaged in agreement with schematic numbering of Figures Je souhaiterais également adresser un remerciement particulier à monsieur le Dr. Guido Ori qui a été comme un co-directeur de ce travail de thèse. Son implication dans l'aboutissement de ce travail a été énorme. Merci d'avoir été très bon conseiller et de m'avoir encouragé dans les moments difficiles. Merci de ton amitié et des discussions fantastiques que nous avons eu tous les deux ces dernières années. Tu as toujours été disponible et à l'écoute pour moi, infiniment merci pour tout.

Je souhaiterais exprimer ma profonde reconnaissance à monsieur le Dr. Carlo Massobrio pour les précieux conseils qu'il a pu donner au cours de ces quatre dernières années et d'avoir facilité mon intégration dans l'équipe. Tu as été un très bon leader et un véritable modèle pour moi. Ton esprit critique m'a permis de m'améliorer au fil des années.

Je souhaiterais exprimer ma gratitude à madame la Dr. Evelyne Martin pour les nombreuses initiatives mises en place pour faciliter les échanges instructifs et enrichissants entre les membres de l'équipe théorique. Ces mécanismes ont non seulement renforcé la cohésion de notre équipe, mais ont également permis de développer de nouvelles idées et perspectives.

Je tiens à remercier chaleureusement le programme IDEX PhD International de l'Université de Strasbourg pour son soutien financier durant mes travaux de recherche. Je suis très reconnaissant pour ce financement précieux qui a rendu possible l'avancement de mes découvertes scientifiques. Un remerciement à l'équipe du meso centre de Strasbourg pour la disponibilité et l'assistance qu'ils m'ont témoigné tout au long de cette thèse de doctorat en me permettant de faire tourner mes calculs sur la plateforme et en résolvant tous les problèmes techniques auxquels j'ai fait face.

Je dédicace cette thèse à mes parents, monsieur Essomba André Pascal et madame Mbala Marie Jeanne. Vous êtes sans doute les parents que tout enfant souhaiterait avoir. Vous m'avez soutenu dans les bons et les moments difficiles sans vous, je ne serai pas ce que je suis aujourd'hui. Vous avez toujours trouvé des mots pour m'encourager à aller de l'avant et ne jamais abandonner. Je souhaiterais particulièrement vous remercier du plus profond de mon coeur pour l'éducation que vous m'avez inculquée. Quelques mots ne seraient suffisants pour quantifier tous les bien faits que vous avez apportés dans ma vie. Je ne peux que résumer cela en une seule phrase "infiniment merci pour tout".. J'ai décidé de vous dédier ce travail de thèse, une façon pour moi de vous dire merci et que je vous aime beaucoup.

Je souhaiterais remercier ma fiancée Mlle. Bakongo Micheline qui a toujours été à mes côtés dans les bons et les mauvais moments. Sa présence a eu un impact positif sur ma motivation à aller jusqu'au bout de ce travail. Elle a été ma force ces quatre dernières années. Infiniment merci de ton soutien et de tes précieux conseils, tu as toujours été et tu resteras ma force, une source de motivation et d'inspiration pour moi.

Je souhaiterais également remercier du plus profond de mon coeur toute ma famille pour les différentes aides, les soutiens, les encouragements et précieux conseils qu'ils m'ont prodigués. Vous avez toujours été là pour moi vraiment merci pour tout.

Je souhaiterais exprimer ma profonde reconnaissance et mes remerciements à monsieurs les Dr. Ottou Martin Thierry, Pr. Ndjaka Jean Marie et Pr. Zekeng Serge de l'université de Yaoundé I pour leurs nombreux conseils et leurs implications dans l'obtention de ma bourse. Je ne cesserai de vous dire merci d'avoir cru en moi, vos conseils seront toujours les bienvenus.

Enfin je souhaiterais dire un merci à monsieur Bonne Marc pour avoir été un très bon ami. J'ai passé de bons moments à l'université de Yaoundé I à tes côté, tu as toujours été un très bon conseiller, sincèrement merci.

Chapter 1 Introduction

Scientific context

Two-dimensional (2D) materials emerged as promising candidates for replacing silicon in next-generation electronics and they represent an ever-growing research area that requires a multidisciplinary approach ranging from fundamental physics to advanced applied chemistry and engineering and, as such, spanning from fundamental science to novel device applications [1,2]. This class of systems represents a large family of materials that exhibit unique layered architectures, differing in chemical composition and crystal structure, leading to markedly different properties. For instance, their electronic properties can span from metallic to insulating. Despite their diversity, all these materials share some common features: large spacings between the composing layers interacting via long-range van der Waals (vdW) forces, resulting in peculiar mechanical properties, like interlayer sliding and cleavability. In the forefront field of device electronics based on 2D (nano)materials such as graphene, boron nitride, black phosphorous, transition metal dichalcogenides (TMDC, of type MX 2 , with M a transition-metal atom (Mo, W, etc.) and X a chalcogen atom (S, Se, or Te)) the presence of very high interfacial charge density (HICD) has been shown to be crucial for the arising of exotic phenomena such as superconductivity [1] and quantum interference effects [1,2] as well as to empower emerging memory devices for neuromorphic computing [3,4] and electric double layers field-effect transistors (EDL-FET) [START_REF] Gao | Triboiontronic transistor of mos 2[END_REF]. In the field of post-silicon electronics, great attention was initially devoted to metallic graphene, rather than to 2D vdW semiconductors, for its ultrahigh carrier mobility [START_REF] Dubonos | Electric field effect in atomically thin carbon[END_REF]. Nonetheless, despite this historical predominance of graphene and worldwide efforts, it was realized that it would be extremely difficult to use graphene for any digital application due to the absence of a bandgap. Renewed interest on 2D vdW semiconductors arose in 2011 when Kis et al. [START_REF] Radisavljevic | Single-layer mos 2 transistors[END_REF] reported high carrier mobility in monolayer MoS 2 field-effect transistor. This was subsequently confirmed in 2012 by Strano et al. [START_REF] Qh | Electronics and optoelectronics of two-dimensional transition metal dichalcogenides[END_REF] In this context, ionic liquids (ILs) started becoming a key target to pioneer novel 2D materials-based transistors possessing HICD requirement [START_REF] Jo | Electrostatically induced superconductivity at the surface of ws 2[END_REF]. The interaction of ILs with the TMDC surface can trigger the modulation of electronic states at the liquid/solid interface, thus leading to promising phenomena such as electric-field-induced superconductivity [1,2,[START_REF] Radisavljevic | Single-layer mos 2 transistors[END_REF][START_REF] Costanzo | Gate-induced superconductivity in atomically thin mos 2 crystals[END_REF] triggered by the accumulation of unprecedented high carrier density (Figure 1.1). Yet, several key issues remain to be addressed about the specific interactions occurring between the ILs and the 2D material. Such a lack of fundamental understanding is even more evident if we consider the case of 2D materials with native and/or induced surface defects. The detailed interface structure and the chemical interactions acting in this crucial region of the devices remain to be clarified, thus representing a stumbling block for a detailed understanding of the basic functioning of these interfaces. This, in turn, possesses severe limitations to a rational design of actual devices. The current experimental and computational literature indicate clearly the following major challenges to face, namely two major issues have to be addressed: (i) General lack of fundamental understanding of the interactions between ILs with 2D materials and (ii) the number of ILs used for gating 2D materials proposed in the current literature is rather narrow and their design, so far, based on either pure phenomenology or commercial availability of ILs. Overall, one of the major difficulties is represented by the fact that this type of interface is buried in the inner part of the whole device, thus of particularly difficult access by experimental probes. In this context, computational approaches can offer a unique possibility to complement experiments and to provide detailed information needed to understand and disentangle the different roles played by ILs and the induced charge density at the surface of 2D materials. The theoretical modelling of the ILs at interface with 2D materials (termed ILs@2Dmat interface hereafter) is, to date, only partially covered by classical simulations (mostly classical molecular dynamics, MD) offering qualitative insights. Accurate simulations of the ILs@2Dmat interface's structure able to provide information also on chemical bonding and electronic properties are still lacking or at very best at a pioneering stage, thus still far from being an assessed way of providing the quantitative comprehension and predictive power sought. In the following three subsections, a more detailed account of main physical and chemical characteristics of ILs and TMDC are presented, along with a comprehensive state-of-theart of their atomistic modeling.

Ionic liquids

Ionic liquids (ILs) are salts formed by large ionic molecules (mainly organic) carrying a delocalized electrostatic charge, and characterized by conformational flexibility and asymmetric chemical structure. These intrinsic molecular properties are responsible for a low melting point close to room temperature and hinder crystallisation. They are represented by the imidazolium, pyridinium, quaternary ammonium and phosphonium families (among others) with different substituent groups, which can be combined with a large variety of anions (often inorganic), ranging from small and symmetrical to large and more flexible molecules (Figure 1.2). ILs show a relatively wide spectrum of intra-and intermolecular interactions including electrostatic, van der Waals, π-stacking, H-bonding, etc. Due to their ionic nature, they are often compared to high-temperature molten salts, and with liquid salts they share charge alternation properties and intermediate range ordering. Unlike water and other standard molecular liquids, ILs are single-component systems [START_REF] Plechkova | Applications of ionic liquids in the chemical industry[END_REF] and [START_REF] Zhang | Nanoconfined ionic liquids[END_REF]. in which the cations and anions likely play independent roles in determining the liquid behavior. However, ILs structural patterns may contain different features, such as microphase separation between polar and non-polar domains, associated to charged head groups and long alkyl side-chains. Charge delocalisation and large ion size lead to weaker interactions between ions, attenuating the long range ordering, and possible cation-anion hydrogen bonding results in structural directionality, especially pronounced for ILs in their crystalline phases.

Alkyl-imidazolium ionic liquids (AI-IL)

Structure of AI-IL

Heterocyclic compounds containing nitrogen atoms such as imidazole, pyridine, pyrrolidine, and their derivatives are the most common cations composing different types of ILs when linked to various anions such as basic halides or more complex organic/inorganic moieties. One of the most studied classes of these prototype ILs is the one based on alkyl-imidazolium ILs (AI-ILs) made of a five-member aromatic heterocyclic ring (imidazolium, hydrophilic head) bound to aliphatic chains (hydrophobic) of different lengths. Often such ILs are liquid salt with a melting point below 100 • C and sometimes even below 25 • C (Figure 1.3). The hydrophilic head has two nitrogen atoms in positions 1 and 3, according to the numbering of Figure 1.2, symmetrical with respect to each other, thus allowing for an easy fictionalization of AI-ILs and consequent tuning of their properties by a simple change of the side group. Often, one side group is a methyl group (-CH 3 ) whereas the other side group can be an aliphatic chain of different lengths or carrying other functional groups. In this particular structure, the C atom shared between the two N atoms (C 2 ), is considered to be a donor and its bonded hydrogen (H) atom, being the most acidic of all, is often the only one capable of participating to proton transfer processes or to the formation of hydrogen bonds (H-bonds) with neighbouring acceptors. The hydrophobic part is an alkyl chain of different length: methyl (CH published literature, the following acronyms are used: DMIM (1,3-dimethylimidazolium), EMIM (1-ethyl-3-methylimidazolium), BMIM (1-butyl-3-methylimidazolium), and HMIM (1-hexyl-3-methylimidazolium). Other compounds have also been explored. Here we limit our choice to the ones studied in this research project, being the major target of present studies worldwide. Depending on the length of the alkyl chain of the cation, the physicochemical properties of the ionic liquid vary from one IL to another. On the other hand, the IL anion can be of different nature, whether inorganic or organic, and show different ability to promote different interactions with the IL cation. Imidazolium ILs form an extended H-bond network in the liquid state and are consequently highly structured. For these reasons they are also referred as 'supramolecular fluids'. This structural organization of ILs can be used as the 'entropic driver' for the spontaneous, well-defined, and extended ordering of nanoscale structures [START_REF] Dupont | Physico-chemical processes in imidazolium ionic liquids[END_REF]. In terms of structural properties, X-ray and neutron scattering experiments can be employed to extract a detailed information of the structure of disordered materials from short-to-intermediate range orders. However, to date this type of characterization has been applied only to a few ILs systems. As an example, neutron scattering of 1,3-dimethylimidazolium chlorine [START_REF] Deetlefs | Liquid structure of the ionic liquid 1,3-dimethylimidazolium bis(trifluoromethyl)sulfonylamide[END_REF][START_REF] Hardacre | Structure of molten 1,3-dimethylimidazolium chloride using neutron diffraction[END_REF][START_REF] Hardacre | Liquid structure of 1, 3-dimethylimidazolium salts[END_REF] has been performed by selective deuteration of cations in order to vary the weight of certain contributions, and where the symmetry of the cation and the absence of flexible side chains facilitate the interpretation of the diffraction patterns. The results obtained show a more intense and structured first peak characteristic of the structural order at very short interatomic distance, along with fading oscillations at large interatomic distances averaged around the unit, due to the amorphous nature of ILs.

Properties of AI-IL

Alkyl-imidazolium ILs show physico-chemical peculiar properties remarkably different from ordinary molten salts. Many studies have shown that the melting temperature or the solubility properties of these ILs generally depend on both the cation and anion used (reported in Table 1.1). More importantly, their properties can be tailored to specific applications by a proper combination of the two ionic compounds [START_REF] Bates | Co2 capture by a task-specific ionic liquid[END_REF][START_REF] Giernoth | Task-specific ionic liquids[END_REF]. ILs have melting temperatures below 100 • C [START_REF] Wilkes | Air and water stable 1-ethyl-3-methylimidazolium based ionic liquids[END_REF][START_REF] Fuller | The room-temperature ionic liquid 1ethyl-3-methylimidazolium tetrafluoroborate: electrochemical couples and physical properties[END_REF] which means that often they are liquid at room temperature [START_REF] Jg | Characterization and comparison of hydrophilic and hydrophobic room temperature ionic liquids incorporating the imidazolium cation[END_REF] and therefore good candidates for applications in which they have to be uniformly distributed on a selected area, in particular on nanosize surfaces. For instance, the melting points of the EMIM + cation in combination with different anions, such as BF - 4 and TFSI -(also denoted as Tf 2 N -in literature) are 15 C [START_REF] Jg | Characterization and comparison of hydrophilic and hydrophobic room temperature ionic liquids incorporating the imidazolium cation[END_REF] and -4 • C [START_REF] Kalyanasundaram | Hydrophobic, highly conductive ambient-temperature molten salts[END_REF], respectively. This type of ILs have extremely low vapor pressure and good thermal and chemical stabilities. Often, the thermal stability is limited by the strength of their heteroatom-carbon and their heteroatom H-bonds, respectively [START_REF] Wasserscheid | Ionic liquids-new "solutions" for transition metal catalysis[END_REF]. Wilkes et al. [START_REF] Me | Thermochemistry of ionic liquid heat-transfer fluids[END_REF] reported that the ionic liquids EMIM-BF 4 , BMIM-BF 4 , 1,2-dimethyl-3-propyl imidazolium-TFSI are stable up to temperatures of 445, 423 and 457 • C, respectively. The electrochemical window of ILs is another important property to be kept into account and plays a major role in the use of ionic liquids in a wealth of applications such as FET, energy storage devices and whenever ILs are used in electrodeposition processes. This window determines the electrochemical stability of ILs. Furthermore, it allows to overcome the limited one of water (1.2 V [START_REF] Endres | Air and water stable ionic liquids in physical chemistry[END_REF]) and its typical values are: 4.15 V for BMIM-PF 6 at a platinum electrode [START_REF] Marken | Water-induced accelerated ion diffusion: voltammetric studies in 1-methyl-3-[2,6-(s)-dimethylocten-2-yl]imidazolium tetrafluoroborate,1-butyl-3-methylimidazolium tetrafluoroborate and hexafluorophosphate ionic liquids[END_REF], 4.10 V for BMIM-BF 4 [START_REF] Marken | Water-induced accelerated ion diffusion: voltammetric studies in 1-methyl-3-[2,6-(s)-dimethylocten-2-yl]imidazolium tetrafluoroborate,1-butyl-3-methylimidazolium tetrafluoroborate and hexafluorophosphate ionic liquids[END_REF] and 5.5 V for BMIM-Tf 2 N at glassy carbon electrodes [START_REF] Dr | Pyrrolidinium imides: A new family of molten salts and conductive plastic crystal phases[END_REF]. Moreover, ILs are typically more viscous than molecular organic solvents [START_REF] Crosthwaite | Phase transition and decomposition temperatures, heat capacities and viscosities of pyridinium ionic liquids[END_REF], with a viscosity that typically increases with the hydrophilic part of the cation [START_REF] Endres | Air and water stable ionic liquids in physical chemistry[END_REF]. The unique physico-chemical and viscosity properties of ILs and have shown great potential in many applications with lubrication as one of the latest. While earlier work were dedicated almost exclusively to the exploration of the feasibility of the use of ILs as neat or base lubricants, their use as lubricant additives has become the new focal research topic since the breakthrough in ILs' miscibility in non-polar hydrocarbon oils in early 2012 [START_REF] Zhou | Ionic liquids as lubricant additives: A review[END_REF]. ILs have a higher density than water with values varying from 1.02 to 1.6 g.cm -3 . Typically, ILs densities increase with the length of the alkyl chain of the cation [START_REF] Endres | Air and water stable ionic liquids in physical chemistry[END_REF], and for some of them they show some correlation with temperature [START_REF] Kim | Physical and electrochemical properties of 1-butyl-3-methylimidazolium bromide, 1-butyl-3-methylimidazolium iodide, and 1-butyl-3methylimidazolium tetrafluoroborate[END_REF][32][START_REF] Muhammad | Thermophysical properties of 1-hexyl-3-methyl imidazolium based ionicliquids with tetrafluoroborate, hexafluorophosphate and bis(trifluoromethylsulfonyl)imide anions[END_REF] specifically, a decrease with temperature. They can show very high ionic conductivity (up to ∼10 mS.cm -1 ) [START_REF] Endres | Air and water stable ionic liquids in physical chemistry[END_REF] and are increasingly used as an electrolyte and gating agents as a trigger in certain electronic devices such as electronic double layer field effect transistors [2].

Transition metal dichalcogenides (TMDCs)

Structural properties

Two-dimensional (2D) materials have attracted attention not only because of the rich and fundamental physics due to their intrinsic structural and electronic properties but also because of their potential for nanoscale device applications. Graphene has been and many extent still is the prominent member of the 2D materials family, but its gapless band structure, useful in other contexts, has a non-negligible drawback for realizing switching operation, which is essential for many logic devices. Even though the bandgap of graphene can be engineered by depositing sheets on particular substrates or by fabricating nanoribbons, it deteriorates the carriers mobility. For these reasons, researchers have turned their attention to alternative 2D materials. One of the most promising and so far better studied 2D systems are the transition metal dichalcogenides (TMDCs). In particular, TMDCs having the composition MX 2 , with M being a transition metal atom (Mo, W, Ti, Nb, ...) and X a chalcogen atoms (S, Se, Te), have been the target of major investigations (Figure 1.4). These materials can be easily exfoliated into 2D layers from their stacked crystal structure by using the same method used for graphene production. Most intriguingly, their band structures are layer-thickness-dependent despite the weak interlayer van der Waals interactions, which indicates they are electronically tunable via thickness control. As the number of layer reduces from a bulk system to a monolayer, the bandgap of several TMDC materials changes from indirect to direct. Among those and the most widely studied materials, monolayers of MoS 2 have emerged as a semiconducting alternative to graphene because of a large intrinsic direct bandgap of approximately ∼1.8 eV, which makes it suitable for optoelectronic and nanoelectronic applications. TMDCs consist of Table 1.1. Typical ranges of room-temperature ILs properties. Adapted from [START_REF] Werner | Ionic liquids in chemical engineering[END_REF][START_REF] Fedorov | Ionic liquids at electrified interfaces[END_REF] .

Property

Typical values for most of RTILs Lower Upper limit example melting point/glass transition 0-60 two-dimensional sheets stacked along the crystallographic c axis where each sheet is composed of an M plane inserted between two X planes forming a X-M-X monolayer structure. The M-X bonds inside the same (mono)layer (intralayer bond) are strong chemical bonds of covalent nature, whereas the X-X interactions, responsible for the packing of the layers, are long-range vdW interactions. TMDCs have a hexagonal honeycomb-shaped atomic lattice and they are formed by three atomic X-M-X type planes where each M atom is bonded to six chalcogen atoms. These chalcogen atoms are arranged in prismatic trigonal coordination in the 1H polymorphic phase (Figure 1.5B) and in prismatic octahedral or antiprismatic trigonal coordination in the 1T polymorphic phase (Figure 1.5C). According to the experimental indications, for this project we consider the 2H polymorphic phase of MoS 2 (2H-MoS 2 ). The number 2 in the notation 2H refers to the number of sheets present in a hexagonal honeycomb-shaped lattice stacked as SMoS-SMoS (Figure 1.5D). As mentioned above, the intralayer bonds in MoS 2 are considered to be of covalent type, while the interlayer interactions keeping the stacking are of the vdW type due to the relatively large interlayer distance d MoMo of approximately 6.15 Å. This accounts for the fact that the MoS 2 can be cleaved or exfoliated along the staching direction to obtain individual crystalline layers from the bulk (method pioneered in 1986 [START_REF] Joensen | Single-layer mos 2[END_REF]). This monolayer has a thickness of 3.01 Å (Figure 1.5A). The crystal unit cell has the symmetry of the P63/mmc space group, with lattice parameters equal to a = 3.1612 Å, b = 3.1612 Å and c = 12.2985 Å, α = 90, β = 90 and γ = 120 [START_REF] José-Yacamán | Structure and catalytic properties of hexagonal molybdenum disulfide nanoplates[END_REF][START_REF] Schönfeld | Anisotropic mean-square displacements (msd) in single-crystals of 2h-and 3r-mos 2 s[END_REF].

Electronic and transport properties

2D vdW layered materials like MoX 2 and WX 2 (X = S, Se, Te) can be semiconducting materials and have a bandgap whose width varies from bulk to monolayer and this makes the gap tunable via mechanical manipulations resuming to exfoliation as quoted in the previous paragraph. Figure 1.6A shows the band structure of MoS 2 in the bulk, quadrilayer, bilayer and in monolayer cases with the solid arrows indicating the lowest energy transitions. Bulk 2H-MoS 2 is characterized by an indirect band gap (E g = 1.23 eV) [START_REF] Kk | Detailed photocurrent spectroscopy of the semiconducting group vib transition metal dichalcogenides[END_REF] with the maximum of the valence band located at the Γ-point. Direct excitonic transitions occur at high energies at the K-point (E g = 1.80 eV). Upon reduction of the, the indirect bandgap increases, while the forward excitonic transition hardly changes. For a monolayer of MoS 2 , it becomes a direct bandgap semiconductor because the absolute minimum of the conduction band and the absolute maximum of the valence band become aligned and located in the valley marked by the K-point at the edge of the hexagonal Brillouin zone(BZ) with an energy gap going from E g = 1.23 eV to E g = 1.80 eV [START_REF] Kf | Atomically thin mos2: A new direct-gap semiconductor[END_REF]. This direct transition has also been observed in other TMDC semiconductors (see Table 1.2). This peculiar change in the electronic structure of MoS 2 is mainly due to the suppression of the interactions between the S atoms between the sheets, underscoring the role of the vdW interactions and making possible to sufficiently reduce the maximum of the valence band at the Γ-point, with an important contribution of the 3p z orbitals of S, so that the energy of the indirect transition becomes higher than that of the direct transition [START_REF] Kuc | Influence of quantum confinement on the electronic structure of the transition metal sulfide ts 2[END_REF][START_REF] Cappelluti | Influence of quantum confinement on the electronic structure of the transition metal sulfide ts 2[END_REF]. This transition has been observed experimentally by Mark et al. [START_REF] Kf | Atomically thin mos2: A new direct-gap semiconductor[END_REF] and can be rationalized in terms of a jump observed in the photoluminescence [START_REF] Eda | Photoluminescence from chemically exfoliated mos 2[END_REF][START_REF] Splendiani | Emerging photoluminescence in monolayer mos 2[END_REF] efficiency of MoX 2 and WX 2 (X = S, Se, Te) monolayers, which presents a yield of a factor of 10 4 times greater than that observed for MoS 2 bilayers of micrometric thicknesses (Figure 1.6B). In terms of transport Figure 1.6. A) MoS 2 band structure computed at DFT level [START_REF] Chhowalla | The chemistry of two-dimensional layered transition metal dichalcogenide nanosheets[END_REF][START_REF] Li | Electronic properties of mos2 nanoparticles[END_REF]. Energy dispersion (energy versus wavevector k) in bulk, quadrilayer (4L), bilayer (2L) and monolayer (1L) from left to right. The horizontal dashed line represents the energy of a band maximum at the k point. The red and blue lines represent the conduction and valence band edges, respectively. The lowest energy transition are indicated by the solid arrows. B) Photoluminescent of suspended monolayer and bilayer MoS 2 [START_REF] Kf | Atomically thin mos2: A new direct-gap semiconductor[END_REF].

properties, many studies have been conducted on TMDCs. For instance, by measuring the temperature dependence of the electrical resistivity and the Hall coefficient of semiconductor compounds and by deducing the mobility of the corresponding free carriers. [START_REF] Fivaz | Mobility of charge carriers in semiconducting layer structures[END_REF]. At low temperature it has been shown that TMDCs can have an electron mobility up to 500 cm 2 V -1 s -1 [START_REF] Fivaz | Mobility of charge carriers in semiconducting layer structures[END_REF][START_REF] Toh | Transport properties of monolayer mos 2 grown by chemical vapor deposition[END_REF] for MoS 2 , 480 cm 2 V -1 s -1 for MoSe 2 and WSe 2 [START_REF] Fivaz | Mobility of charge carriers in semiconducting layer structures[END_REF] with a clear signature of metallic conduction at high doping densities. These charge carriers diffusion phenomena in TMDC semiconductors depend on a number of factors, in particular the temperature, the electronic band structure, the charge carrier density, the thickness of the semiconductor material. Although this charge carriers mobility in TMDCs is low compared to that of other materials such as, for instance, graphene (which shows a mobility 200 times greater [START_REF] Bolotin | Ultrahigh electron mobility in suspended graphene[END_REF]) they are nonetheless semiconductors with a non-zero bandgap. This feature allows TMDCs to be much better candidates than graphene for a variety of electronic and optoelectronic applications and also potential candidates to replace silicon in next-generation electronics.

Post-process doping: conventional dielectrics and electrolyte gating

There are currently two processes for doping nanomaterials in electronic devices by gating approaches in FET geometry: using conventional dielectrics such as SiO 2 for FET devices or using electrolytes such as ionic liquids, polyelectrolytes and gels with high ionic conductivity for electronic-double-layer-transistors (EDLT-FET). In the case of a conventional dielectric, this makes possible to control the amount of charge carriers in the nanomaterial (holes and electrons) by applying a gate voltage V sd which generates a perpendicular electric field. In the case of an electrolyte, by applying a gate voltage, the same electric field is created and controls not only the density of charge carriers in the 2D semiconductor material but also the density of electrolyte ions at the solid/liquid interface. This allows to obtain a deeper insight into the transport characteristic (electron transport or/and hole transport) and to evaluate the carrier mobility. The linear part of the transconductance, also termed transfer curve, allows to extract the mobility of the carriers according to the formula:

µ = L CW V sd dI sd dV g (1.1)
where L and W are the length and the width of the channel respectively, C the specific capacitance of the dielectrics or electrolytes used [START_REF] Nguetchuissi | 2D materials, colloidal nanocrystals and mixed dimensional heterostructures for nanoelectronics and optoelectronics in high doping regime[END_REF]. This capacitance is a crucial parameter to assess the performance of the device. Conventional dielectrics are generally used as an insulation layer between the gate and substrate of a FET to prevent current flow between the gate and substrate and to maintain the gate structure. In the typical geometry of FETs, we have three contact electrodes, the source (charge carrier emitter) the drain (charge carrier collector) and the gate allowing to control the doping during scanning by applying a gate voltage V sd which subsequently generates a perpendicular electric field. This electric field, in turn, produces the polarization of the dielectric material by orienting the existing permanent dipoles or by inducing temporary dipole moments which disappear when the electric field is removed (Figure 1.7a and Figure 1.7b) [START_REF] Nguetchuissi | 2D materials, colloidal nanocrystals and mixed dimensional heterostructures for nanoelectronics and optoelectronics in high doping regime[END_REF]. This polarization process depends on the sign (positive or negative) and intensity of the gate voltage, so in the two cases we can have an n-type doping (for electrons) or p-type doping (for holes). Hence, a conduction channel is created and this channel is a semiconductor with an excess of electrons or holes. To prevent charge carriers from crossing into the grid, a dielectric (SiO 2 ) with a permittivity ǫ r (Figure 1.7c) acts as a barrier and the thickness of the latter plays a crucial role in FETs. Typical specific capacitance values for SiO 2 dielectric are in the nF/cm 2 range. The capacitance (C) of such a parallel plate capacitor and the charges (n) accumulated on the surface are given by:

C = ǫ 0 ǫ r A d , n = CV e (1.2)
where ǫ 0 is the dielectric constant of the vacuum, A is the capacitor plate area, d the dielectric thickness, V the apllied voltage and e the elementary charge and ǫ r the relative permittivity of the dielectric material [START_REF] Nguetchuissi | 2D materials, colloidal nanocrystals and mixed dimensional heterostructures for nanoelectronics and optoelectronics in high doping regime[END_REF]. To further miniaturize the electronic compounds and at the same time increase their performance, the thickness (d) of the silicon dioxide grid should be reduced, an operation that at the same time increases its capacitance. However, below a thickness of 2 nm, the leakage current by tunnel effect increases drastically, consequently increasing the energy consumption of the compound and reducing its reliability. The idea would then be to replace the silicon dioxide of the gate with a dielectric with a higher dielectric constant k (HfO 2 , ZrO 2 , Ta 2 O 5 , Al 2 O 3 ), a high-k material, making it possible to increase the capacity of the gate while avoiding undesired leakage effects. However, dielectric integration is a significant challenge for TMDCs with non-dangling MX 2 layers due to their chemical incompatibility with typical oxide dielectrics. Furthermore, charged impurities and traps, as well as structural defects in MX 2 also lead to a much lower mobility than the intrinsic limit. Massive efforts have been devoted over the years to solve these problems. In order to reduce the contact resistance across the metal electrode and semiconductor channel (MX 2 ) interface at the source and drain regions, Appenzeller et al. [START_REF] Shi | Recent advances in controlling syntheses and energy related applications of mx 2 and mx 2 /graphene heterostructures[END_REF] used a low-work-function metal electrode extraction (Sc ∼3.5 eV) for the construction of a 10 nm thick exfoliated MoS 2 -FET and obtained an additional effective mobility of ∼700 cm 2 V -1 s -1 at room temperature. Contrary to FETs, in EDLTs liquid or gel electrolytes are used to replace conventional dielectrics, when a gate voltage is applied, this generates a transverse electric field which makes it possible to polarize the ions by aligning their dipole moments parallel to the electric field thus creating an accumulation of ions at electrolyte/electrode interface (Figure 1.8). The accumulation of negatively charged ions at the positive electrode and of positively charged ions at the negative electrode give formation of a electrical double layer at the interface with nanoscale ordering The EDL thickness makes it possible to have a high specific capacity (of the order of a µF/cm 2 ), allowing the accumulation of high charge carriers (10 13 to 10 14 charges/cm 2 ). As an example, Gao et al. [START_REF] Gao | Triboiontronic transistor of mos 2[END_REF] demonstrated that MoS 2 -based triboiontronic transistors that operate by triboelectric potential modulation allows to achieve highly efficient EDL control, obtaining high current ON/OFF ratio over ∼10 7 , low power consumption and high stability.

Atomic-scale modeling of ILs and TMDCs

Computer simulations are nowadays regarded as a powerful and reliable tool to study atomic-scale properties of many-body interacting systems and has a wide range of applications in chemistry, physics, materials science, and biology. When the considered degrees of freedom are the atoms, computer simulations are referred to as "atomic-scale (or atomistic) simulations", making explicit the fact that the description of a given systems is given at the atomic-scale. Within this general context, the power of computer simulations rests on two factors. First, one can access the microscopic properties by using rigorous statistical mechanics techniques, thus ensuring a well-founded description. Second, the prediction of macroscopic properties becomes accessible, in the range and beyond the range defined by pertinent measurements enriching or complementing the experimental counterpart [START_REF] Bouzid | First-principles investigation of binary and ternary amorphous chalcogenide systems[END_REF]. Molecular modelling grants access to properties that are difficult or impossible to reach experimentally and contributes to the understanding of structural organisation at the nano-scale level, providing direct insight into the nature of physical and chemical processes. [START_REF] Welton | Solvents and sustainable chemistry[END_REF]61]

Classical molecular dynamics: advantages and pitfalls

Classical molecular dynamics (CMD) can be considered one of the most widely used atomistic modeling technique to study ILs and, to a certain extend, 2D layered materials. An all-atom CMD simulation can be regarded as a compromise between level of detail and computational workload, attaining size and time scales of 6-10 nm and 10-100 ns, respectively. Such size and temporal scales have been found to represent, to a certain degree, in a reasonable and reliable way to simulate the typical medium-range structural order (10 1 -10 2 nm) and dynamical properties (10 1 -10 2 ns) of ILs [START_REF] Goloviznina | Development and application of molecular interaction models with explicit polarisation for ionic liquids and eutectic solvents[END_REF]. In the last decades, general non-polarisable force fields (FF), such as the All-Atom Optimized Potentials for Liquid Simulations (OPLS-AA) [START_REF] Wl | Development and testing of the opls all-atom force field on conformational energetics and properties of organic liquids[END_REF], the AMBER FF (AMBER/GAFF) [START_REF] Wang | Development and testing of a general amber force field[END_REF] or the CHARMM FF (CHARMM/CGenFF) [START_REF] Vanommeslaeghe | Charmm general force field: A force field for drug-like molecules compatible with the charmm all-atom additive biological force fields[END_REF], have been widely used to model ILs and related systems. Individual terms the analytical expression of such FFs (e.g. partial charges for ions or torsion energy profiles) were often parametrized to fit ab initio calculations, allowing to achieve semi-quantitative agreement with experimental data in terms of structural properties. Non-polarisable FFs based on fixed charges for ILs started to become popular since the seminal works of Lynden-Bell [START_REF] Cg | Intermolecular potentials for simulations of liquid imidazolium salts[END_REF] and Maginn [START_REF] Morrow | Molecular dynamics study of the ionic liquid 1-n-butyl-3-methylimidazolium hexafluorophosphate[END_REF], and were followed by FF systematic construction by Wang [START_REF] Liu | A refined force field for molecular simulation of imidazolium-based ionic liquids[END_REF], Canongia Lopes and Padua (CL&P FF) [START_REF] Canongia | Modeling ionic liquids using a systematic all-atom force field[END_REF][START_REF] Goloviznina | Development and application of molecular interaction models with explicit polarisation for ionic liquids and eutectic solvents[END_REF][START_REF] Canongia | Molecular force field for ionic liquids composed of triflate or bistriflylimide anions[END_REF][START_REF] Canongia | Molecular force field for ionic liquids iii: Imidazolium, pyridinium, and phosphonium cations, chloride, bromide, and dicyanamide anions[END_REF][START_REF] Canongia | Nanostructural organization in ionic liquids[END_REF][START_REF] Agílio | Molecular force field for ionic liquids iv: Trialkylimidazolium and alkoxycarbonyl-imidazolium cations; alkylsulfonate and alkylsulfate anions[END_REF][START_REF] Canongia | Clp: A generic and systematic force field for ionic liquids modeling[END_REF][START_REF] Shimizu | Molecular force field for ionic liquids v: Hydroxyethylimidazolium, dimethoxy 2-methylimidazolium, and fluoroalkylimidazolium cations and bis(fluorosulfonyl)amide, perfluoroalkanesulfonylamide, and fluoroalkylfluorophosphate anions[END_REF] and Acevedo [START_REF] Doherty | Revisiting opls force field parameters for ionic liquid simulations[END_REF][START_REF] Sv | Development of opls-aa force field parameters for 68 unique ionic liquids[END_REF]. The CL&P FF, which includes imidazolium, pyridinium, ammonium and phosphonium families with different functional groups and a broad series of anions, became one of the most benchmarked and widely-used FFs for ILs due to its compatibility with OPLS-AA, robustness and reliability, and numerous efforts to refine it were done in the last decade. Due to its simplicity and high degree of transferability, CL&P FF has been extensively used to study a wide variety of bulk IL as well as IL at the interface with solid substrates (crystalline and amorphous) and in confined environments [START_REF] Mendonza | Novel ionic lubricants for amorphous carbon surfaces: molecular modeling of the structure and friction[END_REF][START_REF] Coasne | Loading-controlled stiffening in nanoconfined ionic liquids[END_REF][START_REF] Ori | Ionic liquid confined in silica nanopores: molecular dynamics in the isobaric-isothermal ensemble[END_REF][START_REF] Ori | Structure and dynamics of ionic liquids confined in amorphous porous chalcogenides[END_REF]. CL&P provides a good insight into the structure of ILs in solid and liquid states, with deviations from experimental densities not exceeding 3% [START_REF] Canongia | Clp: A generic and systematic force field for ionic liquids modeling[END_REF]. However, they systematically and dramatically fail to accurately reproduce diffusion coefficients, viscosity or conductivity, indicating a misrepresentation of different contributions in the interaction potential. [START_REF] Tsuzuki | Molecular dynamics simulations of ionic liquids: Cation and anion dependence of selfdiffusion coefficients of ions[END_REF][START_REF] Chaban | Polarizability versus mobility: atomistic force field for ionic liquids[END_REF] Improvement of the underestimated dynamics of typical nonpolarisable FF can be achieved though scaling the net charge, [START_REF] Chaban | Polarizability versus mobility: atomistic force field for ionic liquids[END_REF][START_REF] Chaumont | Aqueous interfaces with hydrophobic room-temperature ionic liquids: A molecular dynamics study[END_REF][84] modification of the van der Waals potential [START_REF] Benoît | Molecular dynamics simulations of ionic liquids and electrolytes using polarizable force fields[END_REF][START_REF] Koddermann | Molecular dynamic simulations of ionic liquids: A reliabledescription of structure, thermodynamics and dynamics[END_REF][START_REF] Koddermann | Ionic liquids: Dissecting the enthalpies of vaporization[END_REF] or by introducing explicitly in the potential form polarisation terms to mimic the incorporation of 'electronic effects'. [START_REF] Goloviznina | Development and application of molecular interaction models with explicit polarisation for ionic liquids and eutectic solvents[END_REF][START_REF] Benoît | Molecular dynamics simulations of ionic liquids and electrolytes using polarizable force fields[END_REF][START_REF] Goloviznina | Extension of the clpol polarizable force field to electrolytes, protic ionic liquids and deep eutectic solvents[END_REF][START_REF] Goloviznina | The clpol polarizable force field for the simulation of ionic liquids and eutectic solvents[END_REF] This latter approach is becoming more and more exploited to improve the ability of CMD modeling in reproducing and predicting the dynamical properties of ILs. This is achieved by adding to the FF potential form a simplified representation of electronic clouds of surrounding the atoms, in an attempt at mimicking the presence of an electronic structure. Nevertheless, the major disadvantages of polarisable FF are the typical lack of transferability, poor extendability and considerable computational cost of development. Among the polarizable FF proposed in the last decade [START_REF] Borodin | Polarizable force field development and molecular dynamics simulations of ionic liquids[END_REF][91][START_REF] Jg | Ab initio force fields for organic anions[END_REF], the recent one by Goloviznina et al. [START_REF] Goloviznina | Extension of the clpol polarizable force field to electrolytes, protic ionic liquids and deep eutectic solvents[END_REF], termed CL&Pol, has been developed and benchmarked to be transferable and extendable thanks to a refined parametrization based on polarisation effects represented by Drude induced dipoles [START_REF] Benoît | Molecular dynamics simulations of ionic liquids and electrolytes using polarizable force fields[END_REF][START_REF] Goloviznina | The clpol polarizable force field for the simulation of ionic liquids and eutectic solvents[END_REF]. In the present PhD work both non-polarisable CL&P and polarisable CL&Pol have been used sorting out their advantages and limitations. This exploratory phase has been an essential step in this whole project, being a novel research subject in the team I joined and, more generally, worldwide.

Regarding the atomistic CMD modeling of TMDCs, several FFs have been developed for the specific case of crystalline MoS 2 . The initial motivation behind the first FFs developments have been the study of mechanical and thermal properties of MoS 2 crystalline layers. In spite of these developments, the refinement of these FFs to model a composite system made of MoS 2 and ILs is hampered by the lack of adequate models for the interfacial interactions between the two counterparts. As outlined in a recent review [START_REF] Nicolini | A comparaison of empirical potential for sliding simulations of mos 2[END_REF], many of the atomistic models for bulk and single-layer MoS 2 available in the literature are not suitable for the calculation of Coulomb and vdW interactions with other materials, whether they are solids, liquids, or gases [START_REF] Sresht | Quantitative modeling of mos 2 -solvent interfaces: Predicting contact angles and exfoliation performance using molecular dynamics[END_REF]. The four most tested FFs of this type are the ones proposed by Becker et al. [START_REF] Becker | Metal island growth and dynamics on molybdenite surfaces[END_REF], Morita et al. [START_REF] Morita | Development of a new molecular dynamics method for tribochemical reaction and its application to formation dynamics of mos 2 tribofilm[END_REF] Varshney et al. [START_REF] Varshney | Md simulations of molybdenum disulphide (mos 2 ): Force-field parameterization and thermal transport behavior[END_REF] and Liang et al. [START_REF] Liang | Parametrization of a reactive many-body potential for mo-s systems[END_REF]. Despite these limitations, the parameters of the aforementioned FFs have been used to make quantitative predictions of the friction behavior and wettability of MoS 2 surfaces, [START_REF] Luan | Wettability and friction of water on a mos2 nanosheet[END_REF] the interactions of proteins with MoS 2 nanotubes and slabs, [START_REF] Ling | Structural damage of a β-sheet protein upon adsorption onto molybdenum disulfide nanotubes[END_REF] and the desalination of water using nanoporous MoS 2 [START_REF] Heiranian | Water desalination with a single-layer mos2 nanopore[END_REF][START_REF] Li | Tunable, strain-controlled nanoporous mos2 filter for water desalination[END_REF]. More recently, Sresht et al. [START_REF] Sresht | Quantitative modeling of mos 2 -solvent interfaces: Predicting contact angles and exfoliation performance using molecular dynamics[END_REF] developed a novel FF (denoted as SrFF for the purpose this work) proposing a potential that overcomes the limitations of former FFs. SrFF has been developed to incorporate the decoupling between the intralayer and the interlayer interactions of MoS 2 layers. SrFF has also shwon its ability to simulate mechanical and interfacial systems, and it is transferable to aquaeous and nonaqueous systems [2,[START_REF] Aah | Dominance of dispersion interactions and entropy over electrostatics in determining the wettability and friction of two-dimensional mos 2 surfaces[END_REF]. In particular, SrFF has been recently applied to model systems where an IL is in contact with MoS 2 as pristine surfaces or electrified surfaces (work performed on IL made of a series of imidazolium-based cations and dicyanamide (DCA) as anion) [START_REF] Manish | Molecular transport across the ionic liquid-aqueous electrolyte interface in a mos 2 nanopore[END_REF][START_REF] Gong | Effect of side chain modifications in imidazolium ionic liquids on the properties of the electrical double layer at a molybdenum disulfide electrode[END_REF].

First-principles studies

Static and dynamical quantum mechanical calculations within the framework of the Density Functional Theory (DFT) approach and first-principles (ab initio) molecular dynamics (FPMD or AIMD) [START_REF] Hutter | Car-parrinello molecular dynamics[END_REF][START_REF] Car | Unified approach for molecular dynamics and densityfunctional theory[END_REF] represent a powerful techniques to study the local electronic structure and chemical effects and they have been applied to the study of ILs in the early 2000 [START_REF] Zahn | Real-world predictions from ab initio molecular dynamics simulations[END_REF]. [START_REF] Bl | Intermolecular structure and dynamics in an ionic liquid: A car-parrinello molecular dynamics simulation study of 1,3dimethylimidazolium chloride[END_REF]. The radial pair distribution functions between the most acidic proton of the DMIM cation and Cl anion were deeply different in terms of location and width of the peaks. Furthermore, differences between FPMD and CMD in the spatial distribution of chloride ions around the cation were evidenced. The data were explained in terms of the formation of a Hbond between the most acidic hydrogen of the imidazolium ring and the chloride ion. The cation-anion H-bond present in the melt was observed as a red-shift in the C-H stretching frequency, [START_REF] Bl | Intermolecular structure and dynamics in an ionic liquid: A car-parrinello molecular dynamics simulation study of 1,3dimethylimidazolium chloride[END_REF] keeping into account the fact that using an exchange-correlation functional such as the Becke-Lee-Yan-Parr (BLYP) one [START_REF] Lee | Development of the colle-salvetti correlation-energy formula into a functional of the electron density[END_REF][START_REF] Perdew | Pair-distribution function and its coupling-constant average for the spin-polarized electron gas[END_REF] leads to typical underestimation of vibrational frequencies [START_REF] Mp | Infrared spectroscopy of nmethylacetamide revisited by ab initio molecular dynamics simulations[END_REF]. In general, using generalized gradient approximation (GGA) functionals such as BLYP, [START_REF] Lee | Development of the colle-salvetti correlation-energy formula into a functional of the electron density[END_REF][START_REF] Perdew | Pair-distribution function and its coupling-constant average for the spin-polarized electron gas[END_REF] and PBE [START_REF] Jp | Generalized gradient approximation made simple[END_REF] is not recommended for modeling ionic liquids as none of the GGAs proposed so far in the literature accounts for vdW dispersion forces, which are, instead, crucial in producing reliable geometries of both the liquid and the stacked substrate. With the addition or correction for dispersion forces their performance improves drastically reaching reasonable description of ILs structure and chemical bonding interactions. The protic ionic liquid monomethylammonium nitrate was the first IL simulated by FPMD applying dispersion-corrected functionals [START_REF] Zahn | Real-world predictions from ab initio molecular dynamics simulations[END_REF]. This IL was found to exhibit a fast fluctuating H-bond network. Nowadays, hybrid DFT functionals (such as B3LYP, [START_REF] Ad | Density-functional thermochemistry. iii. the role of exact exchange[END_REF] PBE0, [START_REF] Carlo | Toward reliable density functional methods without adjustable parameters: The pbe0 model[END_REF] and revPBE [START_REF] Zhang | Comment on generalized gradient approximation made simple[END_REF][START_REF] Li. Chun | Structure and dynamics of liquid water from ab initio molecular dynamics-comparison of blyp, pbe, and revpbe density functionals with and without van der waals corrections[END_REF]), meta-GGA functionals, and longrange separated hybrid functionals (such as M05-2X [START_REF] Zhao | Design of density functionals by combining the method of constraint satisfaction with parametrization for thermochemistry, thermochemical kinetics, and noncovalent interactions[END_REF]) including different amounts of Hartree-Fock (HF) exchange are commonly used to study ILs with appreciable accuracy in the predictions of interaction energies and thermodynamic properties [START_REF] Seeger | A systematic study of dft performance for geometry optimizations of ionic liquid clusters[END_REF]. However, this (high) level of computation is affordable only for the study of an individual cationanion pair or, at very best, small clusters of a few IL ions pairs. Regarding the study of TMDC, along the years a strong effort has been done in the computational description and reproduction of their mechanical, electronic and optical properties by using hybrid functionals such as the one proposed by Scuseria and coworkers, HSE06 [START_REF] Jk | The indirect to direct band gap transition in multilayered mos 2 as predicted by screened hybrid density functional theory[END_REF], corrected for dispersion forces or higher level calculations. This degree of computational accuracy has been applied to describe the optoelectronic properties and band structure of crystalline bulk, monolayer and up to a few layers of MoS 2 . To date, the FPMD modeling of systems made of Mo 2 surfaces with a liquid and lasting a sufficient long time to catch their main dynamical features can be afforded only by lowering the computational level at which the first principles driving "core" is treated. For instance, Iannuzzi et al. [START_REF] Bilichenko | Ab initio nanofluidics: disentangling the role of the energy landscape and of density correlations on liquid/solid friction[END_REF] recently simulated by Born-Oppenheimer MD with the electronic structure computed at the DFT level using the optB88-vdW [START_REF] Bowler | vdw-df functional of langreth and lundqvist et al[END_REF] exchange-correlation vdW functional a system made of a few layers of water molecules confined between two layers of MoS 2 with a focus on the frictional properties of such system. This is roughly the present state-of-the art of FPMD simulations for these composite and complex systems.

Objectives and organization of the thesis

The overall aim of this work is to exploit atomic scale modeling to achieve a better understanding of the physical and chemical interactions occurring when a IL is in contact with a 2D TMDC, such as MoS 2 . In particular, a great effort has been done to get an accurate description of the structural properties and chemical interactions (e. g. H-bonding) of different AI-ILs at the interface with pristine or defected MoS 2 slabs. Both CMD and FPMD methods are employed in order, on one hand, to cover size and time scales close to what is expected to be realized in real device and, on the other hand, assessing the fundamental role of the electronic structure of the different components. This latter aspect is achievable only by taking into explicit account the electronic degrees of freedom. The organisation of this manuscript is the following:

• Chapter 2 gives an overview of the theoretical background of the two methods used in this work: classical and first-principles molecular dynamics simulations.

• Chapter 3 assesses the structural, hydrogen bonding and dipolar properties of a variety of AI-IL. For a subset of AI-IL investigated, explicit comparison between CMD (using non-polarisable and polarisable FF) and FPMD data is presented and discussed.

• Chapter 4 is dedicated to the modeling of AI-IL at the interface with MoS 2 slabs with a particular emphasis on the structural behaviour of the IL. The role of AI-IL cation chain length, type of anion and the presence of MoS 2 surface defects is analysed in details. For EMIM-BF 4 , a detailed account of the comparison CMD versus FPMD is presented and discussed.

• Chapter 5 is devoted to a preliminary study of modeling a EMIM-BF 4 in contact with MoS 2 electrified slabs. A particular focus is oriented to the analysis of the structure and organization of the electrical double layer formed by the ILs.

• Chapter 6 outlines the overall conclusions of this work and identifies possible line of future research for this topic.

Chapter 2

Computational methods

Classical molecular dynamics (CMD)

Basic equations of CMD

Classical molecular dynamics is a computational method making use of a numerical solution of the Newton's equations of motion for a many-body system to obtain the trajectory of a set of N particles exploring dynamically the phase space of the system and from which equilibrium properties (structural and dynamical) can be obtained [START_REF] Allen | computer simulations of liquids[END_REF][START_REF] Mp Allen | Computer simulations in Chemical Physics[END_REF][START_REF] Dc Rapaport | The art of the molecular dynamics simulations[END_REF]. In a simulation, the analytically continuous temporal evolution is discretized way into timeframes (t 0 , t 1 ,..., t m ) separated by time-steps (δt) intervals, namely t i+1 -t i = δt. A CMD simulation allows then for a deterministic calculation of position and velocity of the particles at each time step by propagating the results obtained at former steps. By using textbook definition, the velocity of each atom i given by equation (2.1)

v i (t) = d r i (t) dt (2.1)
where r i and v i is the atomic position and the acceleration is equation (2.2).

a i (t) = d v i (t) dt (2.2)
hence, the Newton's equations become

F i (t) = m i d v i (t) dt (2.3) F i (t) = m i d 2 r i (t) dt 2 (2.4)
The major ingredient of CMD is the definition of an appropriate potential U ( r 1 , r 2 , . . . , r N ) able to catch the physical properties sought and from which a gradient can be computed to get the related force field (FF),

F i (t) = -∇U ( r 1 , r 2 , . . . , r N ) (2.5)
The FF is generally calibrated by a set of parameters accounting for intra-and intermolecular interactions. To integrate numerically the equations of motion, several algorithms are available, such as the Euler algorithm [START_REF] Peskin | Molecular dynamics by the backward-euler method[END_REF], Verlet algorithms [START_REF] Martys | Velocity verlet algorithm for dissipative-particledynamics-based models of suspensions[END_REF][START_REF] Ip. Omelyan | Optimized verlet-like algorithms for molecular dynamics simulations[END_REF][START_REF] Grubmüller | Generalized verlet algorithm for efficient molecular dynamics simulations with long-range interactions[END_REF], Verlet LeapFrog [START_REF] Mp Allen | Computer simulations in Chemical Physics[END_REF][START_REF] Van Gunsteren | A leap-frog algorithm for stochastic dynamics[END_REF] or the velocity Verlet integration algorithm [START_REF] Mp Allen | Computer simulations in Chemical Physics[END_REF][START_REF] Martys | Velocity verlet algorithm for dissipative-particledynamics-based models of suspensions[END_REF][START_REF] Wc Swope | A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters[END_REF], etc. . In the codes used for this doctoral project, we resort to the velocity Verlet and Verlet Leapfrog algorithms thus we limit our discussion to a brief overview of these specific methods.

The velocity Verlet algorithm is based on the 1 st order Taylor expansion of v i (t+δt) [START_REF] Wc Swope | A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters[END_REF].

The Verlet positions at the same time t can be then obtained by a 2 nd order expansion of r i (t+δt) [START_REF] Wc Swope | A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters[END_REF] as given by equations (2.7) and (2.6).

r i (t + δt) = r i (t) + v i (t)δt + F i (t) 2m i δt 2 + O(δt 3 ) (2.6) v i (t + δt) = v i (t) + [ F i (t) + F i (t + δt)] 2m i δt + O(δt 2 ) (2.7)
This provides a good approximation of the analytical differential equations on short time scales and has three major features: i) It ensures, over large time scales, the conservation of the constants of motion such as the total energy of the system while maintaining its stability; ii) It ensures temporal reversibility; and iii) It is symplectic, i. e. preserves the manyfold area in the momentum-coordinates space.

The Leapfrog algorithm is a 2 nd order integration algorithm providing a stable discretization for constant time step. In this algorithm, velocities and positions of the atoms are calculated at alternate instants (temporal time half-step) as

v i (t + δ 2 ) = r i (t + δt) -r i (t) δt + O(δt 2 ) (2.8) v i (t - δ 2 ) = r i (t) -r i (t -δt) δt + O(δt 2 ) (2.9)
The position of the particles can be then obtained from

r i (t + δt) = r i (t) + v i (t + δt 2 )δt + O(δt 2 ) (2.10)
By combining the above three equations, the velocity becomes

v i (t + δ 2 ) = v i (t - δ 2 ) + F i (t) m i δt + O(δt 3 ) (2.11)
To avoid the need of over-refinement of the time step for these discretization algorithms, constrains are applied for bonding interactions of high-frequency vibrations (typical of bonds terminating with hydrogen atoms). Among the proposed methods, such as SHAKE- [START_REF] Mp Allen | Computer simulations in Chemical Physics[END_REF][START_REF] Smith | Parallel macromolecular simulations and the replicated data strategy: I. the computation of atomic forces[END_REF][START_REF] Smith | Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes[END_REF], RATTLE [START_REF] Andersen | Rattle: A velocity version of the shake algorithm for molecular dynamics calculations[END_REF], or SETTLE [START_REF] Pa | Settle: An analytical version of the shake and rattle algorithm for rigid water models[END_REF] we limit our discussion to the one used in this work: SHAKE. In this algorithm, the unconstrained atomic positions r ′ i (t δt) are used as a starting guess and adjusted iteratively until the following constraint equations are numerically satisfied [START_REF] Mp Allen | Computer simulations in Chemical Physics[END_REF]

r ′ i (t + δt) = 2 r i (t) -r i (t -δt) + δt 2 m i F i (t) (2.
12)

The corrected positions at the end of the time step are

r i (t + δt) = 2 r i (t) -r i (t -δt) + δt 2 m i ( F i (t) -G i [t, ṙ(t)]) (2.13)
where G[t, ṙ(t)] is the coercion force. Thus, for any atom i we get

r i (t + δt) = r ′ i (t + δt) - δt 2 m i M α=1 λ α ∂σ α ( r(t)) ∂r i
(2.14)

With σ α ( r(t)) being the constrains. Finally, the atomic coordinates at t + δt must also satisfy the condition

σ α ( r(t + δt)) = 0 (2.15)

Force fields for ILs and TMDCs

The quality and accuracy of the FF used are crucial in any CMD simulation. For this reason we give here below a quick outline of the most common analytic forms.

Non-polarizable FFs for ILs

The most common FFs available nowadays for organic-based systems are relatively well standardized. Specifically, AMBER [START_REF] Zhong | Improved classical united-atom force field for imidazolium-based ionic liquids: Tetrafluoroborate, hexafluorophosphate, methylsulfate, trifluoromethylsulfonate, acetate, trifluoroacetate, and bis(trifluoromethylsulfonyl)amide[END_REF][START_REF] Cornell | A second generation force field for the simulation of proteins, nucleic acids, and organic molecules[END_REF], CHARMM [START_REF] Cadena | Why is co 2 so soluble in imidazolium-based ionic liquids?[END_REF][START_REF] Mackerell | Empirical force fields for biological macromolecules: Overview and issues[END_REF], GROMOS [START_REF] Ta | Validation of the gromos force-field parameter set 45a3 against nuclear magnetic resonance data of hen egg lysozyme[END_REF][START_REF] Scott | The gromos biomolecular simulation program package[END_REF][START_REF] Oostenbrink | A biomolecular force field based on the free enthalpy of hydration and solvation: The gromos force-field parameter sets 53a5 and 53a6[END_REF] and GAFF [START_REF] Kg | The general amber force field (gaff) can accurately predict thermodynamic and transport properties of many ionic liquids[END_REF] are the standard ones for modeling biomolecular systems (mainly proteins and nucleic acids), whereas OPLS-AA [START_REF] Wl | Development and testing of the opls all-atom force field on conformational energetics and properties of organic liquids[END_REF][START_REF] José | Modeling ionic liquids using a systematic all-atom force field[END_REF][START_REF] Wl | Optimized intermolecular potential functions for liquid hydrocarbons[END_REF][START_REF] Wl | Transferable intermolecular potential functions for water, alcohols, and ethers. application to liquid water[END_REF] is the best suited for organic liquids and amorphous systems. The CMD simulations using non-polarizable FFs have been found to provide important insights into molecular correlations and structure of ionic systems such as ILs. In particular, the CL&P FF developed by Canongia Lopes and Padua [START_REF] Canongia | Modeling ionic liquids using a systematic all-atom force field[END_REF][START_REF] Canongia | Molecular force field for ionic liquids composed of triflate or bistriflylimide anions[END_REF][START_REF] Canongia | Molecular force field for ionic liquids iii: Imidazolium, pyridinium, and phosphonium cations, chloride, bromide, and dicyanamide anions[END_REF][START_REF] Canongia | Nanostructural organization in ionic liquids[END_REF] has shown its ability in simulating imidazolium, pyridinium, ammonium and phosphonium families with different functional groups and a large series of anions. According to the literature, this is one of the most exploited FF for ILs due to its compatibility with OPLS-AA, by Jorgensen et al. [START_REF] Wl | Development and testing of the opls all-atom force field on conformational energetics and properties of organic liquids[END_REF][START_REF] Ls | Ligpargen web server: an automatic opls-aa parameter generator for organic ligands[END_REF][START_REF] Mj. Robertson | Development and testing of the opls-aa/m force field for rna[END_REF] because of its robustness and reliability. Systematic efforts have been proposed during the last decade to refine this FF. Depending on the model assumptions, the net ionic charge in a non-polarizable FF can be either integer or fractional. The models with unit ionic charges provide a good insight into the structure of ILs in solid and liquid phases, with deviations from experimental densities not exceeding 3% for CL&P. Hence, CL&P can be regarded as a standard FF and will be briefly presented in the ongoing discussion because it is the one used in the present work.

The OPLS-AA and CL&P FF accounts for two types of interactions, bonding and nonbonding, as shown in

Bonding interactions

The bonding intramolecular interactions represent the different internal degrees of freedom of the molecule studied and they consist of the typical two, three and four-body interactions, namely stretching, bending (or rotation) and dihedral (torsional), respectively. The-two body potential has a harmonic form in which the bond between atom i and j is defined as a spring of stiffness k and depends on the (scalar) interatomic distance only. Its general form is the one of equation (2.16) where r 0 ij is the equilibrium distance, r ij the actual distance between atoms i and j and k ij is the related elastic constant:

U (r ij ) = 1 2 k ij (r ij -r 0 ij ) 2 (2.16)
Also the three-body potential has a harmonic form describing the bending angle θ of three atoms i, j and k as written in by equation (2.17) where θ 0 ijk is the equilibrium angle, θ ijk its actual value and K θ is the angular stiffness constant:

U (r ijk ) = 1 2 K θ (θ ijk -θ 0 ijk ) 2 (2.17)
The four-body potential is defined by a dihedral angle potential and it describes torsional forces acting along a chain of four atoms bonded to each other. The dihedral angle φ ijkl describes the rotation of the bonds ij and kl around the axis jk. Its general form is given by equation (2.18) where V ijkl represents the Fourier constant and m runs over each index for the four atoms:

U (r ijkl ) = 1 2 4 m=1 V ijkl [1 + (-1) m cos(mφ ijkl )] (2.18)

Non-bonding interactions

The non-bonding intermolecular interactions include van der Waals (vdW) and long-range electrostatic interactions. The vdW interactions are expressed in terms of the Lennard-Jones potential (LJ 12-6) with the general form given by the equation (2.19).

U (r ij ) = 4ǫ ij   σ ij r ij 12 - σ ij r ij 6   (2.19)
where the parameters σ ij represent the distances between atoms i and j, and ǫ ij the well depth. The LJ parameters of different atoms are usually obtained by using the Lorentz-Berthelot combining rules of like-atom pairs [START_REF] Delhommelle | Inadequacy of the lorentz-berthelot combining rules for accurate predictions of equilibrium properties by molecular simulation[END_REF]:

σ ij = σ ii + σ jj 2 , ǫ ij = √ ǫ ii ǫ jj (2.20)
The σ ij and ǫ ij values used in the present work are the ones proposed by Canongia Lopez et al. [START_REF] José | Modeling ionic liquids using a systematic all-atom force field[END_REF] The LJ potential becomes 0 when σ ij = r ij and has a minimum at the distance r ij = 2 1/6 σ ij . The 6 power term represents the attractive long-range interaction (dominant at large distances), whereas the 12 power one is the short-range repulsive interaction (dominant at short range) accounting for the Pauli's principle as shown in Figure 2.2. In OPLS-AA and in CL&P FF the (short-range) electrostatic interaction is described by a Coulomb potential that takes into account the electrostatic charge of each atom, and it can be attractive or repulsive depending on the values of q i and q j as

U (r ij ) = 1 4πǫ 0 q i q j r ij (2.21)
where q i , and q j are electrostatic charges for atoms i and j respectively, ǫ 0 is the vacuum dielectric constant, and

r ij = | r i -r j | is the interatomic distance.
Nevertheless, longrange electrostatic correction decaying with a power law slower than r -3 are not taken into account by the Coulomb potential. A typical method to describe the long-range electrostatic interactions is the Ewald summation [START_REF] Toukma | Ewald summation techniques in perspective: a survey[END_REF]153]. This is what is usually done to compute the electrostatic interactions in a periodic system [START_REF] Smith | The DL POLY Classic user manual[END_REF]. The total electrostatic energy is expressed as

U c = 1 2V 0 ǫ 0 ǫ ∝ k =0 exp - k 2 4α 2 k 2 N j q j exp(-i k. r j ) 2 + 1 4πǫ 0 N * n<j q j q n r nj erf c(αr nj ) - 1 4πǫ 0 ǫ molecules M * l m q l q m δ lm α √ π + erf (αr lm ) r 1-δ lm lm - 1 4πǫ 0 ǫ    F * l m q l q m δ lm α √ π + erf (α lm ) r 1-δ lm - 1 V 0 α 2    N j q j    2    (2.22)
where N is the number of charged atoms in the system and N * the same number discounting any excluded (intramolecular) interactions. M * is the number of excluded atoms in a given molecule and includes the atomic self correction, F ⋆ represents the number of frozen atoms in the MD cell. The final term is the Fuchs correction for charged systems [START_REF] Fuchs | A quantum mechanical calculation of the elastic constants of monovalent metals[END_REF]. V 0 is the simulation cell volume and k is a reciprocal lattice vector defined by:

k = l u + m v + n w (2.23)
where l m, n are integers and u, v, w are the reciprocal space basis vectors.

V 0 = |a.b × c| (2.24)
In practice, the convergence of the Ewald sum is controlled by three variables: the real space cutoff r cut ; the convergence parameter α and the largest reciprocal space vector k max used in the reciprocal space sum [START_REF] Smith | The DL POLY Classic user manual[END_REF].

The force on an atom j is obtained by direct calculation of the gradient, resulting in

F j = 1 V 0 ǫ 0 ∝ k =0 i k exp(i k. r j ) exp - k 2 4α 2 k 2 N n q n exp(-i k. r n ) + q j 4πǫ 0 N * n q n r 3 nj erf c(αr nj ) 2αr nj √ π exp (-α 2 r 2 nj ) r nj - q j 4πǫ 0 M * l q l r 3 lj erf (αr lj ) - 2αr lj √ π exp (-α 2 r 2 lj ) r lj (2.25)
The Ewald method splits the (electrostatic) sum for an infinite periodic system into a damped real space sum and a reciprocal space sum. The convergence rate of both sums is governed by α. Evaluation of the real space sum is truncated at r = r cut so it is crucial to select α in a way that contributions to the real space sum are negligible for r > r cut [START_REF] Smith | The DL POLY Classic user manual[END_REF].

The relative error ε in the real space sum truncated at r cut is approximated by

ε ≈ 1 r cut exp [-(αr cut ) 2 ] (2.26)
and the relative error in the reciprocal space term is

ε ≈ 1 k 2 max exp - k 2 max 4α , k max = 2π L N max (2.27)
where k max is the largest k-vector considered in reciprocal space, L is the width of the cell in the specified direction and M max is an integer.

Polarizable FFs for ILs

The first attempts to model ionic systems such as ILs using traditional non-polarizable FF evidenced challenges for a quantitative prediction of transport properties. Several studies indicated that special attention should be paid to polarization effects in systems with high ionic concentrations. Non-polarizable FF systematically underestimate diffusion coefficients and fail to predict other dynamical proprieties such as viscosity or conductivity, indicating a misrepresentation of contributions in the interaction potential. As mentioned in Chapter 1, improvement can be achieved trough scaled-net-charges or modification of the vdW potential. In the first approach, denoted as scaled-charge models, the ionic charge is reduced to represent charge transfer and/or polarisation effects observed originally for isolated cation-anion pairs in the gas phase. Scaling factors of a typical range of 0.5-0.9 are obtained from first-principles calculations of charge distributions on a single ion pair or small clusters or chosen empirically [START_REF] Chaban | Polarizability versus mobility: atomistic force field for ionic liquids[END_REF]. For example, Bhargava et al. [START_REF] Bl | Intermolecular structure and dynamics in an ionic liquid: A car-parrinello molecular dynamics simulation study of 1,3dimethylimidazolium chloride[END_REF] proposed a universal value of 0.8 for scaling the charges in the CL&P force field, which increased the individual ionic diffusion coefficients of by a few order of magnitude, improving the agreement with experiments. The enhanced dynamics in a system with reduced charges can be explained by the weakening of cation-anion electrostatic interactions, associated to the decrease in cohesive energy. The second approach to modify non-polarisable FFs consists in keeping integer ionic charges and using empirical adjustment of the vdW potential. With this approach, Ludwig et al. [START_REF] Koddermann | Molecular dynamic simulations of ionic liquids: A reliabledescription of structure, thermodynamics and dynamics[END_REF][START_REF] Koddermann | Ionic liquids: Dissecting the enthalpies of vaporization[END_REF] reparametrized the CL&P FF for imidazolium-based ILs improving the prediction vaporisation heat, self-diffusion coefficients, viscosity and reorientational correlation times and resulting in a better agreement with experimental values. This latter refinement of the (non-polarisable) CL&P FF corresponds to one of the two FFs used in the present work.

A different approach consists in the explicit inclusion of polarization effects to improve FFs for ILs. This research line has been pioneered in the last few years with a recent detailed review summarising the main features [START_REF] Benoît | Molecular dynamics simulations of ionic liquids and electrolytes using polarizable force fields[END_REF]. Polarization can be represented explicitly using different methods, mainly fluctuating charges, induced point dipoles, or Drude-induced dipoles. In particular, the latter one is the most benchmarked nowadays for ILs and ionic systems and it is implemented in several codes. This approach is the one used by Goloviznina et al. [START_REF] Goloviznina | Extension of the clpol polarizable force field to electrolytes, protic ionic liquids and deep eutectic solvents[END_REF][START_REF] Goloviznina | The clpol polarizable force field for the simulation of ionic liquids and eutectic solvents[END_REF][START_REF] Goloviznina | Transferable, polarizable force field for ionic liquids[END_REF] for the development of the polarizable FF CL&Pol, derived from the CL&P fixed-charge FF that describes most families of ILs, in a form compatible with OPLS-AA, one of the major FFs for organic compounds. In the CL&Pol, Drude-induced dipoles are used with parameters determined from atomic polarizabilities. The original CL&P FF is modified with the inclusion of the Drude dipoles to avoid double-counting of polarization effects. This modification makes use of first-principles calculations of the dispersion and includes vdW interactions using symmetry-adapted perturbation theory (SAPT) for a set of dimers composed of positive, negative, and neutral fragments representative of a wide variety of ILs. A fragment approach provides transferability, allowing the representation of a wide spectrum of ILs cation and anion, including different functional groups, without the need to reparametrize (see blocks diagram in Figure 2.3) [START_REF] Goloviznina | Transferable, polarizable force field for ionic liquids[END_REF]. The CL&Pol FF uses the Drude induced-dipole method [START_REF] Pl | The theory of optics[END_REF][START_REF] Lamoureux | A simple polarizable model of water based on classical drude oscillators[END_REF][START_REF] Lamoureux | Modeling induced polarization with classical drude oscillators: Theory and molecular dynamics simulation algorithm[END_REF] to account for polarizable effects and it requires atomic polarisability for the FF set up, although an additional particle per polarisable atom is introduced. A Drude induced dipole is made by two point charges of opposite sign, a positively-charged Drude core (Dc), at the centre the atomic site r, and a negatively-charged Drude particle (Dp), connected to its core by a harmonic spring with an equilibrium distance d of zero. The charges pair induces a dipole under an external electrostatic field E, equation (2.28), and the magnitude is determined by the Drude charges, q D,i , and by the force constant of the spring, κ D , which are related through the polarisability equation (2.29).

µ i = q D,i × d i = q 2 D,i κ D E (2.28
)

α i = q 2 D,i κ D (2.29)
A displacement of a Dp from its Dc should not exceed a typical value of d = 0.10-0.15 1 Å being sufficiently small compared to any interatomic distance, allowing the µ to be consid-ered as a point dipole. The Coulomb interaction of the Dc-Dp pairs with atomic partial charges and induced dipoles then can be evaluated via point charges. The CL&Pol FF is based on the atomic polarisabilities computed by Schröder [START_REF] Heid | Additive polarizabilities of halides in ionic liquids and organic solvents[END_REF][START_REF] Heid | Langevin behavior of the dielectric decrement in ionic liquid water mixtures[END_REF][START_REF] Heid | Quantum mechanical determination of atomic polarizabilities of ionic liquids[END_REF]. In the CL&Pol FF, all heavy atoms are considered as polarizable, while the polarizability of hydrogen atoms is merged into the polarizability of the atoms to which they are bonded. The mass of Dp is set as mD = 0.4 au and the force constant of the harmonic spring between Dc and DP set is κ D = 4184 kJ•mol -1 [START_REF] Goloviznina | Development and application of molecular interaction models with explicit polarisation for ionic liquids and eutectic solvents[END_REF][START_REF] Goloviznina | Transferable, polarizable force field for ionic liquids[END_REF]. Thole damping functions are used in CL&Pol FF to reduce at short range the electrostatic interactions between induced dipoles, to avoid excessive correlation between neighboring Dp. In the present thesis work, the parameters for CL&Pol FF were taken directly from [START_REF] Goloviznina | Extension of the clpol polarizable force field to electrolytes, protic ionic liquids and deep eutectic solvents[END_REF][START_REF] Goloviznina | The clpol polarizable force field for the simulation of ionic liquids and eutectic solvents[END_REF][START_REF] Goloviznina | Transferable, polarizable force field for ionic liquids[END_REF] or obtained following the procedure reported by the developers [163,164].

FFs for TDMCs

CMD simulations offer a general methodology to study 2D layered materials and solid/liquid interfaces involving different types of interactions. In the last decade, several FFs (or the combination of different FFs) have been applied to describe ILs at the interface with solid substrates. However, different performances and degree of accuracy have been reported for a variety of key properties measured experimentally (Table 2.1). For instance, computed surface energies and elastic constants deviate more than 50% with respect to experiments.

Cohesive energies of small molecules on heavy metal surfaces and wettability properties such as water contact angle can have errors amounting even to 100% depending on the type of model used [START_REF] Liu | Interpretable molecular models for molybdenum disulfide and insight into selective peptide recognition[END_REF]. The FF developed by Sresht et al. [START_REF] Sresht | Quantitative modeling of mos 2 -solvent interfaces: Predicting contact angles and exfoliation performance using molecular dynamics[END_REF] in 2017 found several promising applications also to the study of the interactions of MoS 2 with liquids, whether aqueous or organic (IL and polymer). This FF (denoted as srFF hereafter) allows to reproduce the properties of bulk 2H-MoS 2 , with special attention to the distinction between the covalent, intralayer terms and the noncovalent, interlayer Coulombic and van der Waals interactions. The srFF is compatible with MD FFs for liquids and consists of partial charges and 6-12 LJ parameters for S and Mo, as well as of harmonic potentials for bonds and angles. It has been parameterized over DFT atomic charges and experimental lattice parameters, elastic constants, bulk modulus, and Young's modulus of 2H-MoS 2 [START_REF] Sresht | Quantitative modeling of mos 2 -solvent interfaces: Predicting contact angles and exfoliation performance using molecular dynamics[END_REF]. This parametrization has been validated by contact-angle simulations (with water and apolar liquids) to demonstrate the ability of the srFF to simulate both mechanical and interfacial systems, and its transferability to nonaqueous systems.

Thermodynamic ensembles for molecular simulations

A thermodynamic ensemble is a statistical collection of particles in static equilibrium allowing to evaluate the thermodynamic properties of a targeted real system. To understand the process of energy transfer between particles in a microscopic system and its macroscopic effects, the whole system is subdivided into several sub-parts called accessible states and each state is governed by one of the large thermodynamic ensembles: microcanonical ensemble, canonical ensemble, isobaric-isothermal ensemble, grand canonical ensemble. In the following the main aspects of each ensemble are briefly presented: Microcanonical ensemble (NVE): NVE ensemble is a physical system where the total volume V, the number of particles N and the energy E N V E remain constant and, in this respect, can be regarded as an ensemble completely isolated from any external environment in the sense that no exchange of matter or energy occurs outside the ensemble. For this ensemble, the fundamental assumption of statistical mechanics is that every compatible stationary quantum micro-state is equally probable. Properties of a system in a macro-state are calculated by averaging its values over the ensemble micro-states [START_REF] Rf | Microcanonical ensemble[END_REF].

The conservation of the internal energy is a direct consequence of the absence of external

E(N, V, E) = E K + E P (2.30)
where E K and E P represent the kinetic and potential energy, respectively. Canonical ensemble (NVT): In the NVT ensemble the number of particles N, the volume V and the temperature T are constant. The system is in thermal equilibrium and energy can be transferred between a system and its environment, whereas particles are not exchanged. The Canonical ensemble assumes that the system is immersed in a heat bath at a given temperature T (thermostat), that keeps the temperature constant. The Helmholtz free energy F (T, V, N ) is at its minimum at the thermodynamic equilibrium and is defined as

F (N, V, T ) = E(N, V, E) -T × S (2.31)
where T is the absolute temperature, E(N,V,E) is the internal energy (equation (2.30) ) and S is the entropy. Isobaric-Isothermal ensemble (NPT): In the isobaric-isothermal ensemble the number of particles N pressure P and temperature T are constant. A thermostat is used to conserve the temperature, and a barostat to maintain a constant pressure. The volume V is allowed to change, whether isotropically or anisotropically, and the only invariant physical quantity is the Gibbs free energy G(N, P, T ) expressed as

G(N, P, T ) = F (N, V, T ) + P × V (2.32)
where F(N, V, T) is the Helmholtz free energy, P and V are the pressure and volume, respectively (equation (2.31) ).

Thermostats and Barostats

In classical MD to simulate a system at constant temperature a thermostat is applied to control the temperature and a barostat is used to keep the system on average at a given external pressure (stress tensor) during the dynamics. Several thermostats and barostats have been proposed over the years: Nosé-Hoover thermostat [START_REF] Nosé | A molecular dynamics method for simulations in the canonical ensemble[END_REF][START_REF] Nosé | A unified formulation of the constant temperature molecular dynamics methods[END_REF][START_REF] Wg | A molecular dynamics method for simulations in the canonical ensemble[END_REF][START_REF] Wg | Canonical dynamics: Equilibrium phase-space distributions[END_REF], Andersen thermostat [START_REF] Hc | Molecular dynamics simulations at constant pressure and/or temperature[END_REF], Berendsen thermostat [START_REF] Hjc. Berendsen | Molecular dynamics with coupling to an external bath[END_REF], Hoover barostat [START_REF] Wg | Canonical dynamics: Equilibrium phase-space distributions[END_REF][START_REF] Wg | Constant-pressure equations of motion[END_REF], Andersen barostat [START_REF] Hc | Molecular dynamics simulations at constant pressure and/or temperature[END_REF], Parrinello-rhaman [START_REF] Parrinello | Crystal structure and pair potentials: A moleculardynamics study[END_REF][START_REF] Parrinello | Polymorphic transitions in single crystals: A new molecular dynamics method[END_REF][START_REF] Parrinello | Strain fluctuations and elastic constants[END_REF] and Berendsen barostat [START_REF] Hjc. Berendsen | Molecular dynamics with coupling to an external bath[END_REF]. In the present doctoral project we make use of the Nosé-Hoover thermostat and Hoover barostat.

Nosé-Hoover thermostat

In the Nosé-Hoover approach [START_REF] Nosé | A unified formulation of the constant temperature molecular dynamics methods[END_REF][START_REF] Wg | Canonical dynamics: Equilibrium phase-space distributions[END_REF] a single additional (dynamical) variable s is used to complement the Hamiltonian H N V T and to describe the heat bath, thus obtaining

H N V T = i P 2 i 2m i s 2 + 1 2 ij,i =j U (r i -r j ) + p 2 s 2Q + gk b T ln(s) (2.33)
where g is the number of system-independent degrees of freedom, p s is the conjugate momentum, r i and P i are coordinate and momentum of particle i respectively, Q is a fictitious mass , and k b is the Boltzmann constant.

From above Hamiltonian we can get the Nosé equations of motions:

Ṙi = P i m i s 2 , Ṗi = F i , ṡ = p s Q , ṗs = 1 s i P 2 i m i s 2 -gk b T (2.34)
By using this original idea, Hoover has shown that these equations can be simplified by using a scaled momentum and a scaled time:

p ′ s = p s s , ζ = p s Q , t ′ = t 0 1 s dt (2.35)
In the scaled momentum p ′ and scaled time t ′ the equations of motion become

Ṙ′ i = P ′ i m i , Ṗ′ i = F i -ζP ′ i , Q = gk b T ext τ 2 T (2.36)
where τ T is an arbitrary time constant governing the rate at which the equations of motion damp out fluctuations in the kinetic energy (normally in the range [0.5, 2]ps)

The thermostat friction coefficient ζ is controlled by the equation (2.37):

dζ dt ′ = 1 Q i P ′ 2 i m i -gk b T ext =gk b T (t) -T ext Q (2.37)
where T (t) is is the instantaneous temperature of the system at time t and T ext is the external target temperature around which the system oscillates (constant on average). The equations (2.36) and (2.37) are the Nosé-Hover equations for the thermostat. However in the Nosé-Hoover algorithm, a careful choice of fictitious mass Q is crucial, since too large values of Q (loose coupling) can result in a bad temperature control. Although any finite (positive) mass is sufficient to guarantee in principle the generation of a canonical ensemble, excessively large values of Q will be able to reach a canonical distribution only after very long simulation times. Conversely, too small Q values (tight coupling) can cause high frequency temperature oscillations, the variable s would oscillate at a very high frequency, and in general it will be out of resonance with the characteristic frequencies of the real system, thus giving a thermostat decoupled from the particles motion (slow exchange of kinetic energy).

Hoover barostat

In the Hoover approach [START_REF] Wg | Constant-pressure equations of motion[END_REF] the Hamiltonian H N P T is given by :

H N P T = i P 2 i 2m i + 1 2 i =j U (r i -r j ) + 1 2 Qζ(t) 2 + ̺ ext V (t) + 1 2 W η(t) 2 + t 0 Qζ(t ′ ) τ 2 T + k b T ext dt ′ (2.38)
where Q and W are the fictitious mass of the thermostat and the barostat, η is the barostat friction coefficient, ζ is the thermostat friction coefficient, V is the system volume, τ T and τ P are specified time constants for temperature and pressure fluctuation respectively, ̺(t) and T (t) are the instantaneous pressure and temperature, ̺ ext and T ext are the external pressure and temperature. From above Hamiltonian, for isotropic fluctuations the equations of motion are given by:

dP i (t) dt = F i (t) -[ζ(t) + η(t))] P i (t) , dr i (t) dt = P i (t) m i + η[r i (t) + R 0 ] (2.39)
with R 0 the system centre of mass. The thermostat friction coefficient ζ is controlled by the equation (2.40)

dζ(t) dt = gk b T (t) -T ext Q + W η(t) 2 -k b T ext Q , Q = gk b T ext τ 2 T (2.40)
The barostat friction coefficient η is controlled by the equation (2.41)

dη(t) dt = 3V (t) (̺(t) -̺ ext ) W -ζ(t)η(t) , W = gk b T ext τ 2 P (2.41)
and the volume of system by the equation (2.42)

dV (t) dt = 3η(t)V (t) (2.42)
As specified for the Nosé-Hoover thermostat, it is important to make a good choice of the fictitious mass W . The equations (2.39), (2.40), (2.41), (2.42) are the Hoover equations for the barostat.

Density functional theory-based first-principles methods 2.2.1 Density functional theory in a nutshell

Density functional theory (DFT) has become a reliable many-body formulation of quantum mechanics suitable for applications in a wealth of domains, from solid-state physics to materials science and biochemistry. It allows for the calculation of the electronic structure (ground state) of atoms, molecules and extended systems at an affordable computational workload. The original DFT idea derives from the Thomas and Fermi [START_REF] Fermi | Statistical method to determine some properties of atoms[END_REF][START_REF] Lh | The calculation of atomic fields[END_REF] electron "gas" formulation assuming a Fermi-Dirac statistics. Hohenberg, Kohn, and Sham established in the mid-1960s the grounding foundations of DFT [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF] in which electrons are the quantum objects whereas nuclei are treated as classical point-like particles. In the following a brief outline of the DFT formulation is recalled. Hohenberg and Kohn (HK) [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF] demonstrated that the total energy of a multi-electron system under the action of an external potential V ext is determined by the overall electron density of the ground state ρ(r):

E[ρ(r)] = F HK [ρ(r)] + ρ(r)V ext (r) d 3 r (2.43)
where F HK is the HK universal functional independent from the external potential and it entails the quantum complex many-body effects, although its explicit form is unknown. The exact ground state density of the system is given by the specific distribution ρ(r) that minimizes the total energy E[ρ(r)]. Hence, the many-particle wave function can be replaced by a ground state electronic density to simplify the electronic structure problem.

Although the HK theorems constitute a solid basis for DFT, reducing the problem from wavefunction to density approach, they did not provide an explicit way to obtain the ground state density. Subsequently, Kohn and Sham (KS) [START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF] established a procedure to compute the ground state by minimizing the energy functional. In the KS procedure, authors proposed to replace the original interactive problem by a non-interactive one where the instantaneous electrons repulsion is replaced by an effective potential (as in a mean-field approach).This is achieved by the decomposition of the electronic charge density in a set of one-particle orthonormal orbitals (similar to that of the Hartree-Fock approximation) where the only constraint is that their square module sum must be equal to ρ(r). This resumes to use a fictitious non-interacting system that produces the same density of its interacting counterpart. The energy functional is then expressed as:

E[ρ(r)] = T KS [ρ(r)] + E H [ρ(r)] + E XC [ρ(r)] + ρ(r)V ext (r) d 3 r (2.44)
where T KS [ρ(r)] is the kinetic energy functional of non-interacting mono-electronic orbitals:

T KS [ρ(r)] = - 1 2 N i ψ * i (r)∇ 2 ψ i (r) d 3 r (2.45)
with density ρ(r) expressed as the sum of one-electron KS orbitals ψ i and f i being the occupation number of the i th orbital:

ρ(r) = - N i f i |ψ i (r)| 2 (2.46) E H [ρ(r)
] is the electrostatic Hartree interaction:

E H [ρ(r)] = 1 2 ρ(r ′ )ρ(r) |r -r ′ | d 3 r ′ d 3 r (2.47)
The E XC [ρ(r)] is the exchange and correlation energy accounting for the complexity of the unknown part of F HK functional, although its explicit form is unknown.

E XC [ρ(r)] = T [ρ(r)] -T KS [ρ(r)] + E ext [ρ(r)] -E H [ρ(r)] (2.48)
where

E ext [ρ(r)] = -ρ(r) N I=1 Z I |r -R I | d 3 r + N I<J Z I × Z J |R I -R J | (2.49)
The effective potential becomes

V ef f (r) = V ext (r) + V H (r) + V XC (r) (2.50)
where the first term V ext is an external potential due to the ions and, in this respect, the term "external" (currently used in the literature) is not entirely correct, since ions are part of the system and are just "external" to the electronic degrees of freedom. The second term V H is the Coulomb electrostatic repulsion among the electrons (Hartree potential) and the last term V XC is the exchange-correlation potential. These terms are summarized by (2.51), (2.52) and (2.53).

V ext (r) = - N I=1 Z I |r -R I | (2.51) V H (r) = 1 2 ρ(r ′ ) |r -r ′ | d 3 r ′ (2.52) V XC (r) = δE XC [ρ(r)] δρ(r) (2.53) - 1 2 ∇ 2 + V ext (r) + V H (r) + V XC (r) ψ i (r) = ǫ i ψ i (r) - 1 2 ∇ 2 + V ef f (r) ψ i (r) = ǫ i ψ i (r) H KS e ψ i (r) = ǫ i ψ i (r) (2.54)
which are one-electron equations involving an effective one-particle Hamiltonian H KS e with the local effective potential V ef f .

Exchange and correlation approximations

DFT is, by construction, an exact theory for the ground state. Nevertheless, the exact analytic form of E XC [ρ(r)] remains unknown. Many approximations have been proposed to get an (approximated) analytical formulation. Their success depends on the achievable accuracy. In this respect, the functional proposed by Perdew et al. [START_REF] Perdew | Jacob's ladder of density functional approximations for the exchange-correlation energy[END_REF] was a remarkable step toward the so-called chemical accuracy. Depending on their formulation, XC functionals can be either empirical (derived by fitting known data on atoms or molecules) or non-empirical based on some physical properties. In the following, the most common XC functionals used in this work will be briefly presented. Local Density Approximation (LDA): The LDA represents the simplest approximation for E XC [ρ(r)] and consists in assuming that the exchange-correlation energy per each infinitesimal volume of the electron density of the real inhomogeneous system is supposed to be uniform (homogeneous electron gas) [START_REF] Kohn | Quantum density oscillations in an inhomogeneous electron gas[END_REF] as:

E LDA xc [ρ(r)] = ρ(r)ǫ xc (ρ(r)) d 3 r (2.55)
The exchange-correlation energy per particle of a homogeneous electron gas ǫ xc can be separated into an exchange and a correlation contribution (ǫ x and ǫ c respectively). The exchange energy E LDA xc [ρ(r)] for a homogeneous electron gas has an analytical formulation [START_REF] Bloch | Bemerkung zur elektronentheorie des ferromagnetismus und der elektrischen leitfähigkeit[END_REF][START_REF] Dirac | Note on exchange phenomena in the thomas atom[END_REF] whereas the correlation part is not known explicitly but can be approximated as done by Ceperley and Alder via Quantum Monte Carlo calculations [START_REF] Ceperley | Ground state of the electron gas by a stochastic method[END_REF] and then parametrized. [START_REF] Perdew | Self-interaction correction to density-functional approximations for many-electron systems[END_REF] Despite the fact that LDA is intuitively expected to be valid only for systems having slowly varying densities, it was found to have reasonable performance even for strongly inhomogeneous systems (despite some limitations). This is often due to an error compensation between the overestimation of exchange and the underestimation of correlation parts. General Gradient Approximation (GGA): An improvement over the LDA approximation, although not being an analytical expansion, is to include beside the density its gradient. This is what is generally done in Generalized Gradient Approximations functionals (also referred to as semilocal). The accuracy of energy barriers, molecular geometries, atomization energies and total energies improves and the general form reads:

E GGA xc [ρ, ∇ ρ ] = f GGA (ρ(r), ∇ ρ (r)) d 3 r (2.56)
Typically, the contributions to the exchange and correlation terms are calculated separately as in:

E GGA xc [ρ, ∇ ρ ] = E GGA x [ρ, ∇ ρ ] + E GGA c [ρ, ∇ ρ ] (2.57) with the E GGA x [ρ, ∇ ρ ] being defined as E GGA x [ρ, ∇ ρ ] = ρ(r)ǫ x [ρ(r)]F GGA x (s) d 3 r (2.58)
with s being the dimensionless reduced density gradient, or, better, its modulus,

s(r) = |∇ρ(r)| 2k F (r)ρ(r) (2.59)
and with k F (r) = (3π 2 ρ(r)) 1/3 as the the local Fermi wavevector. Among the wide range of GGA functionals such as Perdew-Wang (PW91), PerdewBurke-Ernzerhof (PBE) or Becke-Lee-Yang-Parr (BLYP), the BLYP formulation is the one mostly used in this thesis, composed by the Becke's exchange [START_REF] Ad | Density-functional exchange-energy approximation with correct asymptotic behavior[END_REF] and the Lee-Yang-Parr correlation [START_REF] Perdew | Pair-distribution function and its coupling-constant average for the spin-polarized electron gas[END_REF].

Planewave basis sets and pseudopotentials

To solve explicitly the KS equation a basis set has to be selected to expand the wavefunctions. This choice is crucial for any implementation of the DFT theory and determines the type of code that will be developed. In the specific case of plane waves (PW), used in several packages, and well suited for condensed matter systems, especially when periodic boundary conditions (PBC) are used, this allows to get a set of basis functions no depending on the atomic coordinates. The effective potential V ef f depending on the positions of the ions must satisfy the following periodicity conditions:

V ef f (r + a) = V ef f (r) (2.60)
where a is an arbitrary translation of (an integer number of) the periodicity of the system. Due to the periodicity of the effective potential, the electronic wavefunctions ψ(r) are also periodic and written as:

ψ i,k (r) = e ir•k u i,k (r) (2.61)
where k is the reciprocal space vector and u i,k (r) is the arbitrary function having the periodicity of the system. The PW representation then becomes the Fourier series:

u i,k (r) = G C i,k (G)e ir•G (2.62)
consistently with the Bloch theorem [START_REF] Bloch | Quantum mechanics of electrons in crystal lattices[END_REF] and where C i,k are the Fourier coefficients depending on the G vectors of the reciprocal space. The electronic wavefunctions are then expressed as

ψ i,k (r) = G C i,k (G)e ir•(G+k) (2.63)
Although the vector set k is now discrete and finite, the numerical calculation of the electronic wave functions ψ i,k (r) is infeasible since equation (2.63) would imply an infinite sum, clearly impossible in practice. From a physical and numerical standpoint, such an infinite sum makes however little sense. On one hand the Fourier series has a uniform convergence property making the coefficients C i,k (G) decreasing systematically with increasing G, thus assuming values that can be either below the DFT accuracy or below the computer accuracy, or both. The sum can then be safely truncated at a proper cut-off value G cut . This determines a corresponding cutoff kinetic energy E cut as:

E cut = h 2 8π 2 m e |k + G cut | 2 (2.64)
This cutoff energy is system dependent and must be carefully selected to get accurate results, thus the number of PWs N P W can vary according to the specific system simulated and can be expressed by:

N P W = 1 2π 2 V E 3/2 cut (2.65)
In the representation of the wave functions, a short-scale (yet computationally demanding) part is represented by the core electrons generally inert and do not participate actively in the chemical boding. For this reason, they can be replaced by effective pseudopotentials (PPs) with a considerable reduction in the number of electrons to be explicitly considered.

The PPs [START_REF] Heine | The pseudopotential concept[END_REF] allow to describe the core-valence interaction by using explicitly only the valence states. This approximation amounts at modifying the shape of the external potential V ext (r) by replacing the strong electron-ion potential (external potential) with a softer potential, a pseudo-potential (PP) able to describe all the salient characteristics of a valence electron (Figure 2.4). These pseudo-electronic states have exactly the same potential outside a given core region as the original electrons, but have a much smoother potential inside the core. Analogously, the pseudo-wave functions are smoother and nodeless, matching the all-electron solutions outside the core region, thus making easier and tractable complex extended systems.

Figure 2.4. Schematic illustration of all-electron (solid lines) and pseudo-electron (dashed lines) potentials and their corresponding wave functions. The radius at which all-electron and pseudo-electron values match is designated as r c [START_REF] Bouzid | First-principles investigation of binary and ternary amorphous chalcogenide systems[END_REF].

In norm-conserving pseudopotentials, each eigenstate of the multi-electronic Schrödinger equation is defined by three quantum numbers (n, l, m) where n is the principal quantum number (namely the energy level), l is the secondary quantum number (angular momentum) and m is a magnetic quantum number. The wave function is written analytically as:

ψ n,l,m (r, θ, φ) = R n,l (r)Y l,m (θ, φ) (2.66)
where R n,l (r) is the radial part and Y l,m (θ, ψ) are spherical harmonics. The family of norm-conserved pseudopotentials respects the following conditions:

• Identity of pseudo (PS) and all-electron (AE) eigenvalues for a given configuration

ε AE n,l = ε P S n,l
(2.67)

• The real and pseudo wave functions are equal beyond the chosen cutoff radius r c

R AE n,l (r) = R P S n,l (r) f or r > r c (2.68)
• The integrals of the real and pseudo charge densities agree for each valence state (norm-conserving) [START_REF] Dr. Hamann | Norm-conserving pseudopotentials[END_REF][START_REF] Pickett | Pseudopotential methods in condensed matter applications[END_REF][193]

rc 0 |R AE n,l (r)| 2 r 2 dr = rc 0 |R P S n,l (r)| 2 r 2 dr (2.69)
Many types of norm-conserving PPs have been proposed such as Troullier and Martins [START_REF] Troullier | Efficient pseudopotentials for plane-wave calculations[END_REF][START_REF] Troullier | Efficient pseudopotentials for plane-wave calculations. ii. operators for fast iterative diagonalization[END_REF], Kerker [196], Haman, Schluter and Chiang [START_REF] Dr. Hamann | Norm-conserving pseudopotentials[END_REF][START_REF] Dr | Generalized norm-conserving pseudopotentials[END_REF], Goedecker, Teter and Hutter [START_REF] Goedecker | Separable dual-space gaussian pseudopotentials[END_REF] and even non-norm-conserving (ultrasoft) ones (Vanderbilt[199]). The general form of a PP is composed by a local and a non-local part:

V P P (r, r ′ ) = V P P loc (r) + V P P non-loc (r, r ′ ) = V P P loc (r) + Lmax l=1 l m=-l Y * l,m (θ, φ)∆V l (r)δ(r -r ′ )Y l,m (θ, φ) (2.70)
The local part can be further rewritten as a sum of a V core and a short-range local part ∆V loc :

V P P loc (r) = {V core (r) + ∆V loc (r)} δ(r, r ′ ) (2.71)

Born-Oppenheimer and Car-Parrinello molecular dynamics

Born-Oppenheimer molecular dynamics (BOMD)

The Born-Oppenheimer (BO) molecular dynamics allows to simulate the motion of atoms upon recalculation at each dynamical step of the electronic ground state, generally within the DFT framework. In the BO approximation, the faster electronic motion can be separated from that of the slower nuclei, in the sense that ground state electrons follow in an adiabatic way the change in nuclear configuration. Thus, the electronic structure problem reduces to solving a time-independent stationary Schrödinger-like(actually Kohn-Sham) equations, whereas nuclei are propagated as in classical mechanics with the forces acting on the nuclei computed from the total energy as computed upon the mentioned electronic structure optimization [START_REF] Marx | Ab initio molecular dynamics: basic theory and advanced methods[END_REF]. The standard BOMD protocol can be summarized as follows:

• Fix positions of nuclei; solve the KS equations self-consistently for electronic minimization.

• Compute electrostatic forces on each atom (via the Hellmann-Feynman theorem if satisfied).

• Move atoms via classical mechanics (forward integration) with a given time step and find new positions of atoms.

• Repeat the process until the end of the total simulation time.

The Lagrangian is given by:

L BO = E kin -E pot = 1 2 I M I Ṙ2 I -E pot [ρ(r), {R I })] (2.72)
provided that ρ(r) is the actual ground state. The first term is the kinetic energies of the nuclei, and the second term is the DFT total energy acting as the potential energy for the nuclei. Computationally, this can become a very slow process due to full electronic minimization required at each MD step. The BOMD scheme is summarized by the following coupled equations:

M I RI = -∇ I min Ψ 0 { Ψ 0 |H e |Ψ 0 } H e Ψ 0 = E 0 Ψ 0 (2.73)
where Ψ 0 are the ground state wave functions and M I is the mass of the I th nucleus.

Car-Parrinello molecular dynamics (CPMD)

To avoid the need for optimization of the electronic structure at each simulation step (and related full diagonalizaion of the Hamiltonian), in 1985 Car and Parrinello (CP) [START_REF] Hutter | Car-parrinello molecular dynamics[END_REF] proposed to propagate in a dynamical Lagrangean way the ground state wavefunctions by introducing a fictitious dynamics of the orbitals resulting in the extended Lagrangean:

L CP = 1 2 I M I Ṙ2 I + µ ψ * i (r) ψi (r) d 3 r
Kinetic energy

-Ψ 0 |H KS e |Ψ 0 Potential energy + Constraints orthonormality (2.74) with constraints = ij Λ ij ( ψ i |ψ j -δ ij ) (2.75)
Λ ij are Lagrange multipliers associated with the orthonormalization constraints of the wavefunctions:

ψ i |ψ j = ψ * i (r)ψ j (r) d 3 r = δ ij (2.76)
In the special case µ = 0 the CPMD reduces to the BOMD. In equation (2.74), the first and second terms are the kinetic energies (real) of nuclei and the fictitious one of the orbitals, respectively. The fictitious kinetic energy of the orbitals serves as a measure to control the deviations of the CPMD from the exact BO surface [START_REF] Jg | Computational materials science an introduction[END_REF]. The corresponding equations of motion are obtained from the related Euler-Lagrange equations:

M I RI = - ∂E({ψ i }, {R I }) ∂R I (2.77) µ ψi = - δE({ψ i }, {R I }) δψ * i + j Λ ij ψ j (2.78)
where µ is the fictitious electron mass controlling the adiabaticity of the system. The total conserved CPMD energy is obtained when δ ij = 0 and given as:

E cons = 1 2 I M I Ṙ2 I + µ ψ * i (r) ψi (r) d 3 r + E({ψ i }, {R I }) (2.79)
while the orthonormality constraints on the orbitals ψ i do not introduce any energy dissipation and are computed iteratively along the dynamics.

Thermodynamic control: Nosé-Hoover and Blöchl-Parrinello thermostats

In CPMD, analogously to what is done in CMD, the temperature control is ensured by the use of thermostats. [START_REF] Wg | Canonical dynamics: Equilibrium phase-space distributions[END_REF][START_REF] Gj | Nosé-hoover chains: The canonical ensemble via continuous dynamics[END_REF] The idea of Nosé and Hoover [START_REF] Wg | Canonical dynamics: Equilibrium phase-space distributions[END_REF][START_REF] Gj | Nosé-hoover chains: The canonical ensemble via continuous dynamics[END_REF] for the classical particles (ions) was extended to the electronic degrees of freedom by Blöchl and Parrinello [START_REF] Blöchl | Adiabaticity in first-principles molecular dynamics[END_REF] resulting in the dynamics expressed by:

M I RI = - ∂E({ψ i }, {R I }) ∂R I -M I ṘI ṡ (2.80) µ ψi = - δE({ψ i }, {R I }) δψ * i + j Λ ij ψ j -µ ψi ξ (2.81)
The frictions terms (M I ṘI ṡ and µ ψi ξ) are governed by the two variables s and ξ of the thermostats having their own equations of motion

Q R s = 2 I 1 2 M I Ṙ2 I - 1 2 gk b T (2.82
)

Q e ξ = 2 i µ ψ i |ψ i -E kin e (2.83)
being (1/2)gk b T the average kinetic energy of the ionic sub-system, g = (3N + 1) the number of degrees of freedom for the atomic motion in a system with N atoms, k b the Boltzmann constant and T the physical temperature of the system. E kin e is the average kinetic energy of the fictitious electronic sub-system. The masses Q R and Q e determine the time scale of the thermal fluctuations of the thermostats. In the special case in which all the atoms have the same mass M , the fictitious kinetic energy is simplified as:

E kin ad = 2k b T µ M i ψ i | - 1 2 ∇ 2 |ψ i (2.84)
This E kin ad value is generally used to estimate the numerical value of fictitious kinetic energy E kin e . However, the recommended value is about twice E kin ad [START_REF] Blöchl | Adiabaticity in first-principles molecular dynamics[END_REF] as shown in Figure 2.5. The frequencies, ω T R , of the thermal fluctuations associated to the thermostat controlling the ionic degrees of freedom are given by:

ω T R = 2gk b T Q R (2.85)
and for the frequencies ω T e of the thermal fluctuations associated to the electronic thermostat by:

ω T e = 4E kin e Q e (2.86)

Adiabaticity control

One of the characteristics of the CPMD scheme is that it accounts for the properties of the quantum spectrum of the electronic eigenvalues. By analyzing the harmonic frequency spectrum of the orbitals, we obtain (2.87):

ω ij = 2(ǫ i -ǫ j ) µ (2.87)
Where ǫ i and ǫ j are the eigenvalues of occupied and unoccupied orbitals, respectively. An analytical approximation for the lowest possible electronic frequency can be obtained as:

ω min e ∝ E gap µ 1/2
(2.88)

being E gap the gap energy between the highest occupied orbital and the lowest unoccupied one, and µ the fictitious electronic mass. In this equation we can see that the frequency increases as the square root of the energy gap E gap . Hence, the key parameter controlling the adiabatic separation between nuclear and electronic variables is the fictitious mass µ [START_REF] Marx | Ab initio molecular dynamics: basic theory and implementation[END_REF]. However, a decrease of µ not only causes the shift of the electronic spectrum upwards, but also allows to extend the frequency spectrum according to equation (2.87). This stretching of the frequency spectrum leads to a subsequent increase in the maximum frequency given by

ω max e ∝ E cut µ 1/2 (2.89)
where E cut is the largest kinetic energy in an expansion of the wavefunction in terms of a plane wave basis set. The molecular dynamics time step ∆t max is inversely proportional to the highest frequency of the system, which is ω max e and is given as:

∆t max ∝ µ E cut 1/2
(2.90)

Impact of dispersion forces on the structural properties of liquids

The role of dispersion forces in modeling surfaces, interfaces and liquids has been shown to be crucial for the proper quantitative description of the behaviour of these systems. Disordered phases (organic and inorganic) can play a benchmark role in assessing the sensitivity of structural properties to the inclusion of dispersion forces (vdW) within DFT-based FPMD methods. In this section, two different theoretical frameworks for the inclusion of vdW interactions within DFT-based FPMD methods are briefly discussed and compared for the case study of a liquid chalcogenide system: i. the fully first principles theoretical approach proposed Silvestrelli et al. [START_REF] Silvestrelli | van der waals interactions in dft using wannier functions without empirical parameters[END_REF][START_REF] Pl | Van der waals interactions in dft made easy by wannier functions[END_REF] and Ikeda et al. [START_REF] Ikeda | Role of van der waals corrections in first principles simulations of alkali metal ions in aqueous solutions[END_REF] (making use of maximally localized Wannier functions (MLWF hereafter) and ii. the (semi)empirical correction proposed by Grimme et al. [START_REF] Grimme | Performance of dispersion-corrected density functional theory for the interactions in ionic liquids[END_REF] within the so-called D2 and D3 schemes. The former is based on the update of the dispersion coefficients according to the time evolution of the electronic structure, while the latter features vdW contributions skillfully extracted from a large set of reference calculations but not sensitive to changes in the electronic structure during the dynamics and, as such, not updated along a dynamical simulation. These two approaches consider uniquely the case of dispersion effects treated as a quantity separable from the main body of the Kohn-Sham Hamiltonian and linearly added to it, adopting the expression C 6 /r 6 ij , where r ij is the interatomic distance between a pair of atoms i and j. By following the D2 formulation of Grimme et al. [START_REF] Grimme | Semiempirical gga-type density functional constructed with a long-range dispersion correction[END_REF], the dispersion contribution to be added to the total energy takes the form:

E D2 vdW = -s 6 N i,j>i C ij 6 r 6 ij f damp (r ij ) (2.91)
where C ij 6 is the coefficient for the atom pair i and j parameterized by quantum chemical calculations on a wide set of training systems, s 6 is a scaling factor depending on the type of exchange-correlation functional used in these first-principle-based calculations, and f damp is a suitable damping coefficient. D2 is a largely tested methodology, already applied for the modeling of disordered chalcogenides such glassy GeTe 4 [START_REF] Bouzid | Role of the van der waals interactions and impact of the exchange-correlation functional in determining the structure of glassy gete 4[END_REF] and liquid GeSe 2 [START_REF] Lampin | Impact of dispersion forces on the atomic structure of a prototypical network-forming disordered system: The case of liquid gese 2[END_REF] and nanoporous chalcogenides gels [START_REF] Ziyad | The structure and dipolar properties of co2 adsorbed in a porous glassy chalcogel: Insights from first-principles molecular dynamics[END_REF][START_REF] Essomba | Assessing the versatility of molecular modelling as a strategy for predicting gas adsorption properties of chalcogels[END_REF]. Along the lines pioneered with D2, Grimme and coworkers later proposed a refinement, named D3, including also an additional energy correction complementing equation (2.91) and consisting of a threebody interaction as derived from the third order perturbation theory [START_REF] Grimme | A consistent and accurate ab initio parametrization of density functional dispersion correction (dft-d) for the 94 elements h-pu[END_REF]. This extended vdW correction has been indicated as affected by "less empiricism" and used in connection with several exchange-correlation functionals, but also features improvements in terms of cost and robustness. These ideas are claimed to provide unprecedented performances for the treatment of dispersion forces and have been also tested on a few ILs systems [START_REF] Grimme | Performance of dispersion-corrected density functional theory for the interactions in ionic liquids[END_REF]. However, a common feature of both the D2 and the D3 schemes is that they are not intended to account for on-the-fly evolutions of the electronic structure with time. The approach by Silvestrelli et al. [START_REF] Silvestrelli | van der waals interactions in dft using wannier functions without empirical parameters[END_REF][START_REF] Pl | Van der waals interactions in dft made easy by wannier functions[END_REF] and Ikeda et al. [START_REF] Ikeda | Role of van der waals corrections in first principles simulations of alkali metal ions in aqueous solutions[END_REF] provides this C ij 6 updating by resorting to the use of MLWF. Accordingly, the coefficient C 6 is written as (with indices n, l at the place of i and j):

C nl 6 = 3 32π 3 |r|<rc d 3 r |r ′ |<rc d 3 r ′ w n (r)w l (r ′ ) w n (r) + w l (r ′ ) (2.92)
The advantage of this scheme is its dependence on the localized character of w n (r) since, in principle, the expression based on the local electronic density ρ n (r) is written as

C nl 6 = 3 32π 3 |r|<rc d 3 r |r ′ |<rc d 3 r ′ ρ n (r)ρ l (r ′ ) ρ n (r) + ρ l (r ′ ) (2.93)
but it is affected by the lack of a rigorous local definition for the total electronic density, particularly when non-localized basis sets, such as PWs, do not allow for an unbiased partitioning of the density. An additional simplification in handling the MLWF is that the electronic structure information can be reduced to four numbers, the position (x,y,z) of the MLWF center of mass, and its spread. At note that this implementation based on localized orbitals allows following small electronic effects due to changes in the atomic configurations, readily affecting the values taken by the dispersion coefficients.

In the following, the application and comparison of these different schemes for taking into account the dispersion forces are compared for the case of liquid GeSe 4 . The total neutron structure factor S tot (k) of liquid GeSe 4 , obtained by linear combination of the partial structure factor by accounting for the coherent neutron scattering lengths and the concentration of the two species is shown in Figure 2.6. S tot (k) resulting from the various dispersion schemes agrees well with the experimental data of Maruyama and coworkers [START_REF] Maruyama | Structural investigation of liquid ge-chalcogen mixtures[END_REF] provided one refers to the range of values of k larger than 2 Å-1 in reciprocal space. However, two notable exceptions can be remarked at lower k. The first one is the strong disagreement recorded around 1 Å-1 when using D3 (for both D3(A) and D3(B) trajectories) as if the intermediate range order typically manifesting itself via the presence of a peak around 1 Å-1 were absent. The second, more surprising, is the appearance of an outstanding feature detectable at very low values of k in the D3 cases. A direct visual inspection of the structures obtained for the different dispersion schemes reveals some macroscopic effects consistent with the peculiarity exhibited by D3. Figure 2.6 is insightful in this respect. While the snapshots extracted from No-vdW, D2, and MLWF trajectories exhibit a fully homogeneous arrangement in space, the snapshot relative to D3 is indicative of an inner phase separation. This phenomenon, in addition to being difficult to grasp on physical grounds, bears no resemblance to the phase separation The Ge atoms are in dark green, the Se atoms in brown, and the unit cell is represented with a black line [START_REF] Martin | Impact of dispersion force schemes on liquid systems: Comparing efficiency and drawbacks for well-targeted test cases[END_REF].

between GeSe 4 and Se n subnetworks hypothesized in the past and is not substantiated by any FPMD calculations. The puzzling evolution of the trajectory generated via D3(A) is confirmed when starting from a different initial configuration (final configuration of the D2 trajectory) and performing an additional 30 ps trajectory (D3(B) results). These results call for particular attention on the role of dispersion forces and their implementation within DFT.

Chapter 3

Modelling of ionic liquids in the bulk phase

Context and objectives

Ionic liquids (ILs) bear great promise in green chemistry, environmental science and innovative electronic devices due to their unique properties, such as a tailorable structure, chemical variety, and their mitigated environmental impact [START_REF] Chenlu | Topological engineering of two-dimensional ionic liquid islands for high structural stability and co 2 adsorption selectivity[END_REF][START_REF] Vioux | From ionogels to biredox ionic liquids: some emerging opportunities for electrochemical energy storage and conversion devices[END_REF][START_REF] Welton | Ionic liquids in catalysis[END_REF][START_REF] Jt | Superconducting dome in a gate-tuned band insulator[END_REF]. ILs consist of a mixture of cations and anions and because of their rich interactions and microstructures, they can act as solvents, active moieties in composite materials, and catalysts. This wide spectrum of applications makes ILs useful in the forefront fields of next-generation electronic devices and green chemistry. In particular, alkylimidazolium ionic liquids represent a special family of ILs based on the 1-alkyl-3-methylimidazolium cation in combination with different anions. These ILs are among the most studied systems because of their remarkable structural and electronic properties and the simplicity of synthetic methods for their production. These ILs features strong interactions ranging from bare electrostatic forces to hydrogen bonds (H-bonds) with complex nano-and micro-structural ordering [START_REF] Canongia | Nanostructural organization in ionic liquids[END_REF][START_REF] Oldamur | On the origin of ionicity in ionic liquids. ion pairing versus charge transfer[END_REF][START_REF] Yanlei | Insights into ionic liquids: From z-bonds to quasi-liquids[END_REF]. Specifically, these ILs can realize H-bonds of different nature with interesting networking properties [START_REF] Kun | Understanding the hydrogen bonds in ionic liquids and their roles in properties and reactions[END_REF][START_REF] Bogdan | Crossover from hydrogen to chemical bonding[END_REF][START_REF] Ap | Hydrogen bonding in ionic liquids[END_REF]. As a consequence, this rich variety of interactions promotes the formation of complex nanostructures (including ionic pairs, aggregates, or even ionic clusters) which translates into interesting and unprecedented effects on the performance of ILs in a wealth of applications [START_REF] Yanlei | Insights into ionic liquids: From z-bonds to quasi-liquids[END_REF]. The complex interplay between H-bonds and dipole properties of ILs is directly dependent on the features of the IL used. Beside their tunability in terms of chemical compositions and physicochemical properties, these features have shown to affect the ILs performance whether used as bulk solvent phases for catalysis or as components in gating media for optoelectronic transistors [START_REF] Js | Properties of ionic liquid solvents for catalysis[END_REF][START_REF] Ph | Synergistic hole-doping on ultrathin mote2 for highly stable unipolar field-effect transistor[END_REF][START_REF] Michio | Exceptionally high electric double layer capacitances of oligomeric ionic liquids[END_REF][START_REF] Yogesh | High-performance organic field-effect transistors gated by imidazolium-based ionic liquids[END_REF][START_REF] Young | Dipole-assisted carrier transport in bis (trifluoromethane) sulfonamide-treated o-res 2 field-effect transistor[END_REF][START_REF] Velpula | Graphene meets ionic liquids: Fermi level engineering via electrostatic forces[END_REF]. Atomistic modeling, especially in terms of CMD simulations, has been provided a usefulinsight into the structural properties of ILs in the last decade, improving our understating of their behaviour and complementing the data provided by experiments [START_REF] Canongia | Modeling ionic liquids using a systematic all-atom force field[END_REF][START_REF] Coasne | Loading-controlled stiffening in nanoconfined ionic liquids[END_REF][START_REF] Ori | Ionic liquid confined in silica nanopores: molecular dynamics in the isobaric-isothermal ensemble[END_REF][START_REF] Ori | Structure and dynamics of ionic liquids confined in amorphous porous chalcogenides[END_REF]. Several examples of CMD studies oriented to the comprehension of H-bonds features in ILs are available nowadays [START_REF] Mendez-Morales | Molecular dynamics simulations of the structural and thermodynamic properties of imidazolium-based ionic liquid mixtures[END_REF][START_REF] Agrawal | Structures of binary mixtures of ionic liquid 1-butyl-3-methylimidazolium bis (trifluoromethylsulfonyl) imide with primary alcohols: the role of hydrogen-bonding[END_REF][START_REF] Aa | Effects of water concentration on the structural and diffusion properties of imidazolium-based ionic liquid-water mixtures[END_REF][START_REF] Bhattacharjee | The wetting behavior of aqueous imidazolium based ionic liquids: a molecular dynamics study[END_REF][START_REF] Andreussi | Transport properties of room-temperature ionic liquids from classical molecular dynamics[END_REF]. First-principles molecular dynamics (FPMD) has also been employed to obtain further insight into the structure, H-bonds and electronic properties of a few ILs [START_REF] Bühl | Ab initio molecular dynamics of liquid 1,3-dimethylimidazolium chloride[END_REF][START_REF] Grimme | Performance of dispersion-corrected density functional theory for the interactions in ionic liquids[END_REF][START_REF] Del Popolo | Ab initio molecular dynamics simulation of a room temperature ionic liquid[END_REF][START_REF] Hunt | Characterising the electronic structure of ionic liquids: an examination of the 1-butyl-3-methylimidazolium chloride ion pair[END_REF][START_REF] Cremer | Towards a molecular understanding of cation-anion interactions-probing the electronic structure of imidazolium ionic liquids by nmr spectroscopy, x-ray photoelectron spectroscopy and theoretical calculations[END_REF][START_REF] Kirchner | Ion pairing in ionic liquids[END_REF][START_REF] Kempter | The role of hydrogen atoms in interactions involving imidazolium-based ionic liquids[END_REF][START_REF] Yong-Lei | Hydrogen bonding versus -stacking interactions in imidazolium-oxalatoborate ionic liquid[END_REF]. However these works have often been focused on individual or small cluster of IL pairs, or small series of ILs, mainly because of the computational workload inherent to an explicit treatment of the electronic structure. The purpose of the present chapter is to explore, mainly via CMD methods, structural, hydrogen bonding and dipole moment properties of a wide series of alkyl-imidazolium-based IL, some of which is currently investigated in optoelectronic devices and field-effect transistors. The specific family of ILs on which we focus is constituted by imidazolium-based cations with differ-ent lengths of the alkyl chain and the mostly experimentally exploited inorganic anions, as sketched in Figure 3.1. As mentioned, this computational study makes use mainly of CMD to simulate, within an identical and consistent approach, a large numbers of IL systems and to extract their structural properties and related H-bond network features. This CMD investigation is an essential step to tackle in a computational affordable way a thorough exploration of such a relatively wide variety of ILs. For a selected number of IL systems, first-principles molecular dynamics (FPMD) simulations are used to enrich the information by the inclusion of the electronic structure and to refine our understanding of their interaction without relying on ad-hoc parameterized FFs. This second stage of this doctoral project allows for a better description of the chemical bonding and the electrostatics regulating the H-bond network formation. Furthermore, the explicit inclusion of the electronic degrees of freedom allows to get a deeper insight into the dipolar interactions for a direct comparison with the CMD results. In a sense this is also an indirect assessment of the FFs used in classical approaches.

Computational details

Classical molecular dynamics

The ILs on which this work is focused are made of imidazolium-based cations, specifically we have considered 1-alkyl-3-methylimidazolium cations combined with different inorganic anions having a general formula C x MIM + -A -, where x = 1, 2, 4 and 6 (also termed DMIM, EMIM, BMIM, and HMIM, respectively) and the anion A -is chloride Cl, or tetrafluoro-borate BF 4 , or hexafluoro-phosphate PF 6 or bis(trifluoromethane)-sulfonimide TFSI. As outlined in Chapter 1 and 2, on a general ground, a non-polarizable FF is often retained as sufficient to reproduce with appreciable accuracy the structure of the targeted ILs [START_REF] Andreussi | Transport properties of room-temperature ionic liquids from classical molecular dynamics[END_REF][START_REF] Gehrke | Structure and lifetimes in ionic liquids and their mixtures[END_REF][START_REF] Dommert | Force fields for studying the structure and dynamics of ionic liquids: A critical review of recent developments[END_REF]. The specific FF mainly used in the present work is the CL&P proposed and developed by Canongia Lopez and coworkers [START_REF] Canongia | Modeling ionic liquids using a systematic all-atom force field[END_REF] with the parametrization of Koddermann et al. [246] to better account for the dynamical properties of the bulk ILs simulated. This FF is fitted on ab initio calculations and allows to simulate a wide variety of ILs chemical structures including the specific ones targeted in this project. For the specific case of EMIM-BF 4 , also the polarisable FF CL&Pol was tested to assess the reliability of the results obtained in the absence of polarization. The development of the CL&Pol FF parameters for this IL was obtained by following the procedure reported by the developers [START_REF] Goloviznina | Extension of the clpol polarizable force field to electrolytes, protic ionic liquids and deep eutectic solvents[END_REF][START_REF] Goloviznina | Transferable, polarizable force field for ionic liquids[END_REF]163,[START_REF] Goloviznina | Development and application of molecular interaction models with explicit polarisation for ionic liquids and eutectic solvents[END_REF]. Such a procedure consisted in the parametrization of scaled LJ parameters and point charges for the atoms and in the refining of suitable parameters for the Drude particles allowing to describe the polarisability of this IL. Periodic boundary conditions are applied to the bulk IL simulation cell and the dispersive interactions are neglected beyond a cutoff of 15 Å. A standard Ewald summation, described in the second chapter of this thesis, is used in the calculation of the electrostatic interaction to avoid finite size effects. In the calculations conducted, we verified that the Ewald parameters selected provide an accuracy of 10 -5 atomic units (au). For each IL, we constructed two simulation cells of different sizes, containing, respectively 20 and 200 IL pairs. The reason for these different sizes of the same IL system will be clarified in the next paragraph and in the ongoing discussion. It is anyhow clear that the larger systems accounting for 200 IL pairs are sufficiently large to ensure a statistically significant average of the structural properties. The initial configuration has been realized as a random accommodation of IL pairs inside a cubic box by means of the Packmol package [START_REF] Martínez | Software news and update packmol: A package for building initial configurations for molecular dynamics simulations[END_REF], keeping into account the steric encumbrance of the molecular building blocks, and the initial volume of this box was selected based on the experimental density of the IL. The bulk equilibrated model was then produced following a standard melt-quench procedure heating the systems up to ∼ 1000 K for a simulation time of 2.5 ns, and then by cooling these same systems to room temperature (∼ 300 K) in a sequence of steps in which a temperature decrease of 100 K was imposed and on which the system was allowed to re-equilibrate for 1 ns on each temperature plateau in the canonical ensemble (NVT). The initial 1000 K NVT stage lasted for 2.5 ns to ensure the lose of memory of the initial configuration constructed as reported above. After the NVT thermal cycle, a simulation stage in the isobaricisothermal ensemble (NPT) at T = 300 K and P = 1 atm was performed for 1.5 ns to get well equilibrated and stress-free models, and to check the density of the simulated systems with respect to the assumed experimental values. These preparatory simulations provided the actual systems simulated for about 10 ns in the NVT ensemble. This corresponds to the production stage in which statistics was collected and used to characterize the various ILs. More precisely, the data analyzed and presented in the following paragraphs refer to the final 5 ns in the NVT ensemble where trajectories were sampled each 1000 steps and the system has been assessed as well equilibrated and free from any residual stress or imposed configuration reminiscent of the original construction. Temperature and pressure were controlled with a Nosé-Hoover thermostat [START_REF] Nosé | A molecular dynamics method for simulations in the canonical ensemble[END_REF][START_REF] Nosé | A unified formulation of the constant temperature molecular dynamics methods[END_REF][START_REF] Wg | Canonical dynamics: Equilibrium phase-space distributions[END_REF] and a Hoover barostat [START_REF] Wg | Canonical dynamics: Equilibrium phase-space distributions[END_REF], respectively. The simulations using the non-polarisable CL&P FF were performed with the DLPOLY package [START_REF] Smith | The dl poly molecular dynamics package[END_REF], and a standard Verlet Leapfrog algorithm is used for the numerical integration of the equations of motion with a time step of 1.0 fs. Conversely, the simulations using the polarisable CL&Pol FF were performed with the LAMMPS package [START_REF] Ap | LAMMPS -a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales[END_REF].

First-principles molecular dynamics

The explicit inclusion of the electronic structure is ensured by FPMD simulations done on a selected number of IL models. To this aim, to make the system tractable within this type of approach and to mitigate the computational workload, the smaller IL systems mentioned in the former paragraph were chosen. The dynamical FPMD simulations were performed by resorting to the Car-Parrinello method [START_REF] Car | Unified approach for molecular dynamics and densityfunctional theory[END_REF] as implemented in the developer's version of the CPMD[251] code. The valence-core interaction was described by norm-conserving numerical pseudopotential of the Troullier-Martins [START_REF] Troullier | Efficient pseudopotentials for plane-wave calculations[END_REF][START_REF] Troullier | Efficient pseudopotentials for plane-wave calculations. ii. operators for fast iterative diagonalization[END_REF] type, and the exchange and correlation interactions are expressed according to the functional forms proposed by Becke [START_REF] Ad | Density-functional exchange-energy approximation with correct asymptotic behavior[END_REF] and Lee, Yang, and Parr [START_REF] Lee | Development of the colle-salvetti correlation-energy formula into a functional of the electron density[END_REF] (BLYP), respectively. Long-range van der Waals dispersion interactions, not included neither in the BLYP nor in any GGA functional, are added according to the maximally localized Wannier functions [START_REF] Marzari | Maximally localized wannier functions: Theory and applications[END_REF] and centers approach [START_REF] Pl | Van der waals interactions in dft made easy by wannier functions[END_REF] boosted by a propagation scheme [START_REF] Ikeda | Role of van der waals corrections in first principles simulations of alkali metal ions in aqueous solutions[END_REF] to reduce the computational burden. This vdW implementation is the one summarized in Section 2.2.7. Valence electrons are treated explicitly and their wavefunctions are expanded on a PW basis set with a cut-off of 100 Ry. This value was tuned upon monitoring of the convergence of the stress tensor as a function of the number of PWs in the preparatory stage before moving to the actual production run. We remark that this relatively large cut-off is crucial for an accurate description of the electron density around several chemical species (O, N, S, F, Cl) composing the ILs. Canonical NVT simulations were performed by controlling the temperature via a Nosé-Hover [START_REF] Nosé | A molecular dynamics method for simulations in the canonical ensemble[END_REF][START_REF] Nosé | A unified formulation of the constant temperature molecular dynamics methods[END_REF][START_REF] Wg | Canonical dynamics: Equilibrium phase-space distributions[END_REF] thermostat chain [START_REF] Gj | Nosé-hoover chains: The canonical ensemble via continuous dynamics[END_REF]. For the numerical integration of the Car-Parrinello equations of motion a fictitious electron mass of 400.0 au and a time step of 3.0 au (0.0726 fs) ensured good numerical control of the constants of motion. The ILs simulated with this first principles approach are

[DMIM][Cl], [DMIM][BF 4 ], [DMIM[]PF 6 ], [DMIM][TFSI], and [EMIM][BF 4 ]
. Each simulation started from its own specific liquid phase model amounting to 20 pairs, which translates into a number of atoms between 340 and 620, depending on the specific composition. A pre-equilibration stage was done within our CMD scheme, following the thermal cycle explained above. This allows to generate the initial configurations to be used for the FPMD simulations at 300 K. The stress tensor analysis done during NVT dynamical simulations at different cell sizes was used to assess the actual size of the simulation cell to be used for the production and the resulting theoretical density of the ILs turned out in appreciable agreement with the experimental value. The system displays just a minimal residual pressure, compatible with the DFT accuracy. In the following figure, as an example of the tests performed to optimize each IL model density, we show the trend of the computed stress tensor as a function of PW basis set cut-off for [DMIM][Cl]. At a PW cutoff of 100 Ry and with a volume expansion of ∼5% with respect to the experimental value, the stress tensor converges to a value of about ∼0 GPa.

Trajectories were accumulated in the NVT production run at 300 K for 39 ps for ] is shown. The constant average value of the KS potential energy and the energetically well-separated temperatures of ions and electrons along the simulation time prove the numerical stability of the simulation and the reliability of the computational protocol used [START_REF] Blöchl | Adiabaticity in first-principles molecular dynamics[END_REF]. In particular, for [DMIM][BF 4 ], but also for all the other IL systems investigated, the ratio between the ionic temperature and the fictitious one of the electrons is about ∼15, well in the range of values allowed for a correct behavior of the Car-Parrinello scheme.

Structural properties

The results obtained with both the CMD and FPMD for the IL systems considered are summarized in Table 3.1, along with the data concerning number of atoms and models sizes. The references in square parentheses indicated in the table are the available data published, where the corresponding experimental densities can be found. In terms of theoretical density for the models produced by CMD, as a combination of the simulations done in the NPT and NVT ensembles, the majority of the models have shown small deviations (<3%) from the experimental value (with the exception of [DMIM][Cl] (∼5%)). Concerning the FPMD results, in the second part of this same table, also small deviations (∼3%) affect the [DMIM][BF 4 ,TFSI] and [EMIM][BF 4 ] systems, whereas for [DMIM][Cl,PF 6 ] deviations up to ∼9% were obtained.

In terms of overall structural properties, the first quantification was done in terms of cation-anion pair correlation functions (PCFs) of the center of mass (COM) cation-anion g ij (r) using the standard definition :

g(r) = V 4πN r 2 i i =j δ(r -r ij ) (3.1)
where N is the number of molecules contained in the simulation cell of volume V and r ij is the distance between the COM of each cation and anion constituting the IL system. The result of these COM PCFs for the different ILs within our CMD approach are reported in Figures 3.3 and 3.4.

In particular, in Figure 3.3 we show the comparison between the two models consisting of 20 and 200 IL ion pairs for the DMIM cation with the different anions considered, along with the EMIM-BF 4 system, to assess possible size effects for this particular structural property. As it can be noted, no size effects could be found analysing the first coordination Table 3.1. Box side (L, Å), number of atoms, density (d, g.cm -3 ) and the minimum of the g ij (r) after the first peak used as cutoff for the calculation of n c of the IL models simulated by MD and FPMD. Both the starting and final values of density used in the MD simulations are reported. r c is the distance value used as cutoff for the integration of the first g ij (r) peak for calculating the coordination number n c . CMD, [START_REF] Marx | Ab initio molecular dynamics: basic theory and advanced methods[END_REF] shell the ions COM. Figure 3.4 reports the COM PCF for all the IL systems simulated by CMD. These results give an insightful global picture of the structure of these compounds as a function of both the cation alkyl chain and the type of inorganic anion. For the ILs in which Cl -is the anion, in the first coordination shell of the cation-anion pairs, the shortest alkyl chain (DMIM) shows a broad peak centered at ∼4.5 Å , followed by a smaller peak at ∼5.8 Å . Then, by increasing the length of the alkyl chain from one to six carbon atoms, the first peak is affected by an increasing broadening until the two adjacent peaks nearly merge into a unique feature for the longer HMIM, as an effect of the larger steric volume characterizing this cation. For the BF - 4 series, a similar trend is observed but overall the intensity of the peaks is reduced with respect to the PCFs of the former ILs carrying as an anion Cl -, clearly smaller than BF - 4 . The PCFs of DMIM-and EMIM-BF 4 show a double peak with the two maxima centered at ∼4.5 Å and ∼5.2 Å , a clear indication of two preferred local coordination of the anion with respect of the orientation of the alkyl-imidazolium cation. In the case of BMIM and HMIM, the double peak is replaced by a broad first peak centered at ∼4.4 Å and a lower second peak at ∼6.2 Å. Finally, the PF - 6 and TFSI-systems show analogous trends with a further reduction in the peaks intensity if compared with the system composed by the Cl -and BF - 4 anions. Namely, a double peak present in the case of the short alkyl chains becomes a unique broad peak upon increasing of the chain length, as seen in the former cases.

In terms of first peak positions, good agreement was found with former results published in the literature, such as Liu et al. [START_REF] Liu | A refined force field for molecular simulation of imidazolium-based ionic liquids[END_REF] and Morrow et al. [START_REF] Morrow | Molecular dynamics study of the ionic liquid 1-n-butyl-3-methylimidazolium hexafluorophosphate[END_REF] for the BMIM-PF 6 , EMIM-BF 4 , BMIM-BF 4 and DMIM-Cl pairs.

An additional insightful parameter useful to characterize the structure of ILs is the coordination number n c . This can be obtained by direct integration of the PCF, g ij (r) in the usual way, namely

n c = 4πρ 0 rc 0 r 2 g(r)dr (3.2)
where ρ 0 and r c are the atomic density of each IL pairs and the position of first absolute minimum between two maxima of the PCF, respectively. The values of n c obtained from both classical MD and FPMD are reported in Table 3.2. Using the cut-offs indicated in Table 3.1, the total average coordination numbers calculated for the cation-anion in the first coordination shell lie in the range between 6 and 10. The size and steric hindrance of the anion, which increases along the anion series as Cl -<BF - 4 <PF - 6 <TFSI -, are major contribution to the average value of the coordination number, which turns out to be n c = 6.69, 7.46, 7.92 and 10.11 for DMIM-Cl, BF 4 , PF 6 and TFSI, in that order. The cation-anion coordination number for Cl-based IL shows no particular trend while BF 4 , PF 6 and TFSI show a decrease in n c with increasing alkyl chain length. A figure of merit worthy of note is the fact that the values of n c estimated from the classical MD show a remarkable agreement with those obtained by FPMD Figure 3.5. Overall, the results obtained in terms of PCF of COMs and cation-anion coordination numbers provide a first reliable validation of the FF used in this work for the different IL models at 300 K.

Hydrogen bonding signatures

The H-bond network characterizing the different ILs targeted in this work was inspected to get a clearer picture of their nature and influence on the structure of these systems. A fundamental driving force determining the molecular structuring of IL is the electrostatic interaction occurring between pairs of ILs, responsible for structuring the H-bond network of the system and, therefore, the cohesive properties of the cation and anion moieties. This general observation is supported by the fact that along the evolution observed in the sampled trajectories, the distributions of the anions (Cl -, BF - 4 , PF - 6 and TFSI -) around the cations are mainly governed by electrostatic forces and related induced molecular of the imidazolium ring, in the g H 2 F (r) distribution is located at d H2-F =2.05 Å. This is a clear indication of a strong preference of F atoms to bind to this site. Such a short distance is also characteristic of a tight H-bond [START_REF] Hunt | Quantum chemical modeling of hydrogen bonding in ionic liquids[END_REF][START_REF] Izgorodina | Nature of hydrogen bonding in charged hydrogen-bonded complexes and imidazolium-based ionic liquids[END_REF]. The other H atoms of the cation are characterized by corresponding PCFs peaks at larger distances, namely 2.5-2.7 Å for EMIM-BF 4 , with reduced intensities compared to g H 2 -F (r). These are still active Hbonds, concurring to the overall structure of IL [START_REF] Izgorodina | Nature of hydrogen bonding in charged hydrogen-bonded complexes and imidazolium-based ionic liquids[END_REF], but being somehow weaker, they can be dynamically broken and reformed during the dynamical evolution of the liquid, as observed in the case of EMIM-TFSI [START_REF] Ishisone | Structural, dynamical, and electronic properties of the ionic liquid 1-ethyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide[END_REF]. Analogous results for the other IL targeted in this study are reported in Figure 3.7. However, with respect to EMIM-BF 4 , some differences can observed. In particular, EMIM-Cl PCF shows that together with H 2 atoms also the other H atoms of the imidazolium ring as well as the H atoms of the first C atoms connect to the ring (C 6 and C 7 ) are responsible for clear sharp peaks. These H-bond interactions are found at larger distance with respect to the case of F in EMIM-BF 4 , being at 2.65 Å (H 2 ,H 4,5 ) and 2.85 Å (H 6,7 ), respectively. In the case of EMIM-PF 6 , the interaction between the F atoms and the H 2 atoms again turns out to be the strongest compared to the other H atoms of the cation (2.25 Å versus 2.75-2.95 Å), however the intensities of the peaks of the different H are much less pronounced compared to the case of EMIM-BF 4 and EMIM-Cl. In the case of EMIM-TFSI, the strongest H-bonding interaction occur between the H atoms of the cation and the O atoms of the sulfonyl groups of TFSI -with an analogous behavior observed in the case of EMIM-Cl. In this case, the first H 2 -O PCF peak is found to be the most intense, but all H x -O participate to several extents to the H-bonding interaction at similar distances in the range 2.66-2.70 This effect of alkyl chain length is summarized in Figure 3.8. One feature that stands out is a structuring effect of the H-bond favored by H 2 with the increase in the number of C atoms composing the alkyl chain as mentioned previously, from DMIM + to longer HMIM + . This effect of alkyl chain length is present in the four ILs differing by the anionic moiety. The main difference is a decrease in the intensity of the first peak according to the Cl > BF 4 > PF 6 > TFSI trend. The comparison with the FPMD results for the DMIM and EMIM systems Figure 3.9 confirms the overall trend evidenced by CMD. However, a reduced intensity arises for the first peak for H 2 -Cl and H 2 -O in the cases of DMIM-Cl and DMIM-TFSI with slightly shorter H 2 -Cl and H 2 -O distances and slightly longer H 2 -F distances for DMIM/EMIM-BF 4 and DMIMM-PF 6 , respectively. These data are summarized in Table 3.3. Being the EMIM-BF 4 one of the most studied IL because of its use in many applications (e. g. EDLT-FET) and being also the main focus of Chapter 4 and 5, we selected this IL to make a further comparison with CMD employing the polarisable FF CL&Pol. To this aim, we present in Figure 3.10 the g H 2 -F (r) obtained by FPMD and the CMD ones obtained both with nonpolarisable (CL&P) and with polarisable (CL&Pol) FF. For both CL&P and CL&Pol FFs, no size effects were found comparing models of 20 to 200 ions pairs. Employing a FF that takes explicitly into account the polarizability of the ions, a longer H 2 -F distance is found with respect to the nonpolarisable FF (2.35-2.38 Å versus 2.05 Å ) and with a slightly lower intensity (i.e. less structured) first g H 2 -F (r) peak, which is in better agreement with the value obtained by FPMD (2.28 Å ). Moreover, the position of the second peak of the g H 2 -F (r) is also in better agreement with the outcome of FPMD simulations (4.2 Å ) than the one found using the nonpolarisable FF (3.9 Å ). These results show that the description of the H-bond interactions for this type of IL provided by a nonpolarisable FF is noticeably over-structured with respect to a polarisable FF. Indeed, CL&Pol shows much better agreement with FPMD data, in good accordance with the fact that FPMD takes into account, by construction and inclusion of the electronic degrees of axes represent the distance of the hydrogen atom H 2 of the cation from the F site of the BF - 4 anion participating to the formation of H-bond. The vertical axes of each panel, instead, show the angle defined by the vector oriented from H 2 to the C 2 atom of the cation and the vector connecting H 2 and F as sketched in Figure 3.11. In this configuration the C 2 is a donor and the F is an acceptor, so to validate the hypothesis that we have the formation of H-bonds it is particularly insightful to calculate the angle formed by donor-H-acceptor. In the angle analysis presented in Figure 3.11 value of 180 • indicates that the atoms C 2 , H 2 and F are aligned, thus the H-bond results perfectly linear. An intense maximum arises in the region around 2.0 Å for a range of value between 125 • and 180 • for the IL pair on which we focus. This originates from the H-bond formed by the H 2 atom of the cation and and the acceptor F atom of the BF - 4 . By increasing the length of the cation alkyl chain, a narrowing of the dispersion of H-bonds appears with more localized regions around 2.0 Å and angles varying between 135 • -180 • . This is in perfect agreement with the features displayed by the PCFs in Figure 3. data obtained from FPMD data for DMIM-BF 4 and EMIM-BF 4 systems. The FPMD data show a position of the strong H-bond at values similar to those computed from the classical MD simulations, but the distributions are less localized and affected by a larger broadening. This can be ascribed both to the oscillations of the electronic structure and to the reduced statistics of FPMD simulations in terms of system size and simulation time scales. The Figure 3.12 reports analogous two-dimensional maps of the CDFs in the case of different anions at fixed (DMIM) cation. Note that for BF 4 and PF 6 the acceptor anion atom is F whereas for TFSI is the sulphonyl O atom. In the present case, corresponding to the shortest cation alkyl chain namely DMIM the differences in these distributions can be clearly ascribed to the nature of the anion. CDF peaks are generally well localized. The strongest localization is observed for BF - 4 for a distance of about 2.0 Å and an angle of 150 • -165 • . However for Cl, PF 6 , and TFSI anions a broader distribution is observed corresponding to a distribution of angles oscillating between 100 • -160 • , 85 • -160 • , 120 • -160 • respectively, and this is not entirely surprising. Indeed, the small Cl anion is more mobile and can find stable positions in a slightly wider angular range than BF 4 and PF 6 . On the other hand, the larger and the flexible TFSI anion, the latter having also different conformers (cis and trans) [START_REF] Ori | Ionic liquid confined in silica nanopores: molecular dynamics in the isobaric-isothermal ensemble[END_REF][START_REF] Ishisone | Structural, dynamical, and electronic properties of the ionic liquid 1-ethyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide[END_REF][START_REF] Lassègues | Raman and ab initio study of the conformational isomerism in the 1-ethyl-3-methyl-imidazolium bis (trifluoromethanesulfonyl) imide ionic liquid[END_REF], require more space to be accommodated beside DMIM cations and this accounts for a more dispersed CDF distribution.

Electronic properties and dipole moments

As mentioned, the FPMD approach provides an explicit and direct access to the detailed electronic structure of a system. Beside the occupied states always computed on-the-fly during a FPMD, we included additional twenty unoccupied states in the diagonalization of the Kohn-Sham Hamiltonian on selected configurations specified in the ongoing discussion to compute the electronic density of states (DOS) and to identify the HOMO (higher occupied molecular orbital) and the LUMO (lowest unoccupied molecular orbital). We use the standard notation HOMO and LUMO for these states, although it has to be kept in mind that we are computing Kohn-Sham orbitals and not exactly molecular orbitals. This type of calculation was performed for DMIM-Cl and DMIM-BF 4 by averaging the results over several configurations sampled from the dynamical trajectories at 300 K. These two systems were chosen being the former one of the best characterized in literature, and the latter the main focus of the following Chapters of this work, since it represents a reference system for the realization of EDLT-FET devices. The result is shown in Figure 3.13. Within the computational approach used here and taking into consideration the systematic underestimation of the energy gap affecting DFT based calculations using semilocal XC, a gap of ∼2.1 eV and ∼3.7 eV were found separating the top of the valence band from the bottom of the empty states for the two systems. For the case of DMIM-Cl, we have been able to calculate the energy gap also using the hybrid XC B3LYP [START_REF] Ad | Density-functional thermochemistry. iii. the role of exact exchange[END_REF] by reducing the PW cutoff to 50 Ry, and we computed a gap increased up to ∼3.5 eV. For both ILs, the HOMO wavefunction was found delocalized over the anion (Cl and BF 4 ) and the LUMO wavefunction delocalized over the cation (DMIM and EMIM), which is in fair agreement with previous calculations done with hybrid functionals and the D2/D3 Grimme schemes for the vdW dispersion forces [START_REF] Grimme | Performance of dispersion-corrected density functional theory for the interactions in ionic liquids[END_REF]. To quantify the partial charges of each atom, we exploited the electronic density ρ(r), available on the fly all along the dynamics, and computed for several uncorrelated configurations the value of the Bader charges [START_REF] Henkelman | A fast and robust algorithm for bader decomposition of charge density[END_REF][START_REF] Fw | A quantum theory of molecular structure and its applications[END_REF] of each atom of the cation and the anion composing the IL system. This method, widely applied for many systems, is based on the partitioning of the charge density into atomic volumes referred to as Bader volumes. Typically, there is one charge density maximum at each atomic center and one Bader volume for each atom. The dividing surfaces separating these volumes, also called zero-flux surfaces, lie in the bonding regions between atoms. The result of this analysis is summarized in Table 3.4. The computed Bader charges of each atom do not undergo any significant fluctuation all along the dynamics. The total Bader charge of DMIM and EMIM is equal to +0.95 and +0.94 and that of Cl and BF 4 is -0.94 and -0.95e, confirming the cationic and anionic nature of the two moieties as a whole and the detailed distribution of the partial charges contributing to these overall values. These values are found slightly lower to the point charges used in the original of CL&P FF, which is one of the reasons of the often more-structured IL models obtained by classical MD with such FF with respect to FPMD results. Indeed, to improve the agreement of the description of the dynamical behaviour of IL the atomic charges of this FF are typical down-scaled [START_REF] Yue | Comparison between ab initio molecular dynamics and opls-based force fields for ionic liquid solvent organization[END_REF]. A further approach to char-

µ = N I=1 q I R I - n W F C i=1 f i r i (3.3)
where N is the number of atoms in a IL ion, q I is the valence (in the case of pseudopotentials approaches) atomic charge of atom I, R I are the coordinates of atom I, n WFC is number of WFCs (equal to the number of electrons) and f i is the occupation of the ith Wannier orbital, namely f i = 2 in a spin-restricted calculation and f i = 1 in the spinunrestricted case, and r i are the cartesian coordinates of each WFC (i. e. center of mass of the corresponding Wannier function) belonging to an IL cation or anion. The analysis of the value of the dipole moment µ = |µ| was performed on the basis of WFCs computed over several uncorrelated configurations for each IL system sampled along the dynamical trajectories obtained at 300 K for DMIM-Cl, -BF 4 ,-PF 6 , -TFSI and EMIM-BF 4 . It is of particular interest to recall that the only gauge invariant dipole having a physical meaning in charge neutrality condition of the simulation cell is the one between the cation and the anion [START_REF] Wendler | Locality and fluctuations: Trends in imidazolium-based ionic liquids and beyond[END_REF]. In this analysis, the computed the dipole moments were obtained between each cation-anion pair present in the simulation cells and the statistics is extracted based on both the number of pairs and the different configurations sampled along the trajectories.

The result is shown in Figure 3.15. The distribution of the dipole moment show an average In the case of DMIM-Cl, while the FWHM does not undergo any change, the distribution shifts toward slightly higher values with a peak centered at 28.0 D. The most remarkable displacement of the distribution was found for the case of DMIM-TFSI, where the average dipole moment takes the value of 33.0 D, consistently with the more bulky anion structure. A general picture that can be extracted is the relative stability and insensitivity of the local dipole moment to the type of anion (Cl, BF 4 , PF 6 ), provided that such an anion has not a complex molecular structure, for the same IL cation (DMIM). An analogous minimal effect has been found for different alkyl chains in the cases of DMIM and EMIM when BF 4 is the cation. The larger dipole moment found in the case of DMIM-TFSI, as opposed to the other cases presented in Figure 3.15 is mainly due to the more complex molecular structure of the TFSI with respect to the other smaller and roughly spherically isotropic anions (Cl, BF 4 , PF 6 ). The bigger TFSI molecule has a larger steric hindrance and is accommodated at larger distances from the cation, as shown in the analysis of the PCFs of Figure 3.8. As a consequence, the relative distance between the center of mass of the positive and negative charges increases, thus accounting for an increased dipole moment, in line with the results obtained for EMIM-TFSI [START_REF] Ishisone | Structural, dynamical, and electronic properties of the ionic liquid 1-ethyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide[END_REF]. As expected, the presence of an extended H-bond network induces a significant enhancement of the dipole moment of a corresponding dimer in the gas phase (14.3 D, 13.6 D, 13.9 D, 14.1 D, 17.1 D, for DMIM-Cl, DMIM-BF 4 , DMIM-PF 6 , EMIM-BF 4 and DMIM-TFSI respectively). This enhancement, typical of any polar liquid, is again an indicator of the presence of H-bonds and a fingerprint of the crucial electrostatic interactions between cation and anion. Since the dipole moment of an ionic liquid is experimentally difficult to measure, this data can provide a useful quantitative analysis to complement experimental probes. This type of WFCs data has been often used in the past to quantify sort of "individual" dipole moments for the IL cation and anion [START_REF] Wendler | Locality and fluctuations: Trends in imidazolium-based ionic liquids and beyond[END_REF] and, despite the fact that such WFCs exploitation has no physical meaning whatsoever, it might be useful in the future for further refining of the dipole and polarizability parameters (Drude particles) used in polarisable FF such as Cl&Pol [START_REF] András | Comparison between ab initio and polarizable molecular dynamics simulations of 1-butyl-3-methylimidazolium tetrafluoroborate and chloride in water[END_REF].

Conclusions

In this chapter a series of alkyl-imidazolium-based ILs targeted in a variety of electronic applications was studied with the scope of unraveling their atomistic structure, specific interactions and microscopic features. By resorting to a combination of classical and firstprinciples dynamical simulations, it was shown that the composition of the ILs influences the structure, the nature and topology of the H-bond network and, in turn, the electronic properties, quantified in terms of dipole moment distribution of the different ILs. This study makes use of an identical simulation protocol, based on both classical model potentials and DFT-based first-principles molecular dynamics, of a series of ILs on which forefront research is focused. This comparative study, done within an identical and consistent paradigm is prone to give a useful guidelines in the selection of the best suited IL systems according to the specific application targeted where the ability of H-bond may play a significant role. Moreover, the results obtained here from both CMD and FPMD simulations shall be useful for future refining of CMD FFs that account explicitly for the polarizability features of ILs, with the aim to further improve the quantitative prediction of the outcome of CMD simulations.

Chapter 4

Modelling ionic liquids at the interface with MoS 2 layers

Context and objectives

Knowledge of how the molecular structures of ionic liquids affect their properties at the interface with 2D layered materials is essential for a rational design of ILs for optoelectronic devices and catalysts. In the present chapter, CMD and FPMD simulations are used to investigate the structural, bonding and dynamical properties of a series of IL at the interface with a MoS 2 layers, widely used in a wealth of experiments. This IL-substrate system will be referred to as IL@MoS 2 hereafter. In particular, we pay special attention to the IL 1-ethyl-3-methylimidazolium tetrafluoroborate (EMIM-BF 4 ) being nowadays one of the most tested IL for EDLT-FET applications and for which a detailed analysis of its structural properties (ions density profiles, molecular orientation, pairs correlation functions), electrical properties (charge density profiles) and dynamical properties (mean square displacements, residence times) has been done in the present doctoral project. These analyses are discussed in the case of EMIM-BF 4 at the interface with a pristine MoS 2 bilayer (i.e. without any surface defect) and with defected MoS 2 bilayers. Subsequently, a detailed comparison between CMD and FPMD results is presented and critically analyzed, providing a detailed atomistic insight into this interface, especially when the electronic structure is explicitly taken into account as in the case of FPMD. In the final part of this chapter, based exclusively in CMD simulations, the effects of the type of IL inorganic anion (Cl, BF 4 , TFSI) and of the chain length of the IL cation (DMIM, EMIM, HMIM) on the ions density and charge density profiles at the interface are systematically discussed.

Computational details 4.2.1 Classical molecular dynamics

The interface models simulated by CMD and studied in this chapter were obtained by the following procedure:

1. Production of a representative bulk model made of 500 ions pairs of the selected IL with hexagonal geometry and equilibrated at 300 K, following the same thermal cycle detailed in Chapter 3, the hexagonal geometry being dictated by the crystallographic symmetry of the substrate;

2. Adding the obtained bulk IL model on top of a MoS 2 bilayer (at a distance of 5 Å from the surface) with a vacuum volume on top of the IL phase. The equilibration of the IL phase was achieved through a mild thermal cycle (from 500 K to 300 K with steps of 100 K equilibrated for 2 ns at each temperature);

3. A final equilibration period of 8 ns at 300 K, of which the last 5 ns were used for the production phase and data analysis.

The MoS 2 bilayer studied here was the hexagonal conformation corresponding to the allotropic form of 2H-MoS 2 in the crystallographic unit cell (space group P63/mmc, Crystallography file CIF ID 9007660) [START_REF] José-Yacamán | Structure and catalytic properties of hexagonal molybdenum disulfide nanoplates[END_REF][START_REF] Schönfeld | Anisotropic mean-square displacements (msd) in single-crystals of 2h-and 3r-mos 2 s[END_REF]. This unit cell was replicated by 15x×15y ×1z times for the CMD simulations in order to build a supercell for a total of 1350 MoS 2 atoms (900 S and 450 Mo). Such model corresponds to a cell of a=47.415 Å, b=47.415 Å, c=δ Å and α=90, β =90, γ=120, with δ that depends on the IL system (EMIM-Cl: 200.0 Å ; EMIM-BF 4 : 200.0 Å ; EMIM-TFSI: 500.0 Å ; DMIM-TFSI: 400.0 Å ; HMIM-TFSI: 500.0 Å ). For these interfaces models, it has been used primarily the CL&P force field proposed and developed by Canongia Lopez and coworkers [START_REF] Canongia | Modeling ionic liquids using a systematic all-atom force field[END_REF] with the parametrization of Koddermann et al. [246] to better account for the dynamical properties of the bulk ILs simulated.

As previously mentioned, this specific FF was fitted on ab initio calculations and allows to simulate a wide variety of ILs chemical structures including the specific ones targeted in this chapter and it has been already tested over a variety of layered materials. For the MoS 2 bilayer, the FF of Sresht and coworkers [START_REF] Sresht | Quantitative modeling of mos 2 -solvent interfaces: Predicting contact angles and exfoliation performance using molecular dynamics[END_REF] have been used. In these simulations, the interaction between MoS 2 and the IL were simulated by a Lennard-Jones (LJ) potential with LJ parameters (σ ij and ǫ ij ) between unlike sites calculated using the Lorentz-Berthelot mixing rules [START_REF] Delhommelle | Inadequacy of the lorentz-berthelot combining rules for accurate predictions of equilibrium properties by molecular simulation[END_REF]. Periodic boundary conditions were used in the simulations. The dispersive interactions were neglected beyond a cutoff of 1.2 nm and the electrostatic interactions were calculated using the Ewald sum technique to avoid finite size effects (the Ewald parameters were selected to lead to an accuracy 10 -5 ). The simulations were integrated using the Verlet algorithm with a timestep of 1 fs and all the CMD simulations were performed with the MD code DLPOLY [START_REF] Smith | The dl poly molecular dynamics package[END_REF].

First-principles molecular dynamics

The dynamical FPMD simulations presented and discussed in this Chapter were carried out by resorting to the Car-Parrinello method [START_REF] Car | Unified approach for molecular dynamics and densityfunctional theory[END_REF] and used to model two interface systems made of 20 ions pairs of EMIM-BF 4 (480 atoms) at the interface with a MoS 2 bilayer without any surface defect (150 atoms) and with two S surface vacancies (147 atoms), respectively The valence-core interaction was described by norm-conserving numerical pseudopotential of the Troullier-Martins [START_REF] Troullier | Efficient pseudopotentials for plane-wave calculations[END_REF][START_REF] Troullier | Efficient pseudopotentials for plane-wave calculations. ii. operators for fast iterative diagonalization[END_REF] type, and the exchange and correlation interactions are expressed according to the functional proposed by Becke [START_REF] Ad | Density-functional exchange-energy approximation with correct asymptotic behavior[END_REF] and Lee, Yang, and Parr [START_REF] Lee | Development of the colle-salvetti correlation-energy formula into a functional of the electron density[END_REF] (BLYP), respectively. In the case of Mo, semicore states were included to ensure a good description of the energetics and to account for electronic states that could be sensitive to polarization effects. Valence electrons are treated explicitly and represented on a plane-wave basis set with an energy cutoff of 100 Ry, with the sampling of the Brillouin restricted to the Γ point. A fictitious electron mass of 400 a.u. and a time step of 3.0 au (0.0726 fs) were used to obtain optimal conservation of the constants of motion. A spin-unrestricted approach (LSD) was adopted in all the FPMD simulations where MoS 2 was involved. FPMD simulations were performed in the NVT canonical ensemble with the ionic temperature controlled with a Nosé-Hoover [START_REF] Nosé | A molecular dynamics method for simulations in the canonical ensemble[END_REF][START_REF] Nosé | A unified formulation of the constant temperature molecular dynamics methods[END_REF][START_REF] Wg | A molecular dynamics method for simulations in the canonical ensemble[END_REF][START_REF] Wg | Canonical dynamics: Equilibrium phase-space distributions[END_REF] thermostat chain [START_REF] Gj | Nosé-hoover chains: The canonical ensemble via continuous dynamics[END_REF]. Periodic boundary conditions were applied throughout the simulations. A preliminary study of the MoS 2 bulk phase has been performed, for which the long-range van der Waals dispersion interactions, not included in the BLYP functional, were modelled according to the maximally localized Wannier functions (MLWFs) and centers [START_REF] Pl | Van der waals interactions in dft made easy by wannier functions[END_REF][START_REF] Ikeda | Role of van der waals corrections in first principles simulations of alkali metal ions in aqueous solutions[END_REF][START_REF] Marzari | Maximally localized wannier functions: Theory and applications[END_REF] (WFCs) as previously used for the bulk IL in Chapter 3. WFCs allowed obtaining detailed insights into the bonding features of MoS 2 . However, for the modeling of the interface models, the Grimme formalism with the D2 scheme was employed. The reason behind this choice was based on the lower computational cost needed with D2 (-24% time per step) with respect to the MLWFs-based method. The initial configurations of the EMIM-BF 4 @MoS 2 interfaces were prepared by classical MD as described in the previous subsection. To this purpose, the final configuration equilibrated at 300 K by classical MD was then used as starting point of the FPMD simulation that was performed to equilibrate the interface model up to 300 K. More details about the FPMD simulations are given in section 4.3.3 of this Chapter.

The case study of EMIM-BF 4 @MoS 2 interface

In this section the chemical and physical properties at play at the interface between EMIM-BF 4 and a bilayer of MoS 2 are reported. Firstly, the results obtained via CMD are presented, discussed and compared with the FPMD results. In Figure 4.1 it is reported a series of snapshots showing the temporal evolution of the formation process used to build representative models of IL@MoS 2 interfaces by CMD. What is described in the following for EMIM-BF 4 4@MoS 2 can be extended to all the interface models produced by CMD discussed in this chapter. At the initial time t 0 , the pre-equilibrated bulk EMIM-BF 4 phase is positioned at distance of ∼5 Å from the MoS 2 surface. After a short time, within few ps, the first ions start to move in the direction of the MoS 2 surface, attracted by the electrostatic charge of MoS 2 layers. Within the first ∼100ps the entire EMIM-BF 4 bulk phase, made of 500 ions pairs, moved in the proximity of the MoS 2 surface, however at this stage not all the MoS 2 surface is completely covered by the IL. It is only by ∼500 ps that the MoS 2 surface is completely covered by the EMIM-BF 4 showing a thickness of the IL layer of about ∼8 nm along the c axes. At this stage, the ions at the proximity MoS 2 surface begin to form molecular layers. The relatively mild thermal cycle (from 500 K to 300 K) allows to promote local rearrangements through the first few layers from the surface and within the layers. At the end of this mild thermal process, the IL model if found equilibrated within 2 ns at 300 K. The last 5 ns at 300 K are used as production phase and for the data analysis with respect to the properties discussed in the following sections. The chemical physical properties of this interface models will be discussed for the two cases with MoS 2 surface without any surface defect (i.e. pristine surface) and with surface defects.

EMIM-BF 4 at the interface with a pristine MoS 2 bilayer

Structural properties: ions density profiles, pair correlation functions and ions orientation

To understand the structural organization of the IL ions and their interactions at the interface with MoS 2 a few, yet crucial, structural properties were calculated. The first one is the analysis of the ions density profiles with respect to the distance normal to the MoS 2 surface. The very sharp and intense first density peak of EMIM suggests a well-defined arrangement of the cation ions at the interface. While the very broad first density peak of BF 4 suggests a much more disorganized arrangement of the anion ions at the liquid-solid inter-face. This result, which reveals the propensity of the IL cation to interact in a stronger way with the MoS 2 surface than the anion, can be ascribed to the electronegative character of the S atoms in the top layer of MoS 2 . Beyond the third layer, i.e. >1.8-2 nm, the typical bulk IL density is recovered. The analysis of the composition of the first IL layer in contact with MoS 2 pinpoint to a close charge-neutral first layer with ∼47% and ∼53% of EMIM and BF 4 , respectively. These results suggest that the exact structure of the confined IL is driven by the subtle interplay between overscreening and crowding effects [START_REF] Mz | Double layer in ionic liquids: Overscreening versus crowding[END_REF]. Indeed, given the fact that in this system the IL is both the charged component and the solvent, large ion concentrations are reached so that the corresponding Debye length is very small ∼0.1 nm. Such a small Debye length indicates that the electrostatic interactions are 'overscreened' (despite the fact the dielectric constant remains relatively small, ǫ ∼ 10ǫ 0 ) and that other effects such as the so-called ion specific effects and crowding effects become predominant.

To gain insights into the local structure EMIM ions at the interface, the pair correlation functions gij (r) between the H atoms of EMIM and the S atoms of the MoS 2 surface were analysed. Of notice, the fact that gij (r) have been corrected for the finite size of the sample [START_REF] Ori | Surface of glassy ges 2 : A model based on a first-principles approach[END_REF].

gij (r) = g ij (r) f (r) with f (r) = 1 - r 2h (4.1)
where h is the thickness of the sample. While this correction in general is not needed to compare the position of the peaks in the g ij (r) functions, it allows correcting the peak intensities for finite size of the samples. Two types of gij (r) functions were analysed, the first one between the sulfur atoms of the MoS 2 top layer (H X -S) and the hydrogen atoms of the EMIM cation located in the first layer with respect to the solid surface (Figure 4.3A), i.e. Z < 1 nm, and the second between hydrogen atoms of the EMIM and the fluorine atoms of the BF 4 (H X -F) in the bulk phase, i.e. Z in the range 3-6 nm (Figure 4.3B). The g H-F (r) functions show the typical preferential H-bond promoted by H 2 atom of EMIM, as found in bulk EMIM-BF 4 (Chapter 3), showing the shortest H-F distance with respect to the other H atoms of EMIM (0.21 nm). Whereas the g H-S (r) functions show the H 6,7,8 of the -CH 3 (methyl) and -CH 2 CH 3 (ethyl) alkyl chains of EMIM to be the closest to the S atoms of MoS 2 layer (0.31 nm) with the imidazolium ring H atoms to be farther away (0.45 nm). This difference suggests a preferential flat orientation of the imidazolium ring of EMIM with respect to the MoS 2 surface, also supporting possible vdW interactions promoted by the short alkyl chains of EMIM and MoS 2 . This observation is also supported to the visual inspection of the top and lateral views the EMIM ions within the first layer at the interface (Figure 4.3C).

The structure of EMIM-BF 4 in contact with MoS 2 was further investigated by measuring a structural order parameters. The following order parameter was computed to monitor the orientation of the ring of the alkyl-imidazolium cation in the first layer in contact with MoS 2 :

S ring (θ r ) = 3 2 cos 2 (θ r ) - 1 2 (4.2)
where the brackets mean that the order parameters were averaged over all the configurations recorded during the classical MD simulations. This orientational parameter is defined from the angle θ r between the normal vector passing through the geometry center of the cation ring and the normal vector to the MoS 2 surface. The Figure 4.4 shows the distribution of the order parameter S ring (θ r ) of EMIM in the first interfacial layer. This interfacial layer is defined when the cation center of mass is located at a distance within 1 nm from the surface. Near the surface EMIM ions have a preferential orientation, the ring of the cation tends to be flat i.e. parallel to the crystalline surface showing a S ring (θ r ) ∼ 1.0, and with than 90% of the cations at the interface showing an orientation parameter close to 0.925. This result further supports the previous observation about the fact that the H 2 atom of EMIM did not show any preferential interaction with MoS 2 , being this located and constrained between the two N atom of the flat imidazolium ring. While, the H atoms of the methyl and ethyl groups of the cation having more degrees of freedom, can promote easily a stronger interaction with the S atoms of MoS 2 . Of notice, however, the fact that the latter H atoms show as shortest distance to the S atoms of MoS 2 a value of about ∼ 0.28 nm, which indicates very weak H-bond at maximum.

Electrical properties: charge density profiles

Although in this Chapter it is considered a MoS 2 bilayer that is not electrified, which will be the focus of Chapter 5, a first hint on the electrical properties of the IL phase at the interface could be obtained by analysing its charge density distribution as a function of the distance from the S atoms of the top layer of MoS 2 . Figure 4.5 shows the total charge distribution of EMIM-BF 4 as well as the individual contribution from EMIM and BF 4 . The total charge profile shows a first positive peak at ∼0.28 nm, followed by a pronounced negative peak at ∼0.34 nm, a subsequent pronounced positive peak at ∼0.41 nm and a final broader negative peak at ∼0.5 nm. Farther away from the MoS 2 surface, it is still visible a very small positive contribution at ∼0.8 nm before vanishing to a constant zero value. It is evident that the ions charge balance of the IL phase occurs at a shorter distance (∼0.9 nm) to the MoS 2 layer with respect to the ions density profiles (∼1.6-1.8 nm). By comparing the total profile with those of EMIM and BF 4 , it is clear that the above peaks are the results of distinct charge peaks of the two counter-ions. The first peak can be ascribed to the H atoms of the methyl and ethyl groups of EMIM that are directly interacting with the S atom of MoS 2 , whereas the following peaks are the result of the overlaps between the positively charged imidazolium rings (that brings the majority of the positive charge) and the negative charges of BF 4 ions.

Dynamical properties: mean squared displacements and residence times

To study the dynamical properties of EMIM-BF 4 at the interface, two properties were evaluated: the mean-square-displacement (MSD(t)) of the center of mass of the ions as a function of time and the residence time (R(t)). The 3D MSD(t) was calculated as follows:

M SD(t) = |r(t) -r(0)| 2 = 1 N N i=1 |r i (t) -r i (0)| 2 (4.3)
where N is the number of particles to be averaged, r i (0) is the reference position of the i-th particle at an initial time t = 0, and the r i (t) is the position of the same i-th particle at a later time t. In order to assess the influence of the interaction with the MoS 2 surface could have on the dynamical properties of the IL ions at the interface, the MSD was calculated at 300K in three different regions as a function of the distance from MoS 2 surface. The first region L1 is defined as the one close to the solid surface (with a thickness of 1 nm).

The second region L2 is defined as the next region apart from the solid surface with a thickness of 2 nm. And the third region L3 is defined as the farthest apart from the solid surface with a thickness of 4 nm. This analysis shows a significant difference in the mobility of the two ions in the three considered regions. In region L1, the closest to the MoS 2 surface, the two ions show minimal mobility, reaching a value after 4 ns of about ∼30 Å2 and ∼40 Å2 for EMIM and BF 4 , respectively. In region L2, that corresponds to the region between 1 to 3 nm from the solid surface (meaning farther than the two first ions density peaks), the two ions show a larger MSD than the one found in region L1 reaching a value of about ∼100-110 Å2 after 4 ns. In region L3, the farthest from the solid surface, the two ions show an even larger MSD of the order of ∼400 Å2 after 4 ns. The very slow dynamics of both ions in region L1 is ascribed to the high ions density of this region and further support a finite chemical interaction of the IL with the surface of MoS 2 as a result of electrostatic and vdW interactions. The faster dynamic of the two ions in region L2 is ascribed to the lower density and farther distance to the MoS 2 surface. The very high mobility of the two ions found in region L3 might also be due to the effect of the IL-vacuum interface which increases the spatial degrees of freedom of both ions. Moreover, the dynamics of the cation is always found faster then that of the anion, as it is typically reported for the bulk IL. The faster dynamics of EMIM arises from its homogeneous mild electrostatic field, which leads to a weaker interaction with the its environment. In contrast, the BF 4 , which shows a strong electrostatic field highly localized due to its size, interacts strongly with its environment so that its dynamics is slower. Also, the flexible (short) ethyl alkyl chain and the planar structure of the cation might also contribute to its faster dynamics [START_REF] Ori | Ionic liquid confined in silica nanopores: molecular dynamics in the isobaric-isothermal ensemble[END_REF][START_REF] Ori | Structure and dynamics of ionic liquids confined in amorphous porous chalcogenides[END_REF].

Residence time

To obtain further insight into the dynamical properties of EMIM-BF 4 , the average time spent by the IL at the MoS 2 is a useful parameter. This quantity, generally referred to as residence time R(t) inn the literature, was calculated from the cation and the anion auto-correlation functions,

R(t) = β i (t)β i (0) (4.4)
where β i (t) = 1 if the ion is located at time t within the first layers in contact with the MoS 2 surface and 0 otherwise. A cation (anion) is considered to be within the first layers if its distance form the surface is <1.5 nm. R(t) was computed at 300K, 350K and 400K and the calculated R(t) are reported in Figure 4.7. Two important remarks should be noted, firstly a large variation in the residence time of the ions as a function of the temperature (decreasing in the time) can be observed. The residence time R(t) is greatly influenced by the temperature just as the MSD(t) which makes it possible to indirectly establish a certain correlation between R(t) and the M SD(t). An increase in temperature, however small, leads to an increase in thermal agitation in the liquid medium, and therefore an increase in the kinetic energy of ions. The latter thus see their mobility increase, so it becomes easier for an anion or a cation located in the first interfacial layers to cross the border of the first layer and to move to the upper layers (space jump). These jumps are accompanied by a reduction in the residence time in the first layers, and the more the temperature increases the more this transition becomes frequent. Secondly, it can be noted that the residence time R(t) of the cation is reduced compared to that of the anion. This final result, is apparently in contradiction with the MSD results reported previously. However, this behaviour can be explained considering that EMIM cation, as described above, shows a larger flexibility with respect to the BF 4 anion due to its structure and electrostatic potential. These features allow the EMIM to show a higher vibrational mobility with respect to BF 4 as well as along the xy plane, however, due to its favorable interactions with the S atoms of MoS 2 , its leaning to move farther away from the solid surface (along the z axes) is reduced.

EMIM-BF 4 at the interface with MoS 2 defected surfaces

In this section, the effect of the presence of surface defects on the MoS 2 layer on the degree of structuring of the IL EMIM-BF 4 are evaluated. In particular, for this specific substrate, the most common defects found (or induced by electron beam, plasma or chemical treatment) experimentally are sulfur vacancies [START_REF] Bertolazzi | Molecular chemistry approaches for tuning the properties of two-dimensional transition metal dichalcogenides[END_REF]. For this reason, we have focused specifically on this type of defects. The defected MoS 2 surface models studied in this work were obtained by removing a certain amount of S on the top layer of MoS 2 and to maintain the overall neutral balance of the bilayer the appropriate number of Mo atoms was removed in the bottom MoS 2 layer, far from the liquid-solid interface (balance of the electrostatic charge means two S atoms were removed together with one Mo atom at a time). Three different surface defects densities, with randomly distributed sites, have been studied: 4.0×10 14 sites/cm 2 (corresponding to 10 sites in total), 8.0×10 14 sites/cm 2 (corresponding to 20 sites in total) and 2.0×10 15 sites/cm 2 (corresponding to 50 sites in total). These values are comparable to the ones typically found (or induced) in real MoS 2 samples (7.4×10 14 -1.8×10 16 sites/cm 2 ) [START_REF] Hu | Wafer-scale sulfur vacancy-rich monolayer mos2 for massive hydrogen production[END_REF]. In Figure 4.8 are shown the top views of the surface models of pristine MoS 2 bilayer as well as two of the three models with surface defects (here shown the two with 4.0×10 14 sites/cm 2 and 8.0×10 14 sites/cm 2 ). In the following the structural, electrical and dynamical properties of the interfaces made of MoS 2 with surface defects are presented and discussed in comparison with those of MoS 2 without defects.

Structural and electrical properties

In order to understand the degree of influence and the impact that sulfur surface vacancies can have on the interaction between the MoS 2 bilayer and the IL ions, the ions density profiles of the models with the three different surface defects densities were analyzed and compared with the original MoS 2 pristine model without defects. Figure 4.9 shows the total and the individual ions density profiles of the four systems. The calculated density profiles show clear differences, especially within the region close to the MoS 2 surface. In terms of total ions density profiles, the main observation that can be note is the decrease of the intensity of the shoulder (centered at about ∼0.5 nm) of the first density peak with increasing the defects density. This shoulder completely disappears for the model with the highest surface defect density. The decrease of this shoulder, which it is reminded was mainly due to a density peak of BF 4 ions at this distance from the MoS 2 surface, is replaced by the appearance of a peak at much shorter distance (centered at ∼0.23 nm). This new density peak, entirely due to BF 4 ions increases with the increase of the surface defects density. A consequence of the formation of this new BF 4 ions density peak, is the sharpening of the first density peak of EMIM ions (intensity increased from ∼12.8 ions/nm 3 to ∼15 ions/nm 3 ), which can be ascribed to a further 'flattening' of the cations parallel to the MoS 2 surface. This arrangement picture is supported by the visual inspection reported in Figure 4.10, where it can be observed the clear saturation of all the surface S vacancies by BF 4 ions and a very pronounced flat arrangement of the EMIM cations. Over the surface defects, the BF 4 ions are often found on top of the S vacancies with one B-F bond aligned normally to the MoS 2 surface, showing a F• • • Mo distance of about ∼0.33-34 nm with the three Mo atoms of the layer below and a distance of about ∼0.34-0.35 nm with the six S atoms of the surface hexagonal layer. By integrating the ions density profiles up to the density minimum after the first main peak (0.6 nm), the composition of the first layer has been analysed. The model with a defect density of 4.0×10 14 sites/cm 2 showed a first layer made of about ∼47% of EMIM and ∼53% of BF 4 ions very similar to the case with no surface defects. While the two models with a higher surface defects density showed a slight majority of cations in the first layer with a ∼53% of EMIM and ∼47% of BF 4 . At the same time, we can report that the surface defects of the models with 4.0×10 14 sites/cm 2 and 8.0×10 14 sites/cm 2 have been found all (100%) saturated with BF 4 ions whereas the model with largest surface defects density was found saturated only for 76% of the defect sites. These results can be ascribed to the fact that the model with the largest surface defects density shows S vacancies very close to each other in some cases, which can not be saturated by BF 4 because of steric hindrance and electrostatic repulsion. To get a better insight into the electrical properties of the IL phase in proximity of the defected MoS 2 surface models, the ions charges density profiles provide useful information. This information is reported in Figure 4.9 where these charge density profiles are shown along with the ones computed on the non-defective pristine MoS 2 surface for a direct comparison. These results support the previous analysis based on the ions density profiles. Indeed, with the presence of surface defects and the formation of a first BF 4 ions density peak at short distance to the surface also a first negatively charged peak is found. The intensity of negative peak increases with the increase of surface defects density. Moreover, the following positive peak becomes more and more narrow with the increase of the surface defects intensity in complete accordance with the previous observations made based on the ions density profiles. One last observation is of interest, and related to the decrease of all the following charge density peaks after the second one with vanishing charge peaks farther than ∼0.8 nm.

Dynamical properties

To obtain insights into the effects of the surface MoS 2 defects on the dynamical properties of EMIM-BF 4 ions, the MSD(t) analysed in region closest to the MoS 2 (corresponding to the region labelled L1 in the previous section). Figure 4.11 shows the calculations of the calculated MSD values of the center of masses of the IL cations and anions along 4 ns of CMD simulation at 300K. The main observation that can be pointed out from this analysis is related to the fact that at short time, the increase of surface defects induces a slow down of the ions for both EMIM and BF 4 which is found close to be proportional to the surface defects density. However, at longer times (>3 ns) the mobility of EMIM and BF 4 ions in the case of defected MoS 2 surfaces all tend to the same MSD values(with EMIM cation remaining slightly more mobile than BF 4 anion) suggesting that this mobility dependence on the defects density is lost once the defects sites are well saturated. Overall, the mobility of both cations and anions in the case of the pristine MoS 2 models always show faster dynamics with respect to the defected MoS 2 models. This can be easily be motivated due to the stronger interaction of BF 4 with the defects sites and consequent increase in density of the first layers of EMIM ions which further slow down the ions mobility.

EMIM-BF 4 @MoS 2 interfaces: FPMD results

In this section, the results obtained by FPMD are presented and discussed. On a first instance, we presented some preliminary calculations performed on a MoS 2 bulk phase as a propaedeutic study necessary before moving to the composite EMIM-BF 4 @MoS 2 system. Then we focus on the results obtained on our EMIM-BF 4 @MoS 2 interface model in which the MoS 2 surface is supposed to be defect-free before moving to a surface model containing two defects as detailed in the following discussion.

Modeling 2H-MoS 2 layered structure

To check the applicability of the DFT and FPMD frameworks used so far for the study of the different ILs, we did a series of calculations to study 2H-MoS 2 layers. It is wellknown and widely reported in the literature that to accurately describe the structural, mechanical and electronic properties of MoS 2 crystal phases, the use of hybrid function-als (such as B3LYP and HSE06) combined with an accurate inclusion of vdW dispersion forces (such as the Grimme-D2 scheme) or the use of more accurate approaches such as post-Hartee-Fock and GW approximation methods, along with spin-orbit coupling, are needed [START_REF] Jk | The indirect to direct band gap transition in multilayered mos 2 as predicted by screened hybrid density functional theory[END_REF][START_REF] Je | Nature and evolution of the band-edge states in mos 2: From monolayer to bulk[END_REF][START_REF] He | Stacking effects on the electronic and optical properties of bilayer transition metal dichalcogenides mos 2 , mose 2 , ws 2 , and wse 2[END_REF][START_REF] Peelaers | First-principles study of van der waals interactions in mos 2 and moo 3[END_REF]. These methods, have been shown to be essential for an accurate description of the electronic properties of 2H-MoS 2 crystal phase, especially for a correct estimation of the peculiar (i.e. direct or indirect) band gaps of 2H-MoS 2 monolayer and the bulk phase as already explained in the first Chapter of this doctoral thesis and summarized in Figure 1.6. Nevertheless, these methods can only be applied to systems of a few tens of atoms because of the heavy computational workload due to the inclusion of the exact exchange, spin-orbit coupling and Green Functions calculation in the case of GW approaches. This is in clear contrast with the scope of modeling a realistic system targeted in the present thesis work. In fact, we aim at modeling interface systems made of over 600 atoms to be closer to the actual experimental conditions, a system size clearly beyond the reach of these computationally expensive approaches. Furthermore, doing dynamics with these post-HF or GW methods is unfeasible and simple static calculations are insufficient to describe the complexity of an Il-substrate system. For these reasons, in this work, the FPMD simulations of the interface models are performed employing the semi-local BLYP functional complemented with a suitable vdW scheme for the inclusion of dispersion forces. In the following we present the preliminary calculations used for the assessment of the crystal structure parameters that correspond to the minimum (Kohn-Sham) total energy for the 2H-MoS 2 model within the BLYP [START_REF] Lee | Development of the colle-salvetti correlation-energy formula into a functional of the electron density[END_REF][START_REF] Ad | Density-functional exchange-energy approximation with correct asymptotic behavior[END_REF] exchange-correlation and the vdW interaction computed via the maximally localized Wannier functions (MLWFs) and related Wannier centers (WFCs) [START_REF] Pl | Van der waals interactions in dft made easy by wannier functions[END_REF][START_REF] Ikeda | Role of van der waals corrections in first principles simulations of alkali metal ions in aqueous solutions[END_REF][START_REF] Marzari | Maximally localized wannier functions: Theory and applications[END_REF]. The optimal values of 2H-MoS 2 crystal parameters were found by systematically varying the a and b parameters of a supercell obtained by replicating 5×5×1 the unit cell along a, b and c axes. Such supercell resulted made of 150 total atoms (50 Mo and 100 S). The computational approach used, is based on the energy and wavefunction optimization of the systems trough a preconditioned conjugate gradient method with a convergence criteria of 10 -6 , followed by a constrained FPMD simulation applying a friction force (ion velocities scaled by a factor 0.98 at each step) for a total of 2000 steps so as to optimize the model at T∼0 K and atomic forces <10 -4 . The Kohn-Sham total energy of the different models was always accompanied by the calculation of the corresponding stress tensor to quantify the degree of residual stress on the systems. Figure 4.12 shows both profiles of KS energy and stress tensor as a function of the variation of a and b parameters.

From the analysis of this plot, the optimal a and b parameters at which it corresponded the minimum of KS energy with the chosen computational approach were found, which corresponded to variation of 3% with respect of the experimental values and with a residual stress of about ∼1.5 GPa. In Table 4.1 are reported the corresponding obtained cell parameters (in terms of unit cells) together with the main interatomic distances of interest for 2H-MoS 2 : Mo-S bond distance, the distance between Mo atoms of the two MoS 2 layers (Mo• • • Mo), and the distance between S atoms of the two MoS 2 layers, corresponding to the vdW gap interlayer distance (S• • • S). These data is compared with the experimental values and with other computational approaches reported in literature. In comparison with the experimental values, it can be observed that the obtained model shows a slight overestimation of the Mo-S bond distance and a reduced S• • • S interlayer distance as a result of adapting to the larger a and b cell parameters. As expected, by employing more advanced methods, such as those-based on hybrid functionals, the structural parameters of 2H-MoS 2 cell can be described with higher accuracy in addition to a better description of its electronic properties.

To probe the chemical bonding in the 2H-MoS 2 layers, we analysed the electronic struc- ture through the WFCs positions and associated dispersion of the MLWFs. As outlined in Chapter 2, a WFC indicates the maximum probability for the location of an electron (or electron pair in spin-unrestricted approaches) in a quantum system [START_REF] Pl | Van der waals interactions in dft made easy by wannier functions[END_REF][START_REF] Ikeda | Role of van der waals corrections in first principles simulations of alkali metal ions in aqueous solutions[END_REF][START_REF] Marzari | Maximally localized wannier functions: Theory and applications[END_REF]. The analysis of the WFCs with respect to the nuclear positions provides an insight into the chemical bonding of several types of materials, and in the past it has been applied to the study of liquids such as water [START_REF] Boero | First principles study of propene polymerization in ziegler-natta heterogeneous catalysis[END_REF][START_REF] Sharma | Ab initio molecular dynamics with maximally localized wannier functions[END_REF][START_REF] Resta | Electron localization in the insulating state[END_REF] and amorphous chalcogenides such as selenides and sulphides [START_REF] Ori | Surface of glassy ges 2 : A model based on a first-principles approach[END_REF][START_REF] Le Roux | Structural properties of glassy ge 2 se 3 from first-principles molecular dynamics[END_REF]. The nth MLWF w n (r) and its WFC r WFC n = (x n , y n , z n ) are obtained as unitary transformation (iterative) of the Kohn-Sham orbitals ψ i (r)

w n (r) = Nocc i=1 p e -A p i,n ψ i (r) (4.5)
where N occ is the total number of occupied states (=electrons in a spin-unrestricted approach), A p i,n is a matrix generalization of the Berry phase connector, and p is the order of the iteration [START_REF] Gh | The structure of electronic excitation levels in insulating crystals[END_REF]. Among all the possible (equivalent) unitary transformations, in the case of MLWF we chose the one that minimized the spread of the orbitals. This is obtained by minimizing the spread functional Ω. This quantity is the difference between the Wannier charge density with respect to its own center of charge, thereby representing the spatial extension of the Wannier orbital:

Ω = Nocc i=1 i|r 2 |i -i|r|i 2 (4.6)
The WFC of each Wannier orbital n is given by

x n = - L x 2π Im ln n |e - i2π x L x | n (4.7)
where L x corresponds to the length of the simulation cell along x. Similar equations are used for the coordinates along y and z. 89 Å corresponding to the 56% and 43% of the sum of the radii, respectively [START_REF] Alvarez | A cartography of the van der waals territories[END_REF]. and the covalent radii 1.54 Å , and 1.05 Å ; corresponding to the 60% and 40% of the sum of the radii, respectively, [START_REF] Verónica | Covalent radii revisited[END_REF]), we observe that the computed Mo-S bond distance is very closed to the sum of the reported covalent radii. On the basis of the obtained results and these considerations, the Mo-S bonds found in our models of 2H-MoS 2 corresponds to what in the literature is retained as a covalent bond. Yet, this is a too simplified description of a proper Mo-S bond as explicitly indicated by the location of the WFCs. In fact, our MLWF-based analysis shows rather clearly that a polarization effect exists and translates into WFC MSB closer to S than Mo (for distance corresponding to the 35% of the bond distance versus the 40% as reported in literature) [START_REF] Verónica | Covalent radii revisited[END_REF]. In this respect, we have here a much more precise characterization of nature of a Mo-S bond than the simple one given just in terms of textbook values of vdW and covalent radii.

Chemical interactions at the liquid-solid interface: ions density profiles and EMIM orientation

In this section, the results obtained by the FPMD simulations of the EMIM-BF 4 @MoS 2 interfaces are presented and discussed. The initial configuration of an interface model made of 20 EMIM-BF 4 ions pairs in contact with a bilayer of 2H-MoS 2 was obtained by classical MD. Then, the c parameter of cell was set to 5 nm, which allowed sufficient vacuum space on top of the IL phase to avoid any significant interaction with the replicas boxes due to the applied PBC along the FPMD simulations. The FPMD computational approach has been applied as follows: in a first stage a friction force was applied on the atoms' velocities in order to bring the system close to 0K within 2000 steps. This stage was followed by a short simulation in the NVE ensemble to monitor its stability and then by over more than ∼10ps in the NVT ensemble at a temperature of 300 K, as described in the computational methods section. Along this FPMD simulation, the bottom MoS 2 layer has been kept frozen throughout the simulation whereas the top layer has been let free to move. Figure 4.14 shows the profile of the Kohn-Sham potential energy along the simulation at 300 K for EMIM-BF 4 @MoS 2 interface. The constant average value of the KS potential energy beyond ∼4 ps and the energetically well-separated temperatures of ions and electrons along the simulation time prove the numerical stability of the simulation and the reliability of the computational protocol used [START_REF] Blöchl | Adiabaticity in first-principles molecular dynamics[END_REF]. A small expansion along the z direction of the IL phase was recorded during the ini- The total ions density profiles obtained by FPMD shows significant layering of the IL, overall supporting the results obtained by classical MD although with some appreciable differences. Overall, three layers can be identified by the peaks' maximum of the density distribution. The first one, the closest to the MoS 2 surface, is found more broadened with respect to the the CMD density data and with a peak centered at a distance 0.45 nm, a value that lies in between the first peak position and its shoulder from the CMD data. The second and the third density peaks are found with similar intensity whereas the following density oscillations vanish at larger distance from the MoS 2 surface. The second peak' intensity and position are found very close to those obtained by CMD (0.72 nm), whereas the third peak is slightly more intense in the FPMD data and with broadened maximum over the range 1.1-1.35 nm. By breaking down the total profile into the contributions from EMIM and BF 4 ions, we can infer further observations. The cation density profile shows the first two density peaks very narrow and located at found in CMD. Only the first peak is found at lower intensity with respect to the CMD data. The third peak obtained by FPMD is found at a relatively shorter distance with respect to the CMD data (1.1 nm versus 1.25 nm). The BF 4 ions density profile obtained by FPMD shows the first three peaks with higher intensity and slightly farther shifted at larger distances with respect to the CMD data. In particular, the first peak is found much more structured and centered at 0.5 nm, whereas the second and third peaks are found less structured and centered at 1.1-1.35 nm and 1.75-1.9 nm (from CMD: 0.3-0.6 nm (1 st ), 1.09 nm (2 nd ), and 1.6 nm (3 rd )). Overall, the IL ions density profiles obtained by FPMD support the ions density layering at the interface obtained by CMD, with density distributions that start at the same distance from the MoS 2 surface. However, the FPMD data show density peaks less intense and centered at distances from the interface shifted at larger values. The analysis of the first density peaks obtained by FPMD, support the preferable flat orientation of EMIM with respect to the MoS 2 surface and depicted by the narrow density peak. Whereas it suggests a more organized arrangement of BF 4 ions within the first layer, showing a first peak particularly pronounced with repsect to the broad one found in CMD. In terms of the structural order S ring (θ r ) of EMIM, the FPMD data show a trend very similar to the one obtained by CMD, with the majority of EMIM lying close to parallel with respect to the MoS 2 surface. The only appreciable difference, is related to a S ring (θ r ) distribution peak more broadened in the range 0.7-0.9 whereas in CMD 85% of EMIM at the interface was found showing an orientation parameter close to ∼0.90.

Electronic properties at the liquid-solid interface: Bader charges, charge density profiles and electron density

As mentioned previously, the FPMD approach provides an explicit and direct access to the detailed electronic structure of a system. This information allowed us to quantify the the partial charges of each atom, exploiting the electronic density ρ(r), available on the fly and computed for several uncorrelated interface configurations the Bader charges, as done previously for the case of bulk IL. The result of this analysis is summarized in Table 4. In a similar way, the Bader charges of B and F are almost the double in absolute values with respect to the values used in CL&P FF (B,F: 3.0,-0.99 versus 0.96,-0.49), while the total anion charge is similar. Analysing the first molecular layers region (<1.0 nm) of both the total and individual cation and anion charge density profiles, we can observe a general shift of the negative and positive peaks towards larger distance with respect to the MoS 2 surface and in comparison to the CMD data. This observation is in agreement with our previous analysis based on the ions density profiles. Moreover, the charge density peaks beyond the first one are found more broadened in comparison to CMD data, the farther from the MoS 2 . A further possible contribution to the noisy character of FPMD data could also be originated by the relatively small size of the IL phase. The direct access to the detailed electronic structure of a system provided by FPMD, allows to inspect the electron density distribution of the EMIM-BF 4 interface. Figure 4.20 shows the electron density distribution of the interface at study, computed at ∼9.5ps of the In this section, we present the preliminary results obtained by modeling a EMIM-BF 4 @MoS 2 interface with two S vacancies on the top MoS 2 . We recall that the initial configuration was obtained by classical MD, as previously described, and such preparation allowed the saturation of the two S vacancies with two BF 4 ions in close contact with the S and Mo atoms neighbouring the S vacancies (∼3.4 Å and ∼3.5 Å, respectively). The FPMD computational procedure adopted here is the same as the one of the previous section, however in the present case we simulated up to about ∼6ps of total dynamics, reaching 300K at about ∼4ps. We recall that prior to the NVT simulation, we performed a constrained run in order to optimize the initial interface model (obtained by classical MD) at 0K. This first stage allows minimal local molecular rearrangements in the search of a relative minimum of energy. 4.2). The relatively stronger nature of these chemical interactions is supported by the analysis of the electronic density (Figure 4.20). Indeed, contrary to the interface without any S surface defect, the interface model investigated here shows a non-zero electron density between one F atom of a BF 4 anion and the Mo and S atoms neighbouring the S vacancy site. These results suggest a possible way to induce specific chemical interactions at MoS 2 surface defects that could be exploited for its rational surface chemical doping, and thus a possible way to modify MoS 2 electronic properties. In this section, we undergo a systematic study by means of classical MD simulations of different IL@MoS 2 interfaces with the aim to assess the effect of the IL anion size and chemistry and the role of the length of the alkyl chain of AI-IL cation.

Effect of AI-IL anions size and chemistry

In order to evaluate the size and chemistry effects of the IL anion on the interface properties IL@MoS 2 interfaces, we simulated three interfaces bearing Cl, BF 4 and TFSI as IL anion.

The same computational procedure used previously for EMIM-BF 4 @MoS 2 interface was used here. Figure 4.23 shows typical snapshots of the final configuration equilibrated at 300K. Already inspecting the different final thicknesses of the three IL phases (EMIM-Cl: 7nm; EMIM-BF 4 : 8nm; EMIM-TFSI: 14nm) for the same number ions pairs simulated (500 pairs), it is possible to observe the effect of the different anion size. To gain further insights in the different structure and behaviour of the three IL pairs, we plot in Figure 4.24 their ions density and charge density profiles with respect to the distance from the MoS 2 surface).

Several distinct remarks can be outlined. The first peak of the total ions density profiles is found at the same distance from the MoS 2 surface for all the three pairs, underlining no anion effect whatsoever on the distance interaction at the interface. However, the intensity of this first peak is drastically reduced moving from Cl>BF 4 >TFSI, which can be ascribed directly to a size effect of the anion. The smaller the anion size the more the cation first layer is found structured, showing the first density peak with increasing intensity. The anion size effect induces also the broadening of the first peak of the total density profile. This is ascribed to the volume occupied by each anion, resulting in the decrease of the intensity of the original first peak' shoulder of the BF 4 density profile (located at 0.5 nm) for the case of Cl and an increase of such shoulder for the case of TFSI. The IL density layering oscillations show different extension departing from the MoS 2 surface for the three IL pairs. EMIM-Cl and EMIM-BF 4 show vanishing of the density oscillations at about 1.5nm whereas TFSI shows it at about 2.3-2.5nm. The analysis of the ions charge density profiles show the most striking differences for the three systems. In particular, the small size of Cl ions allows a very large number of both ions within the first layer which result in very intense positive peak close to the solid surface. This first positive peak is also found in the case of BF 4 and TFSI, however its intensity is drastically reduced with increasing the anions size, with also the consequence of showing a second negative peak with increasing intensity as the anions size increases. A clear contrast between the ions density profiles and the ions charge density profiles is related to the extension of the density oscillations layering with respect to the distance from the MoS 2 surface. In fact, if the former shows extension that increases with the size of the anions, the latter show the opposite. The charge layering is found extended up to 1.8nm for EMIM-Cl, up to 1.4nm for EMIM-BF 4 and up to 1.0nm for EMIM-TFSI.

Effect of AI-IL cation alkyl-chain length

In order to inspect the role of the cation size and especially the length of its alkyl chain, in this section we analysed the ions density and ions charge density profiles with respect to the distance from the MoS 2 of three IL made of DMIM-TFSI, EMIM-TFSI and HMIM-TFSI (Figure 4.25). As previously found in the case of different anions, also in this case the position of the first total ions density peak does not change in the three IL models, supporting the picture of promoted electrostatic and dispersion forces as the main interactions at play at the interface. However, in this case we observe a greater structuring of the first density peak of the cation by reducing alkyl chain length. All the three models show the cations ions to be lying parallel to the solid surface. Both DMIM and EMIM promote clearly up to three ions density peaks as far as 1.6nm and 2.3-2.5 nm, respectively. Whereas HMIM shows a low intense first peak followed by very broadened density peaks of average intensity close to the bulk region. In terms of charge profiles, all the three models show a initial sequence made of a first small positive peak, followed by very intense negative peak and a third very intense positive peak. Further peaks of alternate charges are present up to 1.0nm for the three models. The intensity of the second and third peaks decrease with increasing the length of the alkyl chain of the cation. This result can be ascribed to the tendency of AI-IL cation with long alkyl chains to promote separated polar and apolar regions within the IL phase, the latter being mainly constituted by the alkyl chains where mostly dispersion forces are at play.

Conclusions

We studied a series of alkyl imidazolium-based ILs in contact with pristine and defected 2H-MoS 2 layered substrates to rationalize the different interactions occurring at the liquidsolid interface. To achieve, we used simulation protocols based on classical and firstprinciples molecular dynamics, the study was able to provide insights into the different interactions taking place at the liquid-solid interface with some details on electronic structure. The results of this study can be useful in rationalizing the behavior of these types of materials in real-world applications and can potentially inform the design of new ionic liquids with improved properties for use in various fields. Furthermore, the study's findings can also be useful for understanding the behavior of other liquid-solid systems, as the molecular interactions at play are likely to be similar in nature. Additionally, the use of both classical and first-principles molecular dynamics simulations in this study adds to the robustness of the findings, as it allows for a more comprehensive understanding of the interactions taking place at the liquid-solid interface.

Chapter 5

Preliminary study of EMIM-BF 4 at electrified MoS 2 surfaces

Context and objectives

In electric double layer applications such as batteries, double layer field-effect transistors (EDL-FETs), and double layer capacitors (EDLC), electrolytes are required to function jointly with 2D nanomaterials, such as transition metal dichalcogenides of the same type discussed in the previous chapters. At the electrode-electrolyte interfaces, the ILs form ultra-thin electrical double layers (EDLs). EDL is a common acronym used to identify a complex interfacial structure in systems with mobile charges determined by the electrostatic interaction of ions and/or molecules with a charged surface. Despite of the term 'double' in its name, the EDL consists of several distinct molecular or charge layers, which are relevant in highly concentrated electrolytes as in the case of ILs [START_REF] Kornyshev | Three-dimensional double layers[END_REF][START_REF] Ivaništšev | Restructuring of the electrical double layer in ionic liquids upon charging[END_REF]. The structure and dynamics of EDLs are among the key factors determining the capacitance and charging speed of EDLC and the gating speed of EDL-FETs and thus regulating the performance of the devices [START_REF] Zhao | Adding solvent into ionic liquidgated transistor: the anatomy of enhanced gating performance[END_REF][START_REF] Watanabe | Application of ionic liquids to energy storage and conversion materials and devices[END_REF]. A fundamental understanding of the behavior of EDLs from the molecular level is essential for understanding how electrolytes behave at electrified interfaces and, in turn, how to provide a guideline to design EDLs-based devices with improved performance [START_REF] Fujimoto | Electric-double-layer field-effect transistors with ionic liquids[END_REF]. From the viewpoint of fundamental science, ILs represent an interesting case of highly concentrated electrolytes that possess in general a more chemically inert character with respect to other systems and are more resistant to higher voltages than standard diluted electrolytes such as aqueous ionic solutions. An intense research activity worldwide has shown that the behaviour of ILs near charged surfaces is basically different from the behaviour of standard solvent-based electrolytes. Indeed, ILs display a number of phenomena typically absent in standard electrolytes, namely over-screening, charged driven structural transitions and ion crowding at the luquid/solid interface [START_REF] Fedorov | Ionic liquids at electrified interfaces[END_REF][START_REF] Mz | Double layer in ionic liquids: Overscreening versus crowding[END_REF][START_REF] Baldelli | Surface structure at the ionic liquid-electrified metal interface[END_REF]. In this context, MD simulation provide a powerful tool for investigating the properties of EDLs, leading to major insights into the structure of the interfacial layers at an atomic-scale resolution escaping experimental probes. This helps to understand voltage-capacitance relations and the charging mechanisms. Several works focused on how the properties of EDLs are affected by the nature of the electrodes, such as, for instance, the curvature, the roughness, and the presence of nanopores [START_REF] Gong | Effect of side chain modifications in imidazolium ionic liquids on the properties of the electrical double layer at a molybdenum disulfide electrode[END_REF][START_REF] Sr | Theory and practice in constant potential molecular dynamics simulations[END_REF][START_REF] Ms | Mean-field and modified poisson-boltzmann approaches for modeling electrochemical energy storage systems[END_REF][START_REF] Merlet | On the molecular origin of supercapacitance in nanoporous carbon electrodes[END_REF][START_REF] Merlet | Computer simulations of ionic liquids at electrochemical interfaces[END_REF][START_REF] Burt | A review of molecular modelling of electric double layer capacitors[END_REF][START_REF] Vatamanu | Charge storage at the nanoscale: understanding the trends from the molecular scale perspective[END_REF]. Typically, charge-constant MD and constant-potential MD (ccMD and cpMD, respectively) are the two methods routinely used to get an insight into the structural and electrical properties of EDLs-based systems. The former has the advantage of being applicable to a variety of materials and allows to obtain "rapidly" fundamental insights in EDLs properties. The latter has the advantage of giving access to the charging process of the electrode, although it was not conceived to be applied to electrode materials that are not metallic. The latter method has been exploited in multiple application to the study of EDLC systems in the last decade. However, recent works have shown that these two computational schemes give similar interfacial structures and capacitances within a certain range of voltage drops. In the present Chapter, ccMD simulations are used to investigate the structural and electrical properties of EDLs promoted by EMIM-BF 4 at the interface with charged MoS 2 surfaces. In particular, we give a detailed account of structural properties (ions density profiles, molecular orientation) and electrical properties (charge density profiles) of our targeted systems. These analyses are discussed in terms of overscreening and overcrowding models used to characterize ILs-based systems. In the final part, we show how the use of FPMD allows to obtain a better understanding of the electronic density characterizing the EDL and promoted by the EMIM-BF 4 @MoS 2 interface. Altogether, the results presented in this Chapter offer a critical analysis of EDLs for the case of EMIM-BF 4 in contact with charged MoS 2 surfaces, providing a detailed atomistic insight into this interface.

Computational details

The interface models simulated by ccMD and studied in this Chapter were obtained by the following procedure:

1. Production of a representative bulk model made of 500 ions pairs of EMIM-BF 4 with hexagonal geometry and equilibrated at 300 K, following the same thermal cycle detailed in Chapter 3 and 4, the hexagonal geometry being dictated by the crystallographic symmetry of the substrate;

2. Inserting the obtained bulk IL model in-between two MoS 2 bilayers (at a distance defined by obtaining a final IL density close to the exp. value). The equilibration of the IL phase was achieved through a simulation in the NVT canonical ensemble at 400 K for 1 ns.

3. A final equilibration period of 1.5 ns at 400 K was used for the production phase and data analysis. This model corresponds to the zero-charge system, in which the opposite MoS 2 bilayers (i.e. electrodes) show the same neutral total charge.

4. The final configuration of the previous step was used as the initial structure for the ccMD simulations of the charged interfaces, where the the atoms point charges of the Mo and S atoms of each bilayer were specifically parametrized to obtain a given final surface capacitance.

5. Each charged model was equilibrated for 2.5 ns in the NVT canonical ensemble at 400 K, of which the last 1.5 ns were used for the production phase and data analysis.

The MoS 2 bilayer studied here is one described in the previous Chapter, with a unit cell was replicated by 15x × 15y × 1z times for the ccMD simulations in order to build a supercell for a total size of 1350 MoS 2 atoms (900 S and 450 Mo). This specific model corresponds to a cell of a=47.415 Å, b=47.415 Å and α=90, β =90, γ=120. The c parameter was defined as c=96.920 Å in order to obtain a final IL phase, confined in-between the two MoS 2 bilayers, showing the experimental density. For these interfaces models, it has been used the CL&P force field proposed and developed by Canongia Lopez and coworkers [START_REF] Canongia | Modeling ionic liquids using a systematic all-atom force field[END_REF] with the parametrization of Koddermann et al.[246] to better account for the dynamical properties of the bulk ILs simulated. For the MoS 2 bilayer, the FF of Sresht and coworkers [START_REF] Sresht | Quantitative modeling of mos 2 -solvent interfaces: Predicting contact angles and exfoliation performance using molecular dynamics[END_REF] have been used with slightly modified point charges. We employed the Mo and S charges recently proposed by Bing and coworkers [START_REF] Bi | Co-ion desorption as the main charging mechanism in metallic 1t-mos 2 supercapacitors[END_REF], which have been obtained at zero voltage by cpMD simulations of 1T-MoS 2 supercapacitors and correspond to +0.62 and -0.31 for Mo and S atoms, respectively. The electrified MoS 2 surfaces were obtained by defining Mo and S charges with equal in magnitude but opposite in sign surface charge densities, ±Q S . The following surface charge densities were simulated: ±3.0, ±7.0, ±13.0, ±36.0, ±73.0 µC/cm 2 . Figure 5.1 shows a few examples of models simulated with different surface charge densities. The temperature used for the CMD simulations performed in this Chapter was set to 400 K. This value was chosen to enhance the mobility of the IL ions by thermal energy. Preliminary tests performed in the range 300-500 K did not show significant differences in the final results. Following the approach used in recent works, we calculated the electrode potentials (i.e. voltage) corresponding to the charge density distributions of MoS 2 bilayers simulated in this Chapter. In the practice, the potential φ(z) profile along the z direction (electrode-to-electrode) was calculated by double integration of the charge distribution of the electrolytes, ρ chg (z), and of the charge density on the electrodes, Q S , based on the one-dimensional Poisson equation [START_REF] Gong | Effect of side chain modifications in imidazolium ionic liquids on the properties of the electrical double layer at a molybdenum disulfide electrode[END_REF][START_REF] Zhang | Evaluation of static differential capacitance at the [c 4 mim + ][tfsa -]/electrode interface using molecular dynamics simulation combined with electrochemical surface plasmon resonance measurements[END_REF].

d 2 φ(z) dz 2 = - ρ chg (z) ǫ r ǫ 0 (5.1)
where ǫ 0 is the dielectric constant of the vacuum and ǫ r is the relative permittivity. The values φ(z) obtained correspond to: 0.70 V, 1.60 V, 2.96 V, 8.15V and 16.5V for the five surface charge densities simulated here, respectively. Table 5.1 reports the point charges used here for S and Mo at the two electrodes. [START_REF] Bi | Co-ion desorption as the main charging mechanism in metallic 1t-mos 2 supercapacitors[END_REF] for positive and negative electrodes set at ±2 V (2 V: +0.583/-0.263; -2 V: +0.664/-0.361).

The dispersive interactions were neglected beyond a cutoff of 1.8 nm and the electrostatic interactions were computing using the Ewald sum technique to avoid finite size effects, by selecting the Ewald parameters to ensure an accuracy of 10 -5 . The equations of motion were numerically integrated using the Verlet algorithm with a timestep of 1 fs and all the ccMD simulations were performed with the DLPOLY package [START_REF] Smith | The dl poly molecular dynamics package[END_REF]. One interface model made of 20 EMIM-BF 4 ions pairs and one MoS 2 bilayer, upon equilibration, has shown on the top MoS 2 monolayer a surface charge density of -73.0 µC/cm 2 , balancing the opposite charge of the bottom monolayer. A vacuum region was inserted on top of the IL resulting in a total c dimension of 5 nm, identical to the one previously used and reported in chapter 4. This interface model was equilibrated at 400 K for 2.5 ns with CMD in the NVT canonical ensemble. The electronic structure and density of the final configuration of this interface was computed with FPMD, by performing a short constrained dynamics (1000 steps) with the application of a friction force (0.95) on the atomic velocities (damped dynamics). The computational details of the FPMD simulation are identical to those reported in Chapter 4 and coded in the CPMD v4.

Electric double layer (EDL) properties of EMIM-BF 4 at

electrified MoS 2

EDL structure: ions density profiles and ions orientation

To understand the structural organization of EMIM-BF 4 ions and their interactions at the interfaces with the electrified MoS 2 surfaces a few, yet crucial, structural properties were calculated. The first one is the ions density profiles with respect to the distance orthogonal to the MoS 2 surfaces. Figure 5.2 shows the distribution of the center of mass (COM) of the EMIM cations and the BF 4 anions, together with the total ions density profile along the z direction (electrode-to electrode) in the condition of zero-charge electrode. For this model and surfaces condition (zero-charge), significant layering was observed, translating into relatively large density oscillations of the IL anion and cation density profiles in proximity of the MoS 2 surfaces, we observe symmetrical profiles of the density oscillation at the two electrolyte-electrode interfaces, in agreement with the fact that at zero-charge condition the two MoS 2 electrodes are identical. The number and intensity of distinct layers (e.g. three), identified by the distinct maxima of the corresponding distributions, with the first peak showing a small shoulder, in analogy with the results reported in Chapter 4. This result indicates a minimal effect of the small change of the Mo and S point charges used in this Chapter and can be regarded as a proof of the reliability of our set-up. The, first peak is relatively sharp and characterized by a high intensity, the second and following peaks show decreasing intensities when the distance from the solid MoS 2 substrate increases. By disentangling the total profile into the contributions of the IL cation and anion, we can identify three peaks for EMIM located at 0.35, 0.72 and 1.2 nm from. While BF 4 shows a first broaden peak extending over the range 0.3-0.6 nm without a clearly defined maximum, followed by a second peak at 1.09 nm and a third one, barely observable, at 1.6 nm. The large first density peak of BF 4 is responsible for the shoulder of the first peak of the ions density profile. The density peaks of EMIM and BF 4 in part alternate and in part overlap, suggesting a certain degree of interdigitation between the IL cation and anion layers. In what follows, we analysed the charged MoS 2 surface and this analysis is done in terms of the characteristic phenomena attributed at IL in these interface models. Specifically, the ions layering in the direction orthogonal to the charged surface and the overscreening and crowding features [START_REF] Mz | Double layer in ionic liquids: Overscreening versus crowding[END_REF][START_REF] Kornyshev | Three-dimensional double layers[END_REF][START_REF] Ivaništšev | Restructuring of the electrical double layer in ionic liquids upon charging[END_REF]. We recall that the overscreening regime is a distinctive feature of the EDL structure in ILs and results in the formation of a layered structure near the charged electrode surface, where the first layer contains more counter-charge than needed to compensate for the surface charge. The excess of charge is typically due to counterions (ions with opposite charge with respect to the surface) and it is compensated in the subsequent ion layers. The crowding scenario refers to the situation in which the number of counterions required to compensate the surface charge cannot be accommodated in a layer with a thickness of a single ion diameter. Therefore, the first ionic layer becomes wider due to inclusion of several sublayers of counter-ions. The crossover between the overscreening and crowding regimes occurs via formation of a clear monolayer of counterions. Typically, this monolayer structure is characterized by an absence of oscillations in the electrostatic potential φ(z)-profiles and a linear potential drop between the electrode and the monolayer of counter-ions. Figure 6.2 in the Annex shows the total and individual cations and anions ions density distributions in the different conditions of electrode charging. Note that in this Figure, the positive electrode is at z=0 nm (left side) and the negative electrode at z=7.0 nm (right side). When the electrodes are charged, the total ions density profiles do not show anymore the characteristic symmetrical distribution of the first density peaks shown at zero-charge. By deepening the analysis of the cations and anions density distributions in the ranges 0-2.0 nm and 5.0-7.0 nm (Figure 5.3), several observations in terms of EDL structure can be highlighted. When the MoS 2 electrodes are charged up to 7µC/cm 2 (i.e. ∼1.6 V), we observe an increase in intensity of the EMIM density in the first layer at the negative electrode concomitantly with a shift towards larger distances from the surface of the BF 4 ions. Conversely, at the positive electrode, the opposite process occurs but to a smaller extent. When the charge on the electrodes increases to 13µC/cm 2 (i.e. ∼2.96 V), at the negative electrode the density profile shows an alternately oscillating multilayer structure, and the EDL is in the overscreening regime. This regime is present at the positive electrode but is much less evident. At higher electrode charge, we observe a clear formation of EMIM and BF 4 monolayers at the negative and positive electrodes, respectively. At the the highest surface charge, we remark a thickening and structuring of the first monolayer, showing a second (positive electrode) and up to a third (negative electrode) small but well-defined density peaks of the same charge, following the first one. These new peaks indicate the formation of a well structured first monolayer. Farther apart from these first layers, we observe again the alternately oscillating multilayer structure of the EDL. Driven by the strong opposite charge, with the increase of the MoS 2 surface charge, both ions of the first layer approach each electrode surface. At high surface charge, the formation of a clear monolayer of counter-ions, with respect to the charge of the surface, identifies a regime attaining the crowding state of a EDL. However, we could not find a model showing the ideal crowding state with a monolayer and following overlapped density peaks of the two ions. The Figure 5.4 shows the distribution of the order parameter S ring (θ r ) of EMIM in the first interfacial layer, the same used in Chapter 4. At low electrode charge (3µC/cm 2 (i.e. ∼0.7 V)), near the negatively charged surface the EMIM ions have a preferential orientation: the ring of the cation tends to be parallel to the crystalline surface showing a S ring (θ r ) ∼ 1.0, as found in the case of neutral charge discussed in Chapter 4. However, a minor amount of EMIM ions exists, characterized by S ring (θ r ) over the whole range between -0.5-0.8, indicating also a random distribution of a few EMIM ions. At higher charge (7 µC/cm 2 (i.e. ∼1.6 V)), the situation evolves with a peak at about S ring (θ r ) ∼ 0.7 but with the majority of the ions randomly distributed. At 13 µC/cm 2 (i.e. ∼2.96 V), all the EMIM ions are randomly distributed, whereas at higher charge the majority of the ions show a preferential orientation perpendicular to the solid surface. By changing the orientation, a higher number of EMIM ions can move-in in the first interfacial layer to counter-balance the higher negative charge of the MoS 2 surface, in agreement with the higher intensity of the first peak of the ions density distributions. The orientation of the BF 4 ions can be deduced from the ions density profiles. Indeed, at the positive electrode, a double first peak arises within the first interfacial layer with an intensity ratio of ∼1/3, which identifies a BF 4 oriented in a way in which three F atoms come in close contact with the Mo S surface and the fourth F at the BF 4 vertex is oriented towards the IL phase. This orientation can be observed at every electrode charge, however its sharpness and structuring increases upon increase of the positive charge of the electrode. Figure 5.5

shows typical examples of the aforementioned EMIM and BF 4 orientations, showing the different local organization and of EMIM and BF 4 ions described above. 5.7), we can see that by increasing the surface charge a remarkable increase of the intensity of the first peaks arises. This is due to the orientational structuring of the ions at the the interfaces. Moreover we can also observe a change of charge of cations and anions within the first peaks, due again to the strong degree of structuring of the ions. In fact, at the positive electrode, the first and third negative peak identify the F atoms of the BF 4 tetrahedral anions, wheres the positive second peak is due to the positive point charge of the B atoms. The same considerations can be extended to the EMIM cations, where the formation of negative peaks are due to EMIM atoms carrying negative point charges. As discussed previously, a drop of the total charge to zero is reached very rapidly after the first peak, which indicates the reach of the crowding-like regime of the EDL. Altogether, these results promote the strengthening of the chemical interactions at play between EMIM-BF 4 ions and MoS 2 when a voltage is applied. The relatively stronger nature of the electrostatic forces at play is supported by the analysis of the electronic density (Figure 5.8). Indeed, contrary to the interface at zero-charge (Figure 4.20), the interface model investigated here with a surface charge of 73µC/cm 2 shows the electron density of the IL due to EMIM imidazolium-rings and the top layer of MoS 2 much closer, with regions of non-zero density. This observation is also supported by the analysis of the computed Bader charges, that shows an increase in the absolute (negative) values of S atoms at the interface layer with respect to the S atoms in the bottom layer (-0.6 versus -0.5).

EDL electrical properties: charge density profiles

In the following final part, we focus on the effect of S vacancies on the EDL properties when EMIM-BF 4 ions are confined in-between charged MoS 2 surfaces. The reason for the introiduction of this type of defects in our simulated systems has already been given. For this purpose, we selected the surface charge densities of 36.0µC/cm 2 and with of 10 and 20 S vacancies/ Å2 , corresponding to 4.0×10 14 sits/cm 2 and 8.0×10 14 sits/cm 2 , respectively. Figure 5.9 shows the ions density and charge density profiles with respect to the distance from each electrode with surface defects of 10 and 20 S vacancies. The calculated density profiles show clear differences with respect to the models without surface defects, especially within the region close to the MoS 2 surface at the positive electrode. In terms of total ions density profiles, the main observation that can be extracted from these data is the slight decrease of the intensity of the second peak (centered at ∼0.5 nm) by increasing the defects density. The decrease of this second peak is replaced by the arise of a peak at much shorter distance (centered at ∼0.18 nm). This new density peak, entirely due to BF 4 ions increases with the increase of the surface defects density. A consequence of the formation of this new BF 4 ions density peak is the decrease of the intensity of the second positive peak of the charge density profiles, which may suggest a change in the local arrangement of the BF 4 ions in proximity of S vacancies sites. These results further promote the possibility to tune the interactions occurring at the interface between a IL and MoS 2 surfaces by inducing surface defects and thus a possible way to modify MoS 2 electronic properties.

Conclusions

The study of the behavior of EMIM-BF 4 at electrified MoS 2 surfaces using constantcharge MD provided valuable insights into the interactions between IL ions and surfaces at a molecular level unreachable by experiments. The results obtained provide a clear picture of the effects of electrical potential on the structuring and orientation of EMIM and BF 4 ions at electrified MoS 2 surfaces, along with a better insight into the electrostatic forces which are the main driving interactions regulating structural and electronic properties of these interfaces. This work improves our understanding of ionic liquids at electrified surfaces of TMDC and highlights the importance of combining CMD and FPMD to obtain the full picture in terms of structural, bonding and electronic properties for such complex interfaces. Additionally, the findings give a clear indication about the importance of resorting to FPMD-based simulations to complement, correct and improve what can currently be obtained by CMD, whose accuracy depends entirely on the quality of the FF and its parameterization, not always trivial in these complex systems. Furthermore, the results of this study provide a basis for future investigations into the behavior of other ions and materials at electrified surfaces.

Chapter 6

General conclusions and future perspectives 2D semiconducting materials and, in particular, the experimentally widely used transition metal dichalcogenides (TMDCs) are far from having unveiled their full potential and their use in the manufacturing of advanced electronic nanocomponents in various fields demonstrates their effectiveness. In the field of electronics and optoelectronics, TMDCs have become prominent materials and are prone to play a central role in the modernization of the technology of the future. Their compatibility with ionic liquids (ILs) based on alkyl-imidazolium cations have made it possible to form unprecedented IL-TMDCs based devices possessing unexpected and extraordinary properties on which advanced research is focused worldwide. Nonetheless, to date, because of the complexity of these systems, very few studies have been carried out describe at an atomic-level detail, the behavior of ILs when they are in contact with a MoS 2 substrate or when they are used as electrolytes in FET devices. This thesis work is intended as a pioneering effort to fill this lack of knowledge in three main stages: On a first instance, resorting to a combination of classical and first-principles dynamical simulations we have been able, as explained in Chapter 3, to shed some light onto the microscopic structure and specific chemical interactions occurring in a large class of alkylimidazolium ILs. More precisely, it has been shown how the chemical structure and nature of IL ions influence their atomistic structure, the nature and local topology of their Hbond networks and, in turn, the electronic properties, with particular emphasis on dipole moment distributions, Bader charges and electronic density of states. We have shown how non-polarizable FF can be useful to obtain reasonable agreements with FPMD structural results, pointing out how the need for more sophisticated FFs able to take into account explicitly polarization effects essential for a better quantitative agreement with FPMD. In this respect, the FPMD results are the ones intended be used for the development of polarizable FF parameters applicable to a wide range of ILs. The comparative study performed in this Chapter, conducted within an identical and consistent paradigm, is prone to offer useful guidelines in the selection of the best suited IL systems according to the specific application targeted where the structuring and networking ability of H-bonds play a major role. On a second instance, the combined approach of CMD and FPMD allowed in Chapter 4 to improve our fundamental insights into the interactions occurring at the liquid-solid interface between a series of alkyl imidazolium-based ILs and 2H-MoS 2 layered substrates, pristine and defected. The results of this study can be useful in rationalizing the behavior of these types of interface systems in real-world applications and can potentially steer the design of new ionic liquids with enhanced properties for different and specifically oriented application in next-generation electronics. Altogether, Chapter 4 proposes a more comprehensive understanding of the interactions taking place at the targeted liquid-solid interfaces.

On a third instance, Chapter 5 capitalizes the knowledge acquired in the previous Chapters for the atomistic modeling of EMIM-BF 4 at electrified MoS 2 surfaces. In particular, the use of constant-charge MD provided valuable insights into the interactions between this IL and the charged surfaces at the molecular level. The results obtained have shown the effects of electrical potential on the structuring and orientation of EMIM and BF 4 ions at electrified MoS 2 surfaces, pointing out how the electrostatic interactions are the main driving force regulating the behavior and the macroscopic response of these interfaces. This Chapter evidenced the need for improvement of the currently used classical FFs with insights driven by the FPMD approach in order to achieve a full quantitative picture of the structural, bonding and electronic properties for such complex interfaces. From this standpoint, the results presented in Chapter 5 are prone to disclose new perspectives in the design of novel ILs and realization of tailored IL-2D materials interfaces for future nanoelectronics and optoelectronics in addition of possible usefulness also in the field of TDMC-driven surface catalysis. As general perspectives, two main lines for future research can be outlined from this thesis work. First, with the aim of extending the reach of the quantitative degree of prediction power of FPMD, there is the need to refine currently used classical FFs to better account for polarizable effects, which is crucial for ionic liquids-based systems, although sometimes still overlooked. This is not only needed for ILs but also for the study of interfaces made by ILs and 2D semiconducting materials such as MoS 2 . Moreover, the use of these advanced FFs also on the modeling of electrified interfaces through constantcharge and constant potential MD is expected to drastically improve our current, yet still limited, knowledge of such systems which is based mostly on simplified models. These observations are in line with the direction in which a part of the scientific community is going towards, see [START_REF] Gong | Effect of side chain modifications in imidazolium ionic liquids on the properties of the electrical double layer at a molybdenum disulfide electrode[END_REF][START_REF] Bi | Co-ion desorption as the main charging mechanism in metallic 1t-mos 2 supercapacitors[END_REF]. Secondly, due to the exploratory nature of this thesis work, some electronic properties of the 2D TMDC material studied here have been mostly overlooked (i.e. band gap). This is in part due to the primary aim of this study, that was devoted to the set up of a computational approach that assessed the ability of a FPMD scheme (e.g. dispersion scheme) to properly describe the chemical interactions a play in this liquid-solid interface. Indeed, this work set the basis for a future work where also the DFT and FPMD approaches need to be refined in order to directly assessed, in a quantitative way, the effects of the interactions with ionic liquids on the electronic properties of TMDCs.

I wish to reiterate here what I already stated in the first Chapter. Namely that, within the research activities of our institute (IPCMS), this is the first step toward a new research line that is just at its exploratory stage in many respects. My work has been interesting and stimulating, as any scientific area in which one starts moving his (or her) first steps. Yet, unexpected difficulties, intrinsic in any research activity, have been further complicated by the worldwide pandemic crisis that affected all of us and, in particular, students unable to access the instruments and devices, including supercomputers, necessary for their research work. I did my best to overcome also this unforeseen difficulty and I believe that this thesis can fully testify my commitment.

Résumé détaillé de la thèse en français

Liquides ioniques à l'interface avec MoS 2 : propriétés physiques et chimiques à l'interface par dynamique moléculaire classique et premiers principes Les dichalcogénures de métaux de transitions(TMDCs) sont des matériaux 2D qui n'ont pas encore révélé tout leur potentiel. Leur utilisation dans les nanocomposants électroniques avancés et dans divers domaines montre leur efficacité. Dans le domaine de l'électronique et de l'optoélectronique, les TMDC sont considérés comme des matériaux importants et seront sans doute un élément clé dans la modernisation de la technologie future. Leur compatibilité avec les liquides ioniques à base de cations alkyles imidazolium a permis la formation de dispositifs basés sur les interfaces IL-TMDC, montrant des propriétés exceptionnelles et inédites qui sont actuellement au centre de la recherche avancée (voir Figure 6.4). Cependant, jusqu'à présent, en raison de la complexité de ces systèmes, peu d'études ont été menées pour décrire de manière détaillée au niveau atomique le comportement des liquides ioniques (ILs) en particulier les liquides ioniques à base de cations d'imidazolium (Figure 6.5) lorsqu'ils sont en contact avec un substrat tel que MoS 2 ou utilisés comme électrolytes dans des dispositifs FETs. Ce travail de thèse, organisé en 6 chapitres, vise à combler ce manque de connaissances en effectuant une première étude approfondie.

Le Chapitre 1 offre un compte-rendu détaillé des principales caractéristiques physiques et chimiques des ILs et des matériaux TMDCs, ainsi qu'un état complet de leur modélisation atomistique.

Le Chapitre 2 fournit un aperçu du contexte théorique et des deux méthodes utilisées dans ce travail de thèse à savoir: la dynamique moléculaire classique(CMD) et la dynamique moléculaire des premiers principes (FPMD). Les trois chapitres suivants (chapitre 3, 4, 5) présentent les travaux de recherche consacrés à la modélisation des ILs, des TDMCs et de leurs interfaces ILs-TMDCs.

Premièrement, dans le Chapitre 3, la combinaison entre les simulations de dynamiques moléculaires classiques et les simulations de dynamiques moléculaires basées sur les premiers principes de la mécanique quantique a permis d'élucider la structure microscopique et les interactions chimiques spécifiques dans les liquides ioniques. Plus précisément, nous avons montré comment la structure chimique et la nature des ions de ces liquides ioniques influencent leur structure atomique, la nature et la topologie locale de leurs réseaux de liaisons H, en accordant une attention particulière aux propriétés électroniques telles que les distributions de moments dipolaires, les charges de Bader et la densité d'états électroniques (DOS). De plus, l'utilisation d'un champ de force non polarisable nous a permis d'obtenir des accords raisonnables avec les résultats obtenus à partir de la FPMD (voir Figures 6.6 et 6.7) , en mettant en évidence la nécessité de champs de forces plus sophistiqués et complexes, capables de prendre en compte explicitement des effets de la polarisation afin obtenir un meilleur accord quantitatif avec la FPMD. Les résultats de la FPMD obtenus dans ce travail sont destinés à être utilisés pour développer des paramètres de champs de forces polarisables applicables à une large gamme de liquides ioniques.

Abstract

Transition metal dichalcogenides (TMDCs) are 2D materials that have not yet revealed their full potential. Their use in advanced electronic nanocomponents and in various fields shows their effectiveness. In the field of electronics and optoelectronics, TMDCs are considered important materials and will undoubtedly be a key element in the modernization of future technology. Their compatibility with alkylimidazolium cation-based ILs has allowed the formation of IL-TMDC interface-based devices, showing exceptional and unprecedented properties that are the focus of current cutting-edge research. However, to date, due to the complexity of these systems, few studies have been carried out to describe in detail the behavior of ILs at an atomic level when they are in contact with a substrate such as MoS 2 or used as electrolytes in FET devices. This thesis, organized into 6 chapters, aims to fill this knowledge gap by conducting a first comprehensive study. Chapter 1 provides a detailed account of the main physical and chemical characteristics of ILs and TMDC materials, as well as a complete state of their atomic modeling. Chapter 2 provides an overview of the theoretical context and the two methods used in this work, namely: classical MD and first-principles molecular dynamics(FPMD). The next three chapters present the research work dedicated to the modeling of ILs, TDMCs, and their IL-TMDC interfaces. On a first instance, in Chapter 3, the combination of classical MD simulations and FPMD simulations has allowed us to clarify their microscopic structure and specific chemical interactions. Specifically, we showed how the chemical structure and nature of the ions in these ILs influence their atomic structure, the nature and local topology of their H-bonding networks, with a special focus on electronic properties such as dipole moment distributions, Bader charges and electron density of states(DOS). Furthermore, the use of a non-polarizable force field allowed us to obtain reasonable agreement with the results of FPMD, highlighting the need for more sophisticated and complex force fields that can explicitly account for polarization effects to obtain a better quantitative agreement with FPMD. The results from this work are intended to be used to develop polarizable force field parameters applicable to a wide range of ILs. On a second instance, in Chapter 4, the combined approach of CMD and FPMD improved our fundamental understanding of the properties of ILs, particularly the different interactions occurring at the liquid/solid interface between a series of alkyl-imidazolium based ILs and defect-free and defective 2H-MoS 2 substrates. The results of this study can be useful for rationalizing the behavior of these types of interface systems in real-world applications and may potentially lead to the design of new ILs with improved properties for applications in next-generation electronics. In sum, the combined use of classical molecular simulation and first principles has been a valuable tool for understanding the properties and behaviors of ILs at the liquid/solid interface. On a third instance, in Chapter 5, the knowledge gained in previous chapters on the atomic modeling of the EMIM-BF 4 system in contact with charged MoS 2 surfaces was synthesized. The use of constant charge molecular dynamics (ccMD) provided valuable information on the interactions between this ionic liquid and charged surfaces. The results show how the effects of the electric potential influence the structuring and orientation of EMIM and BF 4 ions on charged MoS 2 surfaces, highlighting that electrostatic interactions are the ones that regulate the behavior and macroscopic response of these liquid/solid in-terfaces. This chapter also highlights the need to improve classical force fields by using the knowledge obtained from the FPMD approach to obtain a quantitative and complete description of the structural properties, chemical bonds, and electronic properties of these complex interfaces. From this point of view, the results of Chapter 5 can reveal new perspectives for the design of new ionic liquids and the realization of ILs-TMDC interfaces adapted for future nanoelectronics and optoelectronics. Finally, the Chapter 6 outline the overall conclusion of this work and identifies possible line of future research for this topic.

Keywords : ionic liquids, 2D material, TMDC, MoS 2 , EDL-FETS, molecular dynamics, first-principles, DFT, CMD, FPMD, CPMD, ccMD, modelling

1. 1 a

 1 ) Device characteristics and superconductivity in an IL-gated MoS 2 FET transistor. b) Phase diagram showing the evolution of electronic phases as a function of carrier density for MoS 2 gated with an IL . . . . . . . . . . . . 1.2 Left: Schematic overview of the wide variety of possible ILs used in synthetic chemistry as solvents. Right: A few examples of the most common AI-IL cations and anions used and their abbreviations. . . . . . . . . . . . . . . . 1.3 Schematic representation of the different type of interactions acting in imidazolium-based IL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 The transition metals and the three chalcogen elements that predominantly crystallize in MX 2 layered structure are highlighted in the periodic table of elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 a) Trigonal prismatic coordination of molybdenum in a MoS 2 sheet. b) Side view of two sheets of MoS 2 . c) Top view of two MoS 2 sheets. d) Hexagonal unit cell of MoS 2 in the 2H polymorphic phase . . . . . . . . . . . . . . . . 1.6 A) MoS 2 band structure computed at DFT level. Energy dispersion (energy versus wavevector k) in bulk, quadrilayer (4L), bilayer (2L) and monolayer (1L) from left to right. The horizontal dashed line represents the energy of a band maximum at the k point. The red and blue lines represent the conduction and valence band edges, respectively. The lowest energy transition are indicated by the solid arrows. B) Photoluminescent of suspended monolayer and bilayer MoS 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7 Dielectric gating. a) Random orientation of permanent dipoles when no bias is applied. b) Dipole moment alignment with applying voltage. c) Schematic diagram of the dielectric gating device with both back (SiO 2 ) and top HfO 2 gates on 2D material . . . . . . . . . . . . . . . . . . . . . . . 1.8 Electrolyte gating: (a) Schematic illustration of the EDL formation process and corresponding charge distributions. 2D mobility in FETs versus EDLTs 2.1 OPLS-AA and CL&P force fields potential form and contributions. . . . . . 2.2 Lennard-Jones interatomic potential. . . . . . . . . . . . . . . . . . . . . . . 2.3 Blocks diagram used to develop the CL&Pol force field. . . . . . . . . . . . 2.4 Schematic illustration of all-electron (solid lines) and pseudo-electron (dashed lines) potentials and their corresponding wave functions. The radius at which all-electron and pseudo-electron values match is designated r. . . . . 2.5 Heat transfer as a function of the pre-fixed average fictitious kinetic energy E kin e of the electronic wave functions for solid aluminum at the melting point. The arrow indicates the kinetic energy required for adiabatic motion of the electrons according to equation. (2.84) . . . . . . . . . . . . . . . . .2.6 Left: Total neutron structure factor of liquid GeSe 4 . The experimental data are those of Ref. Haye (1998) stands for Ref. Right: Snapshots of configurations obtained for liquid GeSe 4 at T = 1073 K by using the four different treatments for the inclusion of the dispersion forces detailed in the text. From top to bottom and from left to right: No-vdW (top left), MLWF (top right), D2 (bottom left), D3(A) (bottom centre), D3(B) (bottom right). The Ge atoms are in dark green, the Se atoms in brown, and the unit cell is represented with a black line. . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 A) Chemical structure of the ionic liquid cations and anions simulated in this chapter: as cations 1-methyl(DMIM), 1-ethyl (EMIM), 1-butyl (BMIM) and 1-hexyl-3-methylimidazolium (HMIM) and, as anions chlorine (Cl), tetrafluoro-borate (BF 4 ), hexafluoro-phosphate (PF 6 ) and bis(trifluoromethylsulfonyl)imide(TFSI). B) Example of structure and related atom indexing for 1-ethyl-3-methylimidazolium [EMIM] cation. C) Typical molecular configurations of 200 IL pairs simulated by CMD. The blue spheres are the EMIM whereas green, orange , yellow and purple correspond to Cl, BF 4 , PF 6 and TFSI anions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Upper panel: Computed stress tensor as a function of the number of plane waves of the basis set (PW cut-off) for the [DMIM] + [Cl] -system. The black dots refer to simulations performed at the experimental volume (density) of the IL, blue points refer to simulations done within the identical protocol at a volume increased of 2.5 % with respect to the experimental value, purple point to an increase of 5.0 % and orange to an increase of 7.5 %. Continuous lines are intended only as a guide to the eye. Lower panel: Kohn-Sham potential energy as a function of time for [DMIM] + BF - 4 simulated at 300 K. Inset: Comparison of the temperature of the ions and electrons along the simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Cation-anion center of mass pair correlation functions (g C/A (r)) for DMIM-Cl, -BF 4 , -PF 6 , -TFSI systems simulated by classical MD at 300K. We compare the data computed for 20 and 200 IL pairs models. . . . . . . . . . 3.4 Cation-anion center of mass pair correlation functions (g C/A (r)) for the ILs simulated by classical MD at 300 K (200 IL pairs). The PCFs are shifted. . 3.5 Cation-anion center of mass pair correlation functions (g C/A (r)) for ILs simulated by classical MD and FPMD at 300K. . . . . . . . . . . . . . . . . 3.6 Left panel: Pair correlation functions g ij (r) of the different H atoms of EMIM + with the F atoms of BF - 4 obtained by classical MD at room temperature. Inset: a zoom-in in the 1-6 Å region. Right panel: The atom numbering is the standard one given in the literature and reported explicitly here. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Pair correlation functions of the H atoms constituting the IL cations forming H-bonds with the various anions considered in the present study obtained by CMD at room temperature. Inset: a zoom-in in the 1-6 Å region. . . . . 3.8 Pair correlation functions g ij (r) of the H atom (labeled as H2) shared by the two N atoms of the IL cation and the Cl, F and O atoms atoms of the anion. The data is based on the MD simulations at 300K and is plotted for each anion as a function of the length of the cation alkyl chain. Inset: a zoom-in in the 1-6 Å region. . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Pair correlation functions g ij (r) of the H atom (labeled as H2) shared by the two N atoms of the IL cation and the Cl, F and O atoms atoms of the anion. The data are obtained by classical MD and FPMD at 300K. . . . . . 3.10 Pair correlation functions g H 2 -F (r) obtained by FPMD, the classical MD with nonpolarisable (CL&P) and with polarisable (CL&Pol) FF. . . . . . . 3.11 Top panels: Combined distribution function showing the hydrogen bond geometry between the atom H 2 of the IL cation and the F atoms of of the BF 4 anion as a function of the length of the cation alkyl chain (DMIM + EMIM + , BMIM + and HMIM + ) obtained by classical MD at room temperature. Lower left panels: Combined distribution function showing the hydrogen bond obtained by FPMD. Last panel on the bottom right: Labeling of the main atoms quoted in the main text and vectors defining the angle for which distribution density plots are shown. . . . . . . . . . . . . . 3.12 Combined distribution function showing the hydrogen bond geometry between the atom H2 of the DMIM and the Cl, F and O atoms of the anions obtained by classical MD at room temperature. . . . . . . . . . . . . . . . . 3.13 Computed electronic density of states (EDOS) for DMIM-Cl (computed with BLYP and B3LYP) and EMIM-BF 4 (computed with BLYP) at 300 K. 3.14 Wannier function centers (WFCs) representation of the chemical bonds for the different cations and anions of the ILs simulated in this work. . . . . . . 3.15 Averaged dipole moment distributions of the IL pairs inside the different ionic liquid obtained by FPMD at room temperature . . . . . . . . . . . . . 4.1 Left: Initial stages of the formation process simulated by classical MD of the interface made of 500 pairs of EMIM-BF 4 in contact with a MoS 2 bilayer. Right: typical final EMIM-BF 4 @MoS 2 interface structural model equilibrated at 300 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 EMIM-BF 4 ions COM density profiles with respect to the distance from the MoS 2 pristine surface obtained by classical MD at 300K. Inset: zoomin within the 0-2 nm region. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 A) PCFs g(r) H X -S between the hydrogen atoms of EMIM with respect to the sulfur atoms of the MoS 2 surface. B) PCFs g(r) H X -F between the hydrogen atoms of EMIM with respect to the fluorine atoms of the BF 4 . C,D) Snapshots showing the lateral (C) and top (D) view of the close contact between the IL ions and the MoS 2 layers. For the top view, only the IL ions of the first layer are shown. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Distribution of the order parameter S ring (θ r ) of the cation ring obtained at room-temperature. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 EMIM-BF 4 ions charge density profiles with respect to the distance from the MoS 2 pristine surface obtained by classical MD at 300K. Inset: zoom-in within the 0-2 nm region. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Left: schematic representation of the three different regions defined to analysed the MSD of IL. Right: Mean squared displacement of the COM of EMIM and BF 4 as a function of time for within the three different L1-3 regions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Residence time R(t) of the [EMIM + ][BF - 4 ] at different temperature 300K, 350K and 400K in the first layers of ILs. . . . . . . . . . . . . . . . . . . . . 4.8 Top view of the pristine MoS 2 bilayer (without defects) and two of the three models with surface defects (here shown the two with 4.0×1014 sites/cm 2 and 8.0×10 14 sites/cm 2 , respectively). . . . . . . . . . . . . . . . . . . . . . 4.9 Ion density and ion charge density profiles of EMIM-BF 4 with respect to the distance from the MoS 2 surface by considering three different surface defects densities compared to the pristine one (top). . . . . . . . . . . . . . 4.10 Top (left) and lateral (right) views of the saturation of the surface S vacancies by BF 4 ions for the MoS 2 model with a surface defect density of 4.0×10 14 sites/cm 2 . In the top view, the locations of the defect sites are highlighted with a orange circle. Bottom: (left) typical configuration of a BF 4 ion that is located on top of a S surface vacancy and (right) a model where the F• • • Mo and F• • • S interatomic distances are reported in Å. Color legend: EMIM ions in blue; BF 4 in orange; S atoms in yellow and Mo atoms in cyan. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.11 3D mean-squared displacement of COM of IL cation EMIM and anion BF 4 as a function of time in the IL layer close to the solid MoS 2 surface. Four models are compared; without surface density (solid line), and the three models with increasing surface defects densities: 10 defects (dashed line, 20 defects (dot line), and 50 defects (dashdot line). . . . . . . . . . . . . . . . . 4.12 Computed Kohn-Sham total energy and stress tensor values as a function of a and b % variation of the supercell of 2H-MoS 2 bilayer. . . . . . . . . . . 4.13 Zoom in the lateral and top views of the local environment of the 2H-MoS 2 cell optimized by FPMD in this work, It is shown the typical arrangement of Wannier centers (WFC) around a S atom (S 1) and Mo atom (Mo 1 ). Mo and S atoms appear as cyan and yelwalow spheres, respectively whereas the WFC appear in purple. Both the lone pairs electrons of S and Mo are shown (WFC SLP and WFC MLP , respectively) as well as the bond pairs electrons shared along the Mo-S bond (WFC MSB ). . . . . . . . . . . . . . . . . . . . . 4.14 Computed Kohn-Sham potential energy as a function of time for EMIM-BF 4 @MoS 2 interface simulated at 300 K. Inset: Comparison of the temperature of the ions and electrons along the simulation. . . . . . . . . . . . . . 4.15 Left: Interlayer distance between the top and bottom MoS 2 layers (i.e. vdW gap) as a function of time, calculated from the difference of the z coordinates of the S atoms of the bottom layer of MoS 2 at the interface (mobile layer) and the z coordinates of the S atoms of the top layer of the lower MoS 2 layer (frozen). Errors bars correspond to the 95% confidence interval obtained from averaging over all the S atoms of each layer. Right: xy displacement of S and Mo atoms of the top MoS 2 layer along the >10 ps of FPMD simulation at 300 K. . . . . . . . . . . . . . . . . . . . . . . . . . 4.16 Temporal evolution snapshots of the simulation of EMIM-BF 4 @MoS 2 interface at 300K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.17 EMIM-BF 4 ions COM density profiles with respect to the distance from the MoS 2 pristine surface obtained by FPMD at 300K (z=0 at the position of S atoms of the top layer of the top MoS 2 ). Top: total density profiles. Bottom: Density profiles of EMIM (left) and BF 4 (right) ions. For comparison, the density profiles obtained by classical MD at 300K are also shown. . . . . . . 4.18 Distributions of the order parameter S ring (θ r ) of the EMIM ring computed at 300K. The FPMD data is compared with those obtained by classical MD. 4.19 EMIM-BF 4 ions charge density distributions with respect to the distance from the MoS 2 surface computed by FPMD at 300K. Top: total charge density profiles. Bottom: Charge density profiles of EMIM (left) and BF 4 (right) ions. For comparison, the charge density profiles obtained by classical MD at 300 K are also shown. . . . . . . . . . . . . . . . . . . . . . . . 4.20 Electron density distributions of the EMIM-BF 4 @MoS 2 interfaces made with the pristine MoS 2 surface (left, computed at ∼9.5 ps) and with the MoS 2 with two S vacancies in the top layer (right, computed at ∼5.1 ps) . The electron density is shown at an isosurface value of 0.015 e/ a.u. 3 . Color code is cyan for Mo, yellow for S, organge fo BF 4 ions and blue for EMIM ions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.21 Left: Snapshot of the EMIM-BF 4 @MoS 2 interface with two S vacancies on the top layer at the end of the first stage of constrained FPMD run (T∼0K). Right: Lateral and side views of the interface model where only the two BF 4 ions interacting with the S vacancies sites are shown. . . . . . 4.22 Sequence of snapshots showing the temporal evolution along the ∼6ps of NVT FPMD run of the local rearrangements of the two BF 4 ions at the top layer of MoS 2 surface: initial stage where the two anions are interacting with the two S vacancies (at ∼0.1ps); detachment of one of the anions from the S vacancy site (at ∼2ps); detachment of the second anions (at ∼4ps); re-formation of a close interaction between the second BF 4 anion (after rotation) and the S vacancy site (at ∼5.5ps) while the first anion is still moving within the IL phase. Right: interatomic distances of the F atom of BF 4 at close interaction with the S vacancies site. The closest interatomic distances between the F atom and the three Mo and six S atoms surrounding the vacancies site are reported. . . . . . . . . . . . . . . . . . . . . . . . . . 4.23 Snapshots of the interface models simulated by classical MD and bearing different AI-IL anions: EMIM-Cl, EMIM-BF 4 , EMIM-TFSI. The chemical structure of each IL ions as well as the thickness of the each IL phase made of 500 ions pairs and equilibrated at 300K are also reported. . . . . . . . . . 4.24 EMIM-Cl, EMIM-BF 4 , EMIM-TFSI ion density and ion charge density profiles with respect to the distance from the MoS 2 surface. These results are obtained from the classical MD at 300K. . . . . . . . . . . . . . . . . . . 4.25 DMIM-TFSI, EMIM-TFSI, HMIM-TFSI ions density and ions charge density profiles with respect to the distance from the MoS 2 surface. These results are obtained from the classical MD at 300K. . . . . . . . . . . . . . 5.1 Lateral view of snapshots of the three models equilibrated at 400K and with 0.0, 36.6 and 73.1 µC/cm 2 surface charge densities. . . . . . . . . . . . . . . 5.2 EMIM-BF 4 ions density profiles as a function of the z direction between the two MoS 2 electrode at the condition of zero-charge. The z positions of the uppermost layer of S atoms of the two MoS 2 bilayer are 0.0 nm and 7.0 nm. 5.3 Zoom-in in 0.0-2.0 nm and 5.0-7.0 nm regions of the EMIM-BF 4 ions density profiles as a function of the z direction between the two MoS 2 electrodes at different surface charging. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Distributions of the order parameter S ring (θ r ) of the EMIM ring calculated as a function of the z direction between the two MoS 2 electrodes and at different surface charging. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Snapshots of EMIM and BF 4 ions orientation within the first interfacial layer obtained in the following charged models: majority of EMIM ions lying parallel (a) and perpendicular (b) to the MoS 2 charged surfaces at 3 µC/cm 2 and 73 µC/cm 2 ; c) BF 4 ions structured orientation found with 73 µC/cm 2 surface charge. We show the well organized orientation of BF 4 ions with three F atoms closely interacting with three S atom so the MoS 2 surface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 EMIM-BF 4 ion charge density profiles with respect to the distance from each electrodes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 Zoom-in in 0.0-2.0 nm and 5.0-7.0 nm regions of the EMIM-BF 4 charge density profiles as a function of the z direction between the two MoS 2 electrodes at different surface charging. . . . . . . . . . . . . . . . . . . . . . 5.8 FPMD computed electron density distributions of the EMIM-BF 4 @MoS 2

6. 1

 1 Pair correlation functions of H atoms constituting the IL cations forming H-bonds with the various anions considered in the present study obtained by classical MD at room temperature. Inset: a zoom-in in the 1-6 Å region. 6.2 Total (left) and individual (right) EMIM-BF 4 ions density profiles as a function of the z direction between the two MoS 2 electrodes at different surface charging. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Total (left) and individual (right) EMIM-BF 4 charge density profiles as a function of the z direction between the two MoS 2 electrodes at different surface charging. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Supraconductivité induite par le champ électrique dans les cristaux 2D. a) représentation schématique d'un appareil EDLT-FET. b) Diagramme de phase montrant l'évolution des phases électroniques en fonction de la densité de porteurs pour MoS 2 . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5 Représentation schématique des différents types d'interactions agissant dans les IL à base d'imidazolium. . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.6 Fonctions de corrélation de paires g H 2 -F (r) obtenues par FPMD et par MD classique avec un champ de force non polarisable (CL&P) et un champ de force polarisable (CL&Pol). . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.7 Line supérieure: Fonction de distribution combinée montrant la géom-étrie

3. 1

 1 Box side (L, Å), number of atoms, density (d, g.cm -3 ) and the minimum of the g ij (r) after the first peak used as cutoff for the calculation of n c of the IL models simulated by MD and FPMD. Both the starting and final values of density used in the MD simulations are reported. r c is the distance value used as cutoff for the integration of the first g ij (r) peak for calculating the coordination number n c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Cation-anion coordination numbers n c obtained from classical MD and FPMD simulations at 300 K. Where the values inside brackets represent the relative errors in %. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Position of the first coordination peak H2-anion for all models CL&P, FPMD and CL&Pol for 20 and 200 IL pairs (n). . . . . . . . . . . . . . . . 3.4 Computed Bader charges of the atoms of the DMIM-Cl and EMIM-BF 4 .

Figure 1 .

 1 Figure 1.1. a) Device characteristics and superconductivity in an IL-gated MoS 2 transistor. b) Phase diagram showing the evolution of electronic phases as a function of carrier density for MoS 2 gated with an IL[1].

Figure 1 . 2 .

 12 Figure 1.2. Left: Schematic overview of possible ILs used in synthetic chemistry as solvents. Right: A few examples of the most common AI-IL cations and anions used and their abbreviations. Adapted from[START_REF] Plechkova | Applications of ionic liquids in the chemical industry[END_REF] and[START_REF] Zhang | Nanoconfined ionic liquids[END_REF].

Figure 1 . 3 .

 13 Figure 1.3. Schematic representation of the different type of interactions acting in imidazolium-based IL. Adapted from [13].

Figure 1 . 4 .

 14 Figure1.4. The transition metals and the three chalcogen elements that predominantly crystallize in MX 2 layered structure are highlighted in the periodic table of elements[START_REF] Chhowalla | The chemistry of two-dimensional layered transition metal dichalcogenide nanosheets[END_REF].

Figure 1 .

 1 Figure 1.5. a) Trigonal prismatic coordination of molybdenum in a MoS 2 sheet. b) Side view of two sheets of MoS 2 . c) Top view of two MoS 2 sheets. d) Hexagonal unit cell of MoS 2 in the 2H polymorphic phase[40].

Figure 1 . 7 .

 17 Figure 1.7. Dielectric gating. a) Random orientation of permanent dipoles when no bias is applied. b) Dipole moment alignment with applying voltage. c) Schematic diagram of the dielectric gating device with both back (SiO 2 ) and top HfO 2 gates on 2D material[55].

Figure 1 . 8 .

 18 Figure 1.8. Electrolyte gating: (a) Schematic illustration of the EDL formation process and corresponding charge distributions. 2D mobility in FETs versus EDLTs[57, 58].

Figure 2

 2

Figure 2 . 1 .

 21 Figure 2.1. OPLS-AA and CL&P force fields potential form and contributions. Adapted from [73].

Figure 2 . 2 .

 22 Figure 2.2. Lennard-Jones interatomic potential.

Figure 2 . 3 .

 23 Figure 2.3. Blocks diagram used to develop the CL&Pol force field[156].

Figure 2 . 5 .

 25 Figure 2.5. Heat transfer as a function of the pre-fixed average fictitious kinetic energy E kin e of the electronic wave functions for solid aluminum at the melting point. The arrow indicates the kinetic energy required for adiabatic motion of the electrons according to equation (2.84)[203].

Figure 2 . 6 .

 26 Figure 2.6. Left: Total neutron structure factor of liquid GeSe 4 . The experimental data are those of Ref. [215]. Haye (1998) stands for Ref. [216]. Right: Snapshots of configurations obtained for liquid GeSe 4 at T = 1073 K by using the four different treatments for the inclusion of the dispersion forces detailed in the text. From top to bottom and from left to right: No-vdW (top left), MLWF (top right), D2 (bottom left), D3(A) (bottom centre), D3(B) (bottom right). The Ge atoms are in dark green, the Se atoms in brown, and the unit cell is represented with a black line[217].

Figure 3 .

 3 Figure 3.1. A) Chemical structure of the ionic liquid cations and anions simulated in this chapter: as cations 1-methyl(DMIM), 1-ethyl (EMIM), 1-butyl (BMIM) and 1-hexyl-3-methylimidazolium (HMIM) and, as anions chlorine (Cl), tetrafluoro-borate (BF 4 ), hexafluoro-phosphate (PF 6 ) and bis(trifluoromethyl-sulfonyl)imide(TFSI). B) Example of structure and related atom indexing for 1-ethyl-3-methylimidazolium [EMIM] cation. C) Typical molecular configurations of 200 IL pairs simulated by CMD. The blue spheres are the EMIM whereas green, orange , yellow and purple correspond to Cl, BF 4 , PF 6 and TFSI anions.

  [DMIM][Cl], 34 ps for [DMIM][BF 4 ], 27 ps for [DMIM][PF 6 ]), 15 ps for [DMIM][TFSI] and 20 ps for [EMIM][BF 4 ]. As an example, in Figure 3.2, the profile of the Kohn-Sham (KS) potential energy along the simulation at 300 K for [DMIM][BF 4

Figure 3 . 2 .

 32 Figure 3.2. Upper panel: Computed stress tensor as a function of the number of plane waves of the basis set (PW cut-off) for the [DMIM] + [Cl] -system. The black dots refer to simulations performed at the experimental volume (density) of the IL, blue points refer to simulations done within the identical protocol at a volume increased of 2.5 % with respect to the experimental value, purple point to an increase of 5.0 % and orange to an increase of 7.5 %. Continuous lines are intended only as a guide to the eye. Lower panel: Kohn-Sham potential energy as a function of time for [DMIM] + BF - 4 simulated at 300 K. Inset: Comparison of the temperature of the ions and electrons along the simulation.

Figure 3 . 3 .

 33 Figure 3.3. Cation-anion center of mass pair correlation functions (g C/A (r)) for DMIM-Cl, -BF 4 , -PF 6 , -TFSI systems simulated by classical MD at 300K. We compare the data computed for 20 and 200 IL pairs models.

Figure 3 . 4 .

 34 Figure 3.4. Cation-anion center of mass pair correlation functions (g C/A (r)) for the ILs simulated by classical MD at 300 K (200 IL pairs). The PCFs are shifted.

Figure 3 . 5 .

 35 Figure 3.5. Cation-anion center of mass pair correlation functions (g C/A (r)) for ILs simulated by classical MD and FPMD at 300K.

Figure 3 . 6 .Figure 3 . 7 .

 3637 Figure 3.6. Left panel: Pair correlation functions g ij (r) of the different H atoms of EMIM + with the F atoms of BF - 4 obtained by classical MD at room temperature. Inset: a zoom-in in the 1-6 Å region. Right panel: The atom numbering is the standard one given in the literature and reported explicitly here.

Figure 3 . 8 .

 38 Figure 3.8. Pair correlation functions g ij (r) of the H atom (labeled as H2) shared by the two N atoms of the IL cation and the Cl, F and O atoms atoms of the anion. The data is based on the MD simulations at 300K and is plotted for each anion as a function of the length of the cation alkyl chain. Inset: a zoom-in in the 1-6 Å region.

Figure 3 . 9 .

 39 Figure 3.9. Pair correlation functions g ij (r) of the H atom (labeled as H2) shared by the two N atoms of the IL cation and the Cl, F and O atoms atoms of the anion. The data are obtained by classical MD and FPMD at 300K.

Figure 3 . 10 .

 310 Figure 3.10. Pair correlation functions g H 2 -F (r) obtained by FPMD, the classical MD with nonpolarisable (CL&P) and with polarisable (CL&Pol) FF.

Figure 3 .Figure 3 . 11 .

 3311 Figure 3.11. Top panels: Combined distribution function showing the hydrogen bond geometry between the atom H 2 of the IL cation and the F atoms of of the BF 4 anion as a function of the length of the cation alkyl chain (DMIM + EMIM + , BMIM + and HMIM + ) obtained by classical MD at room temperature. Lower left panels: Combined distribution function showing the hydrogen bond obtained by FPMD. Last panel on the bottom right: Labeling of the main atoms quoted in the main text and vectors defining the angle for which distribution density plots are shown.

  Figure 3.12. Combined distribution function showing the hydrogen bond geometry between the atom H2 of the DMIM and the Cl, F and O atoms of the anions obtained by classical MD at room temperature.

Figure 3 . 13 .

 313 Figure 3.13. Computed electronic density of states (EDOS) for DMIM-Cl (computed with BLYP and B3LYP) and EMIM-BF 4 (computed with BLYP) at 300 K.

  Figure 3.14. Wannier function centers (WFCs) representation of the chemical bonds for the different cations and anions of the ILs simulated in this work.

Figure 3 . 15 .

 315 Figure 3.15. Averaged dipole moment distributions of the IL pairs inside the different ionic liquid obtained by FPMD at room temperature

Figure 4 . 1 .

 41 Figure 4.1. Left: Initial stages of the formation process simulated by classical MD of the interface made of 500 pairs of EMIM-BF 4 in contact with a MoS 2 bilayer. Right: typical final EMIM-BF 4 @MoS 2 interface structural model equilibrated at 300 K.

Figure 4 .

 4 2 shows the distribution of the centers of masses (COM) of the EMIM cations and BF 4 anions together with the total ions density profile along the c direction. For this EMIM-BF 4 @MoS 2 interface, significant layering was observed as the ions density profiles for the IL anion and cation exhibited relatively large density oscillations. Such significant layering, which is characteristic of liquids at interfaces and especially confined liquids, has been already observed in MD simulations of IL at interfaces with crystalline materials, such as graphene and mica, as well as amorphous materials (carbon, silica and GeS 2 ). Because of the atomic charges present at the MoS 2 surfaces, the electrostatic potential plays a crucial role in the structure of the IL at the liquid-solid interface. By looking at the total density profiles of the ions COM, three distinct layers can be identified by the distinct maxima of the corresponding distributions, with the first peak showing a small shoulder. The first peak is relatively sharp and characterized by high intensity, the second and following peaks show decreasing intensities as long as the distance from the solid MoS 2 substrate increases. By disentangling the total profile into the contributions from the IL cation and anion, we can identify three peaks for EMIM located at 0.35, 0.72 and 1.2 nm from the top layer of S atoms of MoS 2 , with the first peak very narrow whereas the second more broadened and the third barely present. While BF 4 is showing a first peak very broadened in the range 0.3-0.6 nm without a well-defined maximum, followed by a second peak at 1.09 nm and a third one, barely observable, at 1.6 nm. The large first density peak of BF 4 is responsible of the shoulder of the first peak of the ions density profile. The density peaks of EMIM and BF 4 in part alternate and in part overlaps, suggesting a certain degree of interdigitation between the IL cation and anion layers.

Figure 4 . 2 .

 42 Figure 4.2. EMIM-BF 4 ions COM density profiles with respect to the distance from the MoS 2 pristine surface obtained by classical MD at 300K. Inset: zoom-in within the 0-2 nm region.

Figure 4 .

 4 Figure 4.3. A) PCFs g(r) H X -S between the hydrogen atoms of EMIM with respect to the sulfur atoms of the MoS 2 surface. B) PCFs g(r) H X -F between the hydrogen atoms of EMIM with respect to the fluorine atoms of the BF 4 . C,D) Snapshots showing the lateral (C) and top (D) view of the close contact between the IL ions and the MoS 2 layers. For the top view, only the IL ions of the first layer are shown.

Figure 4 . 4 .

 44 Figure 4.4. Distribution of the order parameter S ring (θ r ) of the cation ring obtained at room-temperature.

Figure 4 . 5 .

 45 Figure 4.5. EMIM-BF 4 ions charge density profiles with respect to the distance from the MoS 2 pristine surface obtained by classical MD at 300K. Inset: zoom-in within the 0-2 nm region.

Figure 4 .

 4 [START_REF] Dubonos | Electric field effect in atomically thin carbon[END_REF] shows a schematic representation of the locations of the three different regions with respect of the MoS 2 bilayer and the calculated MSD profiles as a function of time.

Figure 4 . 6 .

 46 Figure 4.6. Left: schematic representation of the three different regions defined to analysed the MSD of IL. Right: Mean squared displacement of the COM of EMIM and BF 4 as a function of time for within the three different L1-3 regions.

Figure 4 . 7 .

 47 Figure 4.7. Residence time R(t) of the [EMIM + ][BF - 4 ] at different temperature 300K, 350K and 400K in the first layers of ILs.

Figure 4 . 8 .

 48 Figure 4.8. Top view of the pristine MoS 2 bilayer (without defects) and two of the three models with surface defects (here shown the two with 4.0×10 14 sites/cm 2 and 8.0×10 14 sites/cm 2 , respectively).

Figure 4 . 9 .

 49 Figure 4.9. Ion density and ion charge density profiles of EMIM-BF 4 with respect to the distance from the MoS 2 surface by considering three different surface defects densities compared to the pristine one (top).

Figure 4 .

 4 Figure 4.10. Top (left) and lateral (right) views of the saturation of the surface S vacancies by BF 4 ions for the MoS 2 model with a surface defect density of 4.0×10 14 sites/cm 2 . In the top view, the locations of the defect sites are highlighted with a orange circle. Bottom: (left) typical configuration of a BF 4 ion that is located on top of a S surface vacancy and (right) a model where the F• • • Mo and F• • • S interatomic distances are reported in Å. Color legend: EMIM ions in blue; BF 4 in orange; S atoms in yellow and Mo atoms in cyan.

Figure 4 .

 4 Figure 4.11. 3D mean-squared displacement of COM of IL cation EMIM and anion BF 4 as a function of time in the IL layer close to the solid MoS 2 surface. Four models are compared; without surface density (solid line), and the three models with increasing surface defects densities: 10 defects (dashed line, 20 defects (dot line), and 50 defects (dashdot line).

Figure 4 .

 4 Figure 4.12. Computed Kohn-Sham total energy and stress tensor values as a function of a and b % variation of the supercell of 2H-MoS 2 bilayer.

Figure 4 .

 4 13 shows the positions of the WFCs in a molecular configuration of the MoS 2 layer simulated. For S, the existence of six valence electrons and the threefold coordination of S is at the origin of a specific pattern; one WFC center is localized close to the S atom, representing a lone pair of electrons not involved in chemical bonding. In what follows, these WFCs are referred to as WFC SLP centers. The other three WFCs, which are localized along the Mo-S bonds, reflect the interatomic bonding. In what follows, these are referred to as WFC MSB . The sixfold Mo atoms are characterized by 14 electrons in total, as the sum of the valence and core electrons due to the type of semicore pseudopotential used for Mo. The Mo 1 atom in Figure 4.13 is coordinated by six S atoms and it shows several WFCs very close to the atom site, corresponding the its WFC SLP and core electrons. The S 1 atom in Figure 4.13 is coordinated by three Mo atoms, with bond lengths equal to ∼2.46 Å. The location of the three WFC MSB surrounding S 1 atom is found at a S-WFC MSB distance d S-WFC of ∼0.85 Å and at a distance Mo-WFC MSB distance d Mo-WFC of ∼1.53 Å. By considering the vdW and covalent radii of Mo and S (vdW radii: 2.25 Å , and 1.

Figure 4 . 13 .

 413 Figure 4.13. Zoom in the lateral and top views of the local environment of the 2H-MoS 2 cell optimized by FPMD in this work, It is shown the typical arrangement of Wannier centers (WFC) around a S atom (S 1 ) and Mo atom (Mo 1 ). Mo and S atoms appear as cyan and yelwalow spheres, respectively whereas the WFC appear in purple. Both the lone pairs electrons of S and Mo are shown (WFC SLP and WFC MLP , respectively) as well as the bond pairs electrons shared along the Mo-S bond (WFC MSB ).

Figure 4 .

 4 Figure 4.14. Computed Kohn-Sham potential energy as a function of time for EMIM-BF 4 @MoS 2 interface simulated at 300 K. Inset: Comparison of the temperature of the ions and electrons along the simulation.

Figure 4 .

 4 Figure 4.15 shows the temporal evolution of the interlayer distance between the top and bottom MoS 2 layers (i.e. vdW gap), calculated from the difference of the average z coordinates of the S atoms of the bottom layer of MoS 2 at the interface (mobile layer) and the z coordinates of the S atoms of the top layer of the lower MoS 2 layer (frozen). The constant average value of the vdW gap beyond ∼5 ps is found in good agreement with the experimental value reported for MoS 2 (3.22 Å versus 3.14 Å). The Mo-S bonds distances of the fully mobile top layer restored quantitatively the experimental value (2.37 Å). Moreover, possible sliding along the xy plane of the top MoS 2 layer over the bottom MoS 2 layer was found to be negligible. Indeed, only a minimal displacement (∼0.3 Å ) of S and Mo atoms along the xy plane was recorded during the entire >10 ps of FPMD (Figure 4.15) simulation.Figure 4.16 shows the temporal evolution of the EMIM-BF 4 @MoS 2 interface at 300K, visualized in snapshots of the models taken at 0.1, 4, 6, 8 and 10.5 ps.A small expansion along the z direction of the IL phase was recorded during the ini-

Figure 4 .

 4 Figure 4.15 shows the temporal evolution of the interlayer distance between the top and bottom MoS 2 layers (i.e. vdW gap), calculated from the difference of the average z coordinates of the S atoms of the bottom layer of MoS 2 at the interface (mobile layer) and the z coordinates of the S atoms of the top layer of the lower MoS 2 layer (frozen). The constant average value of the vdW gap beyond ∼5 ps is found in good agreement with the experimental value reported for MoS 2 (3.22 Å versus 3.14 Å). The Mo-S bonds distances of the fully mobile top layer restored quantitatively the experimental value (2.37 Å). Moreover, possible sliding along the xy plane of the top MoS 2 layer over the bottom MoS 2 layer was found to be negligible. Indeed, only a minimal displacement (∼0.3 Å ) of S and Mo atoms along the xy plane was recorded during the entire >10 ps of FPMD (Figure 4.15) simulation.Figure 4.16 shows the temporal evolution of the EMIM-BF 4 @MoS 2 interface at 300K, visualized in snapshots of the models taken at 0.1, 4, 6, 8 and 10.5 ps.A small expansion along the z direction of the IL phase was recorded during the ini-

Figure 4 .

 4 Figure 4.15. Left: Interlayer distance between the top and bottom MoS 2 layers (i.e. vdW gap) as a function of time, calculated from the difference of the z coordinates of the S atoms of the bottom layer of MoS 2 at the interface (mobile layer) and the z coordinates of the S atoms of the top layer of the lower MoS 2 layer (frozen). Errors bars correspond to the 95% confidence interval obtained from averaging over all the S atoms of each layer. Right: xy displacement of S and Mo atoms of the top MoS 2 layer along the >10 ps of FPMD simulation at 300 K.

Figure 4 . 16 .

 416 Figure 4.16. Temporal evolution snapshots of the simulation of EMIM-BF 4 @MoS 2 interface at 300K.

Figure 4 .

 4 [START_REF] Hardacre | Liquid structure of 1, 3-dimethylimidazolium salts[END_REF] shows the distribution of the total ions density and its breakdown into the contribution of EMIM and BF 4 . Being the IL phase simulated by FPMD of 20 IL ions, it is of relevance to discuss the density profiles only for the first layers of IL from the MoS 2 surface (up to 1.8 nm) to avoid any possible influence of the interface between the IL and the upper vacuum.

Figure 4 .

 4 Figure 4.17. EMIM-BF 4 ions COM density profiles with respect to the distance from the MoS 2 pristine surface obtained by FPMD at 300K (z=0 at the position of S atoms of the top layer of the top MoS 2 ). Top: total density profiles. Bottom: Density profiles of EMIM (left) and BF 4 (right) ions. For comparison, the density profiles obtained by classical MD at 300K are also shown.

Figure 4 . 18 .

 418 Figure 4.18. Distributions of the order parameter S ring (θ r ) of the EMIM ring computed at 300K. The FPMD data is compared with those obtained by classical MD.

2 .

 2 The computed Bader charges of each atom of the EMIM-BF 4 @MoS 2 interface do not undergo any significant fluctuation all along the dynamics. The total Bader charges of EMIM and BF 4 are equal to +0.96 and -0.95, very close to the values obtained in the EMIM-BF 4 bulk phase. The MoS 2 bilayer shows Mo and S atoms with an averaged charge of +1.03 and -0.52, respectively. These value are slightly greater, in terms of absolute values, of the values computed in the 2H-MoS 2 bulk phase. Of notice that the charges values computed by FPMD are almost the double of the point charges used in the CL&P force field (+0.50 and -0.25). The availability of the Bader charges for each atoms allows us to determine the ions charge density profiles, by integration of the sum of the charges at a given distance from the MoS 2 surface.

Figure 4 .

 4 [START_REF] Giernoth | Task-specific ionic liquids[END_REF] shows the total charge distribution of EMIM-BF 4 as well as the individual contribution from EMIM and BF 4 . This analysis allows to make a direct comparison with the charge density profiles obtained by CMD, also reported in Figure4.19.

Figure 4 .

 4 Figure 4.19. EMIM-BF 4 ions charge density distributions with respect to the distance from the MoS 2 surface computed by FPMD at 300K. Top: total charge density profiles. Bottom: Charge density profiles of EMIM (left) and BF 4 (right) ions. For comparison, the charge density profiles obtained by classical MD at 300 K are also shown.

Figure 4 .

 4 21 shows a snapshot of this interface model as obtained at the end of this constrained FPMD run. It can be seen that two BF 4 ions are still located in the close proximity of the two S vacancies sites, however at this point the F• • • Mo and F• • • S interatomic distances (∼3.5 Å and ∼3.6 Å , respectively) already show slightly longer values compared to the values registered with classical MD.

Figure 4 .

 4 [START_REF] Jg | Characterization and comparison of hydrophilic and hydrophobic room temperature ionic liquids incorporating the imidazolium cation[END_REF] shows the temporal evolution along the ∼6ps of NVT FPMD run (heating up to 4ps and then equilibrating at 300K) of the local rearrangements of the two BF 4 ions at the top layer of MoS 2 surface.From the initial stage, where the two anions are interacting with the two S vacancies

Figure 4 . 20 .

 420 Figure 4.20. Electron density distributions of the EMIM-BF 4 @MoS 2 interfaces made with the pristine MoS 2 surface (left, computed at ∼9.5 ps) and with the MoS 2 with two S vacancies in the top layer (right, computed at ∼5.1 ps) . The electron density is shown at an isosurface value of 0.015 e/ a.u. 3 . Color code is cyan for Mo, yellow for S, organge fo BF 4 ions and blue for EMIM ions.

Figure 4 .

 4 Figure 4.21. Left: Snapshot of the EMIM-BF 4 @MoS 2 interface with two S vacancies on the top layer at the end of the first stage of constrained FPMD run (T∼0K). Right: Lateral and side views of the interface model where only the two BF 4 ions interacting with the S vacancies sites are shown.

Figure 4 . 22 . 4 . 4

 42244 Figure 4.22. Sequence of snapshots showing the temporal evolution along the ∼6ps of NVT FPMD run of the local rearrangements of the two BF 4 ions at the top layer of MoS 2 surface: initial stage where the two anions are interacting with the two S vacancies (at ∼0.1ps); detachment of one of the anions from the S vacancy site (at ∼2ps); detachment of the second anions (at ∼4ps); re-formation of a close interaction between the second BF 4 anion (after rotation) and the S vacancy site (at ∼5.5ps) while the first anion is still moving within the IL phase. Right: interatomic distances of the F atom of BF 4 at close interaction with the S vacancies site. The closest interatomic distances between the F atom and the three Mo and six S atoms surrounding the vacancies site are reported.

Figure 4 .

 4 Figure 4.23. Snapshots of the interface models simulated by classical MD and bearing different AI-IL anions: EMIM-Cl, EMIM-BF 4 , EMIM-TFSI. The chemical structure of each IL ions as well as the thickness of the each IL phase made of 500 ions pairs and equilibrated at 300K are also reported.

Figure 4 .

 4 Figure 4.24. EMIM-Cl, EMIM-BF 4 , EMIM-TFSI ion density and ion charge density profiles with respect to the distance from the MoS 2 surface. These results are obtained from the classical MD at 300K.

Figure 4 .

 4 Figure 4.25. DMIM-TFSI, EMIM-TFSI, HMIM-TFSI ions density and ions charge density profiles with respect to the distance from the MoS 2 surface. These results are obtained from the classical MD at 300K.

  3 package[251].

Figure 5 . 1 .

 51 Figure 5.1. Lateral view of snapshots of the three models equilibrated at 400K and with 0.0, 36.6 and 73.1 µC/cm 2 surface charge densities.

Figure 5 . 2 .

 52 Figure 5.2. EMIM-BF 4 ions density profiles as a function of the z direction between the two MoS 2 electrode at the condition of zero-charge. The z positions of the uppermost layer of S atoms of the two MoS 2 bilayer are 0.0 nm and 7.0 nm.

Figure 5 . 3 .

 53 Figure 5.3. Zoom-in in 0.0-2.0 nm and 5.0-7.0 nm regions of the EMIM-BF 4 ions density profiles as a function of the z direction between the two MoS 2 electrodes at different surface charging.

Figure 5 . 4 .

 54 Figure 5.4. Distributions of the order parameter S ring (θ r ) of the EMIM ring calculated as a function of the z direction between the two MoS 2 electrodes and at different surface charging.

Figure 5 . 5 .

 55 Figure 5.5. Snapshots of EMIM and BF 4 ions orientation within the first interfacial layer obtained in the following charged models: majority of EMIM ions lying parallel (a) and perpendicular (b) to the MoS 2 charged surfaces at 3 µC/cm 2 and 73 µC/cm 2 ; c) BF 4 ions structured orientation found with 73 µC/cm 2 surface charge. We show the well organized orientation of BF 4 ions with three F atoms closely interacting with three S atom so the MoS 2 surface.

Figure 5 .Figure 5 . 6 .

 556 Figure 5.6 shows the total charge distribution of EMIM-BF 4 at zero-charge state, as well as the individual contribution from EMIM and BF 4 . This profile is rather symmetrical and very similar to the one reported in Chapter 4, highlighting a negligible effect due to different point charges of S and Mo. Figure 6.3 in the Annex shows the charge profiles as a function of the MoS 2 surface charge (i.e. electrode voltage). By increasing the surface charge a clear sharpening of the initial charge peaks appears and the charge balance is reached at closer distance with respect to the MoS 2 surface in comparison with the ions density distributions. By inspecting the regions close to the electrodes surfaces (Figure

Figure 5 . 7 .

 57 Figure 5.7. Zoom-in in 0.0-2.0 nm and 5.0-7.0 nm regions of the EMIM-BF 4 charge density profiles as a function of the z direction between the two MoS 2 electrodes at different surface charging.

Figure 5 . 8 .

 58 Figure 5.8. FPMD computed electron density distributions of the EMIM-BF 4 @MoS 2 interface produced by classical MD with a surface charge of 73 µC/cm 2 . The electron density is shown at an isosurface value of 0.015 e/ a.u. 3 . Color code is cyan for Mo, yellow for S, organge fo BF 4 ions and blue for EMIM ions.

Figure 5 . 9 .

 59 Figure 5.9. EMIM-BF 4 ions density and charge density profiles with respect to the distance from each electrodes with surface defects of 10 and 20 S vacancies.

Figure 6 . 2 .

 62 Figure 6.2. Total (left) and individual (right) EMIM-BF 4 ions density profiles as a function of the z direction between the two MoS 2 electrodes at different surface charging.

Figure 6 . 3 .

 63 Figure 6.3. Total (left) and individual (right) EMIM-BF 4 charge density profiles as a function of the z direction between the two MoS 2 electrodes at different surface charging.

Figure 6 . 4 .

 64 Figure 6.4. Supraconductivité induite par le champ électrique dans les cristaux 2D. a) représentation schématique d'un appareil EDLT-FET. b) Diagramme de phase montrant l'évolution des phases électroniques en fonction de la densité de porteurs pour MoS 2 [1].

Figure 6 . 5 .

 65 Figure 6.5. Représentation schématique des différents types d'interactions agissant dans les IL à base d'imidazolium. Adapté de [13].

Figure 6 . 6 .

 66 Figure 6.6. Fonctions de corrélation de paires g H 2 -F (r) obtenues par FPMD et par MD classique avec un champ de force non polarisable (CL&P) et un champ de force polarisable (CL&Pol).

Figure 6 . 7 .

 67 Figure 6.7. Line supérieure: Fonction de distribution combinée montrant la géométrie des liaisons hydrogènes entre l'atome H 2 du cation du IL et les atomes F de l'anion BF - 4 en fonction de la longueur de la chaîne alkyle du cation (DMIM + EMIM + , BMIM + et HMIM + ) obtenue par MD classique à température ambiante. Line inférieure à gauche: fonction de distribution combinée montrant la liaison hydrogène obtenue par FPMD. À droite : étiquetage des principaux atomes cités dans le texte principal et des vecteurs définissant l'angle pour lequel les diagrammes de densité de distribution sont affichés.

Figure 6 . 8 .

 68 Figure 6.8. À gauche : étapes initiales par MD classique du processus de formation de l'interface constituée de 500 paires de EMIM-BF 4 en contact avec une bicouche de MoS 2 .A droite : modèle structural final typique de l'interface EMIM-BF 4 @MoS 2 équilibrée à 300 K.

Figure 6 . 9 .

 69 Figure 6.9. Distributions de la densité électronique des interfaces EMIM-BF 4 @MoS 2 réalisées avec la surface sans défauts de MoS 2 (à gauche) et avec deux lacunes S dans la couche supérieure de MoS 2 (à droite). Le code couleur est le cyan pour Mo, le jaune pour S, l'orange pour les ions BF 4 et le bleu pour les ions EMIM.

3 Modelling of ionic liquids in the bulk phase

  3.1 Context and objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Computational details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Classical molecular dynamics . . . . . . . . . . . . . . . . . . . . . . 3.2.2 First-principles molecular dynamics . . . . . . . . . . . . . . . . . . 3.3 Structural properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Hydrogen bonding signatures . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Electronic properties and dipole moments . . . . . . . . . . . . . . . . . . . 3.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  2 (à gauche) et avec deux lacunes S dans la couche supérieure de MoS 2 (à droite). Le code couleur est le cyan pour Mo, le jaune pour S, l'orange pour les ions BF 4 et le bleu pour les ions EMIM. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.10 Instantanés de l'orientation des ions EMIM et BF 4 dans la première couche Typical ranges of room-temperature ILs properties. Adapted form . . . . . 1.2 Bandgaps (eV) and characteristics of a few TMDC materials in the bulk and as monolayer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Comparison of key properties of 2H-MoS 2 according to experiments and computational methods. From left to right: method, lattice parameters, Mo-S bond distance (r 0 ), cleavage energy of the (001) basal plane (E

	List of Tables
	1.1

interfaciale obtenus dans les modèles chargés: la majorité des ions EMIM sont parallèles (a) et perpendiculaires (b) aux surfaces chargées de MoS 2 à 3µC/cm 2 et 73µC/cm 2 ; c) orientation structurée des ions BF 4 trouvée avec une charge de surface de 73µC/cm 2 . d) Organisation des ions BF 4 sur la surface de MoS 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . cleav ), major IR/Raman peak, in-plane Young' modulus (Y), bulk modulus (B), contact angle with liquids at 298 K, atom mobility, compatibility and transferability of the potential. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  Computed Bader charges of the Mo and S atoms of 2H-MoS 2 in the bulk phase and of the two EMIM-BF 4 @MoS 2 interfaces simulated with the pristine and defected MoS 2 bilayers. The charge values per atom types have been averaged over the total number of atoms of the same species. Note that for the defected MoS 2 bilayer also the charges of the surrounding Mo and S atoms around the two S vacancies have been reported (*). . . . . . . Je souhaiterais tout d'abord remercier le bon Dieu pour la santé et la persévérance qu'il m'a accordé tout au long de ce travail de thèse. Par la suite, je souhaiterais adresser un remerciement particulier et spécial à mon directeur thèse monsieur le professeur Mauro Boero pour avoir cru en moi et m'avoir donné l'opportunité d'effectuer des travaux de thèse à l'institut de physique et de chimie des matériaux de Strasbourg (IPCMS) et d'intégrer son équipe de recherche. Ses conseils, sa méthode scientifique et son éthique du travail ont

3.1B and 3.6. . . . . . . 4.1 Comparison of the bulk 2H-MoS 2 structural parameters computed in this work at the BLYP + vdW-W level (vdW-W stands for vdW dispersion forces described by the MLWF-based scheme) with the corresponding experimental data and also other theoretical works reported in literature. . . . 4.2 5.1 Point charges used for S and Mo atoms at the two MoS 2 electrodes in order to define the different surface charges distribution (and voltages). Each electrode is made of a bilayer of MoS 2 of 450 Mo atoms and 900 S atoms. . Remerciements été d'une très grande utilité dans l'accomplissement de ce travail, Je tiens à lui dire grandement merci pour son professionalisme, sa disponibilité et tous les échanges enrichissants que nous avons eus. Je lui dis également infiniment merci d'avoir créé une atmosphère sereine au sein de l'équipe et pour le temps accordé à la relecture et correction de mon manuscrit de thèse.

  Je souhaiterais remercier le responsable du programme QMat monsieur le Pr. Bernard Doudin et le responsable du master Matière condensée et nano physique monsieur le Pr. Mebarek Alouani de m'avoir accordé une bourse d'étude pour la poursuite de mes études à l'université de Strasbourg. Je souhaiterais remercier le directeur de l'IPCMS monsieur le Dr. Pierre Rabu, l'ancienne et le nouveau directeur du DCMI madame la Pr. Nathalie Viart et monsieur le Dr. Guillaume Rogez pour mon intégration rapide au sein du laboratoire et d'avoir facilité des échanges entre doctorants et permanents à travers différents séminaires scientifiques. Mes remerciements vont aussi à tous les membres du département du DCMI, pour toute leur bienveillance et leur sympathie à mon égard, ainsi que pour toute l'assistance et l'aide qu'ils n'ont cessé de témoigner à mon endroit durant ces quatre dernières années. Je remercie très vivement monsieurs les Dr. Donald Mouafo et Ulrich Noumbé d'avoir été là à mes débuts. Merci infiniment de l'aide et des différents conseils que vous avez pu me prodiguer, vous avez été comme deux grands frères pour moi, vous m'avez montré la bonne voie et aviez été deux modèles exemplaires, vraiment merci pour tout. Je souhaite témoigner ma reconnaissance à monsieur le Pr. Michael Moseler et madame la Dr. Kerstin Falk de l'université de Freiburg dans le cadre du programme IDEX PhD International, pour leur accompagnement et les précieux échanges scientifiques que nous avons eus tout au long de mes travaux de recherches. Leur expertise et leur soutien continu ont été des éléments clés dans l'aboutissement de ce travail et je leur suis très reconnaissant pour leur engagement. Je tiens également à remercier monsieur le Pr. Emanuele Orgiu de l'Université du Québec (Centre INRS-EMT) pour ses précieux conseils sur ce projet du point de vue de la contrepartie expérimentale. Je tiens à remercier Kana Ishisone, Steve Dave Wansi, Mohammed Guerboub, Cheick Diarra, Achille Lambrecht pour les nombreux échanges instructifs et les bons moments que nous avons passés ensemble, nous avons appris les uns des autres, et je suis reconnaissant pour leur amitié et leur soutien inébranlable. Je les remercie sincèrement pour leur contribution positive.

  + cation with Cl -, BF - 4 , PF - 6 and TFSI -anions have melting points close to 41 • C[START_REF] Jg | Characterization and comparison of hydrophilic and hydrophobic room temperature ionic liquids incorporating the imidazolium cation[END_REF], -81C •[START_REF] Jg | Characterization and comparison of hydrophilic and hydrophobic room temperature ionic liquids incorporating the imidazolium cation[END_REF], 10

• C

[START_REF] Fuller | The room-temperature ionic liquid 1ethyl-3-methylimidazolium tetrafluoroborate: electrochemical couples and physical properties[END_REF]

and -3 • C

[START_REF] Kalyanasundaram | Hydrophobic, highly conductive ambient-temperature molten salts[END_REF]

;

for BMIM

  Cl 7 ] -,[HSO 4 ] -[PR 4 ] + ,[OAc] -

				• C
	density		1.1-1.6 g.L -1	[HM pyr ][N(CN) 2 ]=0.92 [EMIM]Br-AlBr 3 =2.2 g.L -1
	viscosity		40-800 mPa.s	[EMIM]Cl-AlCl 3 =14 [BMIM]Cl-(supercooled)=40.89 Pa.s
	thermal stability	230-300 • C	[EMIM][OAc]=∼ 200 • C [EMIM][NTf 2 ] = 400 • C
	surface tension	30-50 nN.m -1	[C 12 MIM][PF 6 ]=23.6 [MMIM[MeSO 4 ]=59.8 mN.m -1
	heat capacity	300-400 J.mol -1K -1	[OMIM][NTf 2 ]=654 J.mol -1.K -1
	water miscibility	many ILs do mix with water but also can be extracted from water	[NTf2] -,[P(C 2 F 5 ) 3 F 3 ] -[RSO 3 ] -, [RSO 4 ] -, [R 2 PO 4 ] -
	hydrolytic	heterocyclic cations can	[BF4 ] -,[PF6 ] -[PR 4 ] + , [OAc] -
	stability		hydrolyze under extreme condi-tions
	base stability [Al 2 corrosion all 1,3-dialkylimidazolium ILs are subject to deprotonation many RTILs are corrosive ver-sus Cu; additives to inhibit cor-rosion are available [NTf 2 ] -, [OTf] -|Cl -, HF formed from [MF x ] -hydrolysis
	toxicity		toxicity often increases with in-creasing lipophilicity	cholinium][OAc]| [EMIM][CN]
	static constant	dielectric	10-15	[HMIm][PF 6 ]=8.9±0.9 (2-hydroxyethyl)ammonium lactate=85.6± 3
	electrical		1-10 mS.cm -1	[BMIm][C 4 F 9 SO 3 ]=0.45 [EMIm]-[BF4]=16.3±1.6 mS.cm -1
	conductivity

Table 1 .

 1

	2. Bandgaps (eV) and characteristics of a few TMDC materials
	in the bulk and as monolayer. 1
	Bulk	X = S	X = Se X = Te	Characteristics
	M = Mo	1.23 a	1.09 a	1.0 a	Semiconducting, indirect
	M = W	1.35 a	1.20 a	1.1 a,f Semiconducting, indirect
	Monolayer	X = S	X = Se X = Te	Characteristics
	M = Mo	1.8 a,b	1.64 e , 1.5 a	1.1 a	Semiconducting, direct
	M = W	2.10 a , 2.08 c	1.74 d	1.1 a	Semiconducting, direct

a Ref.

[START_REF] Qh | Electronics and optoelectronics of two-dimensional transition metal dichalcogenides[END_REF][START_REF] Kk | Detailed photocurrent spectroscopy of the semiconducting group vib transition metal dichalcogenides[END_REF]

; b Ref.

[START_REF] Kf | Atomically thin mos2: A new direct-gap semiconductor[END_REF]

; c Ref.

[START_REF] He | Layer-dependent modulation of tunsten disulfide photoluminescence by lateral electric fileds[END_REF]

; d Ref.

[START_REF] Arora | Exciton band structure in layered mose 2 : from monolayer to the bulk limit[END_REF]

; e Ref.

[START_REF] Tongay | Thermally driven crossover from indirect toward direct bandgap in 2d semiconductors : mose 2 versus mos 2[END_REF]

; f Ref.

[START_REF] Ar | The reflectivity spectra of some group va transition metal dichalcogenides[END_REF]

.

  In 2005, FPMD simulations on DMIM-Cl carried out by Del Pópolo et al. and Buhl et al. have been the first evidence of significant differences in the H-bond features with respect to both the classical simulations and neutron diffraction results[109, 110]. One year later, Bhargava et al. found fair agreement between the pair correlation functions obtained from classical MD and FPMD, yet with not-negligible differences

Table 2 .

 2 1. Comparison of key properties of 2H-MoS 2 according to experiments and computational methods. From left to right: method, lattice parameters, Mo-S bond distance (r 0 ), cleavage energy of the (001) basal plane (E c leav), major IR/Raman peak, in-plane Young' modulus (Y), bulk modulus (B), contact angle with liquids at 298K, atom mobility, compatibility and transferability of the potential.[START_REF] Liu | Interpretable molecular models for molybdenum disulfide and insight into selective peptide recognition[END_REF]

	Lattice Parametter

• ) Method a( Å) c( Å) r 0 ( Å) E cleav (mJ/m 2 ) (cm -1 ) Y(Gpa) B(GPa) H 2 O CH 2 I 2 All atoms Applicable

Table 3 .

 3

	2. Cation-anion coordination numbers n c obtained from classical MD and
	FPMD simulations at 300 K. Where the values inside brackets represent the relative
	errors in %.
			C x -Cl	C x -BF 4	C x -PF 6	C x -TFSI
	DMIM	CMD	6.69(3)	7.46(3)	7.92(3)	10.11(2)
		FPMD	7.01(4)	7.44(2)	7.12(3)	10.02(2)
	EMIM	CMD	6.99(2)	7.48(3)	7.99(2)	7.51(3)
		FPMD	-	7.32(3)	-	-
	BMIM	CMD	6.10(3)	6.91(2)	7.17(4)	7.30(2)
	HMIM	CMD	6.77(2)	5.93(3)	6.57(3)	6.36(2)

Table 3 .

 3 3. Position of the first coordination peak H2-anion for all models CL&P, FPMD and CL&Pol for 20 and 200 IL pairs (n). If not mentioned explicitly, hereafter every time CMD data are discussed, we refer to simulations performed with the nonpolarisable FF (CL&P) in order to be consistent with the fact that all the other models apart from EMIM-BF 4 have been simulated with this FF. The limitations discussed here should be anyhow kept in mind. In

		n	R CL&P ( Å) R F P M D ( Å) R CL&P ol ( Å)
	DMIM-Cl	20	2.65	2.38	-
	DMIM-BF 4	20	2.05	2.27	-
	DMIM-PF 6	20	2.25	2.35	-
	DMIM-TFSI 20	2.56	2.26	-
	EMIM-BF 4	20	2.05	2.28	2.38
	EMIM-BF 4	200	2.05	-	2.35
	freedom, polarizability effects.

Table 3 .

 3 4. Computed Bader charges of the atoms of the DMIM-Cl and EMIM-BF 4 . Note that the charge values of some equivalent atoms have been averaged in agreement with schematic numbering of Figures 3.1B and 3.6.

		DMIM-Cl EMIM-BF 4
	N 1,2	-2.49	-2.50
	C 2	2.47	2.46
	C 4,5	0.80	0.79
	C 6	0.31	0.33
	C 7	-	0.40
	C 8	0.32	-0.04
	H 2	0.22	0.23
	H 4,5	0.20	0.18
	H 6	0.11	0.10
	H 7	-	0.10
	H 8	0.10	0.04
	DMIM, EMIM	0.95	0.94
	Cl, BF 4	-0.94	-0.95

Table 4 .

 4 1. Comparison of the bulk 2H-MoS 2 structural parameters computed in this work at the BLYP + vdW-W level (vdW-W stands for vdW dispersion forces described by the MLWF-based scheme) with the corresponding experimental data and also other theoretical works reported in literature.

		Exp.	FPMD (this work)	DFT	DFT
		[38, 39] BLYP + vdW-W B3LYP[273] HSE06 + vdW-D2[270, 272]
	a, b	3.16	3.26	3.17	3.16
	Mo-S	2.37	2.46	NR	NR
	Mo• • • Mo S• • • S	6.15 3.14	6.15 2.98	NR NR	6.19 NR

Table 4 .

 4 2. Computed Bader charges of the Mo and S atoms of 2H-MoS 2 in the bulk phase and of the two EMIM-BF 4 @MoS 2 interfaces simulated with the pristine and defected MoS 2 bilayers. The charge values per atom types have been averaged over the total number of atoms of the same species. Note that for the defected MoS 2 bilayer also the charges of the surrounding Mo and S atoms around the two S vacancies have been reported (*). This analysis allows to point-out minimal electronic density between the upper MoS 2 layer and the first EMIM-BF 4 layer at the interface, highlighting any evident direct chemical interactions between the two counterparts. This result, support the overall picture that the primary chemical interactions at play at the interface between EMIM-BF 4 ions and MoS 2 are electrostatic and dispersion (vdW) forces.

	MoS 2	Bulk phase bilayer	defected bilayer
	Mo	0.97(1)	1.03(1)	1.03(2)/0.97(3)*
	S	-0.49(2)	-0.52(1) -0.51(2)/-0.54(1)*
	EMIM-BF 4
	N 1,2	-	-2.50	-2.49
	C 2	-	2.46	2.38
	C 4,5	-	0.79	0.78
	C 6	-	0.34	0.42
	C 7	-	0.40	0.39
	C 8	-	-0.04	-0.06
	H 2	-	0.23	0.24
	H 4,5	-	0.18	0.17
	H 6	-	0.10	0.11
	H 7	-	0.10	0.08
	H 8	-	0.04	0.05
	EMIM	-	0.96	0.94
	BF 4	-	-0.95	-0.96
	simulations at 300K.

4.3.3.4 EMIM-BF 4 in contact with a defected MoS 2 surface: preliminary results

Table 5 .

 5 1. Point charges used for S and Mo atoms at the two MoS 2 electrodes in order to define the different surface charges distribution (and voltages). Each electrode is made of a bilayer of MoS 2 of 450 Mo atoms and 900 S atoms. Charge values very close to the average values computed by cpMD in

	Q S ( µC/cm 2 ) φ(z) (V) -Positive Mo +0.621 a +0.614 +0.606 +0.592 b +0.546 +0.472 0.00 ±3.00 ±7.00 ±13.71 ±36.56 ±73.12 0.00 ±0.70 ±1.60 ±2.96 ±8.15 ±16.5
	electrode S	-0.310 a -0.302	-0.292	-0.275 b -0.216	-0.122
	Negative Mo +0.621 +0.626 +0.634 +0.648 b +0.694 +0.768
	electrode S	-0.310	-0.318	-0.328	-0.345 b -0.404	-0.498

a MoS 2 point charges for zero-charge electrode taken from

[START_REF] Bi | Co-ion desorption as the main charging mechanism in metallic 1t-mos 2 supercapacitors[END_REF]

(computed by cpMD). b

Appendix

Scientific productions

Scientific publications