Quantum electronic properties of low-dimensional transition-metal systems: modeling, simulation and characterization

QUANTUM ELECTRONIC PROPERTIES OF LOW-DIMENSIONAL TRANSITION-METAL SYSTEMS: MODELING, SIMULATION AND CHARACTERIZATION. Thesis by irián sánchez ramírez for the degree of doctor of philosophy in physics Supervised by fernando de juan sanz and maia garcía vergniory University of the Basque Country (EHU/UPV) Donostia, May 2025. (cc) 2025 Irián Sánchez Ramírezs (cc by-nc-sa 4.0) Irián Sánchez Ramírez: Quantum electronic properties of low-dimensional transitionmetal systems: modeling, simulation and characterization, May 2025. A quien me enseñó a andar y aún sigue acompañándome en cada paso. Sobre todo, no pierdas tu deseo de caminar [...] la salud y la salvación sólo pueden ser encontradas en el movimiento. — Søren Kierkegaard CONTENTS Acronyms viii Publications x Abstract xii Resumen xiv i introduction 1 1 introduction 2 1.1Historical and conceptual background to this work . . . . . . 2 1.1.1First classical steps. . . . . . . . . . . . . . . . . . . . . . 3 1.1.2Advent of quantum mechanics. . . . . . . . . . . . . . . 4 1.2State-of-the art, goals and outline of this work . . . . . . . . . 8 1.3A brief comment on our approach . . . . . . . . . . . . . . . . 11 ii theory &methods 15 2 unraveling the many-body electronic hamiltonian 16 2.1A Hamiltonian to gobern them all . . . . . . . . . . . . . . . . 17 2.2Approximations........................... 19 2.2.1Born-Oppenheimer approximation . . . . . . . . . . . . 20 2.2.2Hartree approximation . . . . . . . . . . . . . . . . . . . 23 2.2.3Hartree-Fock approximation . . . . . . . . . . . . . . . 26 2.2.4Koopman’s theorem: from a computational simplification to a physical insight . . . . . . . . . . . . . . . . 28 2.3Density Functional Theory . . . . . . . . . . . . . . . . . . . . . 29 2.3.1The Hohnenberg-Kohn theorems . . . . . . . . . . . . . 30 2.3.2The Kohn-Sham approach . . . . . . . . . . . . . . . . . 33 2.3.3Local Density Approximation (LDA) . . . . . . . . . . . 35 2.3.4Generalized Gradient Approximation (GGA) . . . . . . 36 2.3.5A brief note on modified Becke-Johnson exchange potential............................. 39 2.3.6Concludingremark..................... 40 2.4Pseudopotentials .......................... 41 2.4.1Early precursors to pseudopotentials. . . . . . . . . . . 42 2.4.2Norm-conserving pseudopotentials . . . . . . . . . . . 42 2.4.3Projector augmented wave method . . . . . . . . . . . . 44 2.4.4Conclusion.......................... 46 2.5Bloch theorem and plane-wavebasis . . . . . . . . . . . . . . . 46 2.6How DFT reduces computational cost . . . . . . . . . . . . . . 49 2.7WannierFunctions ......................... 50 2.7.1Generalized Wannier Functions for Multiple Bands . . 51 iv contents v 2.7.2Non-Uniqueness and Maximally localized Wannier functions ............................. 52 2.7.3Applications and Wannier interpolation . . . . . . . . . 52 2.7.4Why and when using Wannier interpolation . . . . . . 53 iii results 56 3 quasi-1d tmcs 57 3.1Introduction ............................. 57 3.1.1Background ......................... 57 3.1.2Open questions and outline . . . . . . . . . . . . . . . . 58 3.1.3AbriefnoteonCDW.................... 60 3.1.4TMCs (MSe4)nIwith M=Nb,Ta............. 63 3.1.5Phenomenology overview on (MSe4)3I(M=Nb,Ta) TMCs. 64 3.2Methods & Theoretical results for (MSe4)3ITMCs. . . . . . . 66 3.2.1DFT details and methodology . . . . . . . . . . . . . . . 66 3.2.2DFT band structures. . . . . . . . . . . . . . . . . . . . . 67 3.2.3Tight-Binding models . . . . . . . . . . . . . . . . . . . 70 3.2.4Unfolding ab-initio and ARPES . . . . . . . . . . . . . . 74 3.2.5Optical conductivity . . . . . . . . . . . . . . . . . . . . 76 3.2.6Doping holes vs. electrons . . . . . . . . . . . . . . . . . 78 3.2.7Concluding remarks on the theoretical results. . . . . . 80 3.3Experimental realization: (NbSe4)3I. ............... 81 3.3.1Introduction: What are dark states. . . . . . . . . . . . . 82 3.3.2Detailed structural classification. . . . . . . . . . . . . . 84 3.3.3Approximately dark states as seen in ARPES. . . . . . 84 3.3.4Discussion & conclusion on the experimental results. . 87 3.4Discovery of a new quasi-1D TMC: (Nb4Se15I2)I2........ 88 3.4.1Introduction & brief theoretical and experimental details............................... 89 3.4.2Chemical & structural details. . . . . . . . . . . . . . . . 89 3.4.3DFT & theoretical analysis. . . . . . . . . . . . . . . . . 91 3.4.4ConcludingRemarks.................... 93 3.5Concluding remarks and outlook . . . . . . . . . . . . . . . . . 94 4 2d,bulk tmds &their heterostructures 97 4.1Introduction: Get to know TMDs . . . . . . . . . . . . . . . . . 97 4.1.1Theirstructure........................ 98 4.1.2TheirCDWs.......................... 99 4.1.3Their superconductivity. . . . . . . . . . . . . . . . . . . 100 4.1.4Open questions and outline. . . . . . . . . . . . . . . . . 101 4.2Background to our theoretial research . . . . . . . . . . . . . . 102 4.2.1A brief note on Mott & Kondo physics. . . . . . . . . . 104 4.3Charge transfer in T/H heterostructures of monolayer TMDs. 111 4.3.1Ab-initio methods......................114 4.3.2Charge density wave distortions . . . . . . . . . . . . . 115 contents vi 4.3.3Work function analysis . . . . . . . . . . . . . . . . . . . 116 4.3.4Chargetransfer .......................116 4.3.5VanderWaalseffect ....................119 4.3.6Udependence........................120 4.3.7Distance dependence . . . . . . . . . . . . . . . . . . . . 122 4.3.8Results: Bulk 4Hbpolytypes ...............123 4.3.9Conclusion..........................124 4.4Experimental realization: Metallicity through random stacking.127 4.4.1Introduction & Background . . . . . . . . . . . . . . . . 127 4.4.2dI/dV spectroscopy.....................128 4.4.3Why stacking disorder? . . . . . . . . . . . . . . . . . . 129 4.4.4Why not domain walls? . . . . . . . . . . . . . . . . . . 131 4.4.5DFTstudy..........................132 4.4.6Concluding remarks to this Section . . . . . . . . . . . 134 4.5Concludingremarks ........................135 iv concluding remarks 137 5 concluding remarks 138 v appendix 144 a introduction 145 a.1Brief philosophical comment . . . . . . . . . . . . . . . . . . . 145 a.1.1What is our approach? . . . . . . . . . . . . . . . . . . . 146 a.1.2What is true, what is real? . . . . . . . . . . . . . . . . . 149 a.1.3A small note on closure and categories . . . . . . . . . 152 a.1.4On discovery and emergence . . . . . . . . . . . . . . . 153 a.1.5Finalremark.........................154 b theory &methods 155 b.1Variationaltheorem.........................155 b.2Functional derivatives . . . . . . . . . . . . . . . . . . . . . . . 156 c results 158 c.1Unfolding ..............................158 c.1.1Whatisasupercell? ....................158 c.1.2Folding and unfolding wave vectors . . . . . . . . . . . 159 c.1.3Folding and unfolding of states . . . . . . . . . . . . . . 160 c.1.4Spectral weight in plane wave basis . . . . . . . . . . . 161 c.2A brief note on Fermi surface nesting & CDWs. . . . . . . . . 162 c.2.1CDW Formation via Nesting . . . . . . . . . . . . . . . 164 c.2.2Geometrical Interpretation of Nesting . . . . . . . . . . 164 c.2.3Quantifying Nesting: The Lindhard Susceptibility . . . 165 c.2.4Experimental Signatures of Nesting-Driven CDWs . . 165 c.3Angle-Resolved Photoemission Spectroscopy (ARPES) . . . . 165 c.3.1Photoemission and energy conservation . . . . . . . . . 166 c.3.2Momentum-Resolved Electronic Band Structure . . . . 167 contents vii c.3.3ARPES in the Study of CDW Materials . . . . . . . . . 167 c.3.4Self-Energy Effects and Many-Body Interactions . . . . 168 c.3.5Superconductivity and ARPES . . . . . . . . . . . . . . 168 c.4Supporting experimental details to Sec. 3.3............169 c.5Phase transition determination, supporting experimental data for Sec. 3.4...............................170 c.6Displacements of 2D-TMDs studied. . . . . . . . . . . . . . . . 171 c.7U bands of 2D-TMDsstudied....................171 c.8A brief note on STS & dI/dV....................171 c.9BCS Theory of Superconductivity . . . . . . . . . . . . . . . . . 175 c.9.1Electron-Phonon Interaction and Cooper Pair Formation176 c.9.2BCS Ground State and Energy Gap . . . . . . . . . . . 177 c.9.3Implications to our research . . . . . . . . . . . . . . . 179 c.10 Supporting information to Sec. 4.4.................180 c.10.1ACBands ..........................180 c.11 Supporting experimental information. . . . . . . . . . . . . . . 180 c.11.1Experimental superconductivity characterization . . . 181 c.11.2Concluding remarks on the experimental work . . . . 185 vi bibliography 186 bibliography 187 ACRONYMS •1D: One Dimensional •2D: Two Dimensional •ARPES: Angle-Resolved Photo-Emission Spectroscopy •BCS: Bardeen-Cooper-Schrieffer (model for superconductivity) •BO: Born Oppenheimer (Approximation) •CDW: Charge Density Waves •CVT: Chemical Vapour Transport •DFT: Density Functional Theory •DOS: Density of States •GEA: Gradient Expansion Approximation •GGA: General Gradient Approximation •IBZ: First Brillouin Zone •LDA: Local Density Approximation •MLWF: Maximally Localized Wannier Functions •mBJ: Modified Becke-Johnson (exchange potential) •NCPP: Norm-Conserving PseudoPotentials •OPW: Orthogonalized Plane Waves •PAW: Projected Augmented Wave •PBE: Perdew Burke Erzerhof •SG: Space Group •STM: Scanning Tunnel Microscope •SOC: Spin-Orbit Coupling •TCC: Categorial Closure Theory, Teoría del Cierre Categorial •TMC: Transition Metal Chalcogenides viii acronyms ix •TMD: Transition Metal Dichalcohenides •PW91: Perdew Wang 1991 •XRD: X-Ray Diffraction resumen xvi Estas hipótesis fueron posteriormente verificadas experimentalmente gracias a la colaboración con el grupo de N. Schröter, quien llevó a cabo mediciones de ARPES en (NbSe4)3I, confirmando la naturaleza indirecta del bandgap anticipada teóricamente. Posteriormente, se aplicó el mismo marco teórico a la caracterización de un nuevo compuesto sintetizado por el grupo de D. Shoemaker, (Nb4Se15I2)I2. A partir de cálculos de DFT, medidas de transporte y análisis de estados de oxidación, se confirmó su comportamiento aislante, demostrando tanto la aplicabilidad general del enfoque desarrollado como la ubicuidad de los TMCs cuasi-unidimensionales. Además, se propuso un flujo de trabajo orientado a la identificación de nuevos sistemas cuasi-1D con "estados aproximadamente oscuros" y propiedades optoelectrónicas ajustables. El Capítulo 4se dedica al estudio de los TMDs bidimensionales y sus heteroestructuras. Estos materiales han cobrado un protagonismo creciente en los últimos años debido a la fácil modulación de sus capacidades electrónicas , ópticas y magnéticas, así como por su aptitud para albergar fenómenos fuertemente correlacionados como CDWs, superconductividad, estados de Mott y el efecto Kondo. En este capítulo se analizan bicapas formadas por capas alternantes de tipo T(octaédricas) y H(prismáticas trigonales), así como otras fases como las 4Hby6R. Con el objetivo de elaborar un estudio sistemático de estos sistemas, en este Capítulo se desarrolla una metodología basada en DFT para calcular la transferencia de carga entre capas, un parámetro crucial para la descripción de fenómenos de correlación electrónica como el aislamiento tipo Mott o el efecto Kondo. En este estudio se contemplaron variaciones en la distancia intercapas, la composición de calcógenos, y la inclusión de correcciones electrón-electrón mediante el parámetro Ude Hubbard y correcciones de Van der Waals. Los resultados presentados revelan que la distancia entre capas es el principal factor que regula la transferencia de carga, mientras que las correcciones asociadas a Uy a las interacciones de Van der Waals presentan un efecto secundario. Asimismo, se identifican tendencias sistemáticas según la composición química del sistema (M = Nb, Ta; X = Se, S), observándose que los compuestos formados por Ta y Se exhiben una transferencia de carga mayor en comparación a los compuestos formados por Nb y S. Estos resultados permiten contextualizar las observaciones experimentales en materiales como TaS2, TaSe2y NbSe2. En colaboración con el grupo de M. Ugeda, en este Capítulo se caracteriza el compuesto 1T-TaSSe, un análogo metálico isoestructural al aislante 1T-TaSe2. Los cálculos realizados muestran que el desorden en el apilamiento de las CDWs facilita el tunelado intercapas, proporcionando una explicación para la metalicidad de este compuesto. La comparación con resumen xvii resultados experimentales obtenidos mediante técnicas de microscopía y espectroscopía de efecto túnel respalda esta interpretación, confirmando el comportamiento metálico del 1T-TaSSe y ayudando a la interpretación experimental de la coexistencia de CDW y superconductividad en este sistema. Finalmente, el último capítulo de la Tesis reflexiona sobre la madurez y las limitaciones de los marcos teórico-computacionales desarrollados. En el caso de los sistemas cuasi-1D, aunque se logra reconciliar discrepancias reportadas en la literatura, se identifica la necesidad de establecer criterios predictivos sistemáticos para la búsqueda de nuevos materiales. En cuanto a las heteroestructuras bidimensionales, el marco propuesto demuestra ser robusto, aunque se prevee que su aplicación a sistemas de mayor complejidad estructural requerirá ajustes adicionales. En ambos casos, los métodos desarrollados muestran potencial de escalabilidad y ofrecen una base sólida para extender su aplicación a nuevas líneas de investigación en física de materiales. En resumen, esta Tesis plantea un enfoque sistemático para el estudio de las propiedades electrónicas de sistemas formados por metales de transición con dimensionalidad reducida. Este enfoque se lleva a cabo a través de la combinación de cálculos basados en DFT, técnicas de band unfolding, modelos de enlace fuerte y validación experimental colaborativa, construyendo así una plataforma metodológica que no solo contribuye a participar y tratar de esclarecer debates abiertos en el campo, sino que también abre nuevas vías para la exploración de materiales con propiedades electrónicas exóticas y que albergan fenómenos relacionados con la fuertes correlaciones. Part I INTRODUCTION 1 INTRODUCTION The research presented in this Thesis lies at the intersection of theoretical, computational and experimental condensed matter physics. Throughout the following pages, our goal is to develop robust theoretical frameworks that serve as platforms to unify, understand and contextualize past research, while paving the way for future advancements in the field. However, these frameworks cannot exist in a vacuum; they are firmly grounded in a broader scientific and historical context. This introduction aims to establish exactly that: the who,what,where, and when of this work1, setting the stage for the questions it addresses and the research it aims to inspire. First, we will provide a brief quasi-cronological outline of the timeline for the development of condensed matter physics as a field, highlighting the discoveries that are most relevant to our research: the techniques employed and the ongoing debates in the field where our work can contribute. Following this, we will present the state-of-the-art context in which our research takes place, along with a brief outline connecting our work to the current advancements in the field. Finally, we will offer a brief comment on our approach. 1.1 historical and conceptual background to this work Condensed matter physics is the branch of physics that studies the macroscopic and microscopic properties of matter in its solid and liquid states, where high particle densities make interactions between constituent ions and electrons central to understanding their behavior. As a result, condensed matter physics focuses on how large ensembles of interacting particles give rise to emergent2collective phenomena: behaviors that cannot be deduced from studying the individual components in isolation. In what follows, we will outline key discoveries that led to the development of the central techniques and concepts relevant to this Thesis: 1These questions can be understood in terms of the "five Ws" checklist. We have excluded the why, as it pertains to final causes, which fall outside the scope of this Thesis. Additionally, the question of where will be addressed in terms of a conceptual location, rather than a spatial one. 2Emergence occurs when a complex system (composed of multiple entities, such as atoms in the case of condensed matter physics) exhibits properties or behaviors that are not present in its individual components, but arise through their interactions. 2 1.1 historical and conceptual background to this work 3 primarily density functional theory (DFT)3, but also several many-body phenomena that, while not the primary focus, form an essential part of the conceptual background of this work, such as charge density waves (CDWs) (Sec. 3.1.3), superconductivity (App. C.9), and the Kondo and Mott effects (Sec. 4.2.1). To do so, we begin by revisiting the early classical foundations of condensed matter physics and exploring how the advent of quantum mechanics in the early 20th century transformed the field. We illustrate this transformation through a selection of representative discoveries, which highlight the power of the quantum formalism. From there, we shift our focus to the historical and conceptual motivations underlying the research presented in Chapters 3and 4, along with the techniques employed and the current state-of-the-art in the field. 1.1.1First classical steps. To provide proper temporal context for our research, we will trace back to the first known model of electrical conduction4, the Drude model. The choice of Drude’s model as the starting point for our sui generis historical review stems from our view that it represents one of the earliest attempts to understand the macroscopic properties of solids from a formal theoretical perspective, moving beyond purely phenomenological explanations. In 1900, Paul Drude proposed a theoretical model for electronic conduction based on classical physics [2]. In simple terms, under Drude’s model, metals are described as a gas of free electrons moving through a periodic lattice of ions, scattering off them in a manner analogous to a pinball bouncing off obstacles, thereby giving rise to resistivity. Despite its limitations, this model was remarkably successful for its time, offering the first microscopic interpretation of macroscopic transport laws in solids. It accounted for sev3In short, DFT is a quantum-mechanical framework for describing many-electron systems based on their ground-state electron density rather than their full many-body wavefunction. Founded on the Hohenberg–Kohn theorems and practically implemented through the Kohn–Sham scheme, DFT provides the basis for widely used computational methods that can predict structural, electronic, and energetic properties of materials with remarkable accuracy. We will motivate and present DFT in Chapter 2. 4Historically, other studies can be related to condensed matter physics, such as the investigations of condensed states carried out by Humphry Davy or Michael Faraday during the nineteenth century. However, we shall begin in 1900 with Drude’s model of conduction, as it provides a natural starting point for building up towards quantum treatments and highly correlated many-body phenomena. It is also the starting point from "Solid state physics" by Aschroft & Mermin [1], one of the go-to book for introductory courses in solid-state physics in the undergraduate level. 1.1 historical and conceptual background to this work 4 eral empirical observations, including the Wiedemann–Franz law5[3]. However, it was not without shortcomings: it significantly overestimated the electronic heat capacity of metals and failed to account for the temperature dependence of resistivity at low temperatures, discrepancies that would later motivate the incorporation of quantum mechanics into solidstate theory. 1.1.2Advent of quantum mechanics. With the advent of quantum mechanics in the 1920s and 1930s, Drude’s model and condensed matter physics as a field was profoundly refined by physicists such as W. Pauli, A. Sommerfeld, and F. Bloch; who employed the emerging quantum formalism to develop more accurate6theories of physical phenomena. For instance, in 1927, Pauli recognized that electrons in metals obey Fermi– Dirac statistics7and used this insight to develop a quantum theory of paramagnetism8[6]. That same year, Sommerfeld incorporated Fermi–Dirac statistics into the classical Drude model, giving rise to the free-electron model, an important conceptual breakthrough that marked the transition from classical to quantum descriptions of electrons in solids. By treating electrons as a degenerate Fermi gas rather than a classical ensemble, this model resolved key discrepancies in the Drude theory and provided a significantly improved account of charge carrier behavior in metals. Notably, it succeeded in explaining several experimental observations, including the markedly lower electronic heat capacity of metals. However, it still neglected electron–electron interactions, limiting its range of applicability and highlighting the need for more sophisticated many-body treatments, which we will explore later. 5The Wiedemann–Franz law states that the ratio between the electronic contribution to the thermal conductivity (κ) and the electrical conductivity (σe) is proportional to the temperature (T), formally expressed as κ/σe=KT, where Kis a constant. 6By more accurate, we mean that these theories offered greater predictive power and showed improved agreement with experimental observations in condensed matter systems. For instance, Sommerfeld’s quantum refinement of Drude’s model successfully accounted for the electronic heat capacity of metals, a key phenomenon that the classical Drude model could not explain, thereby rendering it a more accurate description of metallic behavior. 7Fermi–Dirac statistics, derived independently by E. Fermi and P. Dirac around 1926 [4,5], describe systems of non-interacting, identical fermions with half-integer spin that obey the Pauli exclusion principle. 8Paramagnetism is a magnetic phenomenon in which paramagnetic materials are weakly attracted to externally applied magnetic fields and develop an internal induced field aligned with the applied field. In contrast, diamagnetic materials are repelled by magnetic fields and generate fields in the opposite direction. 1.1 historical and conceptual background to this work 5 Meanwhile, in the late 1920s and early 1930s, another crucial avenue of condensed matter physics emerged: crystallography. Developed by A. Bravais and Y. Fyodorov, crystallography provided essential tools for the structural analysis of crystalline materials9, leading to their classification by symmetry groups and culminating in the publication of the International Tables of Crystallography in 1935. The importance of this development for condensed matter physics cannot be overstated: it offered a rigorous classification of solids based on their symmetry, established a bridge between atomic-scale structure and macroscopic properties, and laid the groundwork for band theory. Around the same time, in 1929, Bloch formulated his eponymous theorem [7] (see Sec. 2.5), which became central to bandstructure theory. In 1930, L. Brillouin introduced the concept of the Brillouin zone, providing a compact framework for describing wave-like solutions in periodic media. According to Bloch’s theorem, all such solutions can be fully characterized within a single Brillouin zone. Crystallography, Bloch’s theorem, and the Brillouin zone remain of paramount importance in condensed matter physics in general, and in this Thesis in particular, as they form the conceptual and computational foundation for bandstructure calculations. Without a well-defined crystal structure, one cannot define the inputs to DFT; and without Bloch’s theorem and Brillouin zones, one cannot compute bandstructures10. DFT is the main computational technique employed in this Thesis. It was first formalized by W. Kohn and P. Hohenberg in 1964 [8], but its conceptual origins trace back to the Thomas–Fermi model (1926), which described electrons in terms of their density and was further refined by Dirac in 1930 by including an analytical expression for the exchange energy[9], key for DFT. The development of DFT built upon a number of key approximations, such as the Born–Oppenheimer approximation (1927) [10], Hartree’s meanfield treatment (1928) [11], and the Hartree–Fock generalization (1930) [12], all of which are discussed in Sec. 2.1. DFT has since transformed condensed matter physics, enabling the discovery of new materials, the prediction of topological phases, and the construction of robust computational frameworks that systematize the study of material properties 11. 9In the context of condensed matter physics, a crystal is a solid whose constituent atoms, ions, or molecules are arranged in a highly ordered and repeating pattern extending in all directions. 10 A bandstructure describes the allowed energy levels of electrons in a solid as a function of their crystal momentum. In a periodic lattice, Bloch’s theorem shows that electron wavefunctions take the form of plane waves modulated by the lattice, leading to the formation of continuous energy bands rather than discrete atomic levels. The structure and filling of these bands determine whether a material behaves as a metal, semiconductor, or insulator. An example of a bandstructure can be found in Fig. 3.3 11 The decision to include a technique such as DFT in the conceptual background of this Thesis is deliberate: we understand theory and technique (whether experimental, theoretical, or 1.1 historical and conceptual background to this work 6 Although DFT is the primary tool used in this Thesis, it is not the only theoretical technique worth contextualizing. In 1928, as part of his doctoral dissertation, F. Bloch introduced the Linear Combination of Atomic Orbitals (LCAO) method for solids12, offering an intuitive approach to constructing electronic states from localized atomic orbitals. While powerful, this method was later complemented by a more practical and widely used approximation: the parametrized tight-binding model, developed by J. C. Slater and G. F. Koster around 1954 [14]. This method is especially effective for studying low-energy phenomena in complex materials and underpins the tight-binding models we use in Sec. 3.2.3, where we describe the quasi-1D behavior of (TaSe4)3I and (NbSe4)3I. Another indispensable technique in modern computational condensed matter is Wannierization (see Sec. 2.7), which relies on the construction of localized Wannier functions, as introduced by G. Wannier in 1937 [15]. While band theory and its computational implementations have successfully explained many aspects of electronic behavior, they remain approximations that do not fully account for electron–electron interactions. For instance, standard DFT includes electron–electron interactions in an average sense through the exchange-correlation functional, but often fails to capture strong correlation effects, such as those arising in systems with partially filled dor forbitals. These many-body interactions give rise to a variety of collective electronic states that go beyond effective single-particle pictures: phenomena such as Mott insulators (see Sec. 4.2.1.1), first proposed by N. F. Mott in 1937; superconductivity, discovered in 1911 by H. K. Onnes and F. Holst and later formalized microscopically by Bardeen, Cooper, and Schrieffer in 1957 (see App. C.9); and the Kondo effect, observed experimentally in the 1960s and explained by J. Kondo in 1964 (see Sec. 4.2.1.2). Another collective phenomenon of particular relevance to this Thesis is the Charge Density Waves (CDW), first described by R. Peierls in 1955 for one-dimensional metals (see Sec. 3.1.3). CDWs were among the earliest motivations for the study of low-dimensional systems, a subject we now turn to. One key aspect of the systems studied in this Thesis is that most of them exhibit effective or real low-dimensionality. Dimensionality has been a topic of debate in both physical sciences and philosophy since the time of Aristotle (384-322 BC). In modern physics and mathematics, we can trace the development of this concept to B. Riemann, who, in 1854, studied highercomputational) to be linked by a cooperative relationship in which one enables and shapes the other, see Sec. 1.3and App. A.1. 12 The method was first introduced for molecules by B. N. Finklestein and G. E. Horowitz, and was heavily influenced by R. Mulliken’s work on molecular orbitals [13]. 1.1 historical and conceptual background to this work 7 dimensional manifolds, influencing both geometry and theoretical physics. Following this, A. Einstein’s General Relativity (1916) combined space and time into a four-dimensional continuum, while T. Kaluza (1921) and F. Klein (1926) extended this framework to five dimensions in an attempt to unify the fundamental forces. More recently, string theories have suggested the possibility of 10 or 11 dimensions13 (for example, M. Green in 1987 and J. Polchinski in 1998). In contrast, condensed matter physics is usually more concerned with the study of reduced dimensionality. The interest in dimensionality can be traced back to the concept of quantum confinement developed in the 1920s and 1930s, which is closely linked with Bloch’s theorem and band theory. By the mid-20th century, theoretical and experimental foundations were established to explore genuine low-dimensional materials. This is the period when, in the 1950s, R. Peierls and H. Fröhlich predicted that one-dimensional metals should spontaneously distort into insulating states, accompanied by periodic modulations of the charge density, an effect closely related to CDWs. Experimental studies on quasi-onedimensional (quasi-1D) conductors, such as transition metal chalcogenides, began in the 1970s and 1980s (carried out by the groups of A. Meerschaut, P. Gressier, and M. Izumi, among others [16–18]), constructing the foundation for the research presented in Chapter 3. This research aimed to identify signatures of Peierls transitions and CDW states. Concurrently, research into layered two-dimensional materials14 particularly graphite intercalation compounds, gained traction, paving the way for future investigations into more complex two-dimensional systems. A major turning point occurred in 2004, when the groundbreaking isolation of graphene (first theorized by P. R. Wallace in 1947 [21]) by A. Geim and K. Novoselov [22] sparked a global wave of interest in two-dimensional materials. This Nobel-winning discovery led to extensive studies on electronic correlations, quantum transport, and topological phenomena in atomically thin layers. Shortly thereafter, attention shifted rapidly toward transition metal dichalcogenides (TMDs), which emerged as versatile, easily exfoliable, and highly tunable materials. Their first discovery dates back to 1923 by L. Pauling [23], and the first reported production of MoS2monolayer suspensions was in 1986. Seminal works from the 2010s demonstrated exceptional electronic, optical and spin-valley physics in TMDs [24–29], revealing phenomena such as indirect-to-direct bandgap transitions, superconductivity, CDWs, and Mott insulating states in the two-dimensional limit [30,31]. More recently, TMDs have evolved in both complexity and 13 In string theories, extra spatial dimensions are typically compactified or curled up at scales smaller than atomic dimensions, rendering them unobservable at macroscopic levels but still influential in the fundamental laws of physics. 14 Initially considered thermodynamically unstable and non-existent by L. Landau [19] and R. Peierls [20]. 1.2 state-of-the art,goals and outline of this work 8 scientific interest, especially with the engineering of TMD heterostructures [32]. In particular, heterobilayers and twisted multilayers have created new platforms for studying strongly-correlated phenomena. This is the research domain explored in Chapter 4. 1.2 state-of-the art,goals and outline of this work As mentioned in the preamble to this introduction, the primary goal of this Thesis is to develop robust theoretical frameworks that serve as platforms to unify, understand, and contextualize both past and future research. These frameworks are embedded within the scientific and sociological landscapes from which they emerge, making this Thesis a product of its time, shaped by contemporary debates in condensed matter physics and informed by the discoveries and contributions of our closest peers. What binds the different chapters of this work together, beyond a shared methodology, is the consistent focus on quantum materials composed of transition metals with varying oxidation states. These variations crucially affect the materials’ electronic structure and collective behavior, offering a rich testing ground for the theoretical tools we develop. In this Section, we aim to contextualize the open questions, answers, and results presented in this Thesis both historically and conceptually, highlighting how our contributions intersect with and advance the state-of-the-art in condensed matter physics. Chapter 2is devoted to a conceptual and methodological review of the many-body Hamiltonian in condensed matter physics, along with the approximations and computational techniques commonly used to address its inherent complexity. We begin by introducing the general form of the Hamiltonian, followed by an exploration of the main approximations that render the problem tractable. These include the Born–Oppenheimer approximation in Section 2.2.1, as well as the Hartree and Hartree–Fock approximations in Sections 2.2.2and 2.2.3, respectively. Subsequently, we review and analyze the advent of DFT in Section 2.3, examining its theoretical foundations, computational cost, and key approximations: such as the Local Density Approximation (LDA) and the Generalized Gradient Approximation (GGA), discussed in Sections 2.3.3and 2.3.4. To conclude the discussion on DFT, we introduce the pseudopotential technique in Section 2.4, which significantly reduces computational cost by replacing the explicit treatment of core electrons with effective ionic potentials. This simplification is essential for efficiently studying the systems addressed in this Thesis. Finally, we introduce Bloch’s theorem in Section 2.5, along with Wannierbased techniques in Section 2.7, which provide a powerful bridge between first-principles calculations and effective tight-binding models. Part II THEORY & METHODS 2 UNRAVELING THE MANY-BODY ELECTRONIC HAMILTONIAN During the introduction to this Thesis, we contextualized our research within the broader landscape of scientific knowledge to which it belongs. Beyond its defining coordinates within this landscape, the cohesion of this manuscript is not solely derived from the similarities among the studied systems but rather from the approach taken to investigate them. As stated in the introduction, the work presented in the following pages aspires to build a robust and systematic understanding of these materials. To that end, we aim to develop theoretical frameworks that not only unify past research but also lay a solid foundation for future theoretical and experimental advancements. To achieve this, we primarily employ ab initio methods, with a particular focus on DFT. DFT serves as an ideal tool for this systematic approach due to its theoretical rigor and computational reliability. The majority of the theoretical results presented in this Thesis are derived within the DFT framework and, therefore, will be considered correct1as long as the approximations inherent to DFT remain applicable. This Chapter is dedicated precisely to understanding the foundations, limitations, and computational techniques associated with DFT. We begin by formulating the most general many-body Hamiltonian for a system of interacting electrons and ions within quantum mechanics, gradually developing the core principles and approximations that define DFT. In doing so, we will highlight its computational scaling, physical approximations, and the techniques that refine and extend its applicability such as Wannier functions. This structured approach will serve as the methodological foundation for the results presented in the subsequent Chapters. As we anticipated during the introduction, the research presented in this Thesis is mainly located within the field of condensed matter physics. In said field, the low-energy behavior of materials is mainly governed by the interactions between electrons and ions in a many-body quantum system. While, theoretically, such systems can be fully described by a general manybody Hamiltonian, the deep complexity arising from the vast number of interacting particles makes an exact solution computationally intractable. Consequently, a hierarchy of approximations is required to make progress in understanding these systems. This Chapter begins by formulating the many-body Hamiltonian for elec1Aligning with comment presented in Sec. 1.3and App. A.1 16 2.1 a hamiltonian to gobern them all 17 tronic systems and systematically introducing the approximations necessary to make the problem computationally tractable2. Starting with fundamental quantum mechanical principles, such as the Born-Oppenheimer approximation, we progressively refine the treatment of electronic interactions through the Hartree and Hartree-Fock methods, ultimately arriving at DFT and the pseudopotential approximation. Each step in this methodological progression reduces the computational complexity while maintaining accuracy by incorporating physically meaningful approximations tailored to the characteristics of the studied systems. The culmination of these developments is DFT, which provides an efficient and accurate framework for describing the electronic structure of condensed matter systems. Its balance between computational feasibility and predictive power has made it the cornerstone of modern computational condensed matter physics. As such, DFT serves as the principal theoretical tool employed throughout this Thesis, underpinning the systematic approach taken to investigate the materials discussed in the following Chapters. 2.1 a hamiltonian to gobern them all In the context of the Quantum Mechanics, the Hamiltonian His an operator that represents the total energy of a quantum system. The time evolution of such system is governed by the Schrödinger equation, which can be expressed as i h∂ ∂t |ΨS(t)⟩=H |ΨS(t)⟩, (2.1) where |ΨS(t)⟩is the quantum state of the system in Dirac notation. In condensed matter systems, which are the primary focus of this Thesis, |ΨS⟩describes the behaviour of all particles in the system: basically electrons and ions. Starting from the Hamiltonian formulation, we can express the Hamiltonian operator as H=ˆ T+ˆ V, (2.2) where ˆ Tand ˆ Vare the kinetic and potential energy operators of the system, respectively. In this manuscript, all the systems of study are constituted 2In light of the ideas discussed in Sec. 1.3and App. A.1, when we say a problem is "computationally tractable," we don’t mean there’s some ultimate truth we’re failing to reach because our computers aren’t fast enough. Rather, we recognize that the tools and techniques available to us define the boundaries of what kinds of operations we can meaningfully carry out in condensed matter physics today. Approximations, then, are not flaws or shortcuts in our theories, but part of the legitimate set of operations we rely on to generate meaningful and reproducible results. They make it possible to construct the results presented in this Thesis, which is why we consider them not only necessary, but fundamental to the work included in this Chapter. 2.1 a hamiltonian to gobern them all 18 by ions and electrons interacting with each other3. In order to capture all relevant interactions in these systems, we can rewrite 2.2as H=ˆ TI+ˆ Te+ˆ VI−I+ˆ Ve−e+ˆ Ve−I, (2.3) where the kinetic operators from ions and electrons are denoted by ˆ TIand ˆ Te, respectively. In 2.3, the interaction terms within ions and within electrons are captured by ˆ VI−Iand ˆ Ve−ewhile ˆ Ve−Irepresents the interaction between electrons and ions. We can further expand 2.3as H=− h2 2meX i∇2 i−X i,I ZIe2 |ri−RI|+1 2X i=j e2 |ri−rj| X I h2 2MI∇2 I+X I=J 1 2 ZIZJe2 |RI−RJ|, (2.4) where his the Plank’s constant and the terms related to electrons are denoted by lower case subscripts and ions, with charge ZIand mass MI, are denoted by upper case subscripts. In 2.4,iruns over the number of electrons Neand Iruns over the number of ions NI. The expression in 2.4represents the simplest form of the Hamiltonian for a condensed matter system, accounting solely for interactions between the particles within the system itself, without external magnetic or electric fields. It might seem intuitive that, given all the necessary information, one could predict the system’s properties and time evolution by solving the Schrödinger equation for the entire system. HΨS({ri,σi;RI};t) = ESΨS({ri,σi;RI};t), (2.5) where ΨSand ESrepresent the wavefunction4and total energy for the entire condensed matter system. However, there are practical reasons why this approach is not feasible, highlighting the need for approximations to solve such Hamiltonians, especially those involving more than two particles. In the next few paragraphs we address the two primary sources of these limitations5: 1. It is a many-body problem of, at least, EXP computational complexity [61]. 3i.e. electrons with electrons, electrons with ions and ions with ions. 4Just a brief note on notation: previously, we were using the Dirac notation to denote the quantum state of the whole system |ΨS⟩and making explicit its time dependence. In Dirac notation, the state is just a vector in the Hilbert space spanned by the eigenstates from H. Now we will use the wavefunction for convenience, which is just the representation of the state in the space-basis. Usually we can relate both by ψ(x) = ⟨x|ψ⟩. 5Actually, both sources of limitations are very similar and could be treated as one. Nonetheless, the choice to make distinction intentionally leaves room for how quantum computation could adress correlation effects in a much more efficient fashion through the secod point. 2.2 approximations 19 The Hamiltonian in 2.4describes a many-body system with interactions among electrons, among ions, and between electrons and ions. The wavefunction for this system depends on the positions of all electrons and ions, so for N=Ne+NIparticles, the wavefunction is defined over a 3N-dimensional continuous space (without considering additional internal degrees of freedom such as spin), Ψ:R3N →C. (2.6) To treat the problem numerically, one must discretize this space using Mgrid points per spatial coordinate, leading to a total of M3N possible configurations. As a result, both storing the wavefunction and performing exact diagonalization of the Hamiltonian require computational resources that scale exponentially with N. This places the manybody problem, in its full generality, within the class of problems of at least EXP computational complexity 6, rendering it intractable for large N. 2.Correlation effects must be taken into account. The third term in the Hamiltonian 1 2Pi=je2 |ri−rj|introduces strong correlations among electrons. Each electron experiences the influence of all other electrons, leading to a complex many-body problem. Electrons obey the Pauli exclusion principle and the effects of exchange interaction, which arises from their indistinguishability and the requeriment of their total wavefunction to be anti-symmetric under particle exchange. The correct treatment of this interactions cannot be attained simply by applying pairwise forces as in a classical context. Instead, the solution would involve non-local interactions and entanglement, where the state of one particle depends on the state of all others. These quantum effects cannot be efficiently simulated by a classical computer due to the exponential scaling problem coming from quantum phenomena. 2.2 approximations While the challenges presented in the previous Section might seem discouraging, they motivate our need to develope approximations that make the 6A problem is EXP-complete if it can be solved in deterministic exponential time and is EXP-hard, meaning that any other problem in the class EXP can be reduced to it using polynomial-time transformations. In the case of many-body systems, the exponential growth of the Hilbert space with particle number implies that both storage and spectral queries (e.g., whether the ground state energy lies below a given threshold) require exponential time, and are among the hardest problems in EXP. 2.2 approximations 20 problem tractable. We will begin by reducing the effective particle count, N, inside the big O notation O(M3N). 2.2.1Born-Oppenheimer approximation The first step in simplifying our problem is to decouple the ionic and electronic degrees of freedom. Kinetic energy terms in the Hamiltonian are linearly proportional to particle mass and quadratically proportional to velocity. According to equipartition theorem7both of these terms will be also proportional to 3/2kBT, implying that the ratio of their velocities will scale as vI ve =rme mI . (2.7) From 2.7and by substituting the mass of the lightest nucleus we can find, i.e. hydrogen with MI≈1836 me, we can conclude that ions move at least around 40 times slower than electrons. Consequently, on the time-scale of nuclear (ionic) motion, electrons are expected to rapidly relax to an instantaneous ground-state configuration. This will allow us to approximate the ions as stationary particles when obtaining the electronic ground-state by solving the time-independent Schrödinger equation, an approximation which is known as the Born-Oppenheimer (BO) approximation [10]. In order to properly introduce the BO approximation, lets start by stating the Adiabatic theorem, Theorem 2.2.1A physical system remains in its instantaneous eigenstate if a given perturbation is acting on it slowly enough and if there is a gap between the eigenvalue and the rest of the Hamiltonian’s spectrum [62]. Then, the BO approximation assumes that electronic eigenstates are not affected by perturbations or movements from the ions. To validate this assumption, we follow [63] and assume that the wavefunction of the entire system, including electrons and ions, can be expressed as ˜ Ψ({ri},{RI}) = Ψ({ri};{RI})Φ({RI}), (2.8) In expression 2.8, the wavefunction of the entire system is divided into two components: Ψ({ri};{RI}), which represents the electronic part influenced by 7The equipartition theorem states that, at thermal equilibrium, each degree of freedom contributes 1 2kBTto the system’s average energy (in the main text we are using the average, Ekin =m 2(v2 x+v2 y+v2 z)→ ⟨Ekin⟩=3/2kBT). Eventhough this theorem is particularly suited for classical systems, since the quantum kinetic enery operator scales inversely with mass − h2 2m ∇2∝1 m, we can use it to state that for particles in motion in a quantum system, their velocity v∝q1 m. 2.2 approximations 21 both the ionic and electronic positions, and Φ({RI}), which represents the ionic part and depends only on the ionic positions. This separation assumes that the two components are independent8. Furthermore, we require that Ψ({ri};{RI})satisfies time-independent Schrödinger equation, 1 2 X iX j=i 1 |ri−rj|−X i∇2 i−X i,I 2ZI |ri−RI|+ Ψ({ri};{RI}) =Ee(RI)Ψ({ri};{RI}), (2.9) where we aknowledge the dependence of the electronic energy eigenvalue (or adiabatic energy contribution) Ee(RI), on the fixed ionic positions9. Applying the full Hamiltonian 2.4to the entire wavefunction and setting h=me=e2=4πϵ0=1, (2.10) we arrive at H˜ Ψ({ri};{RI}) =  −X j 1 2MJ∇2 I+Ee({RI}) +1 2X JX K=J ZJZK |RJ−RK| ˜ Ψ({ri},{RI}) =Ψ({ri};{RI}) −X J 1 2MJ∇2 J+Ee({RI} +1 2X JX K=J ZJZK |RJ−RK| Φ({RI}) −X J 1 2MJ [2∇JΦ({RI})·∇JΨ({ri};{RI}) + Φ({RI})∇2 JΨ({ri};{RI})i, (2.11) where Ee({RI})is called the adiabatic energetic contribution. The remaining non-adiabatic terms contribute neglibly to the energy. To show this, we apply time-independent perturbation theory. The first order correction from non-adiabatic terms has the form 10 8Here we are considering two independent subsystems whose wavefunctions live in a Hilbert space that is the tensor product of the individual Hilbert spaces for each subsystem, since the systems are independent, the tensor product state can be (and is) represented by the product of states. This can also be understood as the joint probability of two independent events. 9Since the electrons are expected to rapidly adapt to a ground-state depending on a set of fixed atomic positions, the electornic energy eigenvalue Ee(or what we will call adiabatic energetic contribution later) only depends on said ionic positions. 10 Here ×is only used to indicate a line break that separates a scalar product, not the vectorial product. 2.2 approximations 22 −ZY j drjY β dRJΨ∗({ri};{RI})Φ∗({RJ})× X K 1 MK [∇KΦ({RI})·∇KΨ({ri};{RI})] = − X KZY J dRJΦ∗({RI})∇KΦ({RJ})×(2.12)  ZY j drjΨ∗({ri};{RI})∇KΨ({ri};{RI}) , (2.13) where the term in brackets vanishes, ZY j drjΨ∗({ri};{RI})∇KΨ({ri};{RI}) =1 2∇KZY j drj|Ψ({ri};{RI})|2=1 2∇K(1) = 0, (2.14) since the normalisation of the electornic wave-function does not change when the ions move, making the first-order contribution to vanish. The second-order shift gives rise to transitions in the electrons when the ions move, this phenomenom its known as electron-phonon interaction, and has an impact on the total energy11. The second non-adiabatic term becomes significant when the electrons iare tightly bound to ions I. In that case, we can express the wavefunction as depending on the distance α(i,I)between electrons and ions, Ψ({ri};{RI}) = Ψ({α(i,I)}), (2.15) where α(i,I)=ri−RI. The first-order correction of this term is −ZY j drjY J drJ{α(i,I)}Φ∗({RI})X K 1 2MKhΦ({rI})∇2 K{α(i,I)}i =X K 1 2"ZY K drJ|Φ({RI}|2# ZY j,J dα(j,J)Ψ∗{α(i,I)}∇2 KΨ{d(i,I)}  =X (k,K) 1 MKZY (j,J) dα(j,J)Ψ∗{α(i,I)}1 2∇2 (k,K)Ψ{α(i,I)}, (2.16) which is of the order of the electornic kinetic energy multiplied by the ratio me/mI, allowing us to treat it as negligible. Consequently, we can then 11 This electron-phonon part is crucial for all the phenomena surrounding the systems studied in this Thesis, more specifically those in Chapter 4. 2.2 approximations 23 neglect all non-adiabatic terms and decouple the Schrödinger equation for ions as  −X J 1 2MJ∇2 J+Ee({RI}) + 1 2X JX K=J ZJZK |rJ−rK| Φ({rI})=EΦ({rI}), (2.17) which allows us to treat the non-adiabatic electronic terms independently. This approximation will drastically reduce Nin O(N3M)through the decoupling of ionic and electornic degrees of freedom, providing a crucial simplification in the endeavor to solve the many-body Hamiltonian. 2.2.2Hartree approximation Thanks to the BO approximation, we have succesfully simplified the ionic + electronic problem into a electronic one by decoupling the motion of electrons and ions. However, even though this has reduced the effective number of particles in the problem12, we are still dealing with a problem of, at least, EXP-complete computational complexity. The next logical step further in our approach to solving the Hamiltonian is the method proposed by Douglas Hartree [11] in 1928, which treats the electrons as independent particles that move independently within an averaged potential created by all other electrons. This approach effectively reduces the complexity of the problem into a more manageable polynomial scaling of O(N). The first step in this method is what is known as the independent particle approximation: In a many-body Hamiltonian for Nparticles of spin σgiven by H(r1,r2,...,rN;t) = N X i=1 ˆ P2 i 2mi +Vσ(r1,r2,...,rN;t), (2.18) where ˆ Pand miare the momentum operator and mass for particle i. If particles in 2.18 do not interact with each other, we can express the potential term as a sum of independent terms Vσ(r1,r2,...,rN;t) = N X i=1 Vσ(i)(ri;t), (2.19) and hence, rewrite the Hamiltonian as H(r1,r2,...,rN;t) = N X i=1 H(ri;t). (2.20) 12 Nin O(N3M). 2.2 approximations 24 Since the particles are non-interacting, we can represent the system’s wavefunction, Ψ, as a product of each particle’s wavefunction ϕσ i: Ψ(r1,r2,...,rN;t) = ϕσ 1(r1,t)ϕσ 2(r2,t)...ϕσ N(rN,t). (2.21) Each of these single-particle wavefunctions ϕievolves over time as ϕi(ri;t) = ϕϵi(ri)e−iϵit h, (2.22) allowing the system’s wavefunction to be expressed as Ψ(r1,r2,...,rn;t) = ΨE(r1,r2,...,rn)e−iEt h, (2.23) where E=PN i=1ϵi. If we combine both independent electron approximation and the BO approximation, we arrive to the Hartree approximation [11]. By expressing 2.21 as ψ(r) = N Y i ϕσ i(ri), (2.24) we can write an effective Hamiltonian for each single-particles wavefunctions ϕσ i(ri)with spin σ, Heffϕσ i(ri) = −1 2∇2+Vσ(i) eff (ri)ϕσ i(ri) = ϵσ iϕσ i(ri). (2.25) In expression 2.25, we can see that the potential has been factorized in an effective potential, Vσ(i) eff (ri), which acts on each electron of spin σat position ri. Due to the fundamental nature of Hartree approximation, this potential must incorporate the effects of all other particles. In Hartree approximation we assume that each electron interacts with every other electorn proportionally to the product of the densities of each pair as follows, ⟨H(e)⟩ψ=X i,σ⟨Heff⟩ϕσ i(2.26) ⟨H(e)⟩ψ=X i,σ⟨ϕσ i|−1 2∇2 i+Ve−I(ri)|ϕσ i⟩ +1 2X i,σX j,σ′=i,σ⟨ϕσ i,ϕσ′ j|1 |ri−rj||ϕσ i,ϕσ′ j⟩. (2.27) Now, lets solve this effective Hamiltonian by using the variational theorem (see App. B.1). To do this, we define a functional F(see App. B.2) with the corresponding Lagrange multipliers ϵi F[{ϕσ i},{ϵi}]=⟨H⟩ψ−X i ϵi X σ⟨ϕσ i|ϕσ i⟩−1!. (2.28) 2.3 density functional theory 31 HK ⇒ ⇒ ⇒ ⇒ Ψ0({r})Ψi({r}) n0(r)Vext(r) Figure 2.1: Schematic representation of Hohenberg-Kohn thorem. Starting from the upper left side, the potential Vext(r)determines (⇒) all the states Ψi({r})of the system, including the ground-state (⇒)Ψ0({r})and hence (⇒) the ground-state density n0(r). The Hohenberg-Kohn theorem states that Vext(r)is uniquely determined by n0(r), closing the circle. If we then add Eqs. 2.47 and 2.48, we arrive to E(1)+E(2)< E(1)+E(2), (2.49) which is a contradiction that renders our hypothesis incorrect. With this proof by contradiction, we can conclude that there is no two different potentials which give rise to the same non-degenerate ground-state charge density n0. Hence, the density uniquely determines the external potential within a constant. Theorem 2.3.2A universal functional for the energy E[n]in terms of the density n(r)can be defined, valid for any external potential Vext(r). For any particular Vext(r), the exact ground-state energy of the system is the global minimum value of this functional, and the density n(r)that minimizes the functional is the exact ground-state density n0(r). Corollary 2.3.2The functional E[n]alone is sufficient to determine the exact ground-state energy and density. Proof. The proof of the second theorem follows a similar structure to that one of the first, with some additional considerations. This proof, as the original from Hohenberg-Kohn, is restricted to V-representable densities, i.e. densities that correspond18 to the ground-state of an electron Hamiltonian with some external potential Vext. Since all properties are uniquely determined by n(r), they can be expressed as functionals of said density, including the total energy functional EHK[n] = T[n] + Eint[n] + ZdrVext(r)n(r) + EII =FHK[n] + ZdrVext(r)n(r) + EII, (2.50) 18 Can be derived as the ground-state density of a physical many-body system. 2.3 density functional theory 32 KS ⇔ {ϕi(r)} VKS(r) n0(r)HK0 ⇒ ⇒ ⇒ ⇒ Φ(r) HK ⇒ ⇒ ⇒ ⇒ Ψ0({r}) Ψi({r}) n0(r)Vext(r) Figure 2.2: Schematic representation of Kohn-Sham approach. HK0denotes Hohenberg-Kohn theorems applied to the non-interacting system. KS bidirectional arrow signals the relation between the many-body and single-particle systems, showing that the solution for the independentparticle Kohn-Sham problem determines the properties from the full many-body system. where T[n]and Eint[n]are the functionals for the kinetic and potential energies respectively and EII is the ion-ion interaction energy. In the second line of Eq. 2.50 we write FHK[n], a universal functional encompassing terms that are independent of the external potential. To prove the theorem, lets consider a system with ground-state density n(1)(r)corresponding to the external potential V(1) ext(r). Using the groundstate wavefunction Ψ(1), the Hohenberg-Kohn functional can be expressed as19 E(1)=EHK hn(1)i=⟨Ψ(1)| H(1)|Ψ(1)⟩. (2.51) Now, consider a different density n(2)(r)associated with a different wavefunction Ψ(2). Following the same logic as in the proof of the first theorem, we can reformulate Eq. 2.45 as E(1)=⟨Ψ(1)| H(1)|Ψ(1)⟩<⟨Ψ(2)| H(2)|Ψ(2)⟩=E(2). (2.52) This inequality shows that any density n(2)(r)different from the groundstate density n(1)(r)yields a higher total energy. Therefore, the global minimum of FHK[n]corresponds to the true ground-state density. Since the Hohnenberg-Kohn functional FHK[n]is universal and independent of Vext(r), the exact ground-state denisty and energy cna be found by minimizing the total energy with respect to variations in n(r). This result supposts Corollary 2.3.2, establishing that E[n]alone is sufficient to determine the exact ground-state energy and density. 19 Here, we omit the explicit dependence of non rto simplify the notation. 2.3 density functional theory 33 2.3.2The Kohn-Sham approach While the Hohenberg-Kohn theorems state that the ground-state electron density n0(r)determines all properties of a system, they do not provide a straightforward method for obtaining n0(r). The Kohn-Sham approach, depicted in Fig.2.2, tackles this problem by introducing an auxiliary system of non-interacting electrons, defined to have the same ground-state density as the interacting system. This recasts the problem into a set of self-consistent equations that are numerically solvable, while retaining the complexity of the real system through an exchange-correlation functional. This approach effectively reduces the many-body problem to an independent-particle one, making it computationally feasible. The construction of the Kohn-Sham auxiliary system is based in two key assumptions, 1. The exact interacting ground-state density, ne 0(r), is representable by the ground-state density of an auxiliary non-interacting system na 0(r). 2. The auxiliary Hamiltonian Haux has the usual kinetic operator and effective local potential Vσ eff(r), acting on an electron of spin σat r. Based on these assumptions, we can formulate the auxiliary Hamiltonian as20 Hσ aux = −1 2∇2+Vσ(r). (2.53) Since the auxiliary Hamiltonian, Hσ aux from Eq. 2.53 is non-interacting, the ground-state of the system is configured by placing one electron in each of the Nσorbitals ψσ i(r)with the lowest eigenvalues ϵσ i. Hence, the electron density of the auxiliary system is defined as n(r) = X σ n(r,σ) = X σ Nσ X i=1 |ψσ i(r)|2. (2.54) The the kinetic energy TSof the auxiliary system is given by TS= −1 2X σ Nσ X i=1⟨ψσ i|∇2|ψσ i⟩=1 2X σ Nσ X i=1Zdr|∇ψσ i(r)|2. (2.55) And the classical Coulomb interaction energy, resulting from the electron density interacting with itself, is EC[n] = 1 2Zdrdr′n(r)n(r′) |r−r′|. (2.56) 20 Using Hartree units, h=me=e=4π/ϵ0=1. 2.3 density functional theory 34 Within the Kohn-Sham approach, the full interacting many-body HohenbergKohn energy functional from Eq. 2.50 can be rewritten as EKS =TS[n] + ZdrVext(r)n(r) + EC[n] + EII +Exc[n], (2.57) where Vext(r)is the external potential owing to ions and external fields, EII is the interaction between ions, ECis the Coulomb energy defined in Eq. 2.56 and Exc is the exchange-correlation energy that encapsulates all many-body effects of exchange and correlation. The exchange-correlation energy can be expressed in terms of the Hohenberg-Kohn functional FHK[n], Exc[n] = FHK[n] − (TS[n] + EC[n]) = ⟨T⟩−TS[n] + ⟨Eint⟩−EC, (2.58) where n(r)is the density and depends upon position rand spin σ. The exchange-correlation energy can be expressed as the difference in kinetic and internal interaction energies between the interacting many-body system and the non-interacting auxiliary system, where electron-electron interactions are approximated by EC[n]21. Once we have an expression for EKS[n]and understand all the components in 2.57 we can develop the Kohn-Sham equations. The Kohn-Sham equations provide a framework to compute the ground-state density and energy of the system by solving an independent particle problem. These equations arise from minimizing EKS[n]with respect to the density n(r,σ)or the effective potential Vσ eff(r). The kinetic energy TSis expressed as a functional of the orbitals, while the remaining terms are functionals of the density. By applying the chain rule, the variational derivative of EKS[n]is ∂EKS ∂ψσ∗ i(r)=∂TS ∂ψσ∗ i(r)+∂Eext ∂n(r,σ)+∂EC ∂n(r,σ)+∂Exc ∂n(r,σ)∂n(r,σ) ∂ψσ∗ i(r). (2.59) From Eqs. 2.54 and 2.55, the partial derivatives for the kinetic energy and electron density are ∂TS ∂ψσ∗ i(r)= −1 2∇2ψσ i(r),∂nσ ∂ψσ∗ i(r)=ψσ i(r). (2.60) Using the method of Lagrange multipliers and imposing orthonormalization constraints for the orbitals, ⟨ψσ i|ψσ′ j⟩=δijδσσ′, we arrive to a Schrödingerlike equation −1 2∇2+Vext +VC+Vσ xcψσ i(r) = ϵσ iψσ i(r), (2.61) 21 Which is, in turn, the Hartree energy, see 3.2from Ref. [66]. 2.3 density functional theory 35 where Vσ KS(r) = Vext(r) + ∂Ec ∂n(r,σ)+∂Exc ∂n(r,σ) =Vext(r) + VC(r) + Vσ xc(r). (2.62) The Hamiltonian can then be expressed as Hσ KS(r)=−1 2∇2+Vσ KS(r), (2.63) and the Kohn-Sham equations take the compact form [Hσ KS −ϵσ i]ψσ i(r) = 0. (2.64) The Kohn-Sham equations reduce the many-body problem to an independent particle problem where the potential must be determined self-consistently with the electron density. If the exact exchange-correlation functional Exc[n] was known, solving these equations would yield the exact ground-state energy and density of the many-body system. However, since Exc[n]encapsulates all many-body correlation effects, which scale exponentially with the number of electrons, the exact Exc[n]cannot be determined explicitly for systems with more than a few particles [69,70]. Indeed, this would defeat our purpose of downscaling the problem’s computational complexity. Instead, approximations for Exc[n]are employed, forming the foundation of modern-day DFT implementations. We will explore these approximations and their application in computational methods in the following Sections. 2.3.3Local Density Approximation (LDA) One of the most simple and intuitive approximations, but also one of the most foundational ones is the LDA, first introduced by W. Kohn and L. J. Sham in 1965 [65]. LDA stands as LDA approximates the exchangecorrelation energy of an electron system by assuming that, at every point in space, the system locally resembles an homogeneous electron gas with density n(r). Along this manuscript, LDA is rarely used, however; we introduce it as the basis foundation for more developed approximations such as the presented in the following Sections. The exchange-correlation energy functional within LDA is expressed as ELDA xc [n] = Zdrϵhom[n](r)n(r), (2.65) 2.3 density functional theory 36 where ϵhom[n]represents the exchange-correlation energy per particle in a homogeneous electron gas with density n(r). We can further decompose it into exchange and correlation contributions, ϵhom xc [n](r) = ϵhom x[n](r) + ϵhom c[n](r), (2.66) where the exchange part can be expressed analytically as ϵhom x[n]=−3 43n π1/3 . (2.67) Meanwhile, due to its nature, the correlation part cannot be expressed analitically since it arises from complex many-body non-local electronic interactions that depend on the instantaneous positions and motions of all other electrons. Nonetheless, ϵc[n]can be approximated with Monte Carlo simulations [71] through forms of ϵcthat are typically fitted as functions of the average electron distance ϵ(rs)where rs=3 4πn1/3. The corresponding correlation potential is given by VC(rs) = ϵc(rs) − rs 3 dϵc(rs) drs . (2.68) Among the plethora of implementations of LDA, probably the best well known one is the one by Perdew and Zunger [72], which expresses ϵc(rs) as ϵPZ c(rs) =    −0.0480 +0.031 ln(rs) − 0.0116rs+0.0020rsln(rs)rs< 1 −0.1423/ 1+1.0529√rs+0.3334rsrs> 1. (2.69) Due to its reliance on the properties of a homogeneous electron gas, LDA performs well for systems where electrons exhibit such behavior, like metals and other bulk solids. However, for inhomogeneous systems like atoms or molecules, where the electron density varies significantly, LDA tends to be less accurate. Despite its limitations, LDA remains foundational in DFT and is widely used due to its simplicity and computational efficiency, serving as a starting point for more sophisticated approximations. 2.3.4Generalized Gradient Approximation (GGA) Building on the success of LDA, the development of GGA has lead to an increase of accuracy an realibility on exchange-correlation functionals. GGA incorporates not only the electron density n(r)but also its gradient ∇n(r) to describe exchange-correlation energy, allowing it to better account for inhomogeneities in the electron density. GGA is a very reliable approximation and it is the approximation of choice of most of the DFT calculations 2.3 density functional theory 37 presented in this Thesis due to its computational cost-accuracy ratio. The first step towards GGA was proposed by Kohn and Sham in their seminal paper [65] and later formalized by Herman et al. in Ref. [73], leading to the development of the Gradient Expansion Approximation (GEA). However, GEA did not represent a significant improvement over the LDA, as it introduced severe issues, such as sum rule violation. The GGA was developed to overcome these limitations by incorporating more sophisticated functionals. We can define the GGA functional as EGGA xc [n] = Zdrn(r)ϵxc(n(r),∇n(r)(2.70) ≡Zdrn(r)ϵhom x(n)Fxc(n,∇n), (2.71) where Fxc(n,∇n)is a dimensionless enhancement factor that depends only on nand its gradient, while ϵhom xc (n)is the exchange-correlation energy per particle of an homogeneous electron gas as happening in LDA. Among the numerous GGA implementations, the most widely used implementation of GGA (aswell as the most used in the work presented in this manuscript) is the one by Perdew Burke and Erzerhof (PBE) [74]. PBE was presented as a GGA improving the previously existing PerdewWang 1991 (PW91) [75] GGA, which was deemed as non-transparent, overparametrized and failing in several situations. PBE is constructed based on a sytematic approach that improves upon LDA by using exact physical constraints rather than fitting to experimental data like purely empirical functionals. Thanks to that, this GGA is more transferable across different systems and satisfies exact constraints imposed upon it as: correct asymptotic behaviour of Ex, correct scaling properties at low and high density regimes and complies with the Coulomb’s energy lower bound proposed by Lieb-Oxford [76], preventing unphysical Exbehaviour and Exchole normalization. The exchange part ot the PBE functional modifies the enhancement factor Fxc in 2.70, that we will denote as Fx, using the dimensionless reduced gradient s, defined as s=|∇n(r)| 2kF(r)n(r), (2.72) where kF(r)is the Fermi wavevector22 and its related to n(r)by n(r) = kF(r)3/(3π)2. (2.73) 22 The Fermi wavevector defines the radius of the Fermi sphere in momentum space. All electron states with wavevector k<kFare occupied at zero temperature while states with k > kFare unoccupied. In the homogeneous free electron gas limit the Fermi sphere coincides with the Fermi surface. 2.3 density functional theory 38 From 2.70,EGGA xis expressed as EGGA x[n] = ZdrFPBE x(s)n(r)ϵhom x[n], (2.74) where FPBE x=1+κ 1−1 1+µs2 κ!, (2.75) with κ=0.804 and µ=0.21951. Meanwhile, the correlation part can be expresed as a correction of LDA correlation as follows, EGGA c[n] = Zdrn(r)hϵhom c[n] + H(rs,t)i, (2.76) where tis a dimensionless gradient and rsis the average inter-electron distance and can be expressed as t≡∇n(r) 2ksn(r), with ks=4kF π2 . (2.77) Meanwhile, H(rs,t)encodes the correction over LDA; in the case of PBE it is HPBE(rs,t) = β2 2α ln1+2α β t2+At4 1+At2+A2t4, (2.78) with                      A=2α βexp−2αVhom c[n] β2−1 , α=0.0716,β=0.066725, t=|∇n(r)| 2ksn(r),ks=4kF π1/2 . Like LDA, GGA is most accurate when exchange-correlation energies are relatively small, such as in weakly correlated systems. This makes the GGA particularly suitable for studying systems like bulk metals, semiconductors, and simple molecules. For the systems explored in this manuscript, the GGA (and most concretely PBE) provides a reliable balance of computational efficiency and accuracy, making it the method of choice for exchangecorrelation approximations. PBE is generally chosen along this manuscript due to its balance between 2.3 density functional theory 39 accuracy and computational cost. It performs well for bulk materials, surfaces and interfaces and has a fairly accurate thermochemical and structural predictions (we also use it for relaxing structures). It has a decent band gap estimation when compared to LDA, however it tends to understimate band gaps due to its self-interaction error; that is why, along this manuscript we use the modified Becke-Johnson (mBJ) exchange potential [77,78] when trying to determine exact band gaps in Chapter 3. Moreover, beyond its applicability, PBE serves as the foundation for many improved functionals, like PBEsol [79], optimized for solids and densely packed systems; HSE06 [80], which is an hybrid functional that corrects PBE’s band gap understimation and vdW-DF [81], correcting its deficiencies in weak interactions by adding van der Waals corrections. 2.3.5A brief note on modified Becke-Johnson exchange potential As already mentioned in the previous Section, in Chapter 3we use mBJ to determine a more accurate bandgap. mBJ is not a GGA in the conventional sense, it is a semilocal exchange potential that goes beyond GGA but does not fuly reach the level of hybrid funtionals like the previously mentioned HSE06 or nonlocal exchange methods like GW approximation [82– 85] (used to obtain self-energies in many-body systems). It is so because it uses a dependence on the kinetic energy density, t(r), and a systemdependent mixing parameter Cthat improves the treatment of Ex, EmBJ x,σ(r) = cEBR x,σ+ (3c −2)1 πs5 6 tσ(r) nσ(r), (2.79) where tσ=1 2PNσ i=1ψ∗ i,σψi,σis the kinetic-energy density and vBR x,σ(r)=− 1 bσ(r)1−e−xσ(r)−1 2xσ(r)e−xσ(r), (2.80) is the Becke-Roussel potential [86], proposed to model the Coulomb potential created by the exchange hole. xσis determined from an equation involving nσ,∇nσ,∇2nσ, and tσ. Then bσis calculated with bσ= [x3 σe−xσ/8πnσ]1/3. cis a parameter chosen to depend linearly on the square root of the average of ∇n(r′) n(r′), c=α+β1 Vcell Zcell dr∇n(r′) n(r′), (2.81) where αand βare two free parameters and Vcell is the unit cell volume. This spatially dependent exchange potential allows mBJ to capture effects similar to those of GW approximation. 2.3 density functional theory 40 The mBJ potential effectively enhances the exchange potential in regions of low electron density and high density gradients. This is crucial because the underestimation of band gaps in LDA and GGA is largely due to an incorrect description of the screened exchange potential. Moreover, the mBJ correction increases the energy separation between occupied and unoccupied states, thereby improving the band gap prediction. This behavior is particularly relevant in semiconductors and insulators, where the conduction band minimum is often too low in LDA and GGA, while the valence band maximum is slightly too high. The mBJ potential appropriately shifts these bands, leading to better agreement with experimental band gaps. Beyond this, unlike hybrid functionals and the GW approximation, which require computationally expensive nonlocal exchange terms, the mBJ potential remains a semilocal functional while incorporating empirical fitting parameters. This allows it to closely match experimental data without a significant computational burden. In conclusion, the mBJ exchange potential significantly improves band gap predictions in semiconductors and insulators by introducing a nonlocal exchange-like correction through a kinetic-energy-dependent term. It offers superior accuracy compared to LDA and GGA and performs competitively with hybrid functionals while maintaining a relatively low computational cost. This balance between accuracy and efficiency is the primary reason we employ it for band accurate band gap calculations. Nonetheless, for general purposes, GGA+PBE is more than enough for an accurate systematic decription of the electronic structures of the systems under study of this Thesis. 2.3.6Concluding remark Until now, we have developed a framework for simulating the electronic structure of matter by employing the Hohenberg-Kohn theorems to derive the Kohn-Sham equations, as well as introducing approximations for the exchange-correlation energy in these equations. However, a practical implementation for solving the Kohn-Sham equations would still require the explicit treatment of all electrons (both core and valence) in the system under study, which can be computationally prohibitive for the systems considered in this Thesis. In the following Section, we will explore how to overcome this challenge by replacing the complex interactions of tightly bound core electrons with an effective potential, enabling a more efficient calculation of valence electronic states, key for our reseach. 2.5 bloch theorem and plane-wavebasis 47 A translation operator TRcan then be defined to shift the system by a lattice vector Ras TRf(r) = f(r+R). (2.98) A crystalline solid is invariant under the action of any primitive translation operator Taior their combinations26. Because these translation operators commute with the Hamiltonian, the Hamiltonian eigenstates can be also labeled by the eigenvalues of the translation operator. These eigenvalues take the form eik·R, where kis the crystal momentum, a triplet of phases. The implications of this effective periodic treatment are captured in the Bloch’s theorem, Theorem 2.5.1The solutions to the Schrödinger equation in a periodic potential can be expressed as plane waves modulated by periodic functions u(r): ψ(r) = eik·ru(r), (2.99) where u(r)is a function periodic with the lattice. When applied to our Kohn-Sham equations, we can write the Kohn-Sham states ϕias ϕσ n,kwhere iis now decoupled into k, associated to a particular plane-waveeik·rand n, associated to a particular periodic modulation un,k(r), ϕnk(r) = eik·run,k(r), (2.100) where kis defined in the first Brillouin zone (IBZ) and n is called band index. Since any periodic function can be expanded in a complete set of Fourier components, we can express the periodic function un,k(r)in terms of reciprocal lattice vectors G un,k(r) = X G cn,k+Gei(G)·r(2.101) where the plane-waveexpansion naturally leads to a Fourier representation for the Kohn-Sham eigenfunctions, ϕnk=X G cn,k+Gei(k+G)·r. (2.102) In order to prevent energy divergences from the infinite number of Gcomponents, we introduce an energy cutoff Ecut, discarding all components with Ecut ⩽1 2|k+G|2. (2.103) 26 Since translation operators commute. 2.5 bloch theorem and plane-wavebasis 48 This truncation reduces the computational cost and ensures the convergence of calculations, usually is refered in a DFT environment as the kinetic cutoff. With this approach, the Kohn-Sham equations become X G′ HG,G′(k)cn,k+G′=ϵi(k)cn,k+G, (2.104) where the Hamiltonian matrix elements are HG,G′(K)≡ ⟨k+G|HKS|k+G′⟩=1 2|k+G|2δG,G′+VKS(G−G′). (2.105) In Eq. 2.104,ihas been split into IBZ wave-vector kand band index n, which allows us to construct bandstructures from solving this Hamiltonian, by labelling each state by the wavevector k. This crystal momentum is conserved analogously to ordinary momentum in free space modulo the addition of any reciprocal lattice vector Gand hence, solutions for the Hamiltonian are periodic in kand unique when expressed in the IBZ of the reciprocal lattice, allowing bands to build a complete description of the electronic behaviour of the system. Furhtermore, the kdescription allows us to integrate properties such as the number of electrons, the total energy etc. For example, for a function fn(k), where ndenotes the discrete band index, the average value is ¯ fn=1 NkX k fn(k)→Ωcell (2π)3ZBZ dkfn(k), (2.106) where Ωcell is the volume of a primitive cell in real space. The DOS per unit cell can be expressed as N(ϵ) = 1 NkX nX k δ(ϵ−ϵn,k) =Ωcell (2π)dZBZ dkδ(ϵ−ϵn,k). (2.107) As DOS, any quantity can be integrated over the IBZ, like optical response tensors27. This Section is one of the most important to our Thesis, thanks to it we can obatin band structures, DOS and a plethora of results from our DFT calculations. In the following, we will ennumerate all the ways in which DFT can reduce computational cost. 27 There is a caveat to this: eventhough the plane-wave basis framework allows for calculation of first and second order dielectric responses, usually said calculations are not performed in a plane-wave basis since the k−grid is rather sparse due to computational limitations. Usually, said calculations are performed in a Wannier basis that allow for denser meshes without a heavy computational cost. We go through this method in Sec. 2.7. 2.6 how dft reduces computational cost 49 2.6 how dft reduces computational cost In the previous Sections of this Chapter we have understood why and how DFT reduces the computational complexity of solving the many-body Schrödinger equation by several interrelated mechanisms. In the following, we will ennumerate and go through them as a form of summary: 1.Reduction to a single-particle framework with the Kohn-Sham approach: The Kohn-Sham formulation of DFT replaces the full manybody interacting electron problem with a system of non-interacting electrons moving in an effective potential, Veff(r). This effective potential encapsulates the effects of electron-electron interactions through the exchange-correlation energy functional, Exc[n], which depends only on the electron density n(r). And, as a result: • The computational cost will scale polynomially with the number of electrons N, typically O(N3); instead of exponentially. • Solving the Kohn-Sham equations involve iterative self-consistent field methods, which are feasible computationally for systems with hundreds to thousands of atoms. 2.Dimensionality reduction: Through the Hohnenberg-Kohn theorems, we can establish that all properties of the electronic systems can be determined uniquely by the ground-state electron density, which is a3D scalar function, independently of the number of electrons in the system. This will dramaticly reduce the dimensitonality from the 3N-dimensional wavefunction into a 3D density that enables tractable calculations for systems with many electrons. 3.Using approximate functionals: The exact exchange-correlation functionals are not known, but using practical approximations like LDA or GGA provide a reasonable accuracy for a wide range of systems while remaining computationally efficient. With these approximations DFT avoids the direct calculation of correlated electron wavefunctions and replace the complex electron-electron interactions with manageable analytically or semi-empirical expressions. 4.Exploting translational symmetry: In periodic systems, DFT uses the Bloch’s theorem to represent electronic states as Bloch functions, reducing the problem to a single unit cell and sampling in reciprocal space. By using translational symmetry we reduce the size of the todiagonalize-Hamiltonian and limit all our computations to the IBZ. 5.Basis set simplifications: By using plane waves or localized basis sets, combining with energy cutoffs, we allow the truncation of highenergy components in the solution, drastically reducing the number 2.7 wannier functions 50 of basis functions needed. This truncation is a middle-ground between accuracy and computational efficiency. 6.Reducing the number of electrons by eliminating core electrons with Pseudopotentials: Pseudopotentials replace the explicit treatment of the bound core electrons with an efectiv epotential acting only on valence electrons. By this approach we reduce the number of wavefunctions to be calculated and simplify tehe potential landscape 28, improving numerical stability and reducing the basis set size requeriments. 7.Efficient algorithms and parallelization: Moder DFT implementations make use of a plethora of advanced numerical methods due to its nature. Between these methods we have sparse matrix methods for large systems, iterative diagonalization schemes and efficient fast fourier transforms for plane-waveexpansions. These methods allow DFT to scale efficiently for large and complex systems. By transforming the many-body problem into a series of single-particle equations and employing approximations, DFT achieves an exceptional balance between accuracy and computational cost. While some challenges remain, such as the exact treatment of exchange-correlation effects and strongly correlated systems, DFT’s efficiency and flexibility have made it the primary tool for investigating the electronic structure of matter. The framework laid out in this Chapter forms the foundation for the methods and approximations explored throughout this Thesis. 2.7 wannier functions In the study of periodic systems, Bloch’s theorem provides a convenient representation of the electronic states as extended plane waves modulated by periodic functions. However, the delocalized nature of Bloch functions blurs a possible real-space interpretation of electronic properties. To this end, Wannier functions, introduced by G. Wannier in 1937 [15], provide a localized real space basis that is orthonormal and spans the same subspace as the Bloch eigenstates for one or several bands. These functions are particularly advantageous for constructing simplified models, analyzing chemical bonding and studing transport phenomena. Moreover, by focusing on specific sets of bands, Wannier function can facilitate the reduction in computational complexity of electronic structure problems, such as the calculation of optical response tensors as in Sec. 3.2.5and a plethora of other 28 The effective potential replacing core electrons is much smoother when compared to the actual Coulomn potential near the nucleus since it does not account for the high-frequency oscillations of core electron. 2.7 wannier functions 51 properties. In the following, we introduce the concept of Wannier functions and their applications, both in our research and in the broader field, highlighting their crucial role in computational condensed matter physics. From the Bloch theorem, we can write the wavefunctions of a periodic Hamiltonian Has ψn,k(r) = eik·run,k(r). (2.108) These Bloch functions are only defined up to a phase factor, allowing for a gauge tranformation of the form ψn,k(r)→˜ ψn,k=eiθn(k)ψn,k(r), (2.109) where θn(k)is a differentiable phase function. The Wannier function for band n, associated with the lattice vector Rm, is defined as the Fourier transform of Bloch states over the Brillouin zone, ωn(r−Rm) = Ωcell (2π)3ZBZ dke−ik·Rmψn,k(r), (2.110) where Ωcell is the unit cell volume and these functions decay exponentially as |r−Rm|→∞, ensuring their localization. The original Bloch functions can be reconstructed from the Wannier functions via the inverse transform, ψn,k(r) = X m eik·Rmωn(r−Rm)(2.111) 2.7.1Generalized Wannier Functions for Multiple Bands Wannier functions can also be constructed for a group of bands by forming linear combinations of Bloch states. For Nbands, the periodic part of the Bloch function, un,k, can be rotated using a unitary matrix U(k), uW i,k=X n Un,i(k)un,k. (2.112) The resulting Wannier functions are given by ωi(r−Tm) = Ωcell (2π)3ZBZ dkeik·(r−Tm)X n Un,i(k)un,k(r)(2.113) These generalized Wannier functions retain orthonormality and provide a flexible basis for capturing electronic properties. 2.7 wannier functions 52 2.7.2Non-Uniqueness and Maximally localized Wannier functions Due to the gauge freedom in choosing θn(k), the Wannier functions are not uniquely defined. Indeed, the choice of gauge can significantly influence their shape and localization. In order to address this, we can construct maximally localized Wannier functions (MLWFs) by minimizing the spread functional Ω=X ih⟨r2⟩i−⟨r⟩2 ii, (2.114) where ⟨r⟩iand ⟨r2⟩iare the expectation values of position operators for the i-th Wannier function. This minimization ensures a consistent and physically meaningful representation of localized states. The pursuit of MLWFs addresses the arbitrariness in gauge choice while offering the additional benefit of simplifying calculations. Their localization reduces the overlap between neighboring functions, resulting in sparser and more computationally efficient Hamiltonians. MLWFs are particularly well-suited for constructing tight-binding models, calculating Berry phases, and determining polarization. By using maximally localized Wannier functions, spurious oscillations and long-range interactions are minimized, enabling precise real-space analyses and providing a clearer physical understanding. This is because MLWFs yield well-defined centers of charge, enhancing their interpretability and making them a valuable tool for exploring electronic structure. For example, the code employed in Ch. 3;Wannier90 [98], is based on finding MLWFs. 2.7.3Applications and Wannier interpolation Wannier functions are key in constructing tight-binding models, calculating transport properties and studying topological invariants. Additionally, their localization enables Wannier interpolation, allowing electronic properties to be computed on finer k−space grids without the need for full diagonalization of the Hamiltonian. The workflow is as follows: 1.Initial DFT calculation: Kohn-Sham equatinos are solved in a coarse k-point grid 29 in order to obtain eigenvalues and eigenfunctions at discrete point in the reciprocal space. 2.Construction of Wannier Functions: Bloch wavefunctions are transformed into Wannier functions using unitary tranformations. After that, the Wannier functions are maximally localized by minimizing the spread funcitonal. 29 Leading these calcualtions to be computationally inexpensive. 2.7 wannier functions 53 3.Interpolation of Hamiltonian in reciprocal space: In reciprocal space, the electronic Hamiltonian is represented in Wannier basis as Hmn(k) = X R eik·R⟨ωm(r)|H|ωm(r+R)⟩, (2.115) where ⟨ωm(r)|H|ωm(r+R)⟩decay rapidly due to MLWFs being spatially localized. This allows the Hamiltonain to be efficiently interpolated from a coarse grid to an arbitrary k-point grid in the BZ 30 4.Applications: This interpolated Hamiltonian can be diagonalized on an arbitrary set of k-points without aditional ab-initio calculations. This has a plethora of applications such as high-resolution bandstructures, Fermi surfaces, DOS, Berry-phase properties and optical responses as in Sec. 3.2.5. Another great application of Wannier functions is that they allow the systematic construction of tight-binding models on a localized basis that accurately reproduce the band structure of the material in a minimal basis set. 2.7.4Why and when using Wannier interpolation Wannier interpolation method and its computational application through Wannier90 [98] is of paramount importance in modern computational condensed matter research. It builds upon DFT as both a theoretical framework and a practical tool, enabling the calculation of a wide range of properties on dense, high-resolution k-grids. For instance, in Ch. 3, we compute the optical conductivity of (TaSe4)3I using a 91 ×91 ×71 k-mesh, whereas the underlying DFT calculations are performed on a significantly coarser 11 ×11 ×5 k-mesh. Directly computing the optical conductivity on such a dense grid using conventional DFT would be computationally prohibitive. To illustrate the efficiency gain, consider a naïve estimation: assuming an ideal parallelization scheme with no overhead and the same kinetic energy cutoff (i.e., the same basis set) for both cases, let us assume that the self-consistent DFT calculation on the 11 ×11 ×5grid takes 1hour to complete. Given the polynomial scaling with Nkand an assumed sublinear scaling efficiency of 80% due to parallelization, a self-consistent calculation on a 91 ×91 ×71 k-grid would take approximately 3days, disregarding practical limitations such as RAM shortages or paralellization overhead. In contrast, in our experience using Wannier interpolation, the optical conductivity calculation required only 1.5times the computational time of the self-consistent DFT calculation, leading to an overall speed-up of a factor of 50. Moreover, this estimate considers only the self-consistent step; computing optical conductivity also requires evaluating the derivatives of the orbitals with respect to k, which in a direct DFT approach would require 30 This process can be furtherly optimized using algorithms as the fast Fourier Transform. 2.7 wannier functions 54 an additional computationally expensive step, further increasing the total runtime. Even though the methods presented in this Section offer immense computational power and significant computational cost reduction, it is important to discuss why, for example, we do not use them to model charge transfer in Chapter 4. The primary challenge when working with Wannier functions lies in the appropriate choice of basis. The selection of a basis in Wannier function calculations is a user-provided input that requires both expertise in the physical system under study and a deep understanding of the computational implementation31. A poorly chosen basis, combined with over-fitting through a long localization process, can lead to a model that lacks meaningful physical insight, ultimately resulting in incorrect predictions. For instance, the process of maximal localization itself can break crucial crystal symmetries, which are essential for understanding phenomena such as optical responses or Berry phase properties32. Assessing whether a Wannierization is well-executed is highly nontrivial and requires deep knowledge of both the system under study and the Wannierization process itself. A particular example of the issues that can arise when using MLWF are the problems for interpreting material properties through Wannierinterpolated tight-binding models. While these models are built using localized Wannier functions (which often originate from initial trial functions that are orbital-like), any orbital-resolved or symmetry-dependent information about the reconstructed bands may become obscured due to the disentangling and Wannierization procedures, particularly when maximal localization is applied. In many cases, to preserve orbital character and symmetry properties, a one-shot Wannierization approach is recommended (see Sec. II.I.1of [99]). Despite these nuances, the potential pitfalls of Wannierization should not be seen as prohibitive but rather as aspects to be carefully managed. A well-chosen basis that respects the symmetries of the system, or the application of symmetry-constrained maximal localization, can mitigate many of the previously mentioned issues. Furthermore, in some cases, the apparent breaking of symmetry may itself reveal hidden physical phenomena, making Wannierization an important diagnostic tool in condensed matter physics. In conclusion, Wannier interpolation and, more specifically ML31 In this text, we consider only Wannier90 [98]. 32 The gauge choice in Wannierization may not always preserve the full symmetry of the original Bloch wavefunctions. For example, if a material exhibits inversion symmetry in its Bloch states but the localization procedure selects an asymmetric Wannier gauge, the final Wannier functions may not respect inversion symmetry. This can significantly affect Berry phase calculations, electronic polarization, and other topological properties. 2.7 wannier functions 55 WFs, constitutes a powerful tool for condensed matter research. However, the development of systematic frameworks based on MLWFs remains a challenge due to the required careful selection of basis functions. To address this, efforts have been made toward automatic Wannierization procedures [100], though we consider them to still be under development. Part III RESULTS 3.1 introduction 63 those found in CDWs. Moreover, understanding CDW transitions is crucial for interpreting charge transport measurements and exploring unconventional ordering phenomena in layered materials. In later Sections, we will revisit these concepts in the context of specific materials and their interplay with other collective quantum states. 3.1.4TMCs (MSe4)nIwith M=Nb,Ta. The series of TMCs (MSe4)nIwith M=Nb,Ta; are made of weakly coupled MSe4chains, each of them hosting a dz2-derived 1D band with fractional filling δ= (n−1)/2n. Reported structures occur at n=2,3,10 3i.e. δ=0.25,0.33,0.35 respectively. These fractionally filled quasi-1D bands are prone to CDW instabilities which have been extensively studied [17,18, 101,103,104,111–119]. Among the Nb-based compounds, the n=2and n=10 3variants exhibit clear incommensurate CDW transitions originating from a metallic parent state. These transitions display well-established CDW phenomenology, including non-linear transport due to CDW depinning, reflecting the strong coupling between the electronic and lattice degrees of freedom. In contrast, the n=3compound, (NbSe4)3I, stands out as the only member of the series that does not exhibit a conventional CDW. Instead, it undergoes a structural transition and shows activated semiconducting transport. A key to understanding this behavior lies in its unusual electronic configuration: the material features a mixed-valence state among the Nb atoms. Assuming a simplified ionic model with Se2−and I−, the average Nb oxidation state is approximately +4.33, implying a periodic coexistence of Nb4+(4d1) and Nb5+(4d0) ions. This leads to a site-selective distribution of d-electrons along the quasi-one-dimensional chains of facesharing NbSe6octahedra. While the Nb4+sites contribute localized d1electrons prone to Peierls instabilities and gap formation, the Nb5+sites act as electronically inert spacers. The resulting periodic modulation of the electronic density mimics the effects of a Peierls distortion, but without requiring a purely electronic instability, as presented in the introduction to this Chapter. Despite this, several puzzles have precluded a consistent understanding of the n=3compound within simple band theory: multiple structural transitions have been reported [103,104,114,115], transport gaps vary significantly across studies [18,112,120–122], and ARPES and optical conductivity experiments [123,124] reveal a gap much larger than those observed in transport, suggesting the presence of hidden electronic structure and low-spectral-weight states not captured by conventional models. For the isoelectronic Ta compounds, only n=2,3variants are reported in the literature [16,120,124]. While (TaSe4)2I has seen a renewed interest [33] in the context of axionic CDW in Weyl semimetals [34,35], the knowledge about (TaSe4)3I is rather scarce, as it was assumed to behave mostly 3.1 introduction 64 like its Nb counterpart [16,120,124]. In a recent development, however, a polytype of (TaSe4)3I has been found to be metallic at room temperature, with a ferromagnetic transition at 8K, and superconductivity coexisting with ferromagnetism at Tc=3K [36]. This coexistence is unusual on its own [37,38], but it is all the more surprising giving the isoelectronic Nb compound semiconducting behaviour. The conflicting transport and ARPES results for (NbSe4)3I, along with the surprising low temperature behavior of (TaSe4)3I reveal that an understanding of these materials in terms of a band structure picture were lacking. Along the next pages of this Section we develop a theoretical and computational framework for understanding the previous literature and further develop this understanding experimentally. Firstly, we combine ab-initio calculations of all the known structures for both compounds along a tightbinding analysis to explain many of the electronic features of these systems. All phases are found to be semiconducting, with gaps that directly correlate with the amount Nb or Ta trimerization distortion6. We discuss how these gaps compare with the transport observations and generally find good agreement. However, we also show that while the gap is formally direct and located at Γ, there is an approximate in-plane translational symmetry that leads to negligible spectral weight for the valence band edge, so ARPES and optical conductivity actually probe higher bands and lead to a misidentification of the true gap. To conclude, we identify a spin-split Van Hove singularity in the valence band edge in the new (TaSe4)3I polytype which crosses the Fermi level for a very small amount of hole doping, and we argue this may provide a consistent explanation for ferromagnetism and potentially fluctuation mediated triplet pairing. 3.1.5Phenomenology overview on (MSe4)3I(M=Nb,Ta) TMCs. For the n=3compounds, (NbSe4)3I is the best studied one. Its high temperature phase has a tetragonal, centrosymmetric crystal structure with space group (SG) P4/mnc (No. 128, point group D4h). As temperature is lowered, it undergoes a structural phase transition at Tc1 ∼274 −280 K with the resulting space group P¯ 421c(No. 114, point group D2d). A second structural transition to a phase with SG P¯ 4(No. 81, point group S4) at Tc2 ∼90 K has also been reported [101,104,114,115], but not in all experiments. In this work, we will refer to these structures only by their point group. The pattern of symmetry breaking is thus D4h →D2d →S4. Resistivity measurements are generally consistent with semiconducting behavior with activated resistivity ρ∝eEg/2kBTand a gap Egthat appears to 6See Eq. 3.5. 3.1 introduction 65 change through the phase transitions. Above Tc1, values of Eg=190-220 meV [16,119] were reported. A resistivity kink was always observed at Tc1, and below it two types of samples were reported to exist [16], initially indistinguishable by structural measurements [16,112]. In type I samples Egis reduced to values in the range 20-70 meV [16,18,112,121] all the way to lowest temperatures measured. In type II samples, a broader kink is observed which leads to low temperature gaps of 110-130 meV [16,119]. However, in some samples ∂lnρ/∂T−1never flattens to a constant Egvalue, but rather continuously decreases after a maximum [112,122], challenging the view of standard activated transport. Differing observations also include a report of Eg=97 meV for T > Tc1 and 222 meV for T << Tc1 [115], or even a an abrupt change at 180 K from 345 meV to 22 meV [113]. Another characteristic feature of type II samples is that switching to a state of lower resistance can be induced at high currents and low temperatures [122]. This behavior disappears at 140 K, and is completely absent in type I crystals. A partial solution to these transport puzzles was offered in Ref. [119] which reported that the second structural transition was only observed in samples assigned to type II in transport. The overall suggested picture would then be that all three structures are semiconductors: the initial gap of 190-220 meV for the D4h structure is reduced to 20-70 the D2d structure, and for samples with the second transition it increases again in the S4structure to 110-130. While outliers to this picture do exist, we consider this to be the average behaviour to compare with our calculations. Early ARPES experiments at 300 K also attempted to determine the spectral gap [123,125]. The valence band maximum was found to be at 750 meV below EF, providing a lower bound on the gap which is much larger than the one obtained from every transport experiment. Optical conductivity also showed a raising edge at ω∼500 meV [123], again too large compared to transport gaps. The fact that the gap derived from ARPES and optical conductivity is much larger than the ones obtained from transport remains an unsolved problem to date, precluding any consistent band-structure understanding of these materials. Finally, much less is known about (TaSe4)3I [16,120,124], which has generally been assumed to behave like its Nb counterpart. The D4h →D2d structural phase transition was measured at at 200 K [124], but the S4phase has not been reported. Very recently, the D2d phase was reported to be metallic with a resistance kink-plateau at 150 K, a ferromagnetic transition at 8K, and a superconducting one at Tc=3K [36], in stark contrast to previous observations. The coexistence of ferromagnetism and superconductivity is a rarely reported phenomenon [37,38], and it is often taken as a hint that pairing could be in a unconventional odd-parity triplet channel. While experiments on superconductivity are at a very early stage, a band structure 3.2 methods &theoretical results for (MSe4)3Itmcs.66 Figure 3.2: (a) Atomic structure of (TaSe4)3I in the D4h phase and (b) Brillouin zones for the atomic structure of in (a) shown in blue, and for the approximate effective structure described in the text which contains one formula unit understanding of the basic properties of (TaSe4)3I is clearly needed as a starting point to understand these unusual behavior. 3.2 methods &theoretical results for (MSe4)3Itmcs. In the present Section we will present our theoretical and computational framework and methodology along the results obtained by applying said methodology. These results will serve as theoretical support to merge the past research presented in Sec. 3.1.5and lay ground for the posterior experimental work presented in Sec. 3.3. 3.2.1DFT details and methodology With the aim of explaining the previously discussed phenomenology, firstprinciples DFT band structure calculations for the different reported structures were performed using VASP [67,68]v.6.2.1with projector-augmented wave pseudopotentials using two approximations: the generalized gradient approximation with PBE parametrization [126] and the mBJ mehtod [78]. The PBE functional serves as a reliable and computationally unexpensive baseline for structural and qualitative electronic properties, while the mBJ potential is known to yield improved band gap estimations7by bet7LDA and GGA significantly underestimate band gaps since they do not fully capture the discontinuity in the exchange-correlation potential. Moreover, they suffer from selfinteraction errors which artificially delocalizes electrons and reduces the effective bandgap. By doing so, mBJ introduces a semilocal correction to the exchange term and empirical parametrization. mBJ effectively mimics the behaviour of the exact exchange potential for 3.2 methods &theoretical results for (MSe4)3Itmcs.67 ter approximating the exact exchange potential. We employ both methods in order to explain the possible metallic behaviour of (TaSe4)3I in the D2d phase and get accurate bandgaps for the different phases. Fig. 3.2shows the crystal structure of (MSe4)3I compounds in the D4h phase. The MSe4 units form one dimensional chains, where each M atom is sandwiched between Se4rectangular units and separated by I−ions. Each Se4adjacent rectangular unit is rotated around 45◦and modulated by small distortions on Se positions. The typical unit cell for these compounds is tetragonal (a=b=c) and contains two metallic chains. The particular case of the Fig. 3.2is (TaSe4)3I in D4h phase, with a=b=9.719 Å and c=19.363 Å. The difference between phases i.e D4h,D2d and S4, is given by the relative distance between M atoms and small modulations in the Se positions. For (NbSe4)3I, a detailed structural characterization is available for the three considered phases D4h [17], D2d [101], and S4[104]. Therefore, all DFT calculations for these compounds were performed using the experimentally measured lattice parameters. Self-consistent calculations considering spin-orbit coupling were found to be well converged for a kinetic energy cutoff of 520 eV and a 9×9×5k-mesh sampling. Conversely, in the case of (TaSe4)3I, no detailed characterization exists. Given the similarity with the Nb compounds [101], we obtained the structures of the T compounds performing a structural relaxation starting with the positions of the experimentally measured Nb-based compounds [18,127]. The relaxation calculations were performed using a conjugate gradient algorithm [128] and were found to be well converged for a kinetic energy cutoff of 520 eV and a 11 ×11 ×5k-mesh sampling. For phases D4h and D2d, the starting point were the pre-relaxed structures in Ref. [127] while, for phase S4, a direct relaxation with the conjugate gradient algorithm was performed from the original positions of the isoelectronic Nb compound in the same phase since no pre-relaxed data was available. This strategy prioritizes the conservation of the space-group symmetries from the original Nb structures. Self-consistent calculations considering spin-orbit coupling were found to be well converged for a 520 eV kinetic energy cutoff and a 11 ×11 ×5kmesh. Density of states calculations were performed using a 15 ×15 ×9 k-mesh and a energy resolution of 0.7meV. 3.2.2DFT band structures. The electronic band structures for (NbSe4)3I and (TaSe4)3I are presented in Fig. 3.3in both PBE and mBJ approximations. The band structures for both Nb and Ta compounds are overall similar for all phases. A set of 8low energy bands is observed in the energy window E∈[−0.5,0.5]eV. This set semiconductors and insulators and it is cheaper computationally than hybrid functionals or GW approximation. 3.2 methods &theoretical results for (MSe4)3Itmcs.68 Structure a(Å)c(Å)1 6−dMi−Mi+1 c2 (×10−4)∆d Gap PBE (eV) Gap mBJ (eV) Nb D4h 9.4891 19.1323 0.11 -0.11 -0.46 -0.11 -0.11 -0.43 1.15 0.09 0.26 D2d 9.4500 19.0799 0.00 -0.47 -0.39 -0.00 -0.47 -0.39 1.32 0.18 0.31 S49.4365 19.0461 1.07 -0.17 -0.12 -0.58 -0.00 -1.01 1.71 0.28 0.40 Ta D4h 9.7192 19.3626 0.09 -0.09 -0.37 -0.09 -0.09 -0.37 1.05 0.11 0.27 D2d 9.4365 19.4365 0.23 -0.00 -0.26 -0.23 -0.00 -0.28 1.00 0.06 0.19 S49.4365 19.0461 0.00 -0.39 -0.26 -0.01 -0.33 -0.30 1.13 0.12 0.24 Table 1: Lattice parameters, metal-metal distances, distortion ∆d and bandgaps for PBE and mBJ approximations and hybrid potential for the different structures considered in the text. 3.2 methods &theoretical results for (MSe4)3Itmcs.69 is found closer to the valence bands in the Nb compounds compared to the Ta compounds. These bands have a dominant orbital weight coming from dz2orbitals of the Ta atoms, as shown in the density of states in Fig. 3.4 and as explained originally by Gressier et. al. [103,120]. The fact that these bands have dominant weight in a single Ta orbital suggests a simple tight binding model should describe these bands correctly, as shown in the following Section. Regarding the band gaps, in the case of (NbSe4)3I, the mBJ method shows a greater band gaps than the PBE parametrization for the three phases. In PBE, the band gap is slightly indirect, with the valence band maximum slightly off Γin the M−Γdirection and the conduction band minimum at Γ. In mBJ, the band gap is also indirect for D2d and S4phases, but it is direct for D4h. Similarly, for (TaSe4)3I the mBJ method also shows a greater band gap than PBE for all phases. Both in mBJ and PBE approximations the band gaps are direct between Γ−Γfor inversion conserving phase D4h and slightly indirect for inversion-broken phases D2d and S4. In order to understand the different values for the obtained from DFT calculations we recall that the relative distances between Nb atoms were reported to be related to the electronic band gap in the past. In the interest of gaining insight of this statement for (MSe4)3I compounds we introduce the concept of trimerization distortion8∆d, ∆d =v u u t 6 X i=11 6−dMi−Mi+1 c2 (3.5) where dMi−Mi+1=ˆ diis the distance between neighbouring M atoms along the cdirection, c=|c|and the sum runs over the 6M atoms in the (MSe4)3I compounds. The quantity ∆d measures how far a 6-atom chain is from being evenly spaced, and vanishes when all distances between M atoms are equal, i.e. ˆ di=ˆ djfor all i,j. A summary of relevant structural data and distortions, along with the band gaps obtained both in the PBE parametrization and mBJ method are presented in Table 1. The results in the table suggest an relation between distortion, cell volume and electronic band gap. In general terms, a greater distortion leads to wider electronic band gaps. Even so, the cell volume also seems to plays a role in enhancing or reducing the band gap, since smaller cell volumes lead to smaller 8Trimerization refers to a structural distortion in the quasi-1D chains of the TMCs, where the spacing between adjacent metal atoms along the chain becomes unequal. This distortion creates a periodic alternation in the distances between neighnoring metal atoms, effectively grouping them into sets of three within the unit cell, hence the name. Trimerization affects the electronic structure by opening a band gap in the material. It is a key feature distinguishing the semiconducting behavior of these compounds and plays a critical role in their transport and optical properties. 3.2 methods &theoretical results for (MSe4)3Itmcs.70 band gaps. The data recollected in Sec. 3.1.5for transport band gaps in Nb-compounds, shows the following trend: ∆ED4h > ∆ES4 > ∆ED2d, while our data suggests that ∆ES4 > ∆ED2d > ∆ED4h. Even though the trends are dissimilar, the magnitudes for the DFT electronic gaps are comparable to the ones obtained with transport measurements. 3.2.3Tight-Binding models Since the orbital-resolved density of states in Fig. 3.4shows that the bands near the Fermi level are dominated by M dz2orbitals, we can gain a deeper understanding of the structural dependence of the band gap from a simple tight-binding model containing only such orbitals in Ta sites. The model will be constrained by the symmetries of each phase: D4h is generated by inversion I, the glide {m110|1 2 1 2 1 2}, the rotoinversion S4and the two-fold screw {2100|1 2 1 2 1 2}. Breaking of inversion then leads to the point group D2d and further breaking the glide and two-fold axis leads finally to point group S4. In the phase of higher symmetry, D4h, there are only two non-equivalent M atoms, as shown in Fig. 3.5a) and b). The simplest model therefore only contains two different on-site potentials ∆1and ∆2, two intrachain hoppings t1 and t2and a single interchain hopping t⊥. As the figure shows, this model actually has an accidental translation symmetry, with a reduced unit cell with 3M sites, representing a single formula unit (MSe4)3I (the true unit cell would have 4formula units). The lattice parameters of this reduced unit cell are c∗=c 2and a∗=a √2, and the corresponding enlarged Brillouin Zone BZ∗is shown in Fig. 3.2, with high symmetry points denoted as A∗, Z∗and so on. The Hamiltonian for this model is H= −     ∆1+t⊥f(k∥)t1eikzc/6 t3e−ikzc/6 t1e−ikzc/6 ∆2+t⊥f(k∥)t2eikzc/6 t3eikzc/6 t2e−ikzc/6 ∆3+t⊥f(k∥)     (3.6) where t1=t2and f(k∥) = 2cos(kx+ky)a 2√2+cos(kx−ky)a 2√2. To account for the trimerization distortion, we set t2=T−δ/2 and t3=T+δ so that the average hopping is tand δparametrizes the trimerization distortion assuming that the hoppings change linearly with the distance between orbitals. The presence of the extra translation symmetries in this model is accidental, as the inclusion of further neighbor intraand interchain hoppings would 3.2 methods &theoretical results for (MSe4)3Itmcs.71 Figure 3.3: Ab-initio bandstructures for the (MSe4)3I compounds in the three different structures D4h,D2d and S4from left to right. Top row (a-c) shows (NbSe4)3I, bottom row (d-f) shows (TaSe4)3I. Solid back lines were obtained using PBE pseudopotentials and dashed blue lines were obtained using mBJ hybrid potentials. 3.2 methods &theoretical results for (MSe4)3Itmcs.72 −0.6−0.4−0.2 0.0 0.2 0.4 0.6 Enegy [eV ] 0 2 4 6 8 10 12 14 DOS a.u. Total DOS Ta d2 z Se px Se py Figure 3.4: Orbital resolved density of states for (TaSe4)3I in the D2d phase computed within the PBE parametrization. This density of states corresponds to the bands shown in Fig. 3.3(e). The bands near the Fermi level have dominant dz2orbital character on the Ta sites. This is true for all computed band structures (not shown). indeed require the use of the full 12 site unit cell. This further neighbor hoppings are however expected to be smaller, so that the model in Eq. 3.6 serves as a good first approximation for the band structure. Physically, this means that the actual positions of the Se4units and I−ions has little effect on the low energy M-derived bands. Fig. 3.6(a) shows the bands obtained from Eq. 3.6for ∆i=0,t=1eV, δ=0.3eV and t⊥=0.05 eV. For comparison, we also show the same plot for δ=0, and we see that a finite distortion δopens a gap a the Fermi level. A finite value of ∆iwould further contribute to the opening of the gap and is not considered for simplicity. We also observe that the presence of small interchain hopping leads to an indirect gap, with the valence band maximum at A∗and the conduction band minimum at Z∗. To compare the tight-binding model with the computed ab-initio band structures, in Fig. 3.6(b) we also plot the bands in the physical 12 site unit cell. Doing so we observe a backfolding of the bands, so that below the Fermi level we now have four bands at Γ, originating from A∗,R∗,M∗ and Γ. The general agreement of the folded bands with the ab-initio band structures in Fig. 3.3suggests that this reduced model is indeed a very 3.2 methods &theoretical results for (MSe4)3Itmcs.79 Figure 3.9: Ab-initio low energy band structure (left) and total DOS (right) obtained using PBE approximation for (a) Ta D4h phase and (b) Ta D2d phase. Note a Van Hove singularity for holes marked with a dotted line in both cases. unintentional doping of the samples. This could have occurred because the samples did show some non-stochiometry [36]. If this is the case, one might ask whether there is any difference between doping electrons or holes, in particular regarding the ferromagnetic and superconducting instabilities. To answer these questions, we have computed the DOS for the (TaSe4)3I D4h and D2d structures, shown in Fig. 3.9. This reveals that while doping with electrons leads to a rather smooth increase of the carriers, doping with holes leads to a faster raise up to a Van Hove singularity at ED2d VH ≈−0.07 eV and ED4h VH ≈−0.09 eV for the two phases. This singularity emerges due to a change of shape of the Fermi surface from convex to concave as the energy is lowered, as shown in Fig. 3.10. 3.2 methods &theoretical results for (MSe4)3Itmcs.80 Figure 3.10: (a) Fermi surface for phase D4h at E= −0.16 eV. The plane kz=0is shown in blue. (b) kz=0Section of the Fermi surface for phase D4h at different energies with respect to the Van Hove singularity: slightly below (dashed), at the singularity (dotted-dashed) and slightly above (straight). The presence of this Van Hove singularity in the valence band suggests that a ferromagnetic instability could be triggered by a small amount of hole doping, as it is predicted to occur in GaSe [136] and similar monochalcogenides [137–139]. Van Hove singularities in general display different types of instabilities due to the enhanced density of states [140], and their generic phase diagrams include both ferromagnetic states and unconventional pairing [141,142] mediated by the repulsive Coulomb interactions [143]. While more work is needed to measure the amount and type of doping, as well as the Fermi surface shape, we believe the existence of this Van Hove singularity is a unique feature of the hole doped system which we conjecture will play a role in the low temperature instabilities. 3.2.7Concluding remarks on the theoretical results. By providing the first detailed band structure characterization of the (MSe4)3I compounds, our analysis has revealed that despite having a very complex lattice structure with 64 atoms in the unit cell, their low energy electronic structure is actually very simple. It can be understood in terms of a quasi 1D effective model of dz2orbitals in a chain with a three site unit cell, and a trimerization that gives rise to a gap [103,120]. Our detailed characterization has explained a number of puzzles in the literature, and will serve to properly interpret future experiments in this class of materials. First, we have provided a quantitative prediction for the transport gaps of the different structures, which are broadly consistent with experimental 3.3 experimental realization: (nbse4)3i.81 observations. Our results show that the gap magnitude is not only correlated with the amount of trimerization in the different structures, but also with the unit cell volume. In the future, our predictions will also be relevant to confirm the proposed distinction between type I and II samples in transport experiments, and to clarify the transport properties of the S4low temperature phase. Second, our work provides a clear resolution to the discrepancy between the gaps reported in transport vs. those in ARPES and optical conductivity. These experiments were carried out in the high temperature D4h phase, whose band structure is essentially that of an indirect gap semiconductor with valence band at A∗and conduction bands at Z∗. A very weak modulation folds both bands to Γ, but observing them in ARPES or probing an optical transition between them is extremely hard due to the very low spectral weight, proportional to the weak modulation. The large gaps quoted in ARPES and optical conductivity actually correspond to the direct gap in the unfolded bands, which is much larger than the true gap. Our predictions can be readily tested in new ARPES experiments probing the A∗ point. In addition, future low temperature measurements in the D2d phase might be more sensitive to detect the true gap in optical conductivity and ARPES, as we have shown. Finally, our work calls for more studies to understand the origin of the metallic behaviour of D2d (TaSe4)3I studied in Ref. [36]. Both ARPES and optical conductivity will be useful to quantify the existence of any extrinsic doping. In addition, if the doping is hole-like, ARPES experiments can directly map the Fermi surface and confirm the existence of a low energy Van Hove singularity, which will be relevant to understand the magnetic and superconducting instabilities. We hope our work will motivate further studies on the subject to explain this unusual coexistence. 3.3 experimental realization: (nbse4)3i. The theoretical insights provided in the previous Section have laid the foundation for the experimental collaboration presented in this Chapter. Building upon these predictions, the experimental work carried out by the group of N. Schröter focuses on verifying the theoretical claims and uncovering additional phenomena linked to the interplay of band structure and optoelectronic properties. By employing advanced ARPES techniques, their experimental work delves into the spectral weight distributions across the three-dimensional Brillouin zone, offering a comprehensive view of the electronic states. This approach not only corroborates the theoretical findings but also extends our understanding of how approximate dark states influence the optical and transport behaviors of (NbSe4)3I and its family. 3.3 experimental realization: (nbse4)3i.82 The importance of this Section lies in its contribution to bridging the gap between theoretical predictions and experimental observations. The experimental results presented herein provide crucial evidence for the impact of approximate translational symmetries and spectral weight modulations on the optoelectronic properties of (NbSe4)3I and its family. These results deepen our understanding of the material’s band structure and pave the way for future studies on related compounds. Ultimately, this work highlights the synergy between theory and experiment in addressing complex problems in condensed matter physics and underscores the significance of (NbSe4)3I as a model system for exploring fundamental concepts in lowdimensional materials. In the following pages we will focus more on understanding the experimental results and their interplay with our theoretical findings rather in the experimental details, for those we refer to the original article, Ref. [144]. As a novel insight, this experimental work identifies the modulations discussed in the previous Section as manifestations of approximately dark states, drawing on the framework introduced in the recent work by the group of K. Su Kim [145]. In the following Section, we introduce the concept of dark states and present a detailed analysis of the experimental results based on this interpretation, thereby extending our theoretical insights and broadening the scope of the developed framework. 3.3.1Introduction: What are dark states. Quantum states that do not allow for optical transitions are often characterized as dark [146–148]. Although no universal framework exists for their origin, a recent publication has proposed that such states commonly arise in crystals where two sublattices are related by a glide symmetry. This symmetry leads to double destructive interference in the optical matrix elements across alternating Brillouin zones, effectively suppressing photoemission signals and preventing these states from being observed in ARPES [145]. Notably, this classification appeared only after the work presented in the previous Sections had been completed. Nonetheless, the low-spectral-weight electronic states we identified (originating from band folding and orbital character associated with the quasi-1D chain geometry), exhibit remarkably similar phenomenology. In retrospect, these features can be naturally interpreted within the framework of approximately dark states, providing a valuable conceptual tool to analyze the experimental results presented below. While Ref. [145] emphasizes the importance of glide symmetries for the appearance of dark states, as discussed in the previous Section, approximate 3.3 experimental realization: (nbse4)3i.83 translation symmetries are often responsible for the suppression of spectral weight and, consequently, optical transitions appearing in ARPES. As explained in Sec. 3.2.3, approximate translational symmetries arise when the lowest-order terms in an effective tight-binding model allow for translational symmetries involving a small unit cell. These symmetries are then broken by higher-order terms that require extending the unit cell to its full primitive size [149]. As shown in Sec. 3.2.4, these higher-order corrections lead to the emergence of folded bands with small spectral weight alongside bright bands with strong spectral weight that follow the symmetry of the smaller unit cell. In the limiting case, the complete suppression of folded bands would only be possible when the approximate translation symmetries become exact symmetries of the crystal structure. In such a case, the system can be described by a smaller unit cell and fewer bands. In the following pages, the experimental work carried out by our collaborators from N. Schröter group confirms experimentally what we already demonstrated theoretically: that what we now call approximate dark states, which result from approximate translation symmetries, can influence the optoelectronic properties by redistributing the spectral weight and causing a direct band gap semiconductor (NbSe4)3I to behave as if it had an indirect band gap. In the previous Section, we concluded that optical transitions in this material are dominated by bright states with strong spectral weight obeying the symmetry of a smaller approximate unit cell, resulting in an enhanced optical band gap compared to the transport gap, in the following, we will test this conclusion by analyzing ARPES results. While previous ARPES experiments failed to capture this behavior properly, due to our theoretical findings, our experimental colleagues observe these approximate dark states by inspecting the A∗point in reciprocal space. By using photon-energy-dependent ARPES measurements, our collaborators found that the electronic structure of (NbSe4)3I deviates from the previously assumed quasi-1D nature and instead displays strongly dispersive bands along all three momentum directions. By obtaining the electronic band structure along all high-symmetry kpaths in a 3D Brillouin zone corresponding to the smaller approximate unit cell containing only 3Nb atoms, we can identify a spectral weight modulation resulting in approximately dark states. Specifically, the data collected by our collaborators supports a relatively small direct transport gap of 0.2eV in (NbSe4)3I that allows for the thermally activated occupation of approximately dark states, which can contribute to transport responses. However, they show that the optical transitions are dominated by the bright states detected in ARPES, which form an indirect band gap. Beyond the specific example of (NbSe4)3I, these experiments underline 3.3 experimental realization: (nbse4)3i.84 the impact of approximate symmetries on spectroscopic measurements, demonstrating how the presence of approximate dark states leads to a direct-to-indirect band gap transition that determines the optoelectronic properties of semiconductors. 3.3.2Detailed structural classification. The main structural details for (MSe4)3I(M=Ta,Nb) TMCs are present in Sec. 3.1.4and Fig. 3.2. The Nb valence states are sequentially 4+,4+,5+, with a shorter bond length for the Nb4+−Nb4+bond (3.02Å vs. 3.25Å for the Nb4+−Nb5+bond ). As shown in the previous Section, there is an approximate symmetry shown in Fig.3.11 (a) which is a unit cell with 1/4 the volume and just 3Nb atoms. The unit cells are related by c∗=c/2 and a∗=a/√2. In line with previous bibliography [123], our collaborators found the cleavage plane10 to be the (110)plane. In this geometry, as shown in Fig. 3.11 (b), we define kch ≡kzas a k-path along the chain, kip ≡kxas the direction accessed by in-plane rotation of the analyser with respect to the sample normal, and k⊥≡kyas the direction normal to the surface. In k-space, the crystallographic 12-Nb BZ is the one marked wiht a dashed blue line in Fig. 3.11 (b), but the BZ corresponding to the approximate 3-Nb unit cell is larger and rotated as shown in Figs. 3.2and 3.11, here marked in red. Each of the Γ∗,Z∗,M∗,A∗points map to Γpoints of the crystallographic unit cell. 3.3.3Approximately dark states as seen in ARPES. The ARPES dispersion along the chain direction, kch is presented in Fig. 3.12 (a). The main notable feature of said dispersion is that the spectral weight is dominated by a bright band whose maximum is located at ≈0.5eV below EFand presents a periodicity of 4π/c. Hence, the brightest features in the data seem to follow the expected periodicity for the approximate 3-Nb unit cell, not the crystallographic 12-Nb one. Moreover, this data is congruent with older literature reports, albeit that the valence band maximum of the main band was found slightly at higher binding energies [123, 125]. Along the bright band, weaker bands above EF−0.5eV are also observed in Fig. 3.12 (a), marked by a black arrow. Eventhough these are true eigenstates of the full system, they are approximately dark since their spectral weight is suppresed due to approximate translation symmetries. This sig10 The cleavage plane is the specific crystallographic plane along which a crystal sample is cleaved to create the surface for measurement. 3.3 experimental realization: (nbse4)3i.85 Nb I Se (a) (b) 3.25Å 3.25Å 3.06Å 4+ 5+ 4+ 4+ 5+ 4+ b a c (110) (001) 3.25Å 3.06Å (c) X* Z* A* Г* M* R* c d (d) R* X e Nb I Se (a) (b) 3.25Å 3.25Å 3.06Å 4+ 5+ 4+ 4+ 5+ 4+ b a c (110) (001) 3.25Å 3.06Å (c) X* Z* A* Г* M* R* c d (d) R* X e Nb I Se (a) (b) 3.25Å 3.25Å 3.06Å 4+ 5+ 4+ 4+ 5+ 4+ b a c (110) (001) 3.25Å 3.06Å (c) X* Z* A* Г* M* R* c d (d) R* X e kz≡kch kx≡kip ky≡k⊥ Figure 3.11: (a) Simplified 12-Nb unit cell (blue line) and approximate 3-Nb unit cell (red line), see Fig. 3.5. (b) Approximate 3-Nb unit cell’s BZ with high simmetry points marked with * to distinguish from initial 12-Nb unit cell’s BZ. gray shaded areas c,d, and e represent the Fermi surface maps shown later. Adapted from Ref. [144]. nals that the 12-Nb unit cell must be considered in order to identify all states within the measured spectral function, and so careful understanding of the interplay of the crystallographic and approximate unit cells is required. Fig. 3.12(b) shows the DFT11 bandstructure for D4h phase. The low energy bands, mainly characterized by Nb’s 4d2 zorbitals have a cosinelike dispersion along kch with a periodicity of 4π/c or 2π/c∗, following the periodicity of the approximate 3-Nb unit cell and not the overall 2π/c periodicity from the crystallographic lattice. There is a multiplicity in these cosine-like dispersions derived from the combination of finite hoppings transverse to the chain direction together with the expansion of the unit cell when compared with the 3-Nb description. However, only one of these branches is bright experimentally, while the three other are approximately dark states. This aligns with the 3-Nb picture, where only a single cosinelike band would be pressent dispersing upwards along kch from a minimum at Γ∗. The rest of them appear due to partial breaking of 3-Nb unit cell translational symmetry. Since ARPES spectral weight seems to follow the 3-Nb approximate cell description (see App. C.4for the details on how Fig. 3.13 (b) was constructed), in Fig. 3.13 we can finally compare our collaborators’ ARPES data with an unfolded DFT calculation obtained by following the methods presented in App. C.1that we already used in Sec. 3.2.4. The distribution of spectral 11 All DFT calculations for this experimental Section were carried out under the same guidelines as those presented in Sec. 3.2.1,i.e. with projector-augmented wave pseudopotentials, GGA approximation with PBE parametrization [74] using a 520eV kinetic energy cutoff and a 9×9×5 k−mesh. 3.3 experimental realization: (nbse4)3i.86 (a) (b) (c) approximately dark state Figure 3.12: Approximately dark states revealed in dispersion of spectral function along chain. (a) Dispersion along the chain direction. The energy cut is measured at 43eV and 90K, red arrows represent the periodicity of 4π/c and the black arrow indicates the approximately dark states. (b) DFT calculation for 12-Nb atom structure unit cell. Blue band correspond to the bright bands shown in (a). Adapted from Ref. [144]. weight agrees well with the unfolded DFT bandstructure and particularly highlights that bright band maximum in the 3D BZ is found at the A∗point, which can only be accessed experimentally by the correct combination of kperp and kch (see App. C.4), a subtlety which was not appreciated in previous ARPES reports. This would explain the incoherent picture presented in Sec. 3.1.4and addressed in Sec. 3.2.4. ARPES measurements were performed at 90 K, which should render the system’s structure to the D2d phase [16]. Due to the issue of thermallyactivated loss of iodine from the surface, as previously studied in (TaSe4)2I in Ref. [150] where sample charging is an issue at low temperatures, a full study on the temperature dependence of band gap was not possible. As a result, other probes should be better suited for studying the temperature evolution of the gap [119]. Nevertheless, our collaborators’ experiment probed bulk-like states in the measurement and hence the main conclusions are valid irrespective of the space group, since the D2h phase has also a12-Nb unit cell. 3.3 experimental realization: (nbse4)3i.87 (a) (b) Г* Z* R* A* M* Г* Z* R* A* M* strong optical transition approximately dark state Figure 3.13: Comparison between unfolded bandstructure and ARPES revealing bright and dark states. (a) Unfolded DFT calculation (following the method presented in App. C.1& Sec. 3.2.4). The physical crystal structure with 12 Nb atoms are interpreted as a periodically repeated 3Nb atoms one. (b) ARPES dispersion along the equivalent kpath. Black arrows indicate approximately dark states. Adapted from Ref. [144]. 3.3.4Discussion & conclusion on the experimental results. Our collaborators’ experimental data plotted in Fig. 3.13 (b) shows that the band structure of (NbSe4)3Iis very close to a direct band gap at the A∗and Z∗points. However, due to the spectral weight modulation of the approximately dark states clearly detected with ARPES, the optical properties are rather dominated by transitions between bands at higher energies witha a larger spectral weight, leading (NbSe4)3Ito behave as a effective indirect band gap semiconductor in optics. This interesting behaviour foreshadowed in Sec. 3.2.5, might be not that uncommon since it can happen as long as both conduction and valence bands derive from the same orbital and Wyckoff position in the unit cell, leading to the approximate translation symmetry to apply to both equally. Another feature worth noticing is the following: even if folding perturbation does not align conduction and valence band edges, it will still be the case that the optical and transport responses will generally be different, although the true gap will not be direct in that case. Moreover, if conduction and valence bands derive from different orbitals and/or Wyckoff positions, then it might be that only one of them shows a dark state modulation, or that such modulation is different for each, leading to the optical matrix elements to be finite. 3.4 discovery of a new quasi-1d tmc:(Nb4Se15 I2)I2.88 In conclusion, the experimental work carried by the group of N. Schröter effectively proves our hypothesis presented in Sec. 3.2and further develops the proposed framework, solving the puzzling literature behind (MSe4)3I (M=Nb,Ta) TMCs in terms of dark states. It does so by showcasing a simple but dramatic example of how approximately dark states may impact the optoelectronic properties of a system and lead to glaring discrepancies with other experimental probes such as transport. These discrepancies can be readily adressed by mapping out the dark bands spectral weight in an ARPES experiment, which stands now as the main probe for understanding such systems and adress similar probelms in other materials beyond (NbSe4)3I. In this last Chapter we conducted a compehensive bibliographic review, effectively proposed a theoretical framework and methods to address discrepancies in existing literature, proposing an experimental realization that would back our hypothesis that was finally proben experimentally. 3.4 discovery of a new quasi-1d tmc:(Nb4Se15 I2)I2. The preceding Sections of this Chapter have explored the intricate interplay between theory and experiment in the study of quasi-1D TMCs. We have seen how discrepancies in the literature concerning charge ordering, electronic structure, and optical responses can be reconciled through a combination of DFT and targeted experimental investigations. This methodological approach allowed us to establish a theoretical framework that was subsequently validated through experimental work from our collaborators. Building upon this perspective, we now shift our focus to the discovery and characterization of a novel quasi-1D TMC, (Nb4Se15I2)I2. While the material was initially uncovered via experimental techniques by the group of D. P. Shoemaker, its electronic and structural properties challenge conventional expectations for quasi-1D systems, making it a fertile ground for further theoretical exploration. Unlike many of its counterparts, (Nb4Se15I2)I2does not exhibit CDW order, despite its quasi-1D structure and elemental composition. Instead, it presents a moderate bandgap of Eg≈0.6eV and a unique chiral stacking, features that distinguish it from other niobium-based TMCs. This Section aims to place the discovery of (Nb4Se15I2)I2within the broader context of TMC diversity, examining how its structural and electronic properties fit within the landscape of previously studied materials. By leveraging DFT calculations, we provide insight into its band structure, orbital character, and transport properties, reinforcing the critical role of first-principles methods in modern materials science. This study not only expands the taxonomy of quasi-1D TMCs but also opens new questions regarding the relationship between dimensionality, 3.5 concluding remarks and outlook 95 Firstly, in Sec. 3.2, we have synthesized and unified conflicting past research on the (MSe4)3ITMC series with M=Nb,Ta. To achieve this, we established a theoretical framework for understanding the electronic properties of quasi-1D TMC compounds, particularly the (MSe4)nIseries with M=Nb,Ta. Through band structure calculations with different approximations we have contextualized and cataloged the different bandgaps for all the family of compounds, identifying trimerization and cell volume as the main factors affecting the bandgap magnitude. By later incorporating band unfolding techniques and symmetry considerations, we have been able to directly compare our results with previous experimental ARPES and transport data, resolving inconsistencies in reported band gaps and Fermi surface features by determining that small structural modulations lead to small spectral weight modulations and effective indirect band gap semiconductor optical behaviour. These results provide a clearer theoretical understanding of past transport and optical experiments and predicts that probing the A∗point with ARPES would unify results on these compounds. This study also calls for new studies to understand the metallic behaviour of D2d (TaSe4)3I. Secondly, in Sec. 3.3, we have successfully applied the theoretical framework presented in Sec. 3.2to understand the experimentally observed indirect bandgap optical behavior of (NbSe4)3I. A compound that DFT predicts to be a direct bandgap semiconductor. Thanks to the measurements of the group from N. Schröter, we determine that said unexpected behaviour comes from spectral weight modulations caused by small structural changes. To determine this, our collaborators probted the A∗and Z∗points with ARPES, successfully identifying what we define as approximately dark states. The framework presented in Sec. 3.2is thus further expanded by recognizing that this behavior is not uncommon when conduction and valence bands originate from the same orbital and Wyckoff position, leading to approximate translation symmetry that applies to both bands equally. These experimental findings open new avenues for studying and understanding quasi-1D TMCs, defining ARPES as one of the primary probes for investigating the observed behavior. Thirdly, in Sec. 3.4we have broadened the exploration of quasi-1D TMCs and opened new avenues for research. We have characterized and studied the discovery of a new quasi-1D TMC, (Nb4Se15I2)I2, discovered and synthesized via CVT by our colleagues and predicted by DFT to be a semiconductor with a moderate band gap of approximately 0.6eV. The transport measurements carried out by the group of D. P. Shoemaker confirmed its semiconducting behavior, with an activation energy of 0.1eV, likely due to shallow defect states. Unlike other materials in the (MSe4)nI 3.5 concluding remarks and outlook 96 family, (Nb4Se15I2)I2does not exhibit CDW order, as confirmed by lowtemperature X-ray diffraction and differential scanning calorimetry. This distinction raises new questions about the role of symmetry and electronic interactions in stabilizing CDWs in closely related systems. Moreover, our work suggests new directions for future research. The discovery of (Nb4Se15I2)I2 opens the possibility of tuning its band structure via chemical doping, external fields, or strain to explore emergent phases beyond semiconducting behavior. Furthermore, given the complexities of magnetotransport in related materials such as (TaSe4)2I, it remains an open question whether (Nb4Se15I2)I2could host nontrivial topological states or unconventional transport properties under specific conditions. Expanding the chemical phase space of quasi-1D TMCs by varying chalcogen and halogen substitutions may lead to the discovery of new materials with exotic electronic phenomena. We hope this work will inspire further experimental and theoretical studies aimed at uncovering novel quantum states in low-dimensional systems.As an intriguing example of the potential applications of quasi1D TMCs in modern technologies, Ref. [158] demonstrates how Nb3Se12I can be used for high-capacity optical processing and encryption. This highlights the rich landscape of functional applications these systems may enable. Finally, the outlook for this Chapter is clear: the systematic application and expansion of the framework presented in Sec. 3.2in order to build sound knowledge that allows the technological application of quasi-1D TMCs. For example, this could be achieved through the direct-to-indirect gap tunning and hence optical response tunning through doping or strain. 4 2D, BULK TMDS & THEIR HETEROSTRUCTURES Building on the insights gained from the study of quasi-1D TMCs, we now shift our focus to TMDs and their heterostructures, marking a simultaneous transition in both dimensionality and electronic configuration. These layered materials have attracted significant attention due to their tunable electronic [24], magnetic [25], and optical [26] properties. They provide a versatile platform for investigating quantum many-body phenomena, including Kondo [30] and Mott [31] physics. Moreover, their monolayer nature, combined with the ability to form heterostructures [32] through stacking or twisting, opens new avenues for exploring phenomena that are difficult to observe in bulk materials. Notably, interlayer coupling [27], moiré superlattices [28], and strong spin-orbit interactions [29] are intrinsic features of these systems that render electronic interactions more complex and interesting. We will begin by introducing the subject of this Chapter: TMDs and their heterostructures. We will present their characteristic phenomenology following Ref. [159] and highlight how they differ from the quasi-1D TMCs studied in the previous Chapter. This includes an overview of their crystal structure, the types of CDW they exhibit, and the mechanisms by which they can host superconductivity. Finally, we will briefly outline the key open questions that this Chapter aims to address, and explain how the subsequent Sections contribute to answering them. 4.1 introduction:get to know tmds TMDs have a long history, dating back to their discovery by Pauling in the 1920s [23]. By the 1960s, approximately 60 TMDs had been identified, with at least 40 of them exhibiting a layered structure. The first reported production of MoS2monolayer suspensions dates back to 1986 [160]. The remarkable success of graphene-related research, beginning in 2004, led to major advancements in techniques for working with layered materials, paving the way for new studies on TMDs and, more specifically, on their ultrathin films. In contrast to quasi-1D TMCs, TMDs such as MoS2, WS2, NbSe2, and TaS2are characterized by a uniform formal oxidation state of the transition metal, typically M4+, which results in a regular distribution of d-electrons throughout the lattice. This uniformity underpins their well-ordered crystal 97 4.1 introduction:get to know tmds 98 Figure 4.1: "Periodic Table" of known layered TMDs, organized based on the transition metal element involved, summarizing their existing structural phases (2H,1T or others). Adapted from Ref. [159] structures, which are typically composed of stacked layers of edge-sharing MX6(M = transition metal, X = chalcogen) units arranged in either trigonal prismatic (2H phase) or octahedral (1T phase) coordination [161] as we see in the following. The homogeneous electronic environment facilitates a broad range of ground states depending on the electron count and lattice symmetry: semiconducting behavior is observed in compounds like MoS2 and WS2, while metallicity and correlated phases such as CDWs and superconductivity emerge in materials like NbSe2and TaS2[133,159]. The absence of mixed valence prevents intrinsic charge modulation, allowing the electronic properties to be more directly tuned by external parameters such as doping, pressure, or dimensional confinement. 4.1.1Their structure TMDs appear in a plethora of structural phases, originating from different coordination environments of the transition metal atoms. The two most common structural phases are characterized by trigonal prismatic (2H) or octahedral (1T) coordination of the metal atoms (see Fig. 4.4). These structural phases can be understood in terms of differences in the stacking order of the three atomic planes chalcogen-metal-chalcogen that form the individual layers of these materials. The 2H phase corresponds to an ABA stacking, where chalcogen atoms in different atomic planes occupy the same position (A) and are aligned vertically. In contrast, the 1T phase is characterized by an ABC stacking order. The thermodynamically stable phase, either 1T or 2H, depends on the specific combination of the transition metal (from groups IV, V, VII, IX, or X) and the chalcogen (S,Se, or Te), see Fig. 4.1. Nevertheless, the other phase can often be obtained as a metastable configuration. 4.1 introduction:get to know tmds 99 The structure of TMDs can be further specified by the stacking configuration of individual layers in multilayer and bulk samples, as well as by possible distortions that reduce periodicity. In the following, we label the two-dimensional polytypes with letters Tand H, while the bulk polytypes are denoted as 1T, 2H, 4Hb, etc., where the number refers to the number of layers in the unit cell. 4.1.2Their CDWs. Group V layered dichalcogenides, in both the 2H and 1T phases of TaS2 and TaSe2, as well as in 2H-NbSe2, exhibit CDW order [162,163]. Their rich and diverse phenomenology strongly depends on the chemical composition and crystalline phase of the material. For instance, 2H-NbSe2undergoes a transition at TCDW =33 K, leading to the formation of a superlattice with approximately a 3×3periodicity. Unlike other CDW materials, this periodicity is almost independent of temperature and does not exhibit any first-order transition to commensurate phases. The atomic displacement is very small, on the order of 0.05 Å, and the CDW has a relatively moderate impact on the physical properties, causing only a slight increase in resistivity [164]. As the temperature is further lowered, the material becomes superconducting [164]. By contrast, 1T-TaSe2undergoes a sequence of incommensurate (T≈600 K), nearly commensurate (T≈350 K), and commensurate (T≈180 K) transitions, the latter exhibiting a √13 ×√13 periodicity with a so-called Star of David (SoD) CDW, see Fig. 4.5. These transitions are accompanied by significant atomic displacements of Ta atoms, reaching up to 0.24 Å, which strongly affects the electronic properties. The transition from the normal state to the incommensurate phase is second-order and is accompanied by a small jump in resistivity. In contrast, the nearly commensurate and commensurate transitions are first-order and lead to a substantial increase in resistivity [165,166]. Following these results on bulk TMDs, the study of CDW phases in monolayers became a focus of research. DFT calculations initially suggested that monolayer H-NbSe2exhibits a CDW phase with a different periodicity, 4×1, compared to the bulk structure, 3×3×1, along with a greater electronic energy gain and a higher transition temperature, increasing from Tbulk CDW =33 K to Tmono CDW =145 K [167]. However, further experimental research found no such enhancement and instead reported a 3×3periodicity for the monolayer [168]. In the case of the isostructural H-TaSe2, first-principles calculations and experimental studies suggest a 3×3CDW structure for both bulk and monolayer [169]. Several studies have investigated the CDW phase in thin films of 1T-TaS2with thicknesses down to a monolayer [170–172]. However, the evolution of both the commensurate 4.1 introduction:get to know tmds 100 and nearly commensurate phases upon reducing dimensionality remains under debate. In this manuscript, we will adopt the bulk √13 ×√13 CDW structure for 1T-TaS2. Although Fermi surface nesting has long been proposed as a mechanism for CDW formation, its applicability to TMDs and other real materials remains an open question. It is now widely recognized that nesting-driven CDW formation is a highly idealized scenario, applicable primarily to strictly one-dimensional systems or in exceptional cases [173,174]. In more complex materials, such as TMDs, the situation is considerably less straightforward: the lack of a clearly defined nesting vector, along with significant electronic structure reconstructions induced by lattice distortions, have led many researchers to argue that Fermi surface nesting plays only a minor role in the CDW formation of these systems. Instead, several studies propose that the observed distortions originate from strongly anisotropic electron–phonon interactions, whose matrix elements are significantly enhanced at the CDW wave vector [167–169,175,176]. In the following discussion, we do not engage in this debate, but rather rely on experimentally reported or theoretically predicted CDW structures. 4.1.3Their superconductivity. Although superconductivity is not the central focus of this Thesis, it emerges naturally within several of the systems under investigation and can be partially understood by the ab initio methods employed. An example of this appears in Sec. 4.4, where we use DFT to understand how the system under study metallices and subsequently hosts superconductivity. To provide the necessary background for its occasional appearance throughout the text, a brief introduction to the fundamentals of superconductivity is included in Appendix C.9. Is common for bulk TMDs that exhibit a CDW state in their phase diagram to also display superconductivity, with exceptions such as 2H-NbS2, which only shows superconductivity without an accompanying CDW phase [177]. In TMDs, superconductivity can be either intrinsic or induced through chemical doping, electrostatic doping, or applied pressure. In 2H-NbSe2, 2H-TaS2, and 2H-TaSe2, superconductivity coexists with the CDW phase at low temperatures [178–180]. In contrast, for 1T-phase TMDs such as 1TTaS2, the superconducting state emerges as the CDW order melts, as happens for the system studied in Sec. 4.4. This transition can be triggered by applying pressure or through chemical doping, such as copper intercalation. Superconductivity in TMDs persists even in the two-dimensional limit. 4.1 introduction:get to know tmds 101 Truly intrinsic superconductivity has been reported in monolayer H-NbSe2 [168], with a critical temperature of T2D C=3K, which is lower than that of the bulk material, TBulk C=7K. The in-plane critical field required to suppress the superconducting state is an order of magnitude larger than that of the bulk material. This phenomenon is understood to arise from the combined effects of reduced dimensionality and strong SOC [181]. In monolayer Hphases, the spin splitting at the Kand K′valleys acts as an effective out-of-plane Zeeman field. Due to spin-valley coupling, this spin splitting remains compatible with Cooper pairing. As a result, the superconducting state in monolayer TMDs exhibits an Ising-like spin polarization in the out-of-plane direction, making it significantly more robust against in-plane magnetic fields. 4.1.4Open questions and outline. Once introduced the key structural and electronic features of TMDs and their heterostructures, we now turn to the central questions this Chapter seeks to address. As we presented, in contrast to the quasi-1D TMCs of the previous Chapter, TMDs offer a versatile platform where dimensionality, interlayer coupling, and charge redistribution can be engineered with remarkable precision. Yet, despite significant experimental and theoretical advances, fundamental aspects of their electronic behavior, particularly regarding charge transfer, Kondo and Mott phenomena along the emergence of metallicity and superconductivity, remain insufficiently understood. One open question, for instance, concerns the electronic behaviour of T/H heterostructures. As we will see in the following, despite both 1T and 2H phases having partially filled metal d-orbital bands, their low-temperature monolayer properties are markedly different: while Hlayers exhibit weak CDWs and remain metallic and often superconducting [162], Tlayers are understood to be Mott insulators in the √13 ×√13 CDW phase [182]. Even moderate Coulomb repulsion in the orbitals from Tlayers is expected to drive the system into a Mott insulating state [182] and potentially a spinliquid phase [183]. However, when a Tlayer is placed on top of an Hlayer, charge transfer from the Hto the Tlayer partially fills the flat band (see Fig. 4.3to see the Tmonolayer flat band when in √13 ×√13 CDW phase), and it remains unclear whether the Mott state survives under these conditions. Additionally, in this partially filled regime, the localized SoD moments can hybridize with itinerant carriers in the Hlayer, giving rise to Kondo screening. Both the suppression of Mott insulating behavior and the emergence of Kondo physics thus depend sensitively on the degree of charge transfer between layers. Motivated by this interplay, the first part of this Chapter, and the frame- 4.2 background to our theoretial research 102 work developed therein, is devoted to understanding charge transfer in T/H heterostructures. In Sec. 4.3, we conduct a systematic ab initio investigation of how various factors, such as interlayer spacing and chalcogen composition, affect charge transfer in TMD heterostructures. We ask: How do these factors influence the redistribution of electronic charge, and how can this be linked to the onset or suppression of correlated electronic phenomena? Our goal is to build a comprehensive picture of how structural and chemical parameters shape the electronic environment in MX2(M=Nb,Ta; X=S,Se) bilayers, and thereby gain insight into the emergent low-energy behavior of Tand Hmonolayers in proximity. Meanwhile, in Sec. 4.4, we combine DFT calculations with experimental results from our collaborators in the group of M. Ugeda to address another open question: Why is bulk 1T-TaSSe metallic? While this question might seem a bit naïve, its all of the more interesting: previous experimental studies based on ARPES and STM claim bulk 1T−TaX2(X=S,Se) to be Mott insulating [59,60]. However, more recent investigations like Ref. [39] point out dimerization as the main mechanism for the insulating behaviour. In Sec. 4.4we identify the pathway by which an a priori insulating system, closely related to the Mott-insulating 1T-TaS2, becomes metallic and superconducting. Through this combined theoretical and experimental approach, we identify stacking disorder in the CDW phase as the key mechanism, providing a microscopic understanding of metallization in bulk 1T-phase TMDs. 4.2 background to our theoretial research Building on the insights collected from the study of TMDs’ phenomenology, we now present our theoretical work on charge transfer in T/H heterostructures. We first motivate our research by discussing the physical phenomena that these systems can host, namely, Mott and Kondo physics. To that end we conduct a thorough DFT study to determine the conditions under which each heterostructure could realize these phases. This study aims to provide a solid theoretical foundation for future experimental and theoretical investigations. Despite both 1T and 2H phases having partially filled metal d-orbital bands, their low-temperature monolayer properties are markedly different: while Hlayers exhibit weak CDWs and remain metallic and often superconducting [162], Tlayers are understood to be Mott insulators [182]. This behavior is attributed to the previously discussed √13 ×√13 SoD CDW reconstruction [184], which gaps out most of the Fermi surface, leaving a half-filled flat band at the Fermi level derived from an isolated orbital centered at the SoD. Even moderate Coulomb repulsion in these orbitals is expected to 4.2 background to our theoretial research 103 drive the system into a Mott insulating state [182] and potentially a spinliquid phase [183]. However, the realization of such states in the bulk 1T polytype, composed of stacked Tlayers, is complicated by interlayer tunneling and the multiple possible CDW stacking patterns [39]. The recent successful synthesis of Tmonolayers [40] has provided stronger evidence supporting the Mott insulator scenario. TMD heterostructures that alternate Tand Hlayers offer an interesting platform to test the Mott insulator hypothesis. In these structures, the localized SoD moments in the Tlayer couple to the metallic electrons in the H layer, naturally realizing a CDW-induced Kondo lattice. Recent studies on synthesized T/H bilayers have indeed revealed prominent zero-bias peaks in TaSe2[40,41], TaS2[42,43], and NbSe2[44,45], which have been interpreted as signatures of the Kondo effect. Beyond isolated bilayers, two naturally occurring bulk polytypes formed by alternating Tand Hlayers also exist: 4Hb[47,48], with an inversion-symmetric stacking of T/H bilayers, and 6R [49,50], with an inversion-breaking rhombohedral stacking. Interest in these compounds has recently been renewed, as they not only exhibit similar Kondo effects [51,52], but also shows unexpected signatures of unconventional superconductivity [53] (specifically in 4Hb-TaS2), including spontaneous time-reversal symmetry breaking at Tc[54], spontaneous vortex formation in the superconducting state [55], superconducting edge modes [56], and transport evidence of a two-component order parameter [57,58]. The interpretation of these experiments, as well as the validity of the Mott insulator picture, critically depends on two key properties of the T/H interface, which remain under debate: the interlayer hybridization Vand the charge transfer ∆C = (CH−CT)/2, where CH=CTin vacuum, with the factor 1/2 accounting for double-counting. The Tlayer can only remain a Mott insulator if the flat band remains nearly half-filled. However, since the work function of the Hlayer is larger than that of the Tlayer, some charge transfer from Tto H is expected [185–187]. Similarly, the Kondo effect can only survive within a limited range of charge transfer and requires a critical hybridization V > Vc. Superconductivity is also expected to behave differently depending on whether the Tlayer contributes magnetic moments. A recent study has proposed a scenario in which the Tlayer behaves as a doped Mott insulator with negligible hybridization V[46]. Importantly, the role of interlayer charge transfer in bulk 4Hbcompounds has not yet been explored, nor has the potential variability among different members of this family, such as TaS2, TaSe2, or NbSe2. 4.2 background to our theoretial research 104 4.2.1A brief note on Mott & Kondo physics. This manuscript is primarily focused on single-particle physics. Nevertheless, a single-particle approach, mainly through ab initio methods, can provide valuable insights into the collective many-body behavior of the systems under study. This is particularly true for the TMD heterostructures that we examine extensively in this Chapter. As mentioned in the previous Section, the interpretation of experiments and the determination of the many-body nature of TMDs critically depend on interlayer hybridization and charge transfer between layers. These key parameters can be extracted from a DFT study, which in turn helps illuminate the physical understanding of these systems, commonly regarded as Mott insulators or Kondo lattices. In the following Sections, we provide a brief introduction to these collective phenomena to establish the groundwork for our research. 4.2.1.1Mott transition As anticipated in the Introduction, reducing TaX2with X=S,Se TMDs to their monolayer structures, T-TaX2, where correlation effects can be more pronounced, raises fundamental questions about the nature of electronic correlations in these systems, questions that we will explore in the following pages. To establish a theoretical foundation, we begin by introducing the metal-to-insulator Mott transition in a basic manner. Condensed matter physics deals with systems composed of an extremely large number of elementary constituents, typically on the order of 1023. This results in a vast number of degrees of freedom, making their understanding and description challenging with a purely reductionist approach. New physical principles are thus required to describe such systems as a whole [188]. Emergent collective phenomena, arising from the interactions among microscopic constituents, are ubiquitous. These phenomena cannot be understood merely by considering individual constituents in isolation; rather, they emerge from the cooperative behavior of the system as a whole. Notable examples include superconductivity and magnetism, both of which are driven by spontaneous symmetry breaking1. In most cases, interacting systems must be approximated using a minimal description in terms of weakly interacting quasiparticles. Examples include Bogoliubov quasiparticles in the Bardeen-Cooper-Schrieffer (BCS) model of superconductivity, phonons in lattices with broken translational symmetry, and spin waves in magnetic systems [189]. In this discussion, we will follow Ref. [190] to introduce and get a big picture on Mott physics and how can it appear on T/H TMD heterostructures. 1It is considered spontaneous as it emerges from the intrinsic interactions of the system without requiring external fields or perturbations. 4.3 charge transfer in T/H heterostructures of monolayer tmds.111 as a characteristic resistivity upturn at low temperatures. This effect is inherently non-perturbative and is best understood within the framework of renormalization group theory, which shows that the effective coupling between the localized moment and conduction electrons grows stronger as temperature decreases. Beyond its fundamental implications in strongly correlated electron systems, Kondo physics plays a crucial role in a wide variety of materials, including heavy fermion compounds, quantum dots, and engineered heterostructures, where localized spins interact with conduction electrons. This last example would be the case on T/H heterostructures. In the context of T/H TMD heterostructures, the potential emergence of Kondo screening depends on the balance between interlayer charge transfer and hybridization, which determines whether the Tlayer retains a localized moment or becomes metallic. All necessary ingredients for the Kondo effect are present: a nearly flat band in the monolayer Tpolytypes (see Fig. 4.3), associated with the √13 ×√13 CDW state, could act as an effective spin-1/2 localized magnetic moment if it remains half-filled. Meanwhile, the metallic monolayer Hpolytype provides a conduction electron sea capable of interacting with a localized moment in the Tlayer, potentially leading to a Kondo state. Charge transfer and interlayer hybridization are thus key parameters determining the regime of T/H heterostructures. A thorough systematic ab-initio charge transfer study will be the guiding theme of the following pages. 4.3 charge transfer in T/H heterostructures of monolayer tmds. In this Section, we systematically investigate charge transfer across all these compounds, analyzing its correlation with work function mismatch, Van der Waals corrections, Hubbard’s Uparameter, and interlayer spacing. Understanding charge transfer is crucial for determining the electronic nature of these systems, as it directly influences whether the Tlayer retains its Mott-insulating character or transitions into a Kondo lattice by interacting with conduction electrons from the Hlayer. Our calculations reveal that the impact of Uand Van der Waals effects is minimal, while charge transfer is predominantly governed by the interlayer spacing and the mismatch in work functions. Moreover, we identify a general trend where Se-based compounds exhibit lower charge transfer than their S-based counterparts, and 4Hbbulk polytypes display stronger charge transfer than isolated bilayers. These findings provide a solid foundation for understanding the degree of charge transfer in T/H heterostructures and its implications for emergent many-body physics, offering a practical cookbook for predicting which systems exhibit greater or lower charge transfer under specific conditions. 4.3 charge transfer in T/H heterostructures of monolayer tmds.112 °M K ° °1.0 °0.5 0.0 0.5 1.0 E-E f(eV) NbS2 °M K ° °1.0 °0.5 0.0 0.5 1.0 E-E f(eV) NbSe2 °M K ° °1.0 °0.5 0.0 0.5 1.0 E-E f(eV) TaS2 °M K ° °1.0 °0.5 0.0 0.5 1.0 E-E f(eV) TaSe2 Figure 4.3: Electronic bandstructures from the Tmonolayers of the systems considered in the √13 ×√13 CDW. Early works [185–187] already anticipated that charge transfer from the Tto the Hlayers must be present in bulk 4HbTMDs. T/H bilayer and bulk 4HbMX2structures are shown in Fig 4.4. At high temperatures where the CDW in the Tlayers is incommensurate, the change in CDW wavevector compared to bulk 1T polytypes was used to estimate a transfer of 0.12 e− per formula unit [185] (1.56 e−per SoD) in 4Hb-TaS2, and a similar estimate leads to 1.20 e−per SoD for TaSe2. Similar theoretical estimates [187] for 4Hb-TaS2similarly ranged between 1.04 to 1.43 e−per SoD. Charge transfer of this magnitude was also reported to be consistent with changes in the optical conductivity in both 4Hband 6R polytypes [186]. More recent ARPES experiments estimated 0.92 e−per SoD in 4Hb-TaS2[201]. All these early estimates are thus consistent with a nearly empty flat band. Recent ab-initio calculations in the high temperature state without CDW [202, 203] also suggest charge transfer from Tto H, but given the strong band reconstruction due to the CDW, it is important to perform these calculations in the CDW state8. Such calculations [52,56] for a T/H bilayer of TaS2 8In the following, in DFT calculations we consider Tlayer √13 ×√13 CDW but not the 3×3 H layer CDW because of the steep computational cost it would imply. Moreover, it 4.3 charge transfer in T/H heterostructures of monolayer tmds.113 M X Hb - MX 4 2 T/H - MX 2 T H T’ H’ T H Figure 4.4: Bulk 4Hb(left) and bilayer T/H (right) MX2(M=Nb,Ta; X=S,Se) structures considered in this Chapter. still report a fully empty flat band, while a bilayer Ton monolayer Hreported 0.31 e−per SoD cell [205]. For TaSe2, a value of 0.32 e−per SoD cell was reported [41], while for NbSe2 0.17 e−was calculated [204]. A recent study has emphasized the importance of the stacking distance on charge transfer [46], revealing that ∆C ranges from 0.4to 1in TaS2as the interlayer distance goes from 7to 5.8Å. Given such variabilty, and the fact that isolated bilayers on substrates may not stack with the same interlayer distance as bulk 4Hbcompounds, it is important to study the distance dependence in detail for TaSe2and NbSe2. Differences may be expected because the flat band in the Se compounds is significantly closer to the CDW valence bands compared to the S compounds, this is shown in Fig. 4.3. For NbSe2one experiment has claimed [45] a charge transfer with an opposite sign to that of the Ta compounds. It is also important to take into account the details of the different ab-initio calculations done for the Ta compounds [206–212] , for example because the exact position of the flat band within the CDW gap is known to depend on the functional used [213], which can influence charge transfer. Including the Hubbard interaction and explicitly accounting for magnetic states with spin-split bands [213–217] can similarly affect the charge transfer. is a common practice since it is not expected to have a huge impact in charge transfer [41, 46,204]. 4.3 charge transfer in T/H heterostructures of monolayer tmds.114 Finally, note in the context of correlated systems the term charge transfer is often used to emphasize the distinction between Mott and charge transfer insulators [218]. In that case, the term refers to charge transfer between correlated d-derived bands and the dispersive p-derived bands of the same compound. In TMDs this phenomenon may also be relevant at least for some compounds like 1T-NbSe2or 1T-TaSexTe1−xwhere the flat band may overlap with the p-derived states [44,213,219]. In our work, unless specified otherwise, charge transfer will rather refer to interlayer charge transfer between the Tand the Hlayers. 4.3.1Ab-initio methods The aim of this chpater is to provide a systematic study of the interlayer charge transfer between Tand HMX2layers. The workflow used to carry out said study is as follows: 1. Single layer in-plane structure relaxation for Hand Tlayers. 2. Single layer workfunction calculation for Hand Tlayers. 3. T-H bilayer relaxation, first relaxing in the ˆzdirection followed by a subsequent in-plane relaxation. 4. Charge transfer calculation for bilayers. In this step we explore different parameters as Hubbard U, Van der Waals corrections and interlayer distance dependence. 5.4Hbcharge transfer calculation using step 2structures and experimental distances. The purpose of this workflow is to identify trends between different factors affecting ab-initio calculations and experimental measurements, aiming to understand the overall behaviour of charge transfer under different conditions. All calculations were performed using Vienna Ab initio Simulation Package (VASP) [67,68] v.6.2.1. with projector-augmented wave pseudopotentials within the Perdew Burke Ernzerhof parametrization [126]. For step 1 and 3, the relaxation was conducted by using the conjugate-gradient algorithm as implemented in VASP, keeping the cell shape and volume fixed while letting atomic positions relax. For steps 1−4, calculations were found to be well converged with a 480 eV kinetic cutoff and a gamma-centered 15 ×15 ×1 k−mesh. Meanwhile, for step 5, the self-consistent calculations were found to be well converged with a 480 eV kinetic cutoff and a gammacentered 13 ×13 ×3 k−mesh. In step 4, when Van der Waals corrections 4.3 charge transfer in T/H heterostructures of monolayer tmds.115 C B A Figure 4.5: Lattice structure of the SoD CDW distortion of the Tlayers. A, B, C inequivalent metal sites are marked in yellow, green and blue respectively. Black arrows show their displacements while a dotted line marks the CDW unit cell. were considered, DFT-D3method [220] with zero damping was used. Also in step 4, in order to study the effect of Coulomb repulsion, the DFT+U rotationally invariant approach [221] was followed by setting different effective on-site UCoulomb interactions in M’s d−orbitals with J=0. These DFT+U calculations are the only collinear spin-polarized ones. The initial magnetization was set to 1.4µBfor the central Ta/Nb atom (A atom, see Sec. 4.3.2), and zero for the remaining atoms. This choice is based on prior studies suggesting that the predominant magnetic moment is concentrated at the center of the SoD [208], and that the total magnetic moment is typically around 1µB[214,216]. The initial value of magnetization is taken slightly larger than the expected result as this is expected to improve convergence [222]. Atom-projected band structures were obtained using PyProcar [223] package for Python. 4.3.2Charge density wave distortions The √13 ×√13 SoD CDW structure is shown in Fig. 4.5. There are three types of symmetry equivalent M sites labeled as A, B and C hereafter, with multiplicities 1,6, and 6respectively. The structure is parameterized by three independent displacements uiwith i=A,B,C, shown in Fig. 4.5 (|uA|=0by symmetry). In table 3we present the displacements for each equivalent metal site in each T-MX2from the bilayers considered. In Appendix C.6we present a graphic depiction of these displacements from the Tlayer of the bilayer for the four compounds (Fig. C.5) along a brief note on how we obtained the CDW positions. 4.3 charge transfer in T/H heterostructures of monolayer tmds.116 4.3.3Work function analysis We define the work function as the absolute Fermi level with respect to vacuum. Different work functions between two structures indicate misaligned Fermi levels and can be used for a qualitative estimate of charge transfer. To compute the work functions, we perform a self-consistent calculation of the monolayer in a unit cell with ≈20 Å of vacuum in the normal direction. From this calculation we extract the local potential V(r)and define the work function as W=Vvac −EF, (4.15) where Vvac is the value of V(r)in vaccuum. We present the work functions for both Tand Hpolytypes of all TMDs considered in table 3.Hwork function is greater than Twork function for all TMDs which anticipates that the charge transfer will occur from Tto Hlayer. In Table 3we can already see some trends: M=Ta compounds have overall smaller work function than M=Nb and so happens with X=Se when compared with X=S ones. In Section 4.3.4we will see how this affect the charge transfer. The work function results in Table 3are consistent with Ref. [205], which reported work function values of W=5.35 eV for T-TaS2in the CDW state, and W=6.07 eV for H-TaS2in the 3×3CDW state. Here we calculated W=5.19 eV for T-TaS2in the CDW state, and W=5.57 eV for H-TaS2 without CDW. We opted not to include the 3×3CDW state in H-TaS2 layers in our calculations due to the excessive computational cost of a unit cell commensurate with both √13 ×√13 and 3×3CDW states. However, our analysis of the charge transfer as function of work function differences suggests this is a good approximation. 4.3.4Charge transfer Charge transfer is the central result of this work. To calculate the charge transfer between the different constituents of a heterostructure we followed the method presented in Refs.[205,224]. This method is based on computing the charge density of the full structure with respect to that of a hypothetical reference structure built from the calculated charge densities of the isolated constituents, positioned in the places they would occupy in the full structure. The choice of this method relays on two main pilars: firstly, it only needs the main output from DFT calculations (i.e. charge density) with small processing, letting our systematic approach to be performed without calculation overhead. And secondly, as showed in Refs. [205,224], it can provide a quantitative measurement for charge transfer. More explicitly, if we consider ˆzas the stacking direction, we can obtain the plane-averaged charge density from the self-consistent calculations ρall(z)and ρi(z)with 4.3 charge transfer in T/H heterostructures of monolayer tmds.117 NbSe2NbS2TaSe2TaS2 |uA|(Å) 0.00 0.00 0.00 0.00 |uB|(Å) 0.26 0.21 0.26 0.20 |uC|(Å) 0.32 0.28 0.31 0.25 d(Å) 7.37 6.89 7.51 7.10 WT(eV) 5.25 5.47 5.02 5.19 WH(eV) 5.57 6.13 5.45 5.96 ∆WF (eV) 0.32 0.66 0.43 0.77 CT (e) 0.12 0.23 0.14 0.26 Table 3: Displacements for each inequivalent metallic position uiwith i=A,B,C, interlayer distance d, work functions WTand WHfor T and Halong ∆WF and charge transfer CT for all four compounds. 4.3 charge transfer in T/H heterostructures of monolayer tmds.118 0.00 2.50 Ω(e/˚ A) 1T 1H 1T/H - TaS2 -0.01 0.00 0.01 Ωdif (e/˚ A) 0 5 10 15 20 25 z(˚ A) 0.00 0.20 q (e) 1.00 2.00 3.00 Ωall (e/˚ A) 1T 1H 1T’ 1H’ 4Hb-TaS 2 0.00 0.03 Ωdif (e/˚ A) 0 5 10 15 20 z(˚ A) 0.00 0.50 q (e) T H Figure 4.6: Electronic density (top), electronic density difference (mid) and total charge change (bottom) for T/H TaS2bilayer. The zero for electronic density difference is marked with a vertical dashed line. i=1,2,...,N, where Nis the number of component, and obtain the overall charge density difference as ρdif(z) = ρall(z) − N X i=1 ρi(z), (4.16) Integrating ρdif(z)we can obtain the total charge difference in a section (z1,z2)as qdif(z1,z2) = Zz2 z1 ρdif(z)dz. (4.17) Determining appropriate values for z1and z2is straightforward for bilayers: as depicted in Fig. 4.6,z1may be positioned anywhere in vacuum, while z2=2marks the point where ρdif(z2) = 0nearest to the interface. q(z1,z2)represents the total charge difference for the first component and the overall absolute charge transfer between components. However, choosing zican be challenging when N > 2 or for periodic systems, such as the 4Hbstructure. A practical approach is to set z1as the point where ρdif(z1) = 0nearest to the interface with the preceding component, and z2 as the point where ρdif(z2) = 0nearest to the interface with the subsequent component. For instance, in Fig. 4.7, the total charge difference for the T layer in 4Hbis obtained by integrating from z1=0Åto z2=6.40 Å (where the vertical and horizontal dashed lines first intersect in the second graph); 4.3 charge transfer in T/H heterostructures of monolayer tmds.119 1.00 2.00 3.00 ρal l (e/˚ A) T H T’ H’ 4HbTaS2 0.00 0.03 ρdif (e/˚ A) 0 5 10 15 20 z(˚ A) 0.00 0.50 q (e) Figure 4.7: Electronic density (top), electronic density difference (mid) and total charge change (bottom) for 4HbTaS2. Zeros for electronic density difference and thus candidates for z1/z2are marked with vertical dashed lines. for H, integration ranges from z′ 1=z2=6.40 Å to z′ 2=11.95 Å , and so forth. These integration limits correspond to the maximum and minimum values of q(0,z)as illustrated in the third graph. In Table 3we present the charge transfer results for the relaxed interlayer distances. These distances are 6.89Å for NbS2,7.37Å for NbSe2,7.10Å for TaS2and 7.51Å for TaSe2. Those values are plotted in Fig. 4.8, where the previously mentioned CT ∝∆WF trend is clear. From these results, we can detect some trends: ab-inito calculations predict Ta compounds to have greater charge transfer than Nb compounds, and S compounds have greater charge transfer than Se compounds. 4.3.5Van der Waals effect In this Section, we compare the calculation of the charge transfer with and without Van der Waals corrections for the T/H structures obtained in Sec. 4.3.4. To do so we incorporated Van der Waals corrections into the electronic self-consistent calculations and repeated the charge transfer calculations in Sec. 4.3.4. Our findings indicate that these corrections have a negligible influence on charge transfer, typically on the order of < 10−3. While Van der Waals corrections therefore do not affect charge transfer directly, the could do so indirectly if we performed a new relaxation of the 4.3 charge transfer in T/H heterostructures of monolayer tmds.120 0.30 0.40 0.50 0.60 0.70 0.80 ∆ Workfunction (eV) 0.10 0.12 0.15 0.18 0.20 0.23 0.25 0.28 Charge Transfer (e) NbSe2 TaSe2 NbS2 TaS2 Figure 4.8: Charge transfer as a function of workfunction difference between layers in T/H bilayer. structure in the presence of such corrections, since the interlayer distance could change upon relaxation. Since our main interest in this work is to establish relative trends in charge transfer, we have rather opted to study charge transfer as a function of interlayer distance without attempting to calculate precisely its equilibrium value, as this is a more complex problem that depends on both calculational details and experimental conditions. 4.3.6U dependence Following the same logic as in the previous Section, we now consider the effect of the Hubbard interaction U, only at the level of electronic selfconsistent calculations keeping the structure fixed 9. We considered only TaS2as an example. The main effect of the Hubbard interaction is to magnetize the flat band near the Fermi level (see Fig. 4.3), producing a spin splitting and pushing one of the spin polarizations above the Fermi level. Since a larger magnetization is produced when the lower spin-split band is closer to half-filling, increasing Ugenerally leads to a reduction in charge 9To verify our decision of discarding Uin the structural relaxation, we conducted a test relaxation using the same procedures described in 4.3for U=2.0eV in TaS2. We found that the effect of including Uon the relaxation only changes ionic positions on the order of 10−4−10−5Å. These results are consistent with the findings of L. Crippa et al. in Ref.[46], which suggest that the overall influence of Uin DFT calculations for these compounds is small. 4.4 experimental realization:metallicity through random stacking.127 charge transfer accross the family of T/H structures, which will lead to a deeper understanding of the unconventional magnetic and superconducting properties in this family of materials. 4.4 experimental realization:metallicity through random stacking. In the previous Section, we explored charge transfer and its potential role in determining the collective electronic state of T/H heterostructures. In what follows, we present DFT calculations that supported our experimental collaborators from the group of M. Ugeda in characterizing and uncovering the microscopic origin of superconductivity in bulk 1T-TaSSe (1T-TaS2−xSex with x=1). To that end, we studied how an a priori insulating system becomes metallic and subsequently hosts superconductivity. 4.4.1Introduction & Background Octahedrally coordinated 1T bulk polytypes, such as 1T-TaS2and 1T-TiSe2, exhibit similar electronic phase diagrams, like the one shown in Fig. 4.13 (b). Notably, neither of them is superconducting in its pristine form; however, superconductivity can be induced in them through doping [170,230], pressure [231], chemical substitution, or intercalation [232,233]. In both compounds, these external perturbations drive the system through a sequence of phase transitions: from a CCDW state to an ICDW state, preceding the emergence of superconductivity. The transition to an ICDW phase results in the formation of CCDW domains separated by sharp, interconnected domain walls, which were hypothesized to be the driving force behind the onset of superconductivity [230,231]. However, despite these insights, exact knowledge on how this transition to a ICDW leads to metallization and the subsequent emergence of superconductivity remains lacking for TaS2. The emergence of superconductivity upon external perturbation is particularly interesting in the case of 1T-TaS2, which exhibits an insulating ground state at low temperatures, characterized by a √13 ×√13 CDW. The origin of this insulating behavior has been largely attributed to a possible Mott state arising from strong electronic correlations within the CDW [234–236]. However, more recent interpretations suggest that the insulating gap originates from interlayer CDW coupling [39,237–239]. The application of the aforementioned perturbations typically leads to metallization of the material, which subsequently enables the emergence of superconductivity, with optimal critical temperatures TCranging from 1to 6K [170,240] as shown in the phase diagram from Fig. 4.13 (b). In the following, we identify the randomization of the CDW stacking or- 4.4 experimental realization:metallicity through random stacking.128 CCDW (Insulating) ICDW (Metallic) SC Normal metal X Temperature a b cde f Uppermost 1T layer Second 1T layer 1.7 Å 01.7 Å 0 2.4 Å 0 10 nm 2 nm 2 nm CCDW (Insulating) ICDW (Metallic) SC Normal metal X Temperature a b cde f Uppermost 1T layer Second 1T layer 1.7 Å 01.7 Å 0 2.4 Å 0 10 nm 2 nm 2 nm a b Figure 4.13: Generic phase diagram of 1T −TMD metals where xvariable in the ˆx axis refers to Se content upon Ssubstitution. der as the main driver of the metallic behavior in 1T-TaSSe, based on a combined analysis of STM/STS experiments performed by our collaborators and our own DFT calculations. This interpretation is motivated by previous studies on isostructural compounds [216,239,241], where stacking disorder has been linked to the suppression of insulating behavior. To investigate this mechanism, STM/STS measurements were performed at different spatial locations, revealing local variations in the electronic structure. By comparing these measurements with DFT simulations of distinct stacking configurations, we were able to associate each local spectrum with a specific stacking order, providing a microscopic understanding of the metallicity in terms of stacking disorder. 4.4.2dI/dV spectroscopy The system under study, 1T-TaSSe, is a related compound to 1T-TaS2. It becomes metallic10 upon selenization [240,242]: a process consistsing in an isovalent substitution of Sby Se. A gradual selenization of 1T-TaS2−xSex leads to the melting of the CCDW for x > 0.8, favoring the ICDW followed by the emergence of superconductivity in the range 0.9<x<1.6, as shown in Fig. 4.13 (b). Throughout this Thesis, we refer to the ICCDW phase observed in 1T-TaSSe as a mosaic phase. This terminology captures the spatially inhomogeneous character of the charge-ordered state, in which locally commensurate SoD clusters coexist with domain boundaries and defects, forming a pattern 10 In the presented references, said compound becomes metalic, hosts and optimizes superconductivity upon selenization. 4.4 experimental realization:metallicity through random stacking.129 reminiscent of a disordered mosaic. This ICCDW phase in 1T-TaSSe exhibits irregular arrangements of SoD clusters due to the competition between different stacking sequences and the chemical substitution of Se for S. This leads to short-range CDW correlations and a lack of global periodicity, characteristics that we deem more accurately conveyed by the term "mosaic" rather than simply "incommensurate." The use of this term emphasizes the structural fragmentation and local order inherent to the ICCDW in this material. This mosaic structure is apparent in the topographical image showed in Fig. 4.14 The first step in the experimental routine carried by our collaborators involves characterizing the large-scale electronic structure of the mosaic phase in order to understand its metallic behaviour. This is achieved through dI/dV spectroscopy and contextualized by using DFT calculations that we present in the following. A brief introduction to the fundamentals of this experimental technique, following Ref. [243] is presented in App.C.8. In short, we can link dI/dV to the local DOS (LDOS) in a region of the sample, this is showcased in Fig. 4.16 where we show the agreement between DFT calculations for the DOS and dI/dV experimental measurements. A representative figure on the topography of the mosaic phase is presented in Fig. 4.14, where the domain walls are marked with a yellow line and the mosaic phase is explicit. The results from the dI/dV spectroscopy are presented in in Fig. 4.15. The orange dI/dV spectrum shown in Fig. 4.15 (a) qualitatively captures the essential features observed in the mosaic phase of 1T-TaSSe: the occupied state region (represented by the negative bias) shows a featureless, increasing DOS and the empty state region (positive bias) presents two pronounced wide peaks C1and C2for energies Vb< 0.6V. These peaks are consistently observed across different regions of the sample and remain approximately equidistant within the range EF< Vb< 0.6V. This behavior is further illustrated by the gray dI/dV spectrum in Fig. 4.15, acquired at a different spatial location. 4.4.3Why stacking disorder? We attribute the variation of the observed empty-state resonances (differences between the orange and gray lines in Fig. 4.15) to changes in the vertical stacking of the CDW between layers, as previously reported in 1TTaS2[241]. We further hypothesize that this random CDW stacking constitutes the main mechanism behind metallization in 1T-TaSSe by a comprenhensive analysis of past literature that we detail in the following. This hypothesis is rooted in previous studies on isostructural compounds. 4.4 experimental realization:metallicity through random stacking.130 a bc 2.4 Å 0 09 µS Topography T = 2 K Conductance, Vs = 0 V (EF) -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 0 2 dI/dV (Arb. Units) Bias voltage (V) 1 -400 mV 5 mV 620 mV Conductance min max 230 mV Figure 4.14: Experimental topography (Vs=0.37 V,I=0.3 nA). Domain walls are marked with a yellow line. Adapted from Ref. [244]. For instance, in 1T-TaS2, vertical CDW stacking forms strong AA dimers (see Fig. 4.5), which render the system insulating both in periodic and nonperiodic dimer stacking sequences. By contrast, 1T-TaSe2adopts a periodic AC stacking sequence, shows no dimerization, and exhibits metallic behavior [239]. While in the previous Section the choice of chalcogen affects the overall charge transfer, we see that here it plays a crucial role in determining the stacking arrangement. This suggests that the stacking sequence in the randomized-calchogen alloy 1T-TaSSe is inherently complex and unlikely to follow a simple pattern. Our hypothesis is further supported by the work of the group from Y. Zhang [216] showing that stacking randomization in the metastable mosaic phase of 1T-TaS2leads to metallization. Moreover, stable patches of the mosaic phase have been observed in equilibrium 1T-TaS2[245], supporting that they can spontaneously happen in 1T-TaSSe. In the case of 1T-TaS2−xSex, the random distribution of Se atoms across different layers likely pins CDW domains independently (each domain separated by the yellow line in Fig. 4.14), naturally leading to stacking disorder. This mechanism is also relevant in related alloys such as 1T-TaxZr1−xSe2, where alloy-induced disorder disrupts CDW stacking [246]. In 1T-TaSSe, this effect is expected to be even stronger, since the system statistically interpolates between the dimerized 1T-TaS2and the metallic 1T-TaSe2, making all stackings nearly equally probable. Given that both ARPES and transport 4.4 experimental realization:metallicity through random stacking.131 a bc 2.4 Å 0 09 µS Topography T = 2 K Conductance, Vs = 0 V (EF) -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 0 2 dI/dV (Arb. Units) Bias voltage (V) 1 -400 mV 5 mV 620 mV Conductance min max 230 mV Figure 4.15: Electronic structure of the mosaic CDW phase. Upper panel, conductance maps taken in the same region of a CDW domain at different bias voltages. Lower panel, dI/dV spectra acquired in different locations in CCDW domains of the mosaic phase (Vac =5mV,T=4.2K). measurements show metallization for Se content x>0.8, and in line with prior theoretical and experimental work, we conclude that metallization in these alloys could primarly arise from random CDW stacking. 4.4.4Why not domain walls? Further support for this interpretation comes from ruling out domain walls themselves as the primary cause of metallicity, since some studies pointed them as the reason of supeconductivity in these compounds [230,231]. To investigate this, our collaborators performed spatially resolved dI/dV mapping of the mosaic phase within a bias range of ±1eV, as shown in the upper panel of Fig. 4.15. These measurements reveal that the conductance is predominantly localized around the central atom and its six nearest neighbors in each SoD cluster, particularly in the energy range −0.35 V< Vb< 0.55 V, as illustrated in the central conductance maps. This orbital texture is qualitatively similar to that observed in 1T-TaS2[216,247] and highlights the dominant role of Ta-derived dz2orbitals near EF. Combined with STS data showing a finite DOS at EF, this observation confirms the metallic character of the system and supports our interpretation. 4.4 experimental realization:metallicity through random stacking.132 a b c -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 0 2 dI/dV (Arb. Units) Energy (eV) 1 0 2 4 a b c 0 2 4 a b c DOS (Arb. Units) DOS (Arb. Units) AB CDW stacking B A BC A A BC AA CDW stacking AA CDW stacking AB CDW stacking Experiment C1C2 C1C2 C1 C2 ΓM K Γ −1.0 −0.5 0.0 0.5 1.0 E - EF(eV) AA A TOP 0.00 0.02 0.04 0.06 0.08 0.10 ΓM K Γ −1.0 −0.5 0.0 0.5 1.0 E - EF(eV) AB A TOP 0.00 0.02 0.04 0.06 0.08 0.10 d a b c -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 0 2 dI/dV (Arb. Units) Energy (eV) 1 0 2 4 a b c 0 2 4 a b c DOS (Arb. Units) DOS (Arb. Units) AB CDW stacking B A BC A A BC AA CDW stacking AA CDW stacking AB CDW stacking Experiment C1C2 C1C2 C1 C2 ΓM K Γ −1.0 −0.5 0.0 0.5 1.0 E - EF(eV) AA A TOP 0.00 0.02 0.04 0.06 0.08 0.10 ΓM K Γ −1.0 −0.5 0.0 0.5 1.0 E - EF(eV) AB A TOP 0.00 0.02 0.04 0.06 0.08 0.10 d d a b c -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 0 2 dI/dV (Arb. Units) Energy (eV) 1 0 2 4 a b c 0 2 4 a b c DOS (Arb. Units) DOS (Arb. Units) AB CDW stacking B A BC A A BC AA CDW stacking AA CDW stacking AB CDW stacking Experiment C1C2 C1C2 C1 C2 ΓM K Γ −1.0 −0.5 0.0 0.5 1.0 E - EF(eV) AA A TOP 0.00 0.02 0.04 0.06 0.08 0.10 ΓM K Γ −1.0 −0.5 0.0 0.5 1.0 E - EF(eV) AB A TOP 0.00 0.02 0.04 0.06 0.08 0.10 d a b c -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 0 2 dI/dV (Arb. Units) Energy (eV) 1 0 2 4 a b c 0 2 4 a b c DOS (Arb. Units) DOS (Arb. Units) AB CDW stacking B A BC A A BC AA CDW stacking AA CDW stacking AB CDW stacking Experiment C1C2 C1C2 C1 C2 ΓM K Γ −1.0 −0.5 0.0 0.5 1.0 E - EF(eV) AA A TOP 0.00 0.02 0.04 0.06 0.08 0.10 ΓM K Γ −1.0 −0.5 0.0 0.5 1.0 E - EF(eV) AB A TOP 0.00 0.02 0.04 0.06 0.08 0.10 d d Figure 4.16: (a) and (b), sketch (left) and calculated band structure of a bilayer T/T TaSSe in the CDW phase with AA and AB stacking respectively (right). (c) Corresponding DOS of the calculated structures in (a) and (b) (upper an middle panel, respectively). The lower panel shows the dI/dV spectra acquired in the mosaic phase. Adapted from Ref. [244] 4.4.5DFT study In order to interpret these experimental results and validate our hypothesis, we performed a DFT calculation of a T/T bilayer of TaSSe with different CDW stackings (namely, AA and AB11 Nonetheless, they are not included in the main text because the experimental measurements from our collaborators were performed in domains with AA and AB stacking.). The goal of these calculations is twofold: firstly, to link the variation of the emptystate resonances to the spatially changing vertical stacking of the CDW. Secondly, to confirm that the DOS is located around the SoD within the CDW domains. All DFT calculations were performed using VASP [67,68]v.6.2.1with projectoraugmented wave pseudopotentials and GGA with PBE parametrization [74]. Self-consistent calculations were found to be well converged for a kinetic cutoff of 420eV and a 15 ×15 ×1 k-mesh sampling. Fat-band orbitalprojected bandstructures were obtained using PyProcar [223] and 2D partial charge density maps were obtained by integrating the charge density re11 We did also perform calculations for the AC stacking that we present in App. C.10.1. 4.4 experimental realization:metallicity through random stacking.133 a b c -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 0 2 dI/dV (Arb. Units) Energy (eV) 1 0 2 4 a b c 0 2 4 a b c DOS (Arb. Units) DOS (Arb. Units) AB CDW stacking B A BC A A BC AA CDW stacking AA CDW stacking AB CDW stacking Experiment C1C2 C1C2 C1 C2 ΓM K Γ −1.0 −0.5 0.0 0.5 1.0 E - EF(eV) AA A TOP 0.00 0.02 0.04 0.06 0.08 0.10 ΓM K Γ −1.0 −0.5 0.0 0.5 1.0 E - EF(eV) AB A TOP 0.00 0.02 0.04 0.06 0.08 0.10 d a b c -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 0 2 dI/dV (Arb. Units) Energy (eV) 1 0 2 4 a b c 0 2 4 a b c DOS (Arb. Units) DOS (Arb. Units) AB CDW stacking B A BC A A BC AA CDW stacking AA CDW stacking AB CDW stacking Experiment C1C2 C1C2 C1 C2 ΓM K Γ −1.0 −0.5 0.0 0.5 1.0 E - EF(eV) AA A TOP 0.00 0.02 0.04 0.06 0.08 0.10 ΓM K Γ −1.0 −0.5 0.0 0.5 1.0 E - EF(eV) AB A TOP 0.00 0.02 0.04 0.06 0.08 0.10 d d a cb d Vs = - 400 mV Vs = +5 mV Vs = +232 mV Vs = +620 mV 2.3 nS0 0.5 nS0 1.3 nS0 2.2 nS0 9 Å a cb d Vs = - 400 mV Vs = +5 mV Vs = +232 mV Vs = +620 mV 2.3 nS0 0.5 nS0 1.3 nS0 2.2 nS0 9 Å a b C1 C2 Figure 4.17: (a) DFT spatial maps of the C1 (left panels) and C2 (right panels) peaks in the AA (upper panels) and AB (lower panels) stacking. (b) Spatial mapping of the orbital texture of the CDW at selected bias, around gray C1 (upper panel) and orange C1 (lower panel). sulting from self-consistent calculations in intervals of 0.05 −0.10eV around the C1and C2peaks of the LDOS in both stackings. The T/T−TaSSe bilayer structure is built from 1T −TaSe2CDW √13 ×√13 layers with 50% randomized Satoms in chalcogen positions [248]. Separation between layers is chosen to be commensurate with bulk and set to 6.29Å for inter-layer Ta −Ta distances [225]. Our results successfully reproduce the main qualitative features of the experimental DOS as well as the peak-position dependence on stacking as shown in Fig. 4.16 (c). Fig. 4.16 (a) and (b) display the A−atom projected electronic band structure for the top 1T layer in both AA and AB stackings respectively. The band strcutres for both configurations exhibit similar features: two almost flat bands, above and below the Fermi level, lying in the gap between two sets of dispersive entangled bands located at greater and lower energies. These bands can be easy to identify when comparing with the bandstructures from the monolayers presented in Fig. 4.3. The atomprojected DOS curves displayed in upper and middle panels in Fig. 4.16. These peaks correspond, from left to right, to the two flat bands and the first set of entangled conduction bands. In order to facilitate the comparison with experimental dI/dV results, the peaks associated with the upper flat band and the first group of conduction bands have been labeled as C1 and C2, respectively. The comparison between DFT, calculated DOS and 4.4 experimental realization:metallicity through random stacking.134 the experimental dI/dV measurements show substantial agreement12, in particular when examining the energy gap between the peaks C1and C2. The separation between the unoccupied flat band and the nearest conduction band in DFT is greater for AA stacking ∆EAA ≈0.30eV, than for AB stacking ∆EAB ≈0.25 eV. This trend helps us identify the orange and gray dI/dV curves in Fig.4.15 to AA and AB stacking respectively. Upon this identification, this trend is corroborated by dI/dV measurements, where ∆EAA < 0.30eV and ∆EAB > 0.30 eV. Moreover, the 2D partial charge maps in Fig. 4.17 (a) furthermore ressemble the conductance maps at Fig. 4.17 (b) further confirming the lozalized nature from C1and C2peaks and showcasing the agreement between DFT and experimental results. Due to the agreement between DFT and experimental data, we can succesfully associate the variation of the observed empty-state resonances (differences between the orange and gray lines in Fig. 4.15) to changes in the vertical stacking of the CDW between layers. 4.4.6Concluding remarks to this Section In summary, through a combined experimental and theoretical investigation, we have demonstrated that the metallic behavior of 1T-TaSSe arises primarily from stacking disorder in the CDW phase. STM/STS spectroscopy revealed local variations in the electronic structure, which we were able to associate with distinct stacking configurations through first-principles DFT calculations. These results rule out domain walls as the main origin of metallicity and instead highlight random CDW stacking (favored by the structural complexity of the alloy and the distribution of chalcogen atoms) as the key mechanism enabling interlayer tunneling and the emergence of a finite DOS at the Fermi level arising from bulk effects. The consistency between the observed spectral features and the calculated band structures and partial charge densities reinforces this interpretation and provides a microscopic understanding of how a nominally insulating layered system becomes metallic (and, under favorable conditions, superconducting) through stacking-driven metallization. The details of superconductivity characterization along some further experimentla details beyond the scope of this Thesis are presented in App.C.10. 12 Even when these calculations reflect the main qualitative features of the surface DOS for a given stacking of the outermost layers, we could not find finite DOS at the Fermi level. This feature is only recovered in a bulk calculation with a distorted stacking arrangement. 4.5 concluding remarks 135 4.5 concluding remarks In conclusion, this Chapter has provided a detailed exploration of the electronic and structural properties of TMDs and their heterostructures, shedding light on their relevance to correlated phenomena such as Mott and Kondo physics, CDW transitions and superconductivity. Although this work is primarily based on ab initio single-particle approaches, we have established a theoretical framework that is instrumental in understanding the complex interactions present in these systems. The systematic study of charge transfer in T/H bilayers of MX2and the 4Hb−TaX2structure (M=Nb,Ta;X=S,Se) stands as the main theoretical contribution of this Chapter. Our results indicate that neither Hubbard U nor van der Waals corrections significantly affect charge transfer, and we rule out interlayer distance as the primary controlling parameter. Instead, we have identified clear trends across different compounds: generally, Sebased systems exhibit higher charge transfer than their S-based counterparts, while Ta-based compounds tend to show greater charge transfer compared to their Nb-based analogs. Moreover, we predict that charge transfer is more pronounced in bulk 4Hbheterostructures than in bilayers. These results provide valuable insights into the broader phenomenology of charge transfer in T/H systems and its implications for electronic correlations. From an experimental perspective, our findings help contextualize previous observations of Kondo behavior in T/H bilayers, particularly in TaSe2[40,41], TaS2[42,43], and NbSe2[44,45], as well as differences in the orbital character of higher-energy states [40]. Rather than indicating inconsistencies between experiments, these variations may reflect intrinsic differences in charge transfer and proximity to a potential Mott insulating state, emphasizing the need for precise experimental verification of charge transfer across individual compounds. Additionally, our study suggests that 4Hb-TaS2exhibits a larger charge transfer, with clearer evidence of a nearly unoccupied flat band [51,52,56], which plays a minor role in superconductivity. However, our calculations indicate that in 4Hb-TaSe2, the flat band may be more populated, potentially having a stronger influence on the superconducting state [227–229]. This finding opens new directions for future research on the role of flat-band physics in T/H heterostructures and its implications for unconventional superconductivity. Furthermore, the combined effort between our experimental collaborators and the DFT calculations, informed by the insights gained in Sec. 4.3, has allowed us to identify the origin of metallicity in the a priori insulating 1T-TaSSe. Our results indicate that the metallic behavior arises from the local randomization of CDW stacking, while the presence of domain walls within individual layers appears to have little impact on the emergence 4.5 concluding remarks 136 of the metallic state. These findings enabled our collaborators to further characterize and investigate the superconducting properties of the system, showcasing the strength of the interplay between theoretical modeling and experimental observation. Taken together, these results highlight the critical role of charge transfer, interlayer interactions, and structural disorder in determining the electronic properties of TMDs and their heterostructures. Future work should aim to systematically map out charge transfer across the entire family of T/H compounds, explore doping and strain engineering as potential tuning parameters, and further investigate the connection between stacking disorder and correlated effects. Expanding the chemical and structural phase space of these materials may lead to the discovery of novel emergent quantum phases, paving the way for future studies at the intersection of strong correlations, topology, and low-dimensional superconductivity.