David Ceperley

Keywords: Quantum Monte Carlo, electronic structure, Ferrite, CEIMC method, extreme

Firstly, I want to express my gratitude to my advisor Dr. Carlo Pierleoni for his advice, support, and patience during the

Synthèse

Le problème de la métallisation de l'hydrogène, posé il y a près de 80 ans, a été désigné comme la troisième question ouverte en physique du XXIe siècle. En effet, en raison de sa légèreté et de sa réactivité, les informations expérimentales sur l'hydrogène à haute pression sont limitées et extrêmement difficiles à obtenir. Il est donc essentiel de mettre au point des méthodes précises pour guider les expériences.

Au début de la thèse, nous présentons la théorie générale des méthodes électroniques de l'état fondamental utilisées dans cette thèse, qui sont principalement la théorie fonctionnelle de la densité (DFT) et la méthode de Monte Carlo quantique (QMC). Une attention particulière est portée à la méthode quantique de Monte Carlo.

Ensuite, dans le Chapitre 2, nous nous concentrons sur l'étude de la structure électronique, y compris les phénomènes d'état excité, en utilisant les techniques de QMC. En particulier, nous développons une nouvelle méthode de calcul pour le gap accompagnée d'un traitement précis de l'erreur induit par la taille finie de la cellule de simulation. Nous établissons un lien formel entre l'erreur de la taille finie et la constante diélectrique du matériau. Avant d'étudier l'hydrogène, la nouvelle méthode est testée sur le silicium cristallin et le carbone de diamant, pour lesquels des informations expérimentales sur le gap sont disponibles. Nos résultats montrent que le biais dû à la supercellule de taille finie peut être corrigé, de sorte que des valeurs précises dans la limite thermodynamique peuvent être obtenues pour les petites supercellules sans avoir besoin d'une extrapolation numérique.

Comme l'hydrogène est un matériau très léger, les effets quantiques nucléaires sont importants. Une description précise des effets nucléaires peut être réalisée par la méthode de Monte Carlo à ions et électrons couplés (CEIMC), une méthode de simulation des premiers principes basée sur le QMC. Dans le Chapitre 4 nous utilisons les résultats de la méthode CEIMC pour discuter des effets quantiques et thermiques des noyaux sur les propriétés électroniques. Nous introduisons une méthode formelle de traitement du gap électronique et de la structure des bandes à température finie dans l'approximation adiabatique et discutons des approximations qui doivent être faites. Nous proposons également une nouvelle méthode pour calculer des propriétés optiques à basse température, qui constituera une amélioration par rapport à l'approximation semi-classique couramment utilisée.

Enfin, nous appliquons l'ensemble du développement méthodologique de cette thèse pour étudier la métallisation de l'hydrogène solide et liquide dans les Chapitres 5 et 6. Nous constatons que pour l'hydrogène moléculaire dans sa structure cristalline parfaite, le gap QMC est en accord avec les calculs précédents de GW. Le traitement des effets quantiques nucléaires entraîne une forte réduction du gap ( 2 eV). Selon la structure, le gap indirect fondamental se ferme entre 380 et 530 GPa pour les cristaux idéaux et 330-380 GPa pour les cristaux quantiques, ce qui dépend moins de la symétrie cristalline. Au-delà de cette pression, le système entre dans une phase de mauvais métal où la densité des états au niveau de Fermi augmente avec la pression jusqu'à 450-500 GPa lorsque le gap direct se ferme. Notre travail confirme partiellement l'interprétation des récentes expériences sur l'hydrogène à haute pression.

Pour l'hydrogène liquide, la principale conclusion est que la fermeture du gap est continue et coïncide avec la transition de dissociation moléculaire. Nous avons été en mesure d'étalonner les fonctions de la théorie fonctionnelle de la densité (DFT) en nous basant sur la densité QMC des états. En utilisant les valeurs propres des calculs Kohn-Sham corrigé par QMC pour calculer les propriétés optiques dans le cadre de la théorie de Kubo-Greenwood , nous avons constaté que l'absorption optique théorique calculée précédemment s'est déplacée vers des énergies plus faibles.

Nous explorons également la possibilité d'utiliser une représentation multidéterminante des états excités pour modéliser les excitations neutres et calculer la conductivité via la formule de Kubo. Nous avons appliqué cette méthodologie à l'hydrogène cristallin idéal et limité au niveau de Monte Carlo variationnel de la théorie, les résultats peuvent être trouvés dans le Chapitre 3. Le développement théorique présenté dans cette thèse n'est pas limité à l'hydrogène et peut être appliqué à différents matériaux, ce qui donne une perspective potentielle pour des travaux futurs.

Quantum Monte Carlo methods for electronic structure calculations: application to hydrogen at extreme conditions

Vitaly GORELOV

This dissertation is submitted for the degree of Doctor of Philosophy

Abstract

The hydrogen metallization problem, posed almost 80 years ago [START_REF] Wigner | On the Possibility of a Metallic Modification of Hydrogen[END_REF], was named as the third open question in physics of the XXI century [START_REF] Vitaly | Nobel Lecture: On superconductivity and superfluidity (what I have and have not managed to do) as well as on the "physical minimum" at the beginning of the XXI century[END_REF]. Indeed, due to its lightness and reactivity, experimental information on high pressure hydrogen is limited and extremely difficult to obtain. Therefore, the development of accurate methods to guide experiments is essential.

In this thesis, we focus on studying the electronic structure, including excited state phenomena, using quantum Monte Carlo (QMC) techniques. In particular, we develop a new method of computing energy gaps accompanied by an accurate treatment of the finite simulation cell error. We formally relate finite size error to the dielectric constant of the material. Before studying hydrogen, the new method is tested on crystalline silicon and carbon diamond, for which experimental information on the gap are available. Although finite-size corrected gap values for carbon and silicon are larger than the experimental ones, our results demonstrate that the bias due to the finite size supercell can be corrected for, so precise values in the thermodynamic limit can be obtained for small supercells without need for numerical extrapolation.

As hydrogen is a very light material, the nuclear quantum effects are important. An accurate capturing of nuclear effects can be done within the Coupled Electron Ion Monte Carlo (CEIMC) method, a QMC-based first-principles simulation method. We use the results of CEIMC to discuss the thermal renormalization of electronic properties. We introduce a formal way of treating the electronic gap and band structure at finite temperature within the adiabatic approximation and discuss the approximations that have to be made. We propose as well a novel way of renormalizing the optical properties at low temperature, which will be an improvement upon the commonly used semiclassical approximation.

Finally, we apply all the methodological development of this thesis to study the metallization of solid and liquid hydrogen. We find that for ideal crystalline molecular hydrogen the QMC gap is in agreement with previous GW calculations [START_REF] Mcminis | Molecular to Atomic Phase Transition in Hydrogen under High Pressure[END_REF]. Treating nuclear zero point effects cause a large reduction in the gap (∼2 eV). Determining the crystalline structure of solid hydrogen is still an open problem. Depending on the structure, the fundamental indirect gap closes between 380 and 530 GPa for ideal crystals and 330-380 GPa for quantum crystals, which depends less on the crystalline symmetry. Beyond this pressure the system enters into a bad metal phase where the density of states at the Fermi level increases with pressure up to ∼450-500 GPa when the direct gap closes. Our work partially supports the interpretation of recent experiments in high pressure hydrogen. However, the scenario where solid hydrogen metallization is accompanied by the structural change, for example a molecular dissociation, can not be disproved.

We also explore the possibility to use a multideterminant representation of excited states to model neutral excitations and compute the conductivity via the Kubo formula [START_REF] Kubo | Statistical-mechanical theory of irreversible processes. i. general theory and simple applications to magnetic and conduction problems[END_REF]. We applied this methodology to ideal crystalline hydrogen and limited to the variational Monte Carlo level of the theory.

For liquid hydrogen the main finding is that the gap closure is continuous and coincides with the molecular dissociation transition. We were able to benchmark density functional theory (DFT) functionals based on QMC density of states. When using the QMC renormalized Kohn-Sham eigenvalues to compute optical properties within the Kubo-Greenwood theory [START_REF] Kubo | Statistical-mechanical theory of irreversible processes. i. general theory and simple applications to magnetic and conduction problems[END_REF][START_REF] Greenwood | The Boltzmann Equation in the Theory of Electrical Conduction in Metals[END_REF], we found that previously calculated theoretical optical absorption [START_REF] Rillo | Optical properties of high-pressure fluid hydrogen across molecular dissociation[END_REF] have a shift towards lower energies.

Introduction

It has been now almost a hundred years since according to the own words of Paul Dirac, "the general theory of quantum mechanics is now almost complete..." [START_REF] Dirac | Quantum mechanics of many-electron systems[END_REF]. Certainly, thanks to quantum mechanics we know the laws of how electrons and nuclei, building blocks of matter, interact with each other. The main problem is that the equations that need to be solved to explain the main properties of matter are too complicated to be solvable and approximations have to be made. However, since that time many important physical phenomena have been discovered that have pushed fundamental science and technological progress towards up to these days.

Year after year in the past few decades we see a dramatic improvement in modern computers. With this, our ability to simulate more and more complex physical systems with better and better accuracy is increasing. Particularly, in quantum mechanics, where the use of analytical methods is limited to only a few simplest cases and, the numerical methods must be used in order to study realistic systems. Nowadays, modern ab-initio methods can be used to accurately model systems comprising a few thousand atoms [START_REF] Steinhauser | Computational Methods on Electronic/Atomistic Scale[END_REF]. Despite all the impressive progress in this field, we are still very far from the ultimate goal, being able to accurately reproduce experimental information based only on the knowledge of atom types.

Nowadays, due to technological progress and algorithmic advances, it is possible to combine electronic and nuclear problems and to accurately obtain the ground-state properties of the full electron-nuclear system from first-principles [START_REF] Pierleoni | The Coupled Electron-Ion Monte Carlo Method[END_REF][START_REF] Marx | Ab Initio Molecular Dynamics[END_REF]. The troubles, however, arise when considering the excited states of the joint system. The treatment of electronic structure properties in the presence of nuclei at finite temperature is among the problems addressed in this work.

Hydrogen under extreme conditions

It is for a single hydrogen atom that the equations of quantum mechanics can be solved exactly. Thus, from the first sight, it seems like bulk hydrogen can be a simple enough material and can be used in testing new theoretical methods. However, being partially true, this assumption also leads to discover that the physics behind bulk hydrogen is much (green squares [START_REF] Loubeyre | Synchrotron infrared spectroscopic evidence of the probable transition to metal hydrogen[END_REF] and red squares [START_REF] Ranga | Observation of the Wigner-Huntington transition to metallic hydrogen[END_REF]) and predicted liquid-liquid transition line (blue dots [START_REF] Pierleoni | Liquid-liquid phase transition in hydrogen by coupled electron-ion Monte Carlo simulations[END_REF]).

more complex than of hydrogen molecule [START_REF] Mcmahon | The properties of hydrogen and helium under extreme conditions[END_REF][START_REF] Goncharov | Phase diagram of hydrogen at extreme pressures and temperatures; Updated through 2019 (Review article)[END_REF]. It is not by chance that hydrogen is involved in the first three questions that, according to the Nobel Laureate Vitaly Ginzburg [START_REF] Vitaly | Nobel Lecture: On superconductivity and superfluidity (what I have and have not managed to do) as well as on the "physical minimum" at the beginning of the XXI century[END_REF], modern physics is facing in the verge of XXI century. In the first problem, controlled nuclear fusion, hydrogen and its isotopes are the main elements, solving this problem would allow to deal with the energy problem for future societies over the next millennia. Second, it was predicted that atomic hydrogen would be an above room temperature superconductor [START_REF] Ashcroft | Metallic hydrogen: A high-temperature superconductor?[END_REF]. Indeed, the highest claimed temperature of superconductivity known today is for hydrogen-lanthanum combination LaH 10 (T c = 260 K at 180 GPa) [START_REF] Somayazulu | Evidence for superconductivity above 260 k in lanthanum superhydride at megabar pressures[END_REF]. The third question is directly formulated as achieving the metallic state of hydrogen, the problem that puzzles condensed matter science since the first prediction by Wigner and Huntington in 1935 [START_REF] Wigner | On the Possibility of a Metallic Modification of Hydrogen[END_REF].

To reach the near metallization pressures hydrogen has to overcome a long path of transitions. From a simple molecular gas phase, the complexity of this material builds up when pressure increases. By arriving at hundreds of GPa, the number of predicted solid-state phases had almost reached six (see fig. 1), however there is uncertainty for the crystal symmetry at each phase. Firmly established experimentally, the first three solid phases and transitions, phase I-III, can be characterized by starting in phase I with an initial rotational symmetry of spherically disordered molecules arranged in a hexagonal close packed (hcp) structure [START_REF] Mcmahon | The properties of hydrogen and helium under extreme conditions[END_REF]. The rotational symmetry then breaks when compressing hydrogen further to phase II, the molecules are now ordered or at least partially ordered, but their exact arrangement and their shape are unknown [START_REF] Goncharenko | Neutron and X-ray diffraction study of the broken symmetry phase transition in solid deuterium[END_REF]. The transition to phase II, also known as the broken symmetry phase (BSP), depends substantially on isotope, which implies an important role of nuclear quantum effects [START_REF] Isaac | The solid molecular hydrogens in the condensed phase: Fundamentals and static properties[END_REF]. The transition to phase III is characterized by the changes in infrared (IR) and Raman spectra [START_REF] Hemley | Phase transition in solid molecular hydrogen at ultrahigh pressures[END_REF][START_REF] Lorenzana | Evidence for a structural phase transition in solid hydrogen at megabar pressures[END_REF]. Also called hydrogen-A, phase III spans from ∼155 GPa at 100 K, over more than 200 GPa at low temperature [START_REF] Hemley | Phase transition in solid molecular hydrogen at ultrahigh pressures[END_REF]. It is from this phase that hydrogen is predicted to transform into a metallic phase, and in this thesis, the main focus will be on phase III. By the spectroscopic measurements, recent experimental investigations predict metallization pressure being between 425 GPa [START_REF] Loubeyre | Synchrotron infrared spectroscopic evidence of the probable transition to metal hydrogen[END_REF] and 495 GPa [START_REF] Ranga | Observation of the Wigner-Huntington transition to metallic hydrogen[END_REF] The main difficulty with achieving high pressures for solid hydrogen, especially above room temperature, is associated with the diffusive and reactive nature of the material in a dense state. Experimentally solid hydrogen is obtained by static compression in a diamond-anvil-cell (DAC), and very often, hydrogen can penetrate inside the diamond, breaking it and interrupting the experiment [START_REF] Goncharov | Probing hydrogen-rich molecular systems at high pressures and temperatures[END_REF].

Going up in the phase diagram by increasing the temperature in the solid phase, the crystal will melt forming a melting line, which was studied extensively over the past 30 years [START_REF] Silvera | Pathways to metallic hydrogen[END_REF]. Interestingly, that having a negative slope (see fig. 1 [START_REF] Deemyad | Melting line of hydrogen at high pressures[END_REF]) the line can be extrapolated to zero temperature, meaning that hydrogen might be a liquid at the ground state. The liquid phase makes up a major fraction of hydrogen in the Universe. Studying this phase can answer many questions in planetary science, regarding the composition of the planets. The predicted first order phase transition [START_REF] Ld Landau | On the relation between the liquid and the gaseous states of metals[END_REF] should occur in hydrogen at temperatures below some critical temperature and should be accompanied by the dissociation of the molecules and metallization. In this regime, besides static compression in DAC with pulse laser heating, hydrogen can be compressed by a dynamic shock wave, which results in higher pressure, however, the pay-off is the uncertainty in temperature [START_REF] Peter M Celliers | Insulator-metal transition in dense fluid deuterium[END_REF][START_REF] Knudson | Direct observation of an abrupt insulator-to-metal transition in dense liquid deuterium[END_REF][START_REF] Nellis | Minimum metallic conductivity of fluid hydrogen at 140 gpa (1.4 mbar)[END_REF]. There is still an open question whether dissociation and metallization occur at the same time. We will try to shed more light on it in chapter [START_REF] Rillo | Optical properties of high-pressure fluid hydrogen across molecular dissociation[END_REF].

The fact that hydrogen can behave as a metal can be easily demonstrated by considering the general Hamiltonian (see eq. 1.2 for example). Its potential energy scales as r -1 s while the kinetic energy goes as r -2 s , where r s is the Wigner-Seitz radius [START_REF] Ashcroft | Solid state physics[END_REF]. Therefore, as pressure increases (r s → 0), the kinetic part will dominate and electrons will favor the free particle regime (simple metal). A simple theoretical description of this transition can be given in terms of band theory. Change of density causes molecular bands to shift which can lead to band-gap closure. The transition to the metallic state does not necessarily involve the dissociation from molecular to atomic hydrogen. We discuss the gap closure of hydrogen in Chapters 5 and 6.

Structure of the thesis

In this thesis, a methodological development to study electronic gaps and excitations within quantum Monte Carlo (QMC) is discussed. The necessary treatment of electronic properties at finite temperature is presented within QMC and density functional theory (DFT) in combination with path integrals used for nuclei. The methodology is then used to discuss the gap closure and the optical properties of liquid and solid hydrogen around the metallization transition.

Chapter 1 presents the general theory of the electronic ground state methods used in this thesis, which are mostly DFT and QMC. Particular attention is drawn to quantum Monte Carlo. In chapter 2 the main methodology for computing energy gaps and excited states within QMC is presented. There we introduce a novel method to compute gaps and an important treatment of finite size effects for electron addition and removal energies, the main problem of QMC simulations of extended systems. The methodology is tested on silicon, carbon, and ideal crystalline hydrogen. A small overview of optical properties calculations within Kubo-Greenwood [START_REF] Kubo | Statistical-mechanical theory of irreversible processes. i. general theory and simple applications to magnetic and conduction problems[END_REF][START_REF] Greenwood | The Boltzmann Equation in the Theory of Electrical Conduction in Metals[END_REF] formalism with application to ideal solid hydrogen is given in chapter 3. In this chapter, within variational Monte Carlo (VMC), the Kubo formula [START_REF] Kubo | Statistical-mechanical theory of irreversible processes. i. general theory and simple applications to magnetic and conduction problems[END_REF] for computing electronic conductivity of solid hydrogen was used. Next, in chapter 4 finite temperature renormalization of electronic properties is discussed. All these new theoretical developments are applied to compute the fundamental gap closure and the optical properties of liquid (chapter 6) and solid (chapter 5) hydrogen at extreme conditions.

Note that the theoretical development presented in this thesis is not limited to hydrogen and can be applied to different materials, which gives a potential perspective for future work. At the end of the thesis, the conclusions, with a summary of the main topics addressed in this work are provided.

Chapter 1 Electronic Ground State Methods

In this chapter, we will lay the foundations of the quantum mechanical simulation techniques used to obtain the electronic ground state properties of systems of interest, which will serve as basis methods throughout the thesis. Firstly, the tremendous problem of electrons in the presence of nuclei as an external field is introduced. Then we discuss the independent particle Hartree-Fock (HF) approximation and the density functional theory (DFT). The major part of the chapter is devoted to the central method of this work: quantum Monte Carlo (QMC). In particular, we introduce variational Monte Carlo (VMC) and reptation quantum Monte Carlo (RQMC) and discuss its advantages and disadvantages. The most important part of any QMC simulation -the wave function, is discussed here in details including its different terms and approximations. Finally, we introduce a technique of treating the finite simulation cell size effects -an important problem of all QMC simulations.

Besides a standard textbook of Martin [START_REF] Martin | Electronic Structure: Basic Theory and Practical Methods[END_REF] which covers the single electron theory, a more recent book by Martin, Reining and Ceperley [START_REF] Martin | Interacting Electrons[END_REF] gives a whole overview on all methods used in treating systems with interacting electrons, including QMC. Moreover, there is a fundamental review on QMC methods by Foulkes et al. [START_REF] Foulkes | Quantum Monte Carlo simulations of solids[END_REF], which was followed by a more recent review by Kolorenč and Mitáš [START_REF] Kolorenč | Applications of quantum Monte Carlo methods in condensed systems[END_REF]. Detailed information on RQMC method can be found in lecture notes of Pierleoni and Ceperley [START_REF] Pierleoni | The Coupled Electron-Ion Monte Carlo Method[END_REF].

Many body problem

Knowing that the matter is a collection of electrons and nuclei interacting via the Coulomb potential, the equation for the stationary states δE/δΨ = 0, with E representing the eigenenergies and Ψ are the eigenfunctions or wave functions,

N i=1 - 1 2 ∇ 2 i + v ext (r i ) + 1 2 N i =j 1 |r i -r j |
Ψ(r 1 , .., r N ) = EΨ(r 1 , .., r N ), (1.1) is the time-independent Schrödinger equation in atomic units ( = 1, electron mass and charge e = m e = 1) for N electrons moving in a static potential v ext created by the presence of N p nuclei, which describes the behavior of electrons in atoms, molecules, and condensed matter. Spin and relativistic effects have been neglected, as they will not be discussed in this thesis. In eq. 1.1, it was assumed that the nuclei are static and placed at their equilibrium positions, the external potential v ext , therefore, depends on particular nuclear geometry {R} = R 1 , .., R Np . What happens if the nuclear effects are included will be discussed in detail in chapter 4.

The eigenvalues of the Schrödinger equation 1.1 are the energies of the electronic system. However, even when the nuclear motion is neglected, solving this equation remains a complex problem. In fact, if the Hamiltonian

Ĥ = Te + Vep + Vee (1.2)
is split into three parts, it is clear that the complexity comes from the last term, which is the Coulomb interaction between electrons. The presence of the Coulomb interaction is the cause of high-dimensionality, indeed, the many-body wave function Ψ(r 1 , .., r N ) is a function of 3N variables. In real systems N is on the order of 10 23 , the Avogadro's number. However, even restricting the number of electrons to a few makes it impossible to solve the Schrödinger equation directly and approximations have to be made. By neglecting the last term in eq. 1.2, one can decouple the Schrödinger equation 1.1 into many one-particle problems of independent electrons. However, despite the success of the independent electron picture, there are many systems, some of them are considered in this work, where this picture is not valid and more correlation has to be incorporated, introducing higher-order approximations, for example as GW [START_REF] Hedin | New method for calculating the one-particle green's function with application to the electron-gas problem[END_REF], Bethe-Salpeter equation [START_REF] Salpeter | A relativistic equation for bound-state problems[END_REF] or quantum Monte Carlo [START_REF] Mcmillan | Ground state of liquid He 4[END_REF]. We will nevertheless start by discussing single electron approximations to compute the electronic ground state.

Hartree-Fock

Starting from the wave function factorized into the single particle states, the first approximation to include electronic correlation is simply based on Pauli exclusion principle [START_REF] Pauli | Über den Zusammenhang des Abschlusses der Elektronengruppen im Atom mit der Komplexstruktur der Spektren[END_REF]. Working directly with electrons one has to take into account its fermionic nature, e.g., the many-body wave function is restricted to be antisymmetric under the particle exchange. It was firstly proposed by Slater in 1929 [START_REF] Slater | The theory of complex spectra[END_REF] that for single particle orbitals this constrain is automatically satisfied when the product wave function is written as a determinant:

Ψ HF = 1 (N !) 1/2 det(ψ σ i (r j )), (1.3)
where ψ σ i (r j ) denotes a normalized single-particle orbital with the spin, σ quantised along z. The energy,

E HF = dr 1 ..dr N Ψ * HF ĤΨ HF , (1.4)
can be calculated exactly and obeys the variational principle, e.g., for any wave function, Ψ, the energy, E, obtained with eq. 1.1, is always an upper bound for a true ground state energy E 0 . The total HF energy for Hamiltonian in eq. 1.

1 is

E HF = -dr i,σ ψ σ * i (r) ∇ 2 2 ψ σ i (r) + drv ext (r)n(r) + E H + E x , (1.5)
where

n(r) = i,σ ψ σ * i (r)ψ σ i (r) (1.6)
the first and the second terms are the kinetic energy and the energy due to the external potential v ext for independent particles,

E H = 1 2 drdr n(r)n(r ) |r -r | (1.7)
is the Hartree contribution, which is the interaction energy if electrons were classical particles,

E x = - 1 2 σ occ i,j drdr ψ σ * j (r )ψ σ i (r ) 1 |r -r | ψ σ j (r)ψ σ * i (r) (1.8)
is the Fock term and is the result of Pauli exclusion principle. Further minimisation of the HF total energy with respect to the orthonormal single-body orbitals ψ σ i leads to the Hartree-Fock equations for single-body orbitals [START_REF] Fock | Näherungsmethode zur lösung des quantenmechanischen mehrkörperproblems[END_REF] -

1 2 ∇ 2 + v ext (r) + v H (r) ψ σ i (r) + dr Σ x,σ (r, r )ψ σ i (r ) = σ i ψ σ i (r), (1.9)
where v H (r) is the Hartree potential that acts locally on a wave function at each point r.

It arises from the charge of all the electrons, including each electron acting on itself,

v H (r) = dr n(r ) |r -r | . (1.10)
The Σ xσ (r, r ) is the non-local Fock operator, contains only like spins

Σ x,σ (r, r ) = - occ j ψ σ * j (r ) 1 |r -r | ψ σ j (r).
(1.11) Equation 1.9 can be solved exactly only in some cases: for spherically symmetric atoms and the homogeneous electron gas. Usually, these equations are written in a finite basis and are solved for the coefficients of the expansion.

The HF ground state wave function is the determinant built from the N lowest-energy single-particle states. Koopmans' theorem [START_REF] Koopmans | Über die zuordnung von wellenfunktionen und eigenwerten zu den einzelnen elektronen eines atoms[END_REF] shows that the eigenvalues of the HF equations correspond to the total energy differences, namely, to the energies to add or subtract electrons,

± N ±1 = E N ±1 -E N , (1.12)
that would result from increasing the size of the matrix by adding an empty orbital or decreasing the size by removing an orbital if all other orbitals are frozen. Excited states can be represented by choosing other combinations of single-particle spin orbitals ψ σ i to build the determinant. However, note that in the absence of relaxation the addition, removal, and excitation energies are often overestimated. Though, performing the separate calculations of N and N ± 1-particle states would include some aspect of correlation, for an infinite system this will be negligible and would merely reproduce the eigenvalues. In general, assuming a single-determinant form for the wave function neglects correlation between electrons, which makes the Hartree-Fock theory insufficient to make accurate quantitative predictions.

Density functional theory (DFT)

In the attempts to solve the many-body problem one can try to reduce the dimensionality of the problem. For example, by working with the electronic density one can go from the 3N dimensional space of the wave function to a function of only 3 coordinates. Indeed, it seems very appealing to be able to solve the quantum mechanical problem just relying on the knowledge of the electronic density. To see how it can be possible, let us first note that the ground state total energy and wave function of an interacting many-electron system can be considered as a functional of the external potential v ext , i.e., it depends on the entire function v ext (r). This potential can be due to the nuclei and other sources (see eq. 1.1). By the use of Legendre transformation, it was shown by Hohenberg and Kohn in 1964 [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF] that the energy can be also written as a functional of the density, which is the main idea of density functional theory (DFT).

E HK [n] = Ψ| T + Vee |Ψ + drv ext (r)n(r) = F HK [n] + drv ext (r)n(r), (1.13)
where F HK [n] is the universal functional of the density. The ground state is determined by minimizing the energy. This functional is unknown and very difficult to approximate in general. However due to the ingenious work by Kohn and Sham in 1965 [START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF] practical use of DFT became possible. In their work, they were able to reformulate the original many body problem and represent it with a set of independent particle problems having an effective potential. Which have led to the set of equations representing an auxiliary system of independent-particles having the same density as the original, interacting one. The properties of this system can be derived from the single-particle equation with an effective potential v ef f (r),

- 1 2 ∇ 2 + v ef f (r) ψ i (r) = i ψ i (r), (1.14)
where all the interactions are placed into an effective potential v ef f (r) which every particle feels independently. The ground state density,

n(r) = occ i ψ * i (r)ψ i (r), (1.15)
then corresponds to the interacting one by construction, ψ i now denote the Kohn-Sham eigenfunctions.

Making the connection to eq. 1.13 the energy functional for the Kohn-Sham system can be rewritten as

E KS [n] = N i=1 - 1 2 |∇ψ i (r)| 2 + drv ext (r)n(r) + E H [n] + E xc [n], (1.16)
where

E xc [n] = Ψ| T + Vee |Ψ -N i=1 -1 2 |∇ψ i (r)| 2 -E H [n]
includes all the contributions of exchange and correlation to the ground state energy. It accounts also for many body corrections to the kinetic energy. The effective potential v ef f (r) is defined by the condition for the minimum energy:

δE KS [n] δn(r) = -v ef f ([n], r) + v ext (r) + δE H [n] δn(r) + δE xc [n] δn(r) = 0 (1.17) v ef f ([n], r) = v ext (r) + v H ([n], r) + v xc ([n], r). (1.18)
The set of eqs. 1.14, 1.15 and 1.17 (the famous Kohn-Sham equations) must be solved self-consistently, as the effective potential, v ef f ([n], r), is itself a functional of the density. Since E xc [n] is a small fraction of the total energy, it is much easier to create approximate expressions for it, than approximate the full energy functional F HK [n]. A first reasonable approximation for E xc [n], called local density approximation (LDA) was originally proposed by Kohn and Sham, is taken from results for the homogeneous electron gas (HEG) and is a local functional of the density:

E LDA xc [n] = dr n(r) HEG xc ([n(r)], r). (1.19)
The exchange energy density HEG xc 3) is the energy per electron of the HEG at point r that depends only upon the density n(r) at r. The correlation energy density has been calculated numerically for a set of densities by Ceperley and Alder (1980) using quantum Monte Carlo [START_REF] Ceperley | Ground state of the electron gas by a stochastic method[END_REF] and parametrized by [START_REF] Perdew | Self-interaction correction to density-functional approximations for many-electron systems[END_REF] [START_REF] Perdew | Self-interaction correction to density-functional approximations for many-electron systems[END_REF]. It is expected that it will be best for solids close to a homogeneous gas (like a nearly-free-electron metal) and worst for very inhomogeneous cases like atoms where the density must go continuously to zero outside the atom. Due to the success of LDA there were developed functionals which include dependence on gradients and higher derivatives of the density. However, the straightforward expansion has not solved the problem and lead to worse results as the gradients in real materials are so large that the expansion breaks down. Thus, the term generalized-gradient approximation (GGA) denotes a variety of parametrization proposed for functions that modify the behavior at large gradients in such a way as to preserve desired properties. For spin unpolarised system the generalization of eq. 1.19 will take the form of:

([n(r)], r) = -3 4 3 π n (1/
E GGA xc [n] = dr n(r) xc ([n], |∇n|, .., r). (1.20)
The parametrization proposed by Perdew, Burke, and Enzerhof (PBE) [START_REF] Perdew | Generalized gradient approximation made simple[END_REF] is used in this thesis most frequently, for example, to generate the starting Slater determinant for QMC calculations.

Generalized Kohn-Sham theory (GKS)

One can obtain a large improvement in accuracy by considering a non-local effective KS potential instead of a local one. This approach is regarded as a generalisation of the original Kohn-Sham idea. Recently, it was shown by Perdew et al. [START_REF] Perdew | Understanding band gaps of solids in generalized Kohn-Sham theory[END_REF] that in generalized KS theory (GKS), the band gap of an extended system equals to the fundamental gap. One of the ways to include non-locality is to mix the orbital-dependent Hartree-Fock and an explicit density functional, introducing a mixing parameter a. Functionals constructed in this way are called "hybrids". One of the "hybrid" functionals used throughout this manuscript is named after Heyd-Scuseria-Ernzerhof (HSE) and based on a screened Coulomb potential for the exchange interaction [START_REF] Heyd | Hybrid functionals based on a screened Coulomb potential[END_REF]. This form of functional was shown to accurately describe energy band gaps and lattice parameters of many semiconductors [START_REF] Heyd | Energy band gaps and lattice parameters evaluated with the Heyd-Scuseria-Ernzerhof screened hybrid functional[END_REF].

Bloch theorem and pseudopotentials

When dealing with periodic structures, it is possible to use relatively small cells and employ periodic boundary conditions. In fact, for crystal structures, the simulation cell can be as small as the primitive cell. Indeed, in the case of single particle orbitals, the translational symmetries of the system result in Bloch's theorem [START_REF] Ashcroft | Solid state physics[END_REF]. The KS orbitals can be then written as,

ψ ki (r) = e ik • r u ki (r), (1.21)
u ki (r + L α ) = u ki (r),
with u ki (r) having periodicity of the crystal cell with lattice vectors (L 1 , L 2 , L 3 ). By introducing further the reciprocal lattice vectors, G m ,

G m = m 1 B 1 + m 2 B 2 + m 3 B 3 , m = (m 1 , m 2 , m 3 ) ∈ N 3 , (1.22)
B 1 = 2π L 2 × L 3 L 1 • (L 2 × L 3 ) , B 2 = 2π L 3 × L 1 L 2 • (L 3 × L 1 ) , (1.23)
B 3 = 2π L 1 × L 2 L 3 • (L 1 × L 2 )
, the periodic part can then be expanded as

u ki (r) = m C kim e iGm • r (1.24)
where the summation over integers is possible thanks to the periodicity of u ki (r). Moreover, for any vector G n ,

ψ (k+Gn)i (r) = e ik • r m C (k+Gn)im e i(Gm+Gn) • r (1.25) = e ik • r m C kim e iGm • (r+Ln) = ψ ki (r),
meaning that the vector k can be confined to the primitive cell of the reciprocal lattice, conventionally called the first Brillouin zone. This plane-wave expansion is used in codes written to deal with periodic systems including those used throughout this work.

The main drawback of plane waves is the fact that they require highly dense grids to be able to represent the core atomic states. Due to orthogonality requirements, these functions are very oscillatory and localized in the core region. In the majority of cases, these states are chemically inert, so it should be possible to remove them from the calculation without major changes to the properties of the system. Pseudopotentials are renormalized electron-nuclei potentials for the valence states of an atom that include both the Coulomb attraction of the nuclei and the screening effects resulting from the presence of core electrons. By employing pseudopotentials, not only do we remove core states from the calculation, but we also obtain valence states which are smooth in the core region, which greatly reduces the computational demands of the method. For the additional details see ref. [START_REF] Martin | Electronic Structure: Basic Theory and Practical Methods[END_REF].

In the DFT calculations presented in this work, pseudopotentials are used for hydrogen, even though it does not possess a core. That is necessary because the Coulomb potential, 1/r, requires a large number of plane waves G m to achieve the convergence, which greatly increases the computational demands of the calculations. The pseudopotentials are built to reproduce the scattering properties of the atom, therefore, valence states should not be significantly affected outside the core region.

Quantum Monte Carlo (QMC)

All the methods discussed above can be classified as deterministic, which boils down to numerically solving an approximate equation to determine the properties of a quantum system. As we have seen they usually have a trade-off problem of correctly describing correlations while keeping numerical simplicity of equations. Alternatively one can try to solve the Schrödinger equation 1.1 directly by designing an appropriate stochastic method and putting correlations directly into the wave function. Such methods, named quantum Monte Carlo (QMC), provide results more accurate than DFT or HF with just an order of magnitude increase in computational cost (for both methods computational effort grows with the third power of system size, but the prefactor can differ by an order of magnitude). Being an exact method for bosons, in many cases for fermionic systems QMC provides the energies and other properties very close to the exact results.

In this section, the discussion will be limited to zero temperature or ground state QMC methods. Giving a proper description of the ground state methods is necessary, as it serves as the starting point for the development of excited states and finite temperature QMC. For a more complete description of the ground state QMC see refs. [START_REF] Martin | Interacting Electrons[END_REF][START_REF] Foulkes | Quantum Monte Carlo simulations of solids[END_REF][START_REF] Morales | First Principles Methods: A Perspective from Quantum Monte Carlo[END_REF]

Variational Monte Carlo (VMC)

The variational Monte Carlo method (VMC) was first introduced by McMillan in 1965 [START_REF] Mcmillan | Ground state of liquid He 4[END_REF], who observed that the calculation of a quantum system using a correlated wave function can be seen as the evaluation of a classical many-body system of atoms interacting with a pair potential and can be calculated with Monte Carlo techniques. The VMC method is based on the variational principle, which tells us that by choosing an arbitrary trial wave function Ψ T (R) for the electronic configuration R = {r 1 , .., r N }, the expectation value of Ĥ evaluated with the trial wave function Ψ T is always greater or equal, than the exact ground state energy

E V = Ψ * T (R) ĤΨ T (R)dR Ψ * T (R)Ψ T (R)dR E 0 . (1.26)
One can rewrite equation (1.26) in the form

E V = |Ψ T (R)| 2 [Ψ -1 T (R) ĤΨ T (R)]dR |Ψ T (R)| 2 dR E 0 , (1.27)
from which it is easy to separate the integral into the probability distribution P(R) and the observable E L (R)

P(R) = |Ψ T (R)| 2 |Ψ T (R)| 2 dR , E L (R) = Ψ -1 T (R) ĤΨ T (R).
(1.28)

Now the variational energy E V will take the form of an average, that can be evaluated by the Monte Carlo technique By applying the Metropolis algorithm one can sample a set of M points {R m : m = 1, M } from the probability distribution P(R m ) to evaluate the local energy E L (R m ) at each of these points.

E V = P(R)E L (R)dR. ( 1
E, Hartree E L (a 1 ) < E L (a 1 ) > 2 E L (a 2 ) < E L (a 2 ) > 2 E 0
E V ≈ 1 M M m=1 E L (R m ). (1.30)
The variance of the local energy has an important property,

σ 2 (R) = dRP(R)(E L (R) -E V ) 2 = dRP(R)E 2 L (R) -E 2 V 0, (1.31)
if the wave function is an exact eigenfunction of the Hamiltoninan, the local energy E L (R) does not depend on R and is equal to the corresponding eigenvalue and the variance vanishes. Therefore the variance vanishes as the trial wave function approaches the exact one. This is known as zero variance principle in Quantum Monte Carlo. These properties indicate the quality of a trial wave function, since the lower the variational energy and the variance, the better the trial function is. Figure 1.1 illustrates the two VMC trajectories computed for different approximations of the wave function. It shows that as the averaged variational energy lowers, so as the fluctuations, indicating that the system aproaches a true ground state. The advantage of the variance over the energy is that for the variance one also knows the lower bound and can estimate how far one is from the true ground state.

Going into the details of VMC, the first step is to design the Metropolis algorithm [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF] (see Appendix A) to estimate the variational energy. The probability to accept or reject the trial move of electron positions from R to R , according to Metropolis, is proportional to the ratio of the probability distributions P(R )/P(R), more specifically,

A(R → R ) = min 1, T (R → R ) T (R → R) Ψ(R ) Ψ(R) 2 , (1.32)
where T (R → R ) is the transition probability of proposing the trial configuration R given the actual configuration R. In case of uniform displacement of one electron within a volume Ω it becomes T (r i → r i ) = T (r i → r i ) = Ω -1 . The algorithm for performing Variational Monte Carlo simulation is the following:

• Pick the trial wave function. The most common form is Slater-Jastrow wave function, which includes pair correlations. The choice of the trial wave function will be discussed in detail in section 1.4.3.

• Initialise electron positions R (0) ≡ {r

(0) 1 , r (0) 2 , ..., r (0)
N }

• Iterate the loop over electrons i, M times:

-Propose the move r i → r i of electron i in a uniform way within the volume Ω -Determine the acceptance probability:

A = min 1, |Ψ(R )/Ψ(R (n) )| 2 (1.33)
-If A > u with u ∈ (0, 1) a uniform random number, accept the move and update the coordinates (R (n+1) = R ). Otherwise, reject the move and keep the old coordinates (R (n+1) = R (n) ).

-Compute and average the local energy E L (R (n+1) ) and other properties.

• Adjust the trial wave function parameters to minimise the average local energy.

The variational method is very powerful and intuitively understandable. The only assumption is the form of the trial wave function and no further uncontrolled approximations. But this fact can also bring problems to the method, for example, the construction and optimization of trial functions for many-body systems can be time-consuming and can bring an element of human bias (see Appendix B for more information). Moreover, the result of VMC calculation is strongly influenced by the form of the trial wave function.

Reptation Quantum Monte Carlo (RQMC)

Despite all the advantages of VMC, it is still desirable to have a method that will be less dependent on the input trial wave function. That can be achieved based on the idea that starting from a trial wave function and applying a suitable projection operator one can project onto the ground state.

Suppose the exact eigenfunctions and eigenvalues of the Hamiltonian Ĥ in 1.2 are φ i and E i . It is true that any trial state can be decomposed in the eigenstate basis,

|Ψ T = i c i |φ i , (1.34)
where c i is the overlap of the trial state with the i th eigenstate. Let us consider the application of the operator e -β Ĥ onto this state,

|Ψ(β) = e -β Ĥ |Ψ(0) = i c i e -βE i |φ i ∝ β→∞ c 0 e -βE 0 |φ 0 , (1.35)
with the initial state |Ψ(0) = |Ψ T . All excited states will be suppressed exponentially fast with increasing β, the imaginary projection time. The rate of convergence to the ground state depends on the energy gap between the ground state and the first excited state, non-orthogonal to the trial function. The total energy as a function of the imaginary time is defined as follows,

E(β) = Ψ T |e -β 2 Ĥ Ĥe -β 2 Ĥ |Ψ T Ψ T |e -β Ĥ |Ψ T . (1.36)
The reason for splitting the projector operator and putting it on two sides is to get a pure (unmixed) estimator. In other words, to make the same wave function appear on both sides. That is the main difference of the method described below and used throughout this thesis, called reptation quantum Monte Carlo (RQMC) [START_REF] Baroni | Reptation quantum monte carlo: A method for unbiased ground-state averages and imaginary-time correlations[END_REF], from the conventional diffusion Monte Carlo (DMC) [START_REF] Ceperley | Ground State of the Electron Gas by a Stochastic Method[END_REF] that was developed first and is more widespread.

Similarly to the thermal partition function, one can define the generating function of the moments of Ĥ as

Z(β) = Ψ T |e -β Ĥ |Ψ T . (1.37)
The total energy at β can be simply expressed as the derivative of the logarithm of Z(β)

E(β) = - ∂ ∂β ln(Z(β)). (1.38)
Now, we will prove that the energy as a function of β decreases and converges to the ground state E 0 at β → ∞. To do so it is convenient to define the variance of energy σ 2 E (β), which is by definition positive,

σ 2 E (β) = Ĥ2 -Ĥ 2 = ∂ 2 ∂ 2 β ln(Z(β)) = - ∂ ∂β E(β) 0, (1.39)
the derivative of the energy with respect to β should then be negative, meaning that it monotonically decreases as β increases and therefore the following relations must hold,

lim β→∞ E(β) → E 0 , (1.40) lim β→∞ σ 2 E (β) → 0, (1.41)
the last relation is also known as zero variance principle Eq. (1.31). Now it is useful to define the operator e -β Ĥ as the density matrix.

ρ(R, R , β) = R|e -β Ĥ |R , (1.42)
where R again represents a set of electronic coordinates. For exact eigenstates φ i (R) the density matrix writes,

ρ(R, R , β) = i φ * i (R)e -βE i φ i (R ). (1.43)
The partition function Z can be then expressed as,

Z(β) = dRdR Ψ T |R ρ(R, R , β) R |Ψ T , (1.44)
which allows us to write the average of an observable Â, not commuting with Hamiltonian Ĥ, as follows,

 (β) = 1 Z(β) d R 1 ...d R 4 Φ T |R 1 ρ(R 1 , R 2 ; β 2 ) R 2 | Â|R 3 ρ(R 3 , R 4 ; β 2 ) R 4 |Φ T .
(1.45) To get the best estimate to the ground state average it is necessary to put the observable at β/2, as this is where the left and right trial functions have been projected equally. Putting observable at any other β will provide a mixed estimator since the observable sandwiches with different wave functions.

To compute an average over the ground state one needs to know the density matrix at large β. For that it is necessary to factorize β at small imaginary time steps τ = β/M with M being the number of steps. If the time step τ is short enough the system approaches its classical limit and one can do approximations. Consider factorizing the density matrix,

ρ(R, R , β) = R|e (-τ Ĥ) M |R = dR 1 ...dR M-1 M k=1 R k-1 |e -τ Ĥ |R k , (1.46)
with the boundary conditions: R 0 = R and R M = R . Now, one has to evaluate the short time propagator. Applying the Trotter [START_REF] Trotter | On the product of semi-groups of operators[END_REF] split-up one gets,

R k-1 | exp[-τ ( T + V )]|R k ≈ e -τ V (R k-1 )/2 R k-1 |e -τ T |R k e -τ V (R k )/2 , (1.47)
where it has been used the fact that the potential energy operator is diagonal in the position representation and its matrix elements are trivial,

R| exp(-τ V )|R = e -τ V (R) δ(R -R). (1.48)
For the kinetic propagator the solution is the Green's function of the free particle diffusion equation, which in configurational space is,

∂ t Φ(R, t) = 1 2 N i=1 ∇ 2 i Φ(R, t), (1.49)
which in the integral form becomes,

Φ(R, t + τ ) = G(R ← R , τ )Φ(R , t)dR , (1.50)
where

G(R ← R , τ ) = R| exp[-τ T ]|R . (1.51)
The solution is well known and has a form of a gaussian with variance τ ,

R k-1 |e -τ T |R k = (2πτ ) -3N/2 exp - (R k -R k-1 ) 2 2τ . (1.52)
Putting everything back together to obtain density matrix ρ(R, R , t) one gets,

ρ(R, R , β) = dR 1 ...dR M-1      M k=1 e - (R k -R k-1 ) 2 2τ (2πτ ) 3N/2      e -τ V (R 0 ) 2 + M -1 k=1 V (R k )+ V (R M ) 2
.

(1.53) In the continuous limit (M → ∞, τ → 0, β = M τ = const.) it becomes the Feynman-Kac formula [START_REF] Kac | On Distributions of Certain Wiener Functionals[END_REF]:

ρ(R, R , β) = exp - β 0 dτ V (R(τ )) RW , (1.54)
where ... RW indicate a path average over gaussian random walks R(τ ) starting at R(0) = R and ending at R(β) = R at time β.

The last formula, in principle, describes the general scheme of how to perform the ground state quantum Monte-Carlo. However, it can be foreseen that the potential energy V (R) will be largely varying in the electronic configurational space R, lowering the efficiency of the method. There exist the solution to that problem, which is analogous to the classical Monte-Carlo importance sampling.

Importance sampling

To substitute the strongly fluctuating potential V (R) from the Feynman-Kac formula in eq. 1.54 by the smooth value of the local energy E L (R), it is convenient to rewrite the Hamiltonian Ĥ as follows,

Ĥ = Ĥ + E L (R), (1.55)
where

Ĥ = 1 2 i (-∇ 2 i + ∇ 2 i Ψ T Ψ T ), (1.56)
E L (R) = V (R) - 1 2 i ∇ 2 i Ψ T Ψ T . (1.57)
Now the short time propagator 1.47 becomes,

R k-1 | exp[-τ ( Ĥ)]|R k ≈ e -τ E L (R k-1 )/2 R k-1 |e -τ Ĥ|R k e -τ E L (R k )/2 . (1.58)
The solution for the short time kinetic propagator can be obtain using the same strategy as for diffusion propagator Eq. (1.52). The difference will be the correction

F k = F (R k ) = ∇ ln Ψ T (R k )
, which can be seen as the drift force,

R k-1 |e -τ Ĥ|R k = Ψ T (R k-1 ) Ψ T (R k ) 1 2πτ -3N 2 exp - (R k -R k-1 -τ F k-1 ) 2 2τ . (1.59)
The problem now is that the density matrix is not symmetric under the exchange of left and right legs R and R . The solution can be achieved by symmetrizing the short time propagator. Then it is convenient to define symmetrized kinetic link action

L s k (R k-1 , R k ; τ ), R k-1 |e -τ Ĥ|R k = e -L s k (R k-1 ,R k ;τ ) (2πτ ) 3N 2 , (1.60)
where

L s k (R k-1 , R k ; τ ) = (R k -R k-1 ) 2 2τ + τ 4 (F 2 k + F 2 k-1 ) + (R k -R k-1 )(F k -F k-1 ) 2 . (1.61)
Putting eqs. 1.58, 1.60 and 1.46 together, the density matrix, ρ(R, R , β), becomes,

ρ(R, R , β) = dR 1 ...dR M-1 M k=1 e -L s k (R k-1 ,R k ;τ ) (2πτ ) 3N/2 e -τ E L (R 0 ) 2 + M -1 k=1 E L (R k )+ E L (R M ) 2
.

(1.62) The final result in the continuum limit is called the generalised Feynman-Kac formula [START_REF] Pierleoni | The Coupled Electron-Ion Monte Carlo Method[END_REF][START_REF] Glimm | Quantum Physics, A Functional Integral Point of View[END_REF][START_REF] Ceperley | Quantum monte carlo for molecules: Greenâs function and nodal release[END_REF]:

ρ(R, R , β) = exp - β 0 dτ E L (R(τ )) DRW , (1.63)
where ... DRW indicates a path averaged over drifted random walk starting at R(0) = R and ending at R(β) = R .

It is important to define the estimator for the total energy. However, to compute the total energy one can put the local energy operator either on the left or on the right, meaning,

E(β) = 1 Z(β) dRdR E L (R)Ψ T (R)ρ(R, R , β)Ψ T (R ) = 1 Z(β) dRdR Ψ T (R)ρ(R, R , β)Ψ T (R )E L (R ) = 1 2 E L (R) + E L (R ) . (1.64)
For the variance calculation it make sense to split the local energy operator in two and insert it on both sides of the density matrix, leading to

σ 2 (β) = dRdR E L (R)Ψ T (R)ρ(R, R , β)Ψ T (R )E L (R ) -E 2 (β). (1.65)
Now we will sketch the algorithm of the importance sampling reptation QMC [START_REF] Baroni | Reptation quantum monte carlo: A method for unbiased ground-state averages and imaginary-time correlations[END_REF]. Let's define a discretised path X = {R 0 , R 1 , ..., R M }, consisting of M electronic configurations, then the probability of this path π(X) will be,

π(X) = Ψ(R 0 )ρ(R 0 , R 1 , τ )...ρ(R M-1 , R M , τ )Ψ(R M ) Z(t) .
(1.66)

The transition probability to grow the path either left or right can be defined as, According to Metropolis algorithm, the acceptance rate will then become,

T +1 = e -(R -R M -τ F M ) 2 /2τ , T -1 = e -(R -R 0 -τ F 0 ) 2 /2τ . (1.67)
A d (X → X ) = min 1, π(X )T -d (X → X) π(X)T +d (X → X ) (1.68)
The procedure of moving the path around is schematically illustrated on figure 1.2 and goes as follows:

• Randomly choosing one of the two path's ends.

• Grow the chosen side by adding to it one slice.

• Accept or reject, according to acceptance probability Eq. (1.68).

• If the move is accepted, discard one slice on the opposite side.

In the way suggested above, in order to update all slices, it takes about M 2 moves, of course, one can move k slices at a time, but it will still take about (M/k) 2 steps. Moreover, adding slices randomly either to the left or to the right side will not be efficient, as in this way the middle slice almost never moves. D. M. Ceperley and C. Pierleoni came up with the so-called bounce algorithm [START_REF] Pierleoni | Computational methods in coupled electronion Monte Carlo simulations[END_REF], where one keeps moving in the same direction until rejection occurs.

As it was stated above, in eq. 1.45, the easy access to the pure distribution makes RQMC ideal for calculations of unbiased observables and correlation functions, doing so in a more efficient manner than simple forward-walking in other methods like Diffusion Monte Carlo, which is based on a similar idea of projecting out the ground state. Estimation of observables over the pure distribution works whenever one can write a meaningful estimator in terms of position space coordinates. Diagonal position space observables, like the average potential energy and pair-correlation function, can be measured directly from the sampled pure distribution. Observables that are not diagonal in position space, like off-diagonal density matrix elements and the momentum distribution, can be estimated from the pure distribution with some suitable additions to the basic algorithm [START_REF] Pierleoni | Electron localization properties in high pressure hydrogen at the liquid-liquid phase transition by Coupled Electron-Ion Monte Carlo[END_REF].

Even though, the projection result guarantees a closer answer to the ground state that the VMC, there is still some dependency on having a good trial wave function. Firstly, for the reasons of convergence efficiency. Secondly, and more importantly, when dealing with Fermi statistics the density matrix between configurations R and R is no longer positive everywhere since the permutation of the final configuration PR , multiplied by the sign of the permutation has to be taken care of. In this form the fermionic density matrix can not serve as the probability distribution. Information on the regions where the wave function changes sign should be then taken from the trial wave function, which introduces a crucial approximation to the projection methods for fermions.

Fermion sign problem

So far, there was only defined the general Boltzmann density matrix 1.46, which serves us as a probability distribution in the MC run. No information concerning the bosonic or fermionic nature of the particles was provided up to this point. In order to generalize the density matrix for bosons and fermions one has to take into account permutations and antisymmetry:

ρ B/F (R, R , t) = 1 N ! P (±1) P ρ(R, PR , t) (1. 69
)
where P is one of the N ! permutations of particle labels. One can, therefore, think of a path in the configurational space of distinguishable particles as an object carrying not only a weight (given by the exponential of minus the integral of the local energy along the path eq. 1.63), but also a sign fixed by its boundary conditions in time. Hence, as far as this sign is positive (ρ(R, R , t) ≥ 0 for any R and t at given R, i. e. as far as one is dealing with bosons, the sum over permutations can be easily carried out. However, for fermions the density matrix ρ F can be negative, which gives rise to the fermion sign problem [START_REF] Anderson | Quantum chemistry by random walk[END_REF]. The origin of the problem is the fact that the sign has to be left in the estimator for the averages in order to have a density matrix as a sampling probability. Therefore, for such alternating series, the signal to noise ratio will decay exponentially and it will be very inefficient to use the direct method to sample fermions. The fermion sing problem is the most challenging problem in QMC, which many scientists are trying to solve. However, nowadays pretty robust approximations are available. The most widely used is the so-called restricted path or fixed node method [START_REF] Ceperley | Fermion nodes[END_REF]. Within the fixed node method, one needs to consider the nodal surfaces of the fermion density matrix. For any given configuration R, these are defined by the implicit equation ρ F (R, R , t) = 0, as the set of locations R at which the density matrix at time t vanishes.

The nodal surfaces of the initial configuration R divide the configurational space of R in regions of positive ρ F and regions of negative ρ F . One can define Υ(R, t) as the set of points that can be reached from R in time t without having crossed the nodal surfaces at previous times. Formally it can be written as the restricted path identity by restricting the functional integral to paths inside Υ(R, t) in such a way that Ψ(R)Ψ(R ) ≥ 0,

ρ F (R, R , t) = 1 N ! P (-1) P Y(0)=R,Y(t)=PR DYe -S[Y] Υ(R,t) , (1.70)
where S[Y] represents the action over generic path Y. Therefore, using the restricted path identity it can be shown [START_REF] Pierleoni | The Coupled Electron-Ion Monte Carlo Method[END_REF] that the generating function is a positive function at any time t and can be computed considering only positive paths which do not cross the nodal surfaces. A further very important property of the fixed node method is the existence of a variational theorem: the FN-RQMC energy is an upper bound of the true ground state energy E T (∞) ≥ E 0 , and the equality holds if the trial nodes coincide with the nodes of the exact ground state [START_REF] Ceperley | The stochastic solution of the many-body Schroedinger equation for fermions[END_REF]. Therefore even for fermions, the projection methods such as RQMC are variational with respect to the nodal positions.

In principle, if one knows the exact nodal surface, then the solution will be exact. However, knowing the exact nodes means knowing the density matrix itself, which is the final goal of the simulations. Thus, nodal surface has to be approximated by the nodes of the trial wave function and usually is defined by the Slater determinant. In a real simulation of many-body, the nodal surface can be extremely complex, despite that, comparing to VMC, the fixed node approximation dramatically improves energies, however other properties such as the momentum distribution may not be improved to the same order.

Trial wave function

The choice of trial wave function is critical, especially in VMC calculations. The power of Quantum Monte Carlo methods lies in the flexibility of the form of the trial wave function.

Slater Jastrow

The trial wave function in most of quantum Monte Carlo simulations is based on the Slater determinant of single particle orbitals, det(ψ σ k (r i )), first introduced in eq. 1.3. The next step will be to add the pair correlation by simply multiplying the Slater determinant by symmetric pair function i<j f (r ij ), such that f (r ij ) > 0 everywhere. For mathematical convenience one can use the logarithm u(r) =ln(f (r)). Therefore, the pair product trial wave function can be written as:

Ψ SJ (R) = σ det(ψ σ k (r i )) exp[- i<j u(r ij )], (1.71)
where u(r ij ) can be understood as a "pseudopotential" acting between particles i and j at a distance r apart. Pair-product trial functions are employed because they are quite accurate in studies of solid hydrogen [START_REF] Ceperley | Ground state of solid hydrogen at high pressures[END_REF]. The trial function of Eq. (1.71) contains two pseudopotentials, which act between pairs of electrons, and between electrons and nuclei and depend parametrically on nuclei coordinates, R I = {R 1 , ..., R Np },

Ψ SJ (R; R I ) = σ det(ψ σ k (r i )) exp - Ne i=1 1 2 Ne j =i u ee (r ij ) - Np I=1 u ep (|r i -R I |) . (1.72)
In order to reduce the number of variational parameters to a minimum, the randomphase approximation (RPA) for pseudopotentials [START_REF] Gaskell | The Collective Treatment of a Fermi Gas: II[END_REF] was employed. Derived for the electron gas, this approximation assumes that the Hamiltonian of the system is a sum of a short-range interaction among electrons and a long range part described by collective oscillations (plasmons). As it was shown in [START_REF] Ceperley | Ground state of solid hydrogen at high pressures[END_REF], to provide lower energies for solid hydrogen, the two body pseudopotential can be supplemented by gaussian functions:

u α (r) = u RP A α (r) + λ α 2b exp[-(r/w α 2b ) 2 ], (1.73)
with variational parameters λ α 2b , w α 2b and α = (ee, ep), which makes it four variational parameters to be optimised.

Backflow transformation

The same two body correction can be derived iteratively, using the Feynman-Kac formula. Starting from a determinant of single electron orbitals which will be assumed as an initial ansatz, Ψ T (R), the projection in eq. 1.35 can be expressed via the density matrix in eq. 1.63

Ψ 0 (R) ∝ CΨ T (R) exp - β 0 dτ E L (R(τ )) DRW .
(1.74)

Consider that for any stochastic process, one can write the average of the exponent as the exponential of the cumulant expansion. Then, the first iteration generates a bosonic (symmetric) two body correlation function (Jastrow) while the next iteration naturally provides the backflow transformation of the orbitals and a three-body bosonic correlation term [START_REF] Kwon | Transient-estimate Monte Carlo in the twodimensional electron gas[END_REF][START_REF] Kwon | Effects of backflow correlation in the threedimensional electron gas: Quantum Monte Carlo study[END_REF][START_REF] Holzmann | Backflow correlations for the electron gas and metallic hydrogen[END_REF]. The second iteration suggests the backflow transformation of the orbitals:

x i = r i + j ν ij (|r ij |)r ij , (1.75)
where ν ij are the electron-electron and electron-proton backflow functions that must be parameterized. When the single body-orbitals in the determinant det(ψ σ k (x i )) are expressed in terms of the quasiparticle (QP) coordinates x i , the nodal surfaces of the trial wave function become explicitly dependent on the backflow functions ν, a crucial characteristic for electron systems which will provide more accurate energy. Similar to the two body term Eq. (1.73), the RPA was employed to get an analytical expression for backflow function plus gaussian parametrization,

x i = r i + Ne j =i [y RP A ee (r ij ) + ν ee (r ij )(r i -r j )] + Np I=1 [y RP A ee (|r i -R I |) + ν ep (|r i -R I |)(r i -R I )], (1.76) where ν α (r) = λ α b exp[-((r -r α b )/w α b ) 2
] and again α = (ee, ep). Besides the backflow correction, the second iteration also gives the three body correlation factor. However, it was shown [START_REF] Holzmann | Backflow correlations for the electron gas and metallic hydrogen[END_REF] that the inclusion of three-body terms does not noticeably affect the energies when using projector methods like RQMC or DMC.

For high pressure hydrogen, the described above trial wave function with the Slater orbitals coming from different exchange-correlation functionals was extensively tested and compared with each other in ref. [START_REF] Pierleoni | Trial wave functions for high-pressure metallic hydrogen[END_REF]. It has been shown by the authors that for all different parts of the wave function described in this section, combined with the Slater determinant obtained from the DFT-PBE functional, the QMC ground state energy is the lowest.

Size effects

The goal of computational physics is to calculate the properties of the real systems, which can contain on the order of 10 23 electrons. From first sight, it seems impossible to simulate such large systems. Indeed, typical QMC calculations are limited to fewer than several thousand electrons. Nonetheless, as discussed in [START_REF] Lin | Twist-averaged boundary conditions in continuum quantum Monte Carlo algorithms[END_REF][START_REF] Chiesa | Finite-size error in many-body simulations with long-range interactions[END_REF][START_REF] Holzmann | Theory of finite size effects for electronic quantum Monte Carlo calculations of liquids and solids[END_REF], there are different methods that one can use to extract properties in thermodynamic limit from the small simulation cell size.

Twist averaged boundary conditions

By analogy with the single-particle theory, where each single-particle orbital can be taken of the Bloch form ψ k (r) = exp[ik • r]u k (r) (see section 1.3.2), one can write the generalization of the Bloch theorem for a many-body wave function. Ψ(r 1 + Lx, r 2 , ..., r 3 ) = e iθ(x) Ψ(r 1 , r 2 , ..., r 3 ), (1.77) if θ = 0 one restore the PBC, if θ = π the boundary conditions are called anti-periodic boundary conditions (ABC) and the general θ = 0 are twisted boundary conditions (TBC).

The "twist", or the phase can be written as a product of supercell momentum and the cell vector θ(x) = K • Lx. θ(x) is a 3D vector, so it has components θ i , i = {1, 2, 3}. It is straightforward that most observables are (triply) periodic in twist A(θ i + 2π) = A(θ i ) (does not hold in the presence of magnetic field), from where it follows that each component of the twist can be restricted to the range:

-π < θ i π.

(1.78)

It has been shown [START_REF] Lin | Twist-averaged boundary conditions in continuum quantum Monte Carlo algorithms[END_REF], that finite-size effects are much reduced if the twist angle is averaged over, which is called twist averaged boundary conditions (TABC).

A T ABC ≡ 1 (2π) 3 π -π d 3 θ(x) Ψ θ | Â|Ψ θ . (1.79)
As can be seen, TABC is equivalent to integrating over the Brillouin zone in single particle representation, thus removing only independent particles finite size effects. Therefore, by doing TABC in QMC one cannot completely eliminate finite size effects, as further corrections are arising from many particle correlations.

Corrections arising from two-particle correlations

In the following we will discuss finite size error and correction schemes for electronic energies arising from the two-particle correlations [START_REF] Chiesa | Finite-size error in many-body simulations with long-range interactions[END_REF][START_REF] Holzmann | Theory of finite size effects for electronic quantum Monte Carlo calculations of liquids and solids[END_REF], we will restrict our discussion to a single component system to simplify the notation. One can understand the potential energy finite size effects writing the electron-electron interaction energy:

V N = 1 2Ω k =0 v k [S N (k) -1], (1.80)
where v k = 4π k 2 is the Fourier transform of the Coulomb potential, S(k) is the charged structure factor,

S(k) = 1 N | j e ir j • k | 2 = 1 N ρ k ρ -k , (1.81)
where the sum j is over all particles and ρ k = i e ikr i are the collective density fluctuations. The dominant finite-size error of the potential energy comes from the fact that the smallest k in eq. 1.80 is limited by the simulation cell size, which can be written as:

[ lim

N →∞ V N ] -V N = 1 2 dk (2π) 3 v k [S ∞ (k) -1] - 1 2Ω k =0 v k [S N (k) -1].
(1.82)

The leading order size correction, assuming S(k → 0) = 0, which holds for any chargeneutral system, is given by the Madelung constant, v M ,

v M 2 = - dk (2π) 3 - 1 Ω k =0 v k 2 (1.83)
with the remaining subleading term of the potential energy corrections being

∆V N = 1 2 dk (2π) 3 v k S ∞ (k) - 1 2Ω k =0 v k S N (k).
(

Practically, subleading corrections from eq. 1.84 may be evaluated by integrating asymptotic expansions of the structure factor around k = 0, effectively taking into account the contributions only from the volume element around the origin.

The kinetic energy expressed as a function of the momentum distribution writes as

T N = 1 N k k 2 2 n k , (1.85)
where the momentum distribution,

n k = 1 Ω dr dr e ik • (r-r ) ρ(r, r ), (1.86)
is written in terms of the single particle density matrix, ρ(r, r ). Again, the finite-size error is due to the discrete underlying mesh in Fourier space

∆T N = T ∞ -T N = 1 2ρ dk (2π) 3 k 2 n ∞ k - 1 Ω k k 2 n N k (1.87)
Note that, since the direct interpolation to the continuum for the momentum distribution is nontrivial. However, as long as one is only interested in the kinetic energy, it is easier to express it in terms of a different purely local estimators and discuss size effects using them.

Considering that the electronic ground state can be described with the Slater-Jastrow form, Ψ T = De -U , where D is a Slater determinant and ensures the antisymmetry of fermions and U is a many-body symmetric function (for electrons). The kinetic energy contribution involving only the Jastrow factor, U = 1 2Ω k u k ρ k ρ -k , with u k being the Fourier transform of Jastrow potential, can be expressed as,

T U = 1 N i 1 2 [∇ i U ] 2 (1.88) = 1 2N Ω 2 i k =0 (ik)u k e ikr i ρ -k k =0 (ik )u k e ik r i ρ -k = - 1 2N Ω 2 k =0,k =0 (k • k )u k u k ρ k+k ρ -k ρ -k 1 Ω k =0 k 2 2 ρu k u -k S N (k),
where the singular contributions involving ρ k=0 ≡ N e were isolated and all the terms with k = -k have been neglected, which correspond to the RPA [START_REF] Gaskell | The Collective Treatment of a Fermi Gas: II[END_REF] approximation. The integration error then, by analogy to the potential energy, becomes,

∆T U = ρ 2 dk (2π) 3 k 2 u k u -k S ∞ (k) - 1 Ω k =0 k 2 u k .u -k S N (k) (1.89)
Again, in practice to analyze the finite size corrections one can interpolate the Jastrow potential, u k , from its values at discrete k points to all values of k, which allows applying the same procedure as for the potential energy correction via the integration of the asymptotic expansion of the structure factor around k = 0. The remaining part of the finite size correction scheme comes from backflow. Significantly improving the accuracy of QMC calculations, backflow is an important part of the QMC wave function. The corrections coming from backflow coordinates can be derived in the same fashion as the corrections presented above. However, being quite lengthy, will not be discussed here. Detailed derivation and discussion of the correction schemes can be found in ref. [START_REF] Holzmann | Theory of finite size effects for electronic quantum Monte Carlo calculations of liquids and solids[END_REF].

Summary

The foundations of the quantum mechanical simulation techniques for the electronic ground state properties were established. We discussed independent particle approxima-tions: Hartree-Fock (HF) and density functional theory (DFT). The main problem of the Hartree-Fock and DFT methods is the treatment of electron correlation. In the HF method, electron correlations beyond a mean-field picture are entirely neglected. In DFT they are included via a functional of the density E xc [n]. DFT methods require careful calibration to establish their accuracy on a case by case basis.

The major part of the chapter was devoted to the quantum Monte Carlo (QMC). We introduced the variational Monte Carlo (VMC) and reptation quantum Monte Carlo (RQMC). A very practical feature of Quantum Monte Carlo methods is that they provide a direct measure of the accuracy obtained via the variance of the local energy of the wave function. The VMC method provides a way to statistically evaluate observables for electronic systems using a given many-body trial wave function. Although, the wave functions must be designed using physical insight, as the VMC method does not provide direct insight into the parameterization of a wave function.

The fixed-node RQMC method is a variational method for finding the fermionic ground state of a Hamiltonian, with the approximation being the use of an approximate nodal surface from a trial wave function. In practice, the optimized VMC wave function is sufficient to obtain a good quality trial wave function for RQMC calculations.

We have presented different pieces of the trial wave function, with the most important being Slater-Jastrow plus backflow. In our simulation, the wave function is to the large extend analytical with only a few optimization parameters. Finally, we have introduced a technique of treating the finite simulation cell size effects for potential and kinetic energy, which is based on extrapolating the structure factor.

The presented methods are used extensively throughout the thesis.

Chapter 2

Excited States

For an interacting many electron system, besides the ground state properties, it is important to know the electronic excited state related to spectroscopic properties. Properties that can be measured in photoemission, transport, and tunneling experiments involve creating an excited particle (electron or hole) above the ground state. It is as an N + 1 particle problem that requires a different theoretical treatment from those of the ground-state approaches [START_REF] Hybertsen | First-principles theory of quasiparticles: Calculation of band gaps in semiconductors and insulators[END_REF]. Theoretical understanding of the optical properties, such as absorption, from first principles is yet a further challenge because it is an N + 2 particle problem [START_REF] Albrecht | Ab initio calculation of excitonic effects in the optical spectra of semiconductors[END_REF]: when an electron is promoted to the conduction state it leaves the hole in the valence band and an electron-hole pair is created. In this case, the electron-hole interaction needs to be involved, which can be quite important in many systems [START_REF] Rohlfing | Electron-hole excitations in semiconductors and insulators[END_REF].

In the present chapter, the different ways to calculate excitation energies will be discussed. We will limit ourselves to two types of excitations: (charged) quasiparticle and neutral. First, we will focus on the fundamental quasiparticle gap, which is defined as the energy difference required to add and subtract an electron. We start by discussing the fundamental gap within the single electron theory, further switching to many body, where we provide a new method to obtain a gap within QMC. A very important result of this chapter is the treatment of finite size effects on addition and removal energies, which we formally demonstrate to be related to the dielectric constant of the material. We apply the new methodology to determine gaps of silicon, carbon diamond, and ideal crystalline hydrogen. In the second part of this chapter, we discuss the neutral or optical excitations. Within the QMC, we construct a wave function to model the excited states, which are orthogonal to the lower lying states. Applying it to compute optical excitations of solid crystalline hydrogen, shows that by using such wave function we can capture an electron-hole correlation and lower the gap.

Having a lot of reviews (such e.g. [START_REF] Onida | Electronic excitations: densityfunctional versus many-body Green's-function approaches[END_REF]) and books (such e.g [START_REF] Louie | Conceptual Foundations of materials[END_REF]) written on determining excited states within the perturbation theory, in QMC this field is still in the stage of development with some information on the excited state provided in the general QMC reviews [START_REF] Foulkes | Quantum Monte Carlo simulations of solids[END_REF][START_REF] Kolorenč | Applications of quantum Monte Carlo methods in condensed systems[END_REF].

Electron addition/removal

The first kind of excitations discussed in this chapter will be the excitations by adding or removing an electron. The process can be illustrated as follows: consider, for example, a missing electron. Positive charge that is left behind will induce other electrons to oscillate, effectively creating an electron-hole pair. The induced oscillations will change the kinetic energy of the outgoing electron. An experiment that measures the kinetic energy and angle (momentum) of an electron emitted by the photon irradiation is called angle-resolved photoemission spectroscopy [START_REF] Hüfner | Photoelectron Spectroscopy. Advanced Texts in Physics[END_REF]. Combination with an inverse process of electron irradiation and photon emission, inverse photoemission (IPES), makes these two methods to be the most direct measurements for determining N ± 1 excitations.

In insulators and semiconductors, it is energetically more favorable to remove an electron than to add one, the resulting difference is a fundamental gap defined as [START_REF] Martin | Interacting Electrons[END_REF],

E g = [E 0 (N + 1) -E 0 (N )] -[E 0 (N ) -E 0 (N -1)], (2.1)
the lowest possible energy to add and electron minus the lowest possible energy to remove one. The energies involved are all at ground states without the possibility of decay. Correct determination of the gap became like some sort of a test for the theoretical methods aiming to correctly describe the excited states. We will see further that for some theories it is by definition not possible to correctly predict the fundamental gap.

Single electron fundamental gap

Within the Hartree-Fock theory there exist the Koopmans' theorem [START_REF] Koopmans | Über die zuordnung von wellenfunktionen und eigenwerten zu den einzelnen elektronen eines atoms[END_REF], introduced in section 1.2, which states that the HF eigenvalues correspond to the total energy differences when an electron is added or subtracted (eq. 1.12) assuming that the orbitals do not relax when an electron is added or removed.

In turn, when the Kohn-Sham system is considered, the Koopmans' theorem holds only for the highest occupied eigenvalue N in a finite system. Generalizing the DFT to fractional particle number [START_REF] Perdew | Density-functional theory for fractional particle number: Derivative discontinuities of the energy[END_REF], the concept of fractional particle number with occupations f i can be introduced. The Kohn-Sham eigenvalues can be then written according to Janak's theorem [START_REF] Janak | Proof that ∂e ∂n i = in density-functional theory[END_REF] as,

KS i = ∂E ∂f i . (2.2)
Thus the ionization energy E N -E N -1 is then,

-I = E N -E N -1 = N N -1 df N ∂E ∂f N = N N -1 df N N = N , (2.3)
which is true if N does not depend on f N . For the exact functional, E(N ) is linear within the integration region, with the slope discontinuity at integer N [START_REF] Perdew | Density-functional theory for fractional particle number: Derivative discontinuities of the energy[END_REF]. This reasoning, however, does not apply to the electron addition energy, N +1 = E N +1 -E N as N +1 and N are not eigenvalues of the same KS hamiltonian. The fundamental gap, expressed in terms of Kohn-Sham eigenvalues, therefore is,

E g = N +1 (N + 1) -N (N ), (2.4)
where the number of electrons is indicated in parenthesis. However, the KS gap, KS g = N +1 (N ) -N (N ), calculated for N -electron KS hamiltonian, is just a part of the gap defined in eq. 2.4,

E g = KS g + N +1 (N + 1) -N +1 (N ), (2.5)
where the correction ∆ = N +1 (N + 1) -N +1 (N ) is called the derivative discontinuity.

Knowing that the density changes by an infinitesimal amount when adding an electron, the discontinuity should be directly related to the rigid shift in the KS potential. To illustrate this consider the exchange correlation (XC) potential v xc ([n], r) within LDA, which is the functional derivative of E xc with respect to density n(r) (see eq. 1.19 and eq. 1.18),

v xc ([n], r) = HEG xc ([n], r) + n(r) δ HEG xc ([n(r)], r) δn(r) . (2.6)
In an insulator, the functional derivative in the second term of eq. 2.6 should be discontinuous at a band gap, where the nature of the states changes discontinuously as a function of n. This is a "derivative discontinuity" where the Kohn-Sham potential for all the electrons changes by a constant amount when a single electron is added. Therefore, if one thinks of using directly eq. 2.4 to determine the gap, then when computing the total energy of the system with an extra electron there should be a derivative discontinuity artificially incorporated into the XC potential [START_REF] Grüning | Density functionals from manybody perturbation theory: The band gap for semiconductors and insulators[END_REF]. Standard functionals like LDA do not have this property, thus, even within exact KS, the electron addition-removal gap will differ from the experimental one.

Therefore, it is strictly speaking inadequate to use KS eigenvalues as electron addition and removal energies. Nonetheless, KS eigenvalues and eigenstates are well defined within the theory and can be used to construct physically meaningful quantities, which is done in quantum Monte Carlo and many-body perturbation approaches, where KS eigenstates are used as an input.

Although, hybrid functionals, which form a so-called generalized Kohn-Sham (GKS) theory, involve a part of non-local HF potential and obey the Koopmans' theorem and can cure the derivative discontinuity problem.

Quasiparticles

Returning to the HF theory, where the eigenvalues have a physical meaning and the gap can be computed straightforwardly with eq. 2.1, in reality, it is found that the gap is always much larger than the experimental gap. For example, for carbon diamond, the HF gap is predicted to be twice the experimental value (12 eV and 5.5 eV [START_REF] Brener | Correlated Hartree-Fock energy bands for diamond[END_REF]). Indeed, in HF an extra electron feels the bare Coulomb interaction, which tends to localize it excessively. Although, in reality, the interaction should be screened by the relaxation of the other electrons. In this picture, we can define the quasiparticles (QP), that interact through the screened Coulomb interaction, which is weaker and more short-ranged, than the bare one. Each quasiparticle can be seen as an electron plus its screening cloud, which moves around when an electron is added into the system.

Quasiparticle excitations in QMC

Within QMC the quasiparticle gap in eq. 2.1 is immediately accessible because electrons are treated explicitly contrary to mean field methods like DFT and HF. One of the first applications of QMC to compute a fundamental gap was by Ceperley and Alder (1987) [START_REF] Ceperley | Ground state of solid hydrogen at high pressures[END_REF] who calculated the energy gap closure with pressure of P a3 solid molecular hydrogen crystal. In this pioneering work, the authors observed significant size effects, however, even considering all the approximations, i.e. small systems, zero momentum excitations, and fixed node approximation, the resulting value was not very far from the present GW calculations [START_REF] Liang | Band gaps and effective oscillator models for solid hydrogen and h 2 o ice at high pressure[END_REF]. Although, in the literature, the prevailing number of QMC calculations is addressing charge neutral excitations, where faster convergence with respect to the size of the supercell is expected [START_REF] Mitáš | Quantum monte carlo of nitrogen: Atom, dimer, atomic, and molecular solids[END_REF][START_REF] Williamson | Diffusion quantum Monte Carlo calculations of the excited states of silicon[END_REF][START_REF] Towler | Minimum principles and level splitting in quantum Monte Carlo excitation energies: Application to diamond[END_REF][START_REF] Kolorenč | Quantum monte carlo calculations of structural properties of feo under pressure[END_REF][START_REF] Ma | Excited state calculations in solids by auxiliary-field quantum monte carlo[END_REF][START_REF] Wagner | Effect of electron correlation on the electronic structure and spin-lattice coupling of high-T c cuprates: Quantum monte carlo calculations[END_REF][START_REF] Yu | Systematic assessment of errors in diffusion Monte Carlo calculations of semiconductors: case study of zinc selenide and zinc oxide[END_REF][START_REF] Zheng | Computation of the correlated metal-insulator transition in vanadium dioxide from first principles[END_REF][START_REF] Frank | Many-Body Quantum Monte Carlo Study of 2D Materials: Cohesion and Band Gap in Single-Layer Phosphorene[END_REF], in the following sections we will focus mainly on the QP excitations within the QMC and will introduce the scheme to correct for the size effects.

Being a good quantum number, the twist θ (introduced in section 1.4.4.1) can be used to resolve the QP energies into the momentum space. We can define the energy differences,

∆ N (θ) = E N (θ) -E N -1 (θ), (2.7)
the quasiparticle gap is then,

E QP (θ, θ ) = ∆ N +1 (θ) -∆ N (θ ), (2.8)
and the fundamental gap, being the minimal energy gap, is therefore,

E g (θ, θ ) = min φ ∆ N +1 (φ) -max φ ∆ N (φ ), (2.9)
where θ and θ correspond to the twists at which the minimum and the maximum are realized. If θ = θ , the gap is indirect, while when θ = θ the gap can be direct if the simulation cell corresponds to the primitive cell. If the simulation cell is the multiple of the primitive cells, see Appendix C on how to determine the crystal momentum. The vertical or direct gap can be accessed by,

E dir g (θ) = min φ [∆ N +1 (φ) -∆ N (φ)] (2.10)
if band folding in the supercell can be excluded, e.g. if Bloch orbitals of the primitive unit cell are used to build the wave functions.

Band structure

In a periodic solid, where Bloch's theorem is applicable, plotting the energy of the quasiparticle as a function of its quantum number, e.g. momentum and band index, will produce a band structure plot. Note, however, that as the relaxation of the screening cloud around the quasiparticle is finite, the QP energy broadens and becomes a spectrum. Such spectrum is called a spectral function and can be directly measured via the photoemission spectroscopy, discussed above. In the traditional band structure plot, only the first quasiparticle peak is represented, which corresponds to the ground state of added or removed electron.

Figure 2.1 illustrates an example of the QP addition/removal RQMC energies for C2/c hydrogen at approximately 234 GPa plotted on top of the DFT band structure computed with vdW-DF XC approximation. The DFT band structure was shifted in order to match the QMC gap. After the shift, we notice that the QP energies lie on top of the DFT eigenvalues. Due to the limitations of our QMC code and finite size effects we are not able to run simulations of hydrogen in the primitive cell and supercells have to be used. Therefore, to match supercell twists and crystal momenta, we need to unfold the bands (see Appendix C). This is the reason for sparse QMC points on the band structure for crystalline hydrogen. In the previous section 1.4.4, I introduced a key quantity in understanding size effects: the static structure factor and its long wavelength behavior, S Ne (k) = ρ -k ρ k /N e where ρ k = j e ik • r j is the Fourier transform of the instantaneous electron density. In the famous work of Ceperley and Alder in 1987 [START_REF] Martin | Interacting Electrons[END_REF][START_REF] Ceperley | Ground state of solid hydrogen at high pressures[END_REF] it was demonstrated that the structure factor for a homogenous system obeys the bound,

S Ne (k) ≤ k 2 2ω p 1 - 1 k 1/2 , (2.11)
where ω p = 4πn e is the plasma frequency and k the static dielectric function at wavevector k (To simplify the notations, we will suppress the dependence on the wave vector in the following). This inequality is derived from the sum rules of the dynamic structure factor S(k, ω) . It implies that the structure factor must vanish quadratically as k → 0 [START_REF] Pines | The Theory of Quantum Liquids[END_REF].

Equality will be obtained if S(k, ω) reduces to a sharp peak at a single frequency, ω at small k. The 1/N e finite-size corrections of the energy per electron is a direct consequence of this behavior of S Ne (k) [START_REF] Chiesa | Finite-size error in many-body simulations with long-range interactions[END_REF]. However, these leading order corrections are not sufficient for excitation energies, since the energy gap is of the same order as finite size corrections to the total energy. As we will show below, the key to understanding the size effects of energy differences is encoded in the change of S Ne (k) as electrons are added or removed. In particular, the limiting behavior of S Ne±1 (k) as k → 0 will provide the dominant finite size correction.

For concreteness, we will assume a Slater-Jastrow form (see eq. 1.71) for the ground state wave function Ψ 0 ∝ D exp[-U ]. The determinant, D, is built out of Bloch orbitals, ψ qn (r) with q inside the first Brillouin zone, n is the band index, and U is, as defined in eq. 1.72, a general symmetric n-body correlation factor [START_REF] Holzmann | Theory of finite size effects for electronic quantum Monte Carlo calculations of liquids and solids[END_REF]. For simplicity we assume it is pairwise additive: U = i<j u(r i , r j ). Let us consider the action of e ik • r j on a single particle orbital ψ qn (r j ) in the Slater determinant of the ground state. In the limit of small k, this can be approximately written as ψ (q+k)n (r j ). According to Jacobi's formula, the determinant can be expanded in terms of its cofactors δD δψqn(r j ) as,

D = q,n δD δψ qn (r j ) ψ qn (r j ) (2.12)
Formally making the excitation means,

ρ k Ψ 0 ∝ j e ik • r j q,n δD δψ qn (r j ) ψ qn (r j )e -U (2.13) = j q,n δD δψ qn (r j ) e ik • r j ψ qn (r j )e -U ,
where the determinant was expanded for each j in the first summation. Notice that after the summation over j the resulting determinant vanishes for small k if the Bloch orbital (q + k, n) is already occupied in the ground state determinant. Considering N e ± 1 electron wave functions, Ψ 0 (N e ± 1; ±q, m) where N e corresponds to the insulating state with fully occupied bands in the Slater determinant, and qm denotes the additional particle/hole orbital, based on equation 2.13 we get,

lim k→0 ρ k Ψ 0 (N e ± 1; q, m) = lim k→0 j q,n δD δψ qn (r j ) e ik • r j ψ qn (r j )e -U = lim k→0 j δD δψ qn (r j ) ψ q+kn (r j )e -U ∼ ±Ψ 0 (N e ± 1; q + k, m), (2.14)
for k = 0 where different sign for particle or hole excitations on the r.h.s. is chosen to match the most common sign convention, e.g. of ref. [START_REF] Kohn | Interaction of Charged Particles in a Dielectric[END_REF]. The limit k → 0 is discontinuous since ρ k=0 Ψ 0 (N e ± 1; q, m) ≡ (N e ± 1)Ψ 0 (N e ± 1; q, m). Kohn [START_REF] Kohn | Interaction of Charged Particles in a Dielectric[END_REF][START_REF] Kohn | Efective mass theory in solids from a many-particle standpoint[END_REF] has pointed out that in the insulating state the matrix elements, lim

q →q Ψ 0 (N e ± 1; q, m)|ρ q-q |Ψ 0 (N e ± 1; q , m) = ± 1 , (2.15)
approach the inverse dielectric constant, -1 , up to a sign. Substituting eq. 2.14 into eq. 2.15, suggests the following,

± 1 ∼ lim q →q Ψ 0 (N e ± 1; q, m)|ρ q-q lim k→0 ρ k |Ψ 0 (N e ± 1; q, m) (2.16) = lim k→0 Ψ 0 (N e ± 1; q, m)|ρ -k ρ k |Ψ 0 (N e ± 1; q, m) = lim k→0 (N e ± 1)S Ne±1 (k),
where we put q = q + k. Permitting us to define a finite size behavior of the static structure factor of insulators,

lim k→0 S ± k = α ± + O(k 2 ), S ± k ≡ (N e ± 1)S Ne±1 (k) -N e S Ne (k), (2.17)
where the N e S Ne (k) term was added to get the size correction for the addition/removal energies and α ± is proportional to -1 . However, α ± in general differs from -1 unless eq. (2.14) is an exact equality. Figure 2.2 shows the behavior of S ± k for carbon and silicon crystals. Note that these functions extrapolate to a nonzero value as k → 0.

The long wavelength behavior of the structure factor, eq. ( 2.17), similarly to the total energy correction (see eq. 1.84), then gives rise to size corrections to excitation energies through the potential energy term,

d 3 k (2π) 3 - 1 Ω k =0 v k 2 S ± k α ± |v M | 2 , (2.18)
where we have defined the Madelung constant as in eq. 1.83,

v M = 1 Ω k =0 - d 3 k (2π) 3 v k ∼ L -1 ∼ N -1/3 e . (2.19)
For the Coulomb potential, v M is proportional to L -1 , the inverse linear extension of the simulation cell. The negative proportionality constant depends on the boundary conditions, e.g. cell geometry, and can be calculated by the Ewald image technique [START_REF] Ewald | Die berechnung optischer und elektrostatischer gitterpotentiale[END_REF].

Kinetic energy

Following ref. [START_REF] Holzmann | Theory of finite size effects for electronic quantum Monte Carlo calculations of liquids and solids[END_REF] and section 1.4.4, we now discuss the kinetic energy contribution [∇U ] 2 /2 which arises from electron correlation. For a two-body Jastrow U = vectors of the crystal, G. Isolating the singular contributions involving ρ k=0 ≡ N e in the spirit of the rotating (random) phase approximation (RPA) [START_REF] Gaskell | The Collective Treatment of a Fermi Gas: II[END_REF] and as in eq. 1.88 we have,

S ± k S + k S - k N=8 N=64 -1 /(1 + 0.41) (a) carbon 0.0 0.2 0.4 0.6 0.8 1.0 1.2 k 2 (a.u.) -0.4 -0.2 0.0 0.2 0.4 0.6 S ± k S + k S - k N=8 N=64 N=216 -1 /(1 + 0.57) (b) silicon
[∇U ] 2 = - 1 Ω 2 k =0,k =0 (k • k )u k u k ρ k+k ρ -k ρ -k 1 Ω 2 k =0 N e k 2 u 2 k ρ k ρ -k . (2.20)
Therefore, based on the asymptotic behavior of static structure factor S ± k , for systems with explicit long-range correlations u k ∼ k -2 , the kinetic energy will contribute also to the leading order size corrections with

d 3 k (2π) 3 - 1 Ω k =0 n e k 2 u 2 k 2 S ± k α ± c |v M | 2 , (2.21)
where c = lim

k→0 n e k 2 u 2 k /(v k ) (2.22)
is approximately given by the ratio of the 1/N e finite-size corrections of the kinetic to potential energy of the ground state energy per particle due to two-body correlations [START_REF] Holzmann | Theory of finite size effects for electronic quantum Monte Carlo calculations of liquids and solids[END_REF].

Total gap corrections from Coulomb singularity

Up to now, we have shown how the long range behavior of the structure factor and Jastrow factor can give rise to a 1/L correction to the excitation gap with a proportionality . Lines are fits to the low-k data. The horizontal lines mark experimental inverse dielectric constants.

factor determined by the structure factor changes. In the following, we will further demonstrate that, given the trial wave functions coincide with the exact ground state wave function for N e and N e ± 1 electrons, this proportionality factor is indeed given by the dielectric constant

∆ ∞ -∆ V = |v M | + O 1 Ω , (2.23)
as phenomenologicaly assumed in previous works [START_REF] Engel | Quasiparticle bands in a two-dimensional crystal found by gw and quantum monte carlo calculations[END_REF][START_REF] Hunt | Quantum Monte Carlo calculations of energy gaps from first principles[END_REF].

In the following this will be proved by an independent argument based on commutation relations. Let us denote the exact insulating ground state of the N e electron system as |Ψ Ne 0 , its energy as E Ne 0 , and the exact excited state of the N e ±1 electron system as |Ψ Ne±1 k with energy E Ne±1 k ; k indicates that the additional/subtracted electron adds/subtracts the crystal momentum k. Considering the creation and annihilation operators a † k and a k for plane wave states of wave vector k acting in the following way,

|Ψ Ne+1 k ∝ a † k |Ψ Ne 0 , (2.24)
|Ψ Ne-1 k ∝ a k |Ψ Ne 0 .
We have for electron addition and removal energies,

E Ne+1 k -E Ne 0 = Ψ Ne+1 k | Ĥ, a † k |Ψ Ne 0 Ψ Ne+1 k |a † k |Ψ Ne 0 (2.25)
for particle and

E Ne-1 k -E Ne 0 = Ψ Ne-1 k | Ĥ, a k |Ψ Ne 0 Ψ Ne-1 k |a k |Ψ Ne 0 (2.26)
for hole excitations.

In second quantization , the Hamiltonian, Ĥ = Te + Vee , is given by

Te = k k 2 2 a † k a k + G u(G)a † k+G a k , (2.27) Vee = 1 2V q =0 v q [ρ q ρ -q -N e ] , (2.28)
u(G) the periodic crystal potential, and v q is the Coulomb potential between electrons, ρ q = k a † k+q a k , and N e = k a † k a k . Considering the commutator relations for fermions,

{a i , a † j } ≡ a i a † j + a † j a i = δ ij , (2.29
)

{a † i , a † j } = {a i , a j } = 0,
the commutator involving the single-particle energy term becomes

Te , a † k = k 2 2 a † k + G u(G)a † G+k . (2.30)
There are corresponding terms for hole excitations, but none of these terms involve singular contributions responsible for anomalous size effects, so that these terms do not contribute at leading order. However,

V ee , a † k = 1 V q =0 v q ρ q a † k-q -1 (2.31) and [V ee , a k ] = - 1 V q =0 v q ρ q a k+q (2.32)
involve terms approaching the Coulomb singularity, v q ∼ q -2 → ∞ for q → 0.

From these terms we get the leading order size corrections by noting that lim k,q→0

Ψ Ne+1 k |ρ q a † k-q |Ψ Ne 0 Ψ Ne+1 k |a † k |Ψ Ne 0 = 1 2 1 + 1 (2.33) and lim k,q→0 Ψ Ne-1 k |ρ q a k+q |Ψ Ne 0 Ψ Ne-1 k |a k |Ψ Ne 0 = - 1 2 1 + 1 . (2.34)
Both relations can be obtained by extending Kohn's diagrammatic approach [START_REF] Kohn | Interaction of Charged Particles in a Dielectric[END_REF] (see supplementary information of [START_REF] Yang | Electronic band gaps from Quantum Monte Carlo methods[END_REF]). Integrating around the v q singularity for small q in eq. ( 2.31), we obtain the leading order finite size corrections. As before, this involves the Madelung constant, Eq (2.19). In the particle channel we get |v M | Previous, heuristic approaches [START_REF] Hunt | Quantum Monte Carlo calculations of energy gaps from first principles[END_REF] have suggested that one can use experimental or DFT values of the dielectric constant for finite-size extrapolation. Our approach further suggests that this value can be determined from the QMC structure factor extrapolated to zero wave vector 2

≡ (1 + c) lim Ne→∞ lim k→0 S + k + S - k , (2.35)
with the singular behavior of the Jastrow factor determining c. We emphasize that the order of the limits involved above is crucial. An independent estimate of the dielectric constant is based on the inequality of eq. 2.11. We can bound and estimate the value of dielectric constant using the structure factor of the insulating ground state. By extrapolating 1 -Γ 2 k vs. k to k = 0 we obtain an upper bound to the inverse dielectric constant, where Γ k ≡ 2ω p S Ne (k)/k 2 . This involves only the extensive part of the density-density correlations, thus, it is less sensitive to noise and has much smaller statistical uncertainty. In fig. 2.3, we show that for C and Si, this upper bound gives accurate values of the dielectric constant (see section 2.1.7 for a more quantitative comparison).

Twist correction of two particle correlations

The above size effects explain the leading order 1/L correction to the single particle gap. However, as we will see in our results, the asymptotic region, where this law can be reliably applied, may still be difficult to reach for currently used system sizes and next-to-leading order effects are important. Here, we show that an important part can be corrected for, by further restoring the full symmetry properties in the contribution of the direct Coulomb interaction.

For non homogeneous systems, it is convenient to separate the mean density from its fluctuating components in the static structure factor [START_REF] Holzmann | Backflow correlations for the electron gas and metallic hydrogen[END_REF], i.e.

S Ne (k) = 1 N e ρ k ρ -k + δS Ne (k) (2.36) δS Ne (k) = 1 N e (ρ k -ρ k ) (ρ -k -ρ -k ) (2.37)
For crystals with periodic density distributions, the Fourier components of the mean density, ρ k , only contribute for reciprocal lattice vectors, k ∈ G. The long wavelength behavior of the structure factor is entirely due to the fluctuating part δS Ne (k), which therefore contains the leading order size effects [START_REF] Holzmann | Theory of finite size effects for electronic quantum Monte Carlo calculations of liquids and solids[END_REF]. However, the mean single particle density, ρ(r) = V -1 k ρ k e ik • r , of the finite system may significantly differ from the infinite one, particularly in cases where the supercell is not compatible with the full symmetry group of the crystal.

Averaging over twisted boundary conditions is designed to restore the symmetry of the crystal and thus accelerate the convergence of single particle densities to the thermodynamic limit. In the following, we denote the twist averaged expectation value by

O ≡ 1 M θ θ O Ne,θ (2.38)
where we have explicitly indicated the N e and θ dependence on the expectation value on the r.h.s. For any single particle theory, ρ(r) approaches its thermodynamic limit for calculations at fixed N e by averaging over a dense grid of twist angles (M θ → ∞). Within many-body calculations, twist-averaging [START_REF] Lin | Twist-averaged boundary conditions in continuum quantum Monte Carlo algorithms[END_REF] takes over a large part of the size effects to any observable linear in the density. Here, we extend this approach to correct also the quadratic expression entering the two-body contributions of the total energy. For the potential energy, this correction to the twist converged QMC calculation is

δV s Ne = 1 2V k v k δC(k) δC(k) = ρ k ρ -k -ρ k ρ -k . (2.39)
For the ground state energies, this correction provides only a small improvement over our previous correction [START_REF] Chiesa | Finite-size error in many-body simulations with long-range interactions[END_REF][START_REF] Holzmann | Theory of finite size effects for electronic quantum Monte Carlo calculations of liquids and solids[END_REF].

For the gap, many terms entering eq. (2.39) cancel, and the expression can be simplified. Let us consider the case of adding/removing one electron at twist ψ to the insulating ground state, denoting Π ± k the difference of the respective densities

Π ± k ≡ ρ k Ne±1,ψ -ρ k Ne,ψ (2.40)
In the thermodynamic limit, the density of the ground state system with N e electrons coincides with the twist averaged ground state density ρ k , whereas we obtain ρ k + Π ± k for the density of the N e ± 1 electron system. Inserting into eq. (2.39), we obtain the correction for the δV s Ne±1,ψ ,

δV s Ne±1,ψ = 1 2V kθ v k (ρ k + Π ± k )(ρ -k + Π ± -k ) -( ρ k Ne,θ + Π ± k δ ψθ )( ρ -k Ne,θ + Π ± -k δ ψθ ) = 1 2V k v k Π ± k ρ -k -ρ -k Ne,ψ + (ρ k -ρ k Ne,ψ ) Π ± -k + δC(k) , (2.41)
which results into the difference between the two states,

δV s Ne±1,ψ -δV s Ne = 1 V k∈G v k Re (ρ k -ρ k Ne,ψ ) Π ± -k (2.42)
where only wave vectors of the reciprocal crystal lattice contribute to the sum. The corresponding finite size correction for the gap, denoted by δ∆ s in the following, is order 1/N e or smaller, mainly determined by the changes of the ground state densities at the first Bragg-peaks due to twist averaging. Equation (2.42) can be understood as follows: it corrects the direct Coulomb interaction between the electron/hole in the excited state (Π ± ) with the unexcited electrons. The density of those electrons is expected to change by ρ kρ k Ne,ψ in the thermodynamic limit.

Converged ground state densities are naturally calculated within Grand-Canonical twist averaged boundary conditions. It is straightforward to apply the correction eq. (2.42) to all excitation energies. Alternatively, the corresponding DFT densities may be used. This removes the stochastic error at the cost of introducing a small bias in the next-to-leading order size correction.

Grand-Canonical twist averaged boundary condition (GCTABC)

In the following, we consider N e electrons in a perfect crystal, neglecting both, zero point motion of the ions and temperature effects. A uniform background charge (depending on N e ) is added to assure global charge neutrality when adding or substracting electrons, without introducing defects in the ionic crystal. The fundamental gap, eq. 2.1, is unaffected by the background energy, because the background charge needed when adding an electron cancels against the one needed when removing an electron. Periodic boundary conditions of the charge densities are used to eliminate surface effects.

The energetic cost of adding an electron to the system at a fixed volume, V = L 3 , defines the chemical potential

µ + Ne = E 0 (N e + 1) -E 0 (N e ). ( 2

.43)

A non-vanishing gap implies a discontinuity in the chemical potential from eq. 2.1.

It is convenient to work in the grand-canonical (GC) ensemble. There, the chemical potential µ is treated as an independent variable. At the ground state and fixed volume in the GC ensemble, the number of electrons assumes the value Ne (µ) that minimize the grand potential, Ω(µ) = min

Ne

[E 0 (N e ) -µN e ] . (2.44)
Insulators then represent an incompressible electronic state within the gap where ∂N e /∂µ = 0 in the thermodynamic limit, i.e. when V → ∞. Periodic boundary conditions of the N e -body density are guaranteed by imposing twisted boundary conditions on many-body wave function where the twist angle θ is applied to the phase of the many body wave function as an electron is moved across the supercell, as defined in eq. 1.77 Different twist angles modify the ground state energy E 0 (N e , θ) and twist averaging can significantly accelerate the convergence to the thermodynamic limit [START_REF] Lin | Twist-averaged boundary conditions in continuum quantum Monte Carlo algorithms[END_REF].

Within the grand-canonical ensemble [START_REF] Chiesa | Finite-size error in many-body simulations with long-range interactions[END_REF][START_REF] Holzmann | Theory of finite size effects for electronic quantum Monte Carlo calculations of liquids and solids[END_REF], the number of electrons Ne (µ, θ) will depend on θ for given chemical potential µ. Single particle finite size effects are reduced by averaging over the twists. For the total number of independent twist angles M θ , the following quantities can be defined,

ω(µ) = 1 M θ V θ e 0 (n e (µ)) -µn e (µ), (2.45
)

n e (µ) = (M θ V ) -1 θ Ne (µ, θ), (2.46)
e 0 (n e (µ)) = (M θ V ) -1 θ E 0 ( Ne (µ), θ), (2.47)
where ω(µ) is the free energy density, e 0 (µ) the energy density and n e (µ) the electron density. By eliminating µ from the last two equations we obtain the discontinuity in the energy derivative at n e corresponding to the energy gap

∆ gc = µ + -µ -= de dn e n + p - de dn e n - p , (2.48)
where the derivatives are computed at n e = n p = N p /V . For any single electron theory the electronic density, n e (µ), as well as the ground state energy density, e 0 (n e (µ)), coincide exactly with the corresponding thermodynamic limit values, where the sum over twists becomes an integral over the Brillouin zone. Size effects remaining after twist averaging are due to electron-electron correlations. Figure 2.4 illustrates e 0 (µ) and n e (µ) with result for solid hydrogen, computed from HSE functional and from QMC (see section 2.1.8 for details). The value of the band gap can be directly extracted from the width of the incompressible region. Alternatively, eliminating µ in favor of n e , and plotting e 0 as a function of n e , the fundamental gap is obtained by the discontinuity of the derivative, according to eq. (2.48). It is straightforward to generalize the fundamental gap for different symmetry sectors. For a perfect ionic crystal, the total momentum of the electrons modulo reciprocal lattice vectors is conserved. Imposing the total momentum of the electrons e.g. using Bloch type orbitals in the Slater-determinant, the full band structure in the Brillouin zone can be mapped out. For a spin-independent Hamiltonian, one can also impose the total spin to determine the fundamental gap in each spin sector (different twist angles for a spin up and down electrons may also be used in the grand-canonical sampling). In practice, the charge gap in the spinless sector can be determined by restricting to even changes in the number of electrons, equal in each spin component. The extensions of our definitions and formulas to this case are straightforward, e.g. ∆ Ne = [E 0 (N e + 2) + E 0 (N e -2) -2E 0 (N e )]/2, and not explicitly detailed in the following.

Results: silicon and carbon diamond

In this section, two insulating solids will be discussed: silicon and carbon in the diamond structure at zero pressure. Since we are interested in the spin-neutral charge gap, we used an equal number of spin up and spin down electrons. We used a Slater-Jastrow trial wave function with backflow corrections [START_REF] Holzmann | Backflow correlations for the electron gas and metallic hydrogen[END_REF][START_REF] Pierleoni | Trial wave functions for high-pressure metallic hydrogen[END_REF]. The Jastrow and backflow functions were fully optimized within variational Monte Carlo including the long-range (reciprocal lattice) contributions. The orbitals in the Slater determinant were taken from DFT-LDA calculations using Quantum Espresso [START_REF] Giannozzi | QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials[END_REF][START_REF] Giannozzi | Advanced capabilities for materials modelling with quantum ESPRESSO[END_REF]. Diffusion Monte Carlo calculations have been performed with the QMCPACK code [START_REF] Kent | Qmcpack: Advances in the development, efficiency, and application of auxiliary field and real-space variational and diffusion quantum monte carlo[END_REF] at the experimentally measured zero pressure valence densities, r s = 1.318 for carbon and r s = 2.005 for silicon. Pseudopotentials were used to remove the core electrons: carbon ions modeled by the Burkatzki-Filippi-Dolg (BFD) pseudopotential [START_REF] Burkatzki | Energy-consistent pseudopotentials for quantum monte carlo calculations[END_REF], and silicon ions by the Trail-Needs

∆ N ∆ N + v M / ∆ N + δ∆ N s ∆ N + v M / + δ∆ N s (a) carbon 1/8 1/3 1/64 1/3 1/216 1/3 1/N 1/3 0.6 0.8 1.0 1.2 1.4 1.6 1.8 ∆ (eV) ∆ N ∆ N + v M / ∆ N + δ∆ N s ∆ N + v M / + δ∆ N s (a) silicon
Figure 2.5: Fundamental gap before and after finite-size corrections. ∆ N is the DMC gap from a simulation with N atoms in the supercell without any finite-size correction, V M / is the leadingorder Madelung correction using the experimental value of -1 , ∆ N s is the next-to-leading-order density correction, which is related to the static part of the structure factor. The line is a fit to ∆ N + δ∆ N s .

(TN) pseudopotential [START_REF] Trail | Smooth relativistic hartreeâfock pseudopotentials for h to ba and lu to hg[END_REF]. For carbon, we used two system sizes: the cubic cell containing 8 atoms and a 2 × 2 × 2 supercell containing 64 atoms. For silicon, in addition to these systems, we used a 3 × 3 × 3 supercell containing 216 atoms. For C and Si, the twist grid density decreases with increasing system size.

The LDA band gaps of carbon and silicon in the diamond structure, are indirect, and lie along the ΓX direction where Γ is the origin of the Brillouin zone and X the Brillouin zone boundary in the (100) direction. By looking directly at the HOMO and LUMO states with LDA it is found that the carbon gap is 3.89 eV for our pseudopotential LDA calculation involves the hole state at (Γ, 13.30 eV) and the particle excitation at (0.75[Γ → X], 17.19 eV), whereas the silicon gap is 0.34 eV between (Γ, 6.36 eV) and (0.85[Γ → X], 6.70 eV), respectively.

The QMC computed gaps for the different sizes of the supercell are summarized in table 2.1. The results from different supercells clearly show the important bias on the gap introduced by the finite size of the supercell. In Figure 2.5, we show the bare gap, ∆ N , the Madelung-corrected one, ∆ N + |v M |/ , and our best correction, ∆ ∞ = ∆ N + |v M |/ + δ∆ s , for both systems against the linear size of the supercell where N is the number of atoms in the supercell. We see that the next-to-leading-order corrections are comparable to the leading-order one, in particular for the 8-atom supercell of Si, whereas they rapidly decay for the larger sizes.

The finite size corrected values, ∆ ∞ , of all different sizes C and Si supercells agree with each other within the statistical uncertainty, yielding the DMC-SJ values ∆ ∞ = 6.8(1) and ∆ ∞ = 1.8(1) for the C and Si gap, respectively. We further note, that these values also Table 2.1: Energy gaps obtained from GCTABC QMC calculations in eV. The bare gap, ∆ N , was calculated from eq. (2.48) for a finite supercell containing N atoms. The leading-order finite-size corrections are given by the screened Madelung constants |v M |/ , the next-to-leading order by the twist correction of two particle density correlations, δ∆ s . We used the experimental value of for C and Si (5.7 and 11.7, respectively) and 18.8 for H 2 extraced from S(k). Finite size corrections were applied also to the band edges, µ ± . The estimate of the gap in the thermodynamic limit is ∆ ∞ = ∆ Ne + |v M |/ + δ∆ s . From our LDA analysis, we estimate a systematic bias of ∼ 0.1 eV from the finite twist grid. This bias is larger than the statistical error. SJ indicates Slater-Jastrow trial wave function, BF backflow. [START_REF] Slater | The theory of complex spectra[END_REF] agree with a numerical N -1/3 extrapolation of the gap values corrected by δ∆ s . For any numerical N -1/3 extrapolation, it is very important to reduce any bias due to higher order corrections as much as possible, since the outcome of a fit is sensitive to the smallest system sizes since they have the smallest statistical uncertainty. For Si, a N -1/3 extrapolation of the bare ∆ N values yields an overestimation of 0.3 eV compared to ∆ ∞ . Since our finite-size corrected gaps show size-convergence for the smallest system size, it is now feasible to address the systematic error due to the fixed node approximation. In order to reduce this bias we have added backflow (BF) correlations in the Slater orbitals. Our backflow correlations lower the SJ gap by 0.1 eV for both, C and Si. Previous BF calculations [START_REF] Hunt | Quantum Monte Carlo calculations of energy gaps from first principles[END_REF] on Si have reported a 0.2 eV lowering compared to SJ. The difference might be due to a different functional form or optimization procedure. A systematic study on the bias of the fixed-node approximation such as done with more general backflow correlations [START_REF] Ruggeri | Nonlinear network description for many-body quantum systems in continuous space[END_REF][START_REF] Holzmann | Orbital-dependent backflow wave functions for real-space quantum monte carlo[END_REF] or multi-determinant trial wave functions [START_REF] Zhao | Variational Excitations in Real Solids: Optical Gaps and Insights into Many-Body Perturbation Theory[END_REF], possible for small supercells, could be done in the future.

r s N ∆ N |v M |/ δ∆ s µ - ∞ µ + ∞ ∆ ∞ H 2 (BF) 1.
So far, in our analysis of C and Si, we have imposed the experimentally known dielectric constant in the leading order Madelung correction. As described in Sec. 2.1.5, there is no need for any external knowledge to perform the size extrapolation as the value of the Madelung correction can be obtained from the behavior of the static structure factor, that can be computed within the same QMC run, see Figs 2.2 and 2.3. However, since the extrapolation involved introduces an additional uncertainty, we have preferred to use the experimental values to benchmark our theory and better distinguish leading from next-to-leading order size effects.

Using the dielectric bound eq. (2.11) on the ground-state structure factor to determine , we get 0 = 6.2 ± 0.4 for C and 0 = 10.3 ± 1.3 for Si, which are compatible with the experimental values of 5.7 and 11.7. The corresponding leading-order finite-size corrections on the gap of the 64-atom system are then 0.92 ± 0.06 eV for C and 0.36 ± 0.14 eV for Si using the ab initio -1 , as opposed to 1.00 eV for C and 0.32 eV for Si based on the experimental values of -1 .

As shown in Fig. 2.2, the asymptotic values of the finite sized structure factors, S ± k , are affected by a much larger uncertainty, introducing larger systematic bias when used for ab-initio size corrections. Still, already the extrapolation to a non-zero value fixes the leading order size corrections to decay as 1/L. This information alone can be crucial as calculations for only two different supercell sizes will be sufficient to determine size effects, whereas more supercell sizes would be needed if the asymptotic form was not known.

Our best values for the fundamental electronic gap (BF-DMC) significantly overestimate the experimentally measured values for C and Si by 1.1 and 0.5 eV, respectively as shown in Table 2.2. There are two main sources of systematic errors that need to be taken into account: the use of pseudo-potentials and the neglect of electron-phonon coupling.

The QMC values for C and Si presented above are based on pseudo-potentials to replace the core electrons of the atoms. Pseudo-potentials are usually designed for accurate prediction of static structural quantities. Excitation spectra, in particular the single particle excitation gap, may be less well described. This has been found in many-body perturbation theory calculations within the GW framework where all-electron calculations have been shown to lower the gap of C and Si by ∼ -0.25 eV [START_REF] Gómez-Abal | Influence of the core-valence interaction and of the pseudopotential approximation on the electron self-energy in semiconductors[END_REF] with respect to pseudopotentials calculations. Although the actual pseudo-potentials of our QMC simulations differ from those used in the GW calculations, we expect that our QMC values will be shifted by a similar amount; we can roughly transfer the all-electron correction of GW to our QMC results.

For lighter atoms, electron-phonon coupling leads to a further reduction of the gap values, even at zero temperature, due to the presence of zero point motion of the ions in the crystal. For C, GW predicts a significant lowering of the gap by -0.6 eV [START_REF] Giustino | Electron-phonon renormalization of the direct band gap of diamond[END_REF], whereas a smaller shift between -60 meV [START_REF] Monserrat | Comparing electron-phonon coupling strength in diamond, silicon, and silicon carbide: First-principles study[END_REF] and -0.1 eV [START_REF] Lautenschlager | Temperature dependence of band gaps in si and ge[END_REF] is expected from DFT for Si.

Considering both the bias due to the pseudo-potential approximation and the neglect of electron-phonon coupling, our BF-DMC for C and Si overestimate the gap by ∼ 0.1-0.2 eV (see table 2.2), larger than our statistical uncertainty. This remaining offset to experiment may either be due to residual bias of the fixed-node approximation, or due to effects in pseudo-potential and e-ph coupling beyond our simple estimations based on GW and DFT.

Table 2.2: Extrapolated band-gap of Si and C from backflow DMC calculations, ∆ BF compared to the experimental values (exp). We tabulated two main corrections: the difference between the gap of an all-electron (AE) and the pseudo-potential (PP) calculation within GW calculations, and the neglect of electron-phonon coupling (e-ph).

∆ BF

AE -PP e-ph exp C 6.6(2) -0.26 (G 0 W 0 ) [106] -0.6 (GW ) [START_REF] Giustino | Electron-phonon renormalization of the direct band gap of diamond[END_REF] 5.48 [START_REF] Madelung | Intrinsic Properties of Group IV Elements and III-V, II-VI and I-VII Compounds[END_REF] Si 1.

Results: hydrogen

Here we will discuss the molecular hydrogen placed in the C2/c-24 structure [START_REF] Chris | Structure of phase III of solid hydrogen[END_REF] at two different densities (r s = 1.38 and r s = 1.34), roughly corresponding to pressures of 234GPa and 285 GPa, respectively. Energies and structure factors were obtained from Reptation Quantum Monte Carlo calculations using the BOPIMC code [START_REF] Pierleoni | The Coupled Electron-Ion Monte Carlo Method[END_REF]. Orbitals were generated using the PBE functional, which has been shown to provide a good trial wave function [START_REF] Pierleoni | Liquid-liquid phase transition in hydrogen by coupled electron-ion Monte Carlo simulations[END_REF][START_REF] Morales | First Principles Methods: A Perspective from Quantum Monte Carlo[END_REF]. We used a supercell with 2 × 2 × 1 primitive cells so that the supercell is nearly cubic and contained 96 protons. the twist convergence has been achieved using a 8 × 8 × 8 twist grid. The crystalline structures have been optimized by variable cell structural relaxation with DFT vdW-DF1 at pressures of 250GPa and 300GPa, respectively. (QMC and DFT-vdW-DF1 pressures differ by ∼10-20GPa [START_REF] Raymond C Clay | Benchmarking exchange-correlation functionals for hydrogen at high pressures using quantum Monte Carlo[END_REF][START_REF] Gorelov | Benchmarking vdW-DF firstprinciples predictions against Coupled Electron-Ion Monte Carlo for high-pressure liquid hydrogen[END_REF]).

We have computed the band gap of solid hydrogen using GCTABC in BF-RQMC The results, in table 2.1 show that the gap and size effects decrease with increasing pressure. For these calculations, we use calculations for one supercell and use its structure . The line is the fit to obtain asymptotic values at k → 0. Black line indicates the inverse dielectric constant extracted from asymptotic behavior of S ± k (Fig. 2.6). Gray shaded area is the error of the fit of S ± k .

factor to estimate the dielectric constant. From Fig. 2.4, we see that HSE DFT slightly underestimates the gap, however, the deviations from the plateau on both sides are quite similar.

For hydrogen, it is difficult to access experimental information on the static dielectric constant. Therefore, we use the asymptotic values of the structure factors S ± k (fig. 2.6) as ab-initio size corrections of the excitation energies. From table 2.1 it can be noticed that the size effects are decreasing with increasing pressure, which can be explained by the fact that the system becomes more "metal-like" with the increased dielectric constant. Figure 2.7 shows the upper bound on the inverse dielectric constant based on eq. 2.11. We see that for hydrogen the above inequality holds for both densities, however, the values of the inverse dielectric constant determined from the asymptotic behavior of S ± k are considerably lower than its upper bound. As the S ± k explicitly contain information on the excitation energies, we used it to correct the gap values. More details on the gap of liquid and solid hydrogen will be given in chapters 5 and 6 respectively.

Neutral excitations

While the previous section was devoted to the electron addition and removal problem, in optical absorption experiments the system gets excited promoting an electron into an empty conduction state. In contrast to (charged) quasiparticle excitations, optical excitations are charge neutral: when the electron is promoted to the conduction state it leaves the hole in the valence band. Since the electron-hole interaction is attractive one may expect a lowering of the excitation energy. The resulting energy gap to the ground state is usually called excitonic or optical gap. Electron-hole interaction effect can also play an important role when calculating the absorption spectrum, because the localization of the electron-hole pair also affects the oscillator strength of the corresponding transitions. Possible bound state are called excitons. Note that usually in solids, the electron-hole pair is delocalized and the exciton is weakly bound, however there are many exceptions [START_REF] Rohlfing | Electron-hole excitations in semiconductors and insulators[END_REF].

Single particle excitations

To zeroth order one might use the KS eigenvalues to describe neutral excitations. Indeed, the unoccupied KS orbitals are one-electron states resulting from exactly the same (local) potential as the occupied orbitals. The excited electron "sees" the same Coulomb field of N -1 electrons as do the ground state electrons. However, excitations from an independent particle picture, like Hartree-Fock, typically overestimate excitation energies, as we have already discussed in section 2.1.1.

QMC excitations

Within the QMC, the trial wave function for neutral excitations is formed by replacing an occupied orbital by an unoccupied in the ground-state determinant of the Slater-Jastrow wave function. However, note that the variational principle only provides an energy which is greater or equal to the lowest eigenstate within the given symmetry class. Therefore one has to construct an excited state trial wave function that will be orthogonal by symmetry to the lower-lying states. Even if the different determinants are orthogonal within the DFT, the Jastrow factor in general causes the wave functions not to be any more strictly orthogonal to the ground state. Running projector Monte Carlo for such states will only guarantee an upper bound of the ground state, not the desired excited one. In practice, one of the ways to achieve orthogonality is by substituting an orbital in the Slater determinant by an unoccupied orbital which has a different crystal momentum k. Note that here we refer to the crystal momentum, which coincides with k-points of the primitive cell, when working in a supercell, to determine the crystal momentum it is necessary to unfold k-point to the Brillouin zone of the primitive cell (see Appendix C).

Assuming that the wave function that models an excited state is appropriately constructed, i.e. is orthogonal to the lower lying excitations, the excitonic/optical gap will, therefore, be the difference of total energies between the system with an electron promoted to the excited state and the ground state,

E ex (k i , k f ) = E + N (k i , k f ) -E N , (2.49)
where E + N (k i , k f ) is the lowest energy of the system with an electron promoted from the initial occupied valence-band orbital at crystal momentum k i to the final unoccupied conduction-band orbital at k f . Since momentum transfer can be negletced in photon absorption, only vertical transitions are possible with k i = k f as long as nuclear motion is not considered. The wave function of an excited state describes a correlated state of an excited electron and a remaining hole, which creates an exciton. The excitonic effect is defined as the difference between the quasiparticle gap, defined in eq. 2.8 and an excitonic gap from eq. 2.49,

∆ ex (k i , k f ) = E QP (k i , k f ) -E ex (k i , k f ) (2.50) = E N +1 (k f ) + E N -1 (k i ) -E + N (k i , k f ) -E N ,
the QMC total energies entering the equation above are statistically independent and can be calculated in parallel.

Excited states expansion in VMC

In order to construct an excited state wave function orthogonal to the lower lying excitations, we will use the generalized variational principle (i.e., the variational principle for the energies of a set of orthogonal trial functions). In QMC this method was first used by Ceperley and Bernu in 1988 [START_REF] Ceperley | The calculation of excited state properties with quantum Monte Carlo[END_REF]. Using the DMC algorithm, authors derived a method for calculating the eigenvalues of several different excited states simultaneously. In the following, we will construct a basis, consisting of M linearly independent excited states Slater-Jastrow wave functions using different orbitals in the Slater determinant as excitations,

Φ i (R) = D i (R)e -U (R) , (2.51)
where R represents electronic coordinate, D i is a Slater determinant with an excitation i and e -U is a Jastrow pair correlation. An approximation to the eigenstates I of the hamiltonian Ĥ can be constructed as a linear combinations of the basis functions,

|Ψ I = M i c I i |Φ i . (2.52)
Upper bound to the exact eigenvalues is then found by minimizing the Rayleigh quotient with respect to c I i [START_REF] Ceperley | The calculation of excited state properties with quantum Monte Carlo[END_REF],

Λ I = M i,j c I * i c I j Φ i | Ĥ|Φ j M i,j c I * i c I j Φ i |Φ j ≥ E I . (2.53)
Minimizing Λ I with respect to c I i , one obtains the many-body generalized eigenvalue equation,

∂Λ I ∂c I i = j c I * j (H ij -Λ I S ij ) = 0, (2.54)
which in matrix notation will become,

Hc I -Λ I Sc I = 0.
(2.55)

There will be M independent solutions (c I , Λ I ) to the equation above, unless S and H are not singular. Due to MacDonald's theorem [START_REF] Macdonald | Successive Approximations by the Rayleigh-Ritz Variation Method[END_REF], the eigenvalues Λ I are upper bounds to the exact excited-state energies of the many-body Schrödinger equation,

Λ I ≥ E I for all 0 ≤ I ≤ M -1. (2.56)
Looking closely into the H matrix, its components can be rewritten in a way that they can be sampled within the VMC run,

H i,j = dRφ * i Ĥφ j |Φ T | 2 (2.57) = dRφ * i φ j E Lj |Φ T | 2 = φ * i φ j E Lj Φ T ≡ H r i,j ,
where φ i = Φ i /Φ T and ... Φ T means the VMC sampling is over |Φ T | 2 . In practice, we choose the Φ T to be the ground state trial wave function. As we work with extended systems, we are sure that the trial wave function does not have zeroes and that there is always an overlap between the ground and the first excited states. Of course, by choosing Ψ T as a combination between the ground and excited states will result in a smaller variance. Index r indicates that we applied the Hamiltonian to the right. One can also apply the Hamiltonian to the left getting,

H l i,j ≡ E * Li φ * i φ j Φ T , (2.58)
taking the conjugate transpose of H l i,j we get,

H l † ≡ E Lj φ j φ * i Φ T (2.59) = φ * i φ j E Lj Φ T ≡ H r .
The property of H l † = H r is satisfied by construction, so the H matrix should be selfadjoint. However, note that due to fluctuations, H will not be symmetric and one should not symmetrize it, because, according to [START_REF] Nightingale | Optimization of ground-and excitedstate wave functions and van der waals clusters[END_REF] that will destroy the zero-variance property.

The estimation of the overlap matrix, S, is straightforward,

S i,j = dRφ * i φ j |Φ T | 2 = φ * i φ j Φ T .
(2.60)

The excitonic gap can be then estimated applying equation 2.49 using the I = 1 and I = 0 solutions of eq. 2.55 as excited and ground energies respectively. So far, the method described above can be easily applied for VMC. There is no particular difficulty in applying DMC to the lowest state of a given symmetry by simply employing a trial wave function of the proper spatial and spin symmetry.

Results: hydrogen

We applied the generalized variational principle in combination with VMC to compute excited state energies for ideal hydrogen crystal. The generalized variational theorem within VMC ensures an upper bound for the excited states [START_REF] Macdonald | Successive Approximations by the Rayleigh-Ritz Variation Method[END_REF]. Note that for high pressure hydrogen we have developed a quite accurate ground state trial function [START_REF] Morales | First Principles Methods: A Perspective from Quantum Monte Carlo[END_REF]. With this fact in mind, we will limit ourselves to VMC calculations only. Furthermore, in the case of excited states, we believe that being largely analytical (including RPA Jastrow and backflow, see section 1.4.3) should make our wave function general enough to treat the excited states as well. The few variational parameters of the trial wave function Ψ T , used to sample excited states, were optimized for the ground state only.

Left panel of figure 2.8 shows quasiparticle and optical gap for ideal C2/c-24 crystalline hydrogen at approximately 285 GPa. On the top panel (fig. 2.8a) the optical gap, computed as described in previous section 2.2.3 and with eq. 2.49 using two basis sizes M = 19 and 25, is presented as a function of only those supercell twists θ, that upon the unfolding (see Appendix C) correspond to the direct excitations of crystal momentum k i = k f . We do not see a significant effect on the gap when the basis size is increased, except that the excitonic gap is slightly larger for a larger basis, therefore, within the error-bars the convergence is reached already at M = 19.

Consider now the different set of supercell twists, θ , that upon unfolding result in indirect excitations, e.g.

k i = k f (fig. 2.8b
). Due to the fact that the orthogonal wave function can be now constructed by simply promoting an electron to the first unoccupied band (equivalent to putting basis size M = 1, see the red line), we can compare the two ways of modeling the excited state. We clearly notice that, by taking a larger basis to model the excitation and orthogonalizing it, correlates the electron-hole more, giving the larger reduction of the optical gap.

On the right panel of figure 2.8 is plotted the excitonic effect energy as defined in eq. 2.49. Large statistical uncertainty does not allow for proper quantitative analyses of the excitonic effect with respect to different basis size, however by increasing the basis the excitonic effect energy decreases. The difference between the direct and indirect excitations is also minor, with the excitonic effect for indirect excitations being lower. Averaged over all twists, the values for the excitonic effect energy are 0.37(0.21) eV for M = 19 and 0.30(0.21) eV for M = 25.

Size effects

Unfortunately, within this work we could not afford to consider systems with different sizes in order to properly evaluate the finite size effect on the neutral excitations. Nevertheless, we will still present general ideas on the origin of the size effects.

In the case of neutral excitations, an effective size of the electron hole interaction plays an important role in determining the size effects. For the simulation cell of the linear size exceeding the effective size of an exciton, the bound state can be formed and the FS error should be dominated by the interaction of excitons with its periodic images and should scale as 1/L 3 , with L being the linear cell size [START_REF] Hunt | Quantum Monte Carlo calculations of energy gaps from first principles[END_REF]. However, when the cell size is smaller than the exciton effective size, the exciton consists of weakly interacting electron and hole and the FS error dominated by the Madelung energies of the free electron and hole and should scale as 1/L. Usually, in solids, an electron-hole interaction is delocalized resulting in a single FS scaling similar to the quasiparticle one (see section 2.1.5) with the error being 1/L.

An estimation of an effective exciton size can be obtained considering a hydrogen-like model of an exciton and a hole interacting with each other. An effective length, L ex , associated to the electron hole interaction can be estimated from the kinetic energy as,

∆ ex = 2 2µL 2 ex , (2.61)
where µ = m e m h /(m e + m h ) is the reduced mass of an electron m e and a hole m h . If assuming the mass of a hole m h = m e and considering the ∆ ex = 0.4 eV of solid C2/c-24 crystalline hydrogen discussed above, then the L ex ≈ 8.3 bohr, which is on the order of the simulation cell size. This indicates that the excitonic effect is indeed delocalized and the FS scaling of the optical gap would potentially be as in quasiparticle excitation of the order of 1/L. However, all this is just a model and extensive studies of size effects of neutral excitation should be performed before making any conclusion.

Conclusion

In this chapter, we have focused on computing the excited states with QMC. The first major part of the discussion is on electron addition and removal energies. We have discussed the fundamental gap, which is the difference between electron addition and removal energy, and have defined its meaning within single electron theories and experiments. We have proposed a new way of computing the fundamental gap using QMC only, which relies on grand canonical twist averaging. Most importantly, we have proposed a scheme for correcting finite simulation cell size effects when computing the fundamental gap. We have shown that for charged systems, finite size supercell calculations are necessarily biased by a finite size error decaying as 1/L, where the prefactor is determined by the absolute value of the Madelung constant and the inverse dielectric constant. We have pointed out that the 1/L functional form is encoded in the long wavelength behavior of the finite size structure factor extrapolating to a nonvanishing value at the origin. We have applied this procedure to determine the fundamental gap of molecular hydrogen at high pressure and carbon and silicon in the diamond structure at zero pressure. Our finite-size corrected gap values for carbon and silicon are larger than the experimental ones. Our results for C and Si demonstrate that the bias due to the finite size supercell can be corrected for, so precise values in the thermodynamic limit can be obtained for small supercells without the need for numerical extrapolation.

In the second part of the chapter, we have discussed the neutral excitations, which are made by promoting an electron from the valence to conduction band and characterized by the coupling between the promoted electron and the hole that was left behind. Within the QMC, we have constructed a wave function to model the excited states by considering a basis of excited states Slater determinants multiplied by Jastrow factor. Applying the generalized variational principle, we can assure that the new excited stated are orthogonal to the lower lying states. Using this procedure, we have computed optical excitations of solid crystalline hydrogen with variational Monte Carlo and have found that we can capture an electron-hole correlation and lower the gap. Chapter 3

Optical properties

We shall discuss in this chapter how the previously introduced electronic excited states are closely related to the optical spectra of materials. By the spectra, we refer to an object's response to a probe, for example, light as a function of probe's energy. Different kinds of responses correspond to different ways of probing the material and measuring the outcome. Here we will limit our focus to a linear response and its formalism within the Kubo theory [START_REF] Kubo | Statistical-mechanical theory of irreversible processes. i. general theory and simple applications to magnetic and conduction problems[END_REF]. With all generality, we will start by defining the response functions from Maxwell's equations. Then, within independent particle theory, we will provide the derivations of the Kubo-Greenwood formula for the conductivity [START_REF] Greenwood | The Boltzmann Equation in the Theory of Electrical Conduction in Metals[END_REF] and will compare the single particle optical spectra of ideal crystalline hydrogen with the results provided in the literature. Finally, we will show the conductivity spectra, computed with many-body Kubo theory within VMC for ideal crystalline hydrogen.

An introduction to the optical properties of solids can be found in the book of Wooten [START_REF] Wooten | Optical properties of Solids[END_REF]. An in-depth theoretical description of the linear response theory is provided in the book of Mahan [START_REF] Mahan | Many-Particle Physics[END_REF]. Moreover, detailed derivations of Kubo and Kubo-Greenwood formalism is provided in the book chapter on electron transport by Allen [START_REF] Allen | Chapter 6 Electron Transport[END_REF].

Dielectric response functions

Propagation of electromagnetic waves in materials is described by Maxwell's equations [START_REF] Wooten | Optical properties of Solids[END_REF][START_REF] Clerk | VIII. A dynamical theory of the electromagnetic field[END_REF]. We shall write these as follows,

∇ × H = 1 c ∂D ∂t + 4π c j, ∇ × E = - 1 c ∂B ∂t , (3.1)
∇ • D = 0, ∇ • B = 0,
where E and D are the electric field strength and displacement, H and B are the magnetic field strength and induction and j is the current density. We have assumed that there are no external sources, therefore the external current and electron density are zero. The electric and magnetic fields E and H are related to their derived fields D and B via the polarisation P and magnetisation M,

D = E + 4πP, H = B -4πM. (3.2)
However, when non-linear effects are neglected we can simplify this relations introducing the coefficients,

D(r, t) = dr dt (r, r , t -t )E(r, t), B(r, t) = dr dt µ(r, r , t -t )H(r, t), j(r, t) = dr dt σ(r, r , t -t )E(r, t), (3.3)
P(r, t) = dr dt χ e (r, r , t -t )E(r, t), M(r, t) = dr dt χ m (r, r , t -t )H(r, t),
where , µ, σ, χ e and χ m are complex dielectric tensor, permeability tensor, conductivity tensor, (electric) susceptibility and (magnetic) permeability. From eqs. 3.2 and 3.3 we have,

= 1 + 4πχ e , (3.4)
µ = 1 + 4πχ m .
In the following, we will consider a non-magnetic media by setting µ = 1 and χ m = 0. Eliminating the magnetic field by substituting the first equation into the second in eqs. 3.2 and using relations in eqs. 3.3, one obtains a wave equation for the electric field E,

∇ 2 E = c 2 ∂ 2 E ∂t 2 + 4πσ c 2 ∂E ∂t . (3.5)
For optical fields, we must look for a sinusoidal wave propagated with dissipation, K, at frequency, ω,

E = E 0 exp{i(K • r -ωt)}, (3.6)
the real part of K can be identified as a wave vector, while the imaginary part accounts for attenuation of the wave inside the material. Then our wave equation requires,

-K 2 = - ω 2 c 2 - 4πiσω c 2 . (3.7)
We can define a complex refractive index N such that,

K = ω c N = ω c (n + ik), (3.8)
where n is the refractive index and k is the extinction coefficient. Rewriting the eq. 3.6 as,

E = E 0 exp - ω c k • r exp i ω c n • r -ωt , (3.9)
we get the equation describing the attenuation of wave amplitude with distance. The fractional decrease in intensity with distance refers to the absorption coefficient, defined as,

α = - 1 I dI dr , (3.10)
where I is the intensity that is proportional to the square of the wave amplitude, from eqs. 3.9 and 3.10 we find,

α = 2ωk/c, (3.11)
here and in the following we will consider that the wave propagation vector is just in one direction (plane waves). Combining eqs. 3.7 and 3.8 we obtain expressions for and σ,

= n 2 -k 2 , (3.12) 4πσ/ω = 2nk, (3.13)
together they can be identified as a complex dielectric function,

= 1 + i 2 = N 2 , (3.14)
All these quantities are, in the most general case, frequency dependent. From eqs. 3.14 we can see that the 1 and 2 are related, together with n and k. These quantities depend on each other in a quite fundamental way by means of the Kramers-Kronig dispersion relations [START_REF] Kramers | La diffusion de la lumière par les atomes[END_REF][START_REF] De | On the theory of dispersion of x-rays[END_REF], which for any complex function χ(ω), which is analytic in the closed upper half-plane of ω give,

Re χ(ω) = 1 π P ∞ -∞ dω Im χ(ω ) ω -ω , (3.15)
Im

χ(ω) = - 1 π P ∞ -∞ dω Re χ(ω ) ω -ω ,
where where P denotes the Cauchy principal value. When concering the optical properties of solids, often the normal incidence reflectivity is involved. In this case we want to construct a solution to Maxwell's equations 3.2 matched to an incident and reflected wave outside. Restricting the wave propagation to x direction we have for the wave outside in the vacuum,

E z = E i e i( ωx c -ωt) + E r e -i( ωx c +ωt) , (3.16)
where the E r < E i . Inside the material we have the transmitted part,

E z = E t e i(kx-ωt) , (3.17)
with the amplitude conservation boundary condition E i + E r = E t . There is as well a magnetic field H y associated with these waves,

-∂ x E z = iω c H y , (3.18)
resulting in the condition E i -E r = NE t . Normal incidence reflectivity is defined as,

R = E r E i 2 = 1 -N 1 + N 2 = (n -1) 2 + k 2 (n + 1) 2 + k 2 < 1. (3.19)
Therefore, by measuring the reflectivity and absorption coefficient, all the optical constant can be determined, which will help to infer many properties of the material. In our interest is to compare measured optical constants to the theoretical predictions. We will only consider here the linear response functions. Theoretically, they can be determined by considering the response of the system to an external perturbation.

Linear response theory

A linear response of an expectation value of an operator Ô to an external field E ext can be written via the coefficient χ, Ô(r, t) ext = Ô(r, t) + dr dt χ(r, r , tt )E ext (r , t ), (3.20) which resembles to the macroscopic response coefficients in eqs. 3.3. It was shown by Kubo in 1957 [4] that the linear response to a time dependent perturbation Ĥext (t) writes as

δ Ô(t) = i t t 0 dt Ψ 0 | Ĥext H (t ), ÔH (t) |Ψ 0 , (3.21)
where subscript H denotes the Heisenberg picture and |Ψ 0 is the many-body ground state. Now, consider the external perturbation to a scalar field Ĥext (t) = drn(r, t)V ext (r, t), where n(r, t) is the electron density operator. Such perturbation would induce a charge density n ind in the system and the response can be written according to eq. 3.21 as

n ind (r, t) = i t t 0 dt dr V ext (r , t ) Ψ 0 | [n(r , t ), n(r, t)] |Ψ 0 . (3.22)
The response function χ now can be written, by combining eq. 3.20 and eq. 3.22,

χ(r, r , t -t ) = -iΘ(t -t ) Ψ 0 | [n(r, t), n(r , t )] |Ψ 0 , (3.23)
where Θ(tt ) has been introduced to meet the causality conditions to be able to write the response in the form of eq. 3.20. Elaborating further by opening the commutator, inserting the completeness relation for the N-particle basis states i |Ψ i Ψ i | = 1 and putting τ = tt we get

χ(r, r , τ ) = -iΘ(τ ) i Ψ 0 |n(r, τ )|Ψ i Ψ i |n(r , 0)|Ψ 0 + c.c. (3.24)
Now, recalling from eq. 3.22 that n(r, τ ) = e i Ĥ0 τ n(r, 0)e -i Ĥ0 τ and denoting the matrix element ρ i (r) = Ψ 0 |n(r, 0)|Ψ i , eq. 3.24 transforms to

χ(r, r , τ ) = -iΘ(τ ) i ρ i (r)ρ * i (r )e i(E 0 -E i )τ -c.c , (3.25)
where E i is the energy of the i th state. The Fourier transform (FT) to a frequency domain brings us to the final expression for χ,

χ(r, r , ω) = i ρ i (r)ρ * i (r ) ω -(E 0 -E i ) + iη + ρ i (r )ρ * i (r) ω + (E 0 -E i ) + iη , (3.26)
where infinitesimal η was introduced to have the time integral in FT converge. In the end one should consider η → 0.

Independent particle polarisability

Considering the independent particle electronic system and writing χ in single particle orbitals, φ n (r), we get

χ 0 (r, r , ω) = nn 2(f n -f n )φ * n (r)φ n (r)φ * n (r )φ n (r ) ω -( n -n ) + iη , (3.27)
where χ 0 states for the independent particle polarizability and f n is the occupation number.

All the quantum numbers are contained in the label n. Independent particle irreducible polarizability, χ 0 , is often called the random phase approximation (RPA), the name that is related historically to the RPA used in section 1.4.4. When Fourier transformed, going from r to q space, and considering the single particles orbitals as Bloch functions with n = (i, k), with k being in the first Brillouin zone, the independent particle polarisability obtains the usual form

χ 0 G,G (q, ω) = 2 Ω k,k ,i,j (f k,i -f k ,j ) φ k,i |e i(q+G)r |φ k ,j φ k ,j |e -i(q+G )r |φ k,i ω -( k,i -k ,j ) + iη , (3.28)
where Ω is the cell volume, which comes from performing the Fourier transform and

φ k,i |e i(q+G)r |φ k ,j ≡ Ω dru * k,i e -ikr i e i(q+G)r u k ,j e -ik r j (3.29)
is calculated as an integral over the cell with volume Ω. Vector q is lying within the first Brillouin zone and vectors G and G are the reciprocal lattice vectors. It can be shown that in order to obtain a macroscopic polarisability it is necessary to take G = G = 0, considering only the first Brillouin zone. Therefore, in the following we omit indexes G and G , considering only zero's components. Note further that due to the momentum conservation k = k + q, the sum over k can be omitted,

χ 0 (q, ω) = 2 Ω k,i,j (f k,i -f k+q,j )| φ k,i |e iqr |φ k+q,j | 2 ω -( k,i -k+q,j ) + iη , (3.30)
the dielectric function in Fourier space, (q, ω) = 1 -4πe 2 q 2 χ(q, ω), becomes

0 (q, ω) = 1 - 8πe 2 Ωq 2 k,i,j (f k,i -f k+q,j )| φ k,i |e iqr |φ k+q,j | 2 ω -( k,i -k+q,j ) + iη , (3.31)
in the following we will consider only the long-wavelength dielectric function, i.e. q → 0 limit. The exponent can by then expanded, e iqr 1 + iqr, and the dielectric function for insulators writes as a tensor,

0;α,β (ω) = 1 - 8πe 2 Ω k,i,j (f k,i -f k,j ) φ k,i |r α |φ k,j φ k,j |r β |φ k,i ω -( k,i -k,j ) + iη , (3.32)
For small η the imaginary part of denominator in eq. 3.32 behaves like the Dirac delta function,

lim η→0 Im 1 ω -( k,i -k,j ) + iη = lim η→0 -η [ω -( k,i -k,j )] 2 -η 2 (3.33) = -πδ(ω -( k,i -k,j )).
The imaginary part of the dielectric function is then,

0;α,β 2 (ω) = 8π 2 e 2 Ω k,i,j (f k,i -f k,j ) φ k,i |r α |φ k,j φ k,j |r β |φ k,i (3.34)
×δ(ω -( k,ik,j )).

Kubo-Greenwood electrical conductivity

For systems with periodic boundary conditions it is useful to transform the matrix elements in eq. 3.35 from r α to ∇ α . Using the commutator relation (in SI units, for comparison to the literature and for local pseudopotentials),

i [ Ĥ, r α ] = - i ∇ α m , (3.35)
the matrix elements of r α can be expressed in terms of ∇ α as follows,

φ k,i |r α |φ k,j = 2 m φ k,i |∇ α |φ k,j k,i -k,j . (3.36)
The real part of the independent particle conductivity σ 1 (ω), according to eq. 3.14, can be expressed as,

σ 1 (ω) = ω 2 4π . (3.37)
Explicitly including the directionality dependence, the conductivity then fully described by the Cartesian tensor,

σ α,β 1 (ω) = 2πe 2 2 m 2 Ωω k,i,j (f k,i -f k,j ) φ k,i |∇ α |φ k,j φ k,j |∇ β |φ k,i δ(ω -( k,i -k,j )) (3.38)
the imaginary part is easily accessible via the Kramers-Kronig relations defined in eq. 3.15.

Optical spectra of hydrogen

In the following section, we will discuss the optical spectra of ideal crystalline hydrogen computed with the Kubo-Greenwood (KG) formula. We considered C2/c-24 and Cmca-12 crystalline symmetry at approximately 300 GPa. Calculations of the eigenstates and eigenvalues were done using HSE [START_REF] Heyd | Energy band gaps and lattice parameters evaluated with the Heyd-Scuseria-Ernzerhof screened hybrid functional[END_REF] and vdW-DF2 [START_REF] Lee | Higher-accuracy van der Waals density functional[END_REF] XC approximation using Quantum Espresso [START_REF] Giannozzi | QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials[END_REF][START_REF] Giannozzi | Advanced capabilities for materials modelling with quantum ESPRESSO[END_REF], Kubo-Greenwood conductivity was then evaluated with the KGEC code [START_REF] Calderín | Kubo-Greenwood electrical conductivity formulation and implementation for projector augmented wave datasets[END_REF]. Note that, as in section 2.1.8, the crystalline structures have been optimized by variable cell structural relaxation with DFT vdW-DF1 at pressures of 300 GPa. From the real and imaginary parts of the conductivity, all the other optical properties can be inferred using relations in eqs. 3.14 and 3.19.

First, we discuss the convergence of the spectra with respect to the simulation parameters. Figure 3.1 shows the real and imaginary part of the Kubo-Greenwood conductivity of C2/c-24 hydrogen crystal at 300 GPa, computed for a 2 × 2 × 1 supercell of 96 nuclei. The convergence is tested for three parameters: number of k-points, the total number of bands included in DFT calculations, and smearing used to represent the delta function in the equation for the KG conductivity. For all the KG calculations within this thesis, we used Gaussians to represent the delta function. Changing the number of bands results in the largest changes in the spectrum, middle panel of fig. 3.1. In order to reach the full convergence, we will have to include much more than 70 bands in our DFT calculation with 96 electrons. However, sometimes one is just interested in the beginning of the spectrum (up to ∼8 eV, sometimes called the edge or profile ), for that region, the convergence can be achieved already with 56 bands. On the top panel of fig. 3.1 conductivity with the different number of k-points is presented, the convergence is reached already with 8 × 8 × 8 grid. Finally, the test of the smearing factor in the lower panel of fig. 3.1 done at 8 × 8 × 8 k-point grid shows that the spectra are best represented at 0.2 eV. The top and the bottom panels are connected, indeed, the dependency is that by increasing the size of k-grid, the smearing needed for a smooth representation of the spectra should decrease. Next, we compare the KG spectra to the available theoretical predictions. For that purpose, we consider Cmca-12 hydrogen crystal at 300 GPa. The only available theoretical calculations of the dielectric function were done by Lebegue et al. [START_REF] Lebegue | Semimetallic dense hydrogen above 260 GPa[END_REF] and Dvorak et al. [START_REF] Dvorak | Quasiparticle energies and excitonic effects in dense solid hydrogen near metallization[END_REF]. We will focus here on Lebegue's work, as they provide the dielectric function at pressure considered in this thesis. In their work, the authors presented single electron RPA calculations of the dielectric tensor with PBE XC approximation. Figure 3.2a shows the imaginary part of the dielectric function for the three diagonal components of the tensor. Note that the crystalline hydrogen is predicted to be a planar material and, therefore, the ZZ component is different. It is important to note also that in their work, the authors optimized their structure with PBE functional while we use vdW-DF1. When the structure is optimized with PBE the calculated dielectric function using PBE and Kubo-Greenwood is in perfect agreement to Lebegue et al. However, geometry optimized with vdW-DF, results in different gap values and shifts the dielectric function to higher energies (see fig. 3.2b). This fact explains the difference between GW and QMC gaps, which will be presented further in chapter 5.

Kubo electrical conductivity

In the same spirit of the linear response, one can consider the response of a current to an electric field. The resulting response function is conductivity, σ α,β (q, ω), which is again given by the Kubo formula [START_REF] Kubo | Statistical-mechanical theory of irreversible processes. i. general theory and simple applications to magnetic and conduction problems[END_REF][START_REF] Mahan | Many-Particle Physics[END_REF] (in SI units, for future comparison to the literature),

σ α,β (q, ω) = 1 Ω ω ∞ 0 dte iωt j † α (q, t), j β (q, 0) (3.39) the angular brackets Ô ≡ 1 Z I e -βE I I| Ô|I (3.40)
mean a canonical ensemble average using the equilibrium density matrix, Z = I e -βE I is the partition function and |I is a complete set of many-body eigenstates of the system before the external field is applied, Ĥ|I = E I |I . Considering Heisenberg time-dependence for the current , j α (q, t) = exp(i Ĥt/ )j α (q, 0) exp(-i Ĥt/ ), (3.41) with the time independent Fourier transform of the current operator being,

j α (q) = - 1 2m n e i p n,α e iq • rn + e iq • rn p n,α . (3.42)
The time integral can be solved by assigning to a frequency a small positive imaginary part, ω → ω + iη. In the further we will as well work in q → 0 limit in the current operator. If we know the complete set of many-body eigenstates |I , then eq. 3.39 reads,

σ α,β (ω) = i Ωω IJ e -βE I -e -βE J Z I|j † α |J J|j β |I ω -(E J -E I ) + iη . (3.43)
When considering just the real part of the conductivity, the denominator can be transformed to a delta function as in eq. 3.34

σ α,β 1 (ω) = π Ωω IJ e -βE I -e -βE J Z I|j † α |J J|j β |I δ( ω -(E J -E I )). (3.44)
Consider now low temperature limit, so only the ground state will contribute to the density matrix. Regrouping the temperature factors, e -βE I (1e -β(E J -E I ) ), we get,

σ α,β 1 (ω) = π 1 -e -β ω ΩωZ IJ e -βE I I|j † α |J J|j β |I δ( ω -(E J -E I )) (3.45) → β→∞ π Ωω J 0|j † α |J J|j β |0 δ( ω -(E J -E 0 )).
Rewriting the current operator in the q → 0 limit,

j α = - e m n p n,α = ie m n ∇ n,α , (3.46)
where p n,α is the momentum of nth electron, e is electronic charge and m electronic mass. The final expression for the real part of the conductivity tensor at zero temperature writes,

σ α,β 1 (ω) = πe 2 2 Ωωm 2 J 0| n ∇ n,α |J J| n ∇ n,β |0 δ( ω -(E J -E 0 )). (3.47)
In this final form, the Kubo formula for conductivity can be used within the VMC to obtain the optical properties of the systems at ground state and zero temperature.

Momentum matrix elements within VMC

The momentum matrix elements in the basis introduced in eq. 2.52,

-iP α ij = Φ i | n ∇ n,α |Φ j , (3.48)
where n indicates the sum over the electrons, can be obtained within VMC in the same fashion as described in section 2.2.3. Indeed, the quantity

-iP α ij = dRφ * i φ j n ∇ n,α log(Φ j )|Φ T | 2 = φ * i φ j n ∇ n,α log(Φ j ) Φ T , (3.49)
with φ i = Φ i /Φ T , can be accumulated during the VMC run. However, in eq. 3.47 the momentum matrix element is between the ground state Ψ 0 and excited state Ψ J . In terms of the basis functions, Φ i , it can be written as,

-iP α 0J = i,j,n c 0 * i φ i |∇ n,α |φ j c J j = i,j c 0 * i P α ij c J j = c 0 † P α c J , (3.50)
where vectors c J are the solutions of the many-body generalized eigenvalue equation 2.55 normalized by M i,j c J * i c J j Φ i |Φ j . With this, the real part of the Kubo electrical conductivity tensor can be obtained within VMC as

σ α,β 1 (ω) = πe 2 2 Ωωm 2 J P α 0J P β * J0 δ( ω -(E j -E 0 )). ( 3

Results: solid hydrogen

.3 illustrates an attempt to use variational Monte Carlo to compute the real part of electrical conductivity using the Kubo formula for the ground state, as presented in eq. 3.51. We use again C2/c-24 ideal crystalline hydrogen at 285 GPa. For comparison, we provide Kubo-Greenwood spectra computed using vdW-DF2 XC functional. It is clear that the number of excitations considered (M = 19 and 25) is not sufficient and additional study is required. Further, we discuss the potential steps that we think are important to consider in order to improve the accuracy of computed spectra:

• we need to study the convergence with respect to the number of the basis functions Φ i used to represent excited states • it is important to take into account the fact that, unless the trial wave function is exact, using VMC always gives a statistical uncertainty on the data, which might make the diagonalization procedure unstable and result in non-physical features in the spectra

• the statistical noise on the matrix elements might give a positive contribution even for the forbidden transitions Trying to investigate these points should be the next step of this section. Unfortunately, future investigations could not be afforded within this work.

Conclusion

We have introduced in this chapter the theoretical background for computing the optical properties of materials. In particular, we have presented the introduction into the linear response theory, discussing different optical response coefficients such as dielectric function, conductivity, and reflectivity and have established the relations between them. We have made the derivations necessary to establish the single particle Kubo-Greenwood formula for conductivity, and have applied it to the ideal solid crystalline hydrogen.

Moreover, we have attempted to compute the Kubo many-body conductivity with variational Monte Carlo. Our preliminary results for the C2c-24 ideal solid crystalline hydrogen indicate that further investigations of the convergence of the spectra are necessary.

Chapter 4 Thermal crystals: renormalization of electronic properties

In the previous chapters, we have assumed that the nuclear lattice is immobile and the protons were point particles, fixed at their equilibrium positions. This assumption can only be valid for some materials with the nuclear mass large enough such that the energy contribution due to the lattice vibration is small compared to the electronic excitation energies (such e.g. eq. 2.49) computed for a static lattice,

E static ex << E vib . (4.1)
Moreover, the above assumption makes sense only if the property of interest is the electronic energy. However, when dealing with light materials such as hydrogen, the vibrational energy is no longer small and has to be taken into account when computing the electronic structure properties. The focus of this chapter, therefore, will be on the renormalization of energy gaps, excited states, and optical spectra of materials at finite temperature. We first define the full electron nuclear hamiltonian and will explain the adiabatic splitting and path integral formulation of the problem. Then, we introduce the formal treatment of electron addition energies and gaps in the canonical and semi-grand canonical ensemble at finite temperature and will discuss the approximations to be made. Further, we proceed with an issue of determining crystal momentum in the presence of finite temperature nuclei. At the end of the chapter, we provide a new way of treating the optical spectra of finite temperature materials, an alternative to the common semiclassical method.

A general introduction to the nuclear dynamics and the adiabatic approximation can be found in the book of Marx [START_REF] Marx | Ab Initio Molecular Dynamics[END_REF]. An introduction to the finite temperature methods used for finite temperature nuclei to study matter at extreme conditions can be found in the review of MaMahon [START_REF] Mcmahon | The properties of hydrogen and helium under extreme conditions[END_REF]. There is as well an instructive review of path integral Monte Carlo [START_REF] Ceperley | Path integrals in the theory of condensed helium[END_REF], together with the lecture notes by Pierleoni and Ceperley [START_REF] Pierleoni | The Coupled Electron-Ion Monte Carlo Method[END_REF].

To properly treat electrons and nuclei as a quantum mechanical system, one has to consider the full Hamiltonian (non relativistic, in atomic units):

Ĥ = T R - Np I=1 ∇ 2 R I 2M I + V R pp Np I Np J>I Z I Z J |R I -R J | + Te + Vee + V R ep ĤR e , (4.2)
where for convenience, ĤR e , is the electronic Hamiltonian, defined in eq. 1.2, plus the nuclear-nuclear interaction potential:

ĤR e = V R pp + Te + Vee + V R ep , (4.3)
which now depends parametrically on the nuclear coordinates. This Hamiltonian can be viewed as coupled nuclear and electronic problems. In fact, if the nuclei were fixed in space, the full Hamiltonian in eq. 4.2 would become just the electronic one, ĤR e . The eigenstates, Φ n (r, R), of the full Hamiltonian will depend on the positions of the electrons and of the nuclei. In the following, unless the index is specifically given, r and R denote the set of all electronic and nuclear positions, r = {r 1 , .., r i , .., r Ne } (4.4)

R = R 1 , .., R I , .., R Np ,
being a 3N e and 3N p dimensional vector respectively. Solving the Schrödinger equation analytically with the Hamiltonian defined in eq. 4.2 for all particles together is essentially an impossible task for most systems except a few extremely simple cases such as the H + 2 . Therefore, when dealing with nuclear motion, one has to decouple electrons and nuclei. In the following, we will consider a non-degenerate case, when the eigenvalues of the electronic hamiltonian are well separated from each other.

Born-Oppenheimer approximation

The mass of the electrons is much smaller than that of the nuclei, therefore the ratio of kinetic to potential energy for nuclei is much smaller than in the case of electrons. As long as the temperature is not too high, the time scale of electronic motion is much smaller than the corresponding time scale for nuclei. As a consequence, one can assume that electrons relax infinitely fast to their equilibrium state as the ions move, neglecting any coupling between different electronic states and also retardation effects in the electron-ion interaction. In these conditions the electrons do not exchange energy with the nuclei; therefore their evolution is adiabatic. If further considering that the electronic ground state contribution will be dominant, then this approximation goes under the name of the Born-Oppenheimer approximation (BOA) [START_REF] Bransde | Physics of atoms and molecules[END_REF][START_REF] Born | Zur quantentheorie der molekeln[END_REF] To be more concrete let us expand the total electron nuclei wave function in the eigenstates of the electronic Hamiltonian,

Φ n (r, R) = α Ψ α (r|R)χ αn (R), (4.5)
where the electronic wave function Ψ α (r|R) depends parametrically on the nuclear coordinates, R. Plugging it into the Schrödinger equation with the full Hamiltonian defined in eq. 4.2,

α T R Ψ α (r|R)χ αn (R) + α ĤR e Ψ α (r|R)χ αn (R) = α E αn Ψ α (r|R)χ αn (R), (4.6)
one can see that the electronic Hamiltonian operates only on the electronic wave function, Ψ α (r|R), and not on the coefficients χ αn (R),

α ĤR e Ψ α (r|R)χ αn (R) = α E R α Ψ α (r|R)χ αn (R). (4.7)
Whereas the nuclear kinetic operator acts on both wave functions,

∇ 2 I 2M I Ψ α (r|R)χ αn (R) = - 1 2M I Ψ α (r|R)(∇ 2 R I χ αn (R)) (4.8) - 1 M I (∇ R I Ψ α (r|R))(∇ R I χ αn (R)) - 1 2M I (∇ 2 R I Ψ α (r|R))χ αn (R).
Plugging eq. 4.8 and eq. 4.7 into the eq. 4.6, one can integrate out the electronic degrees of freedom by multiplying Ψ * β (r|R)) on the left,

T R + E R α -E αn χ αn (R) = Np I=1 1 M I α Ψ β |∇ R I |Ψ α (∇ R I χ αn (R)) + 1 2M I α Ψ β |∇ 2 R I |Ψ α χ αn (R) . (4.9)
The adiabatic approximation relies on neglecting the r.h.s of the previous equation, which would lead to an eigenvalue problem for nuclei with potential given by the electronic eigenvalue, E R β , evaluated at fixed nuclei, R. This approximation assumes that the electrons react instantaneously on the change of the nuclear positions. To understand its validity lets use the following notation for the terms on the r.h.s of eq. 4.9:

A I αβ = 1 M I Ψ β |∇ R I |Ψ α ∇ R I , B I αβ = 1 2M I Ψ β .|∇ 2 R I |Ψ α (4.10)
The off-diagonal terms A I αβ and B I αβ are usually regarded as non-adiabatic corrections and are neglected such that the electrons stay in a given state as the ions move. Note that, for normalized wave functions the diagonal part of A I αα is identically zero, which can be shown by integration by parts. Neglecting further the B I αα term, the Schrödinger equation for the nuclei writes:

T R + E R α χ αn (R) = E αn χ αn (R), (4.11)
where the potential energy surface (PES) felt by the nuclei is defined as the eigenvalue of electronic Hamiltonian in the clamped nuclei approximation, E R β , which includes as well the nuclei-nuclei potential energy contribution, V R pp . As in this work, the typical temperatures of the interest are T T F , where T F is the Fermi energy for the electrons, it can be assumed that electrons remain in their ground state, arriving at the modern version of the BOA,

T R + E R BO χ αn (R) = E 0n χ 0n (R), (4.12)
where we define the BO potential energy surface as the ground state energy of the electronic Hamiltonian, ĤR e .

Ab initio path integrals

Alternatively, one can derive the BOA using the path integral (PI) formulation [START_REF] Richard P Feynman | Statistical Mechanics A Set of Lectures[END_REF]. In analogy to [START_REF] Marx | Ab Initio Molecular Dynamics[END_REF], we will work in a mixed basis combining the position representation for the nuclear degrees of freedom with the energy representation for the electrons. In particular, the product basis, |Ψ α (R), R = |Ψ α (R) |R , involving the complete and orthonormal basis set of the eigenstates of the electronic Hamiltonian. The completeness relation,

α |Ψ α (R), R R, Ψ α (R)|dR = 1, (4.13)
if necessary, has to be extended over the full available space in the integration over the nuclear positions, and the summation over the electronic states α must also include all continuum states. We will be interested in computing the expectation value of an operator Ô, that can be then expressed according to Statistical mechanics,

Ô = 1 Z Tr Ô exp -β Ĥ . (4.14)
The exact partition function, Z, of the quantum-statistical canonical ensemble reads, can be invoked in order to decouple nuclear kinetic energy from the rest of the Hamiltonian. Inserting P -1 times completeness relation, this will allow us to rewrite the partition function as the product, Z = lim

Z = Tr exp -β Ĥ = dR α R, Ψ α (R)|e -β[ T R +H R e ] |Ψ α (R), R (4.15) = dR α ρ(R, Ψ α ; R, Ψ α ; β),
P →∞ P s=1 α (s) R (s+1) , Ψ α (s+1) (R (s+1) )|e -β P T R e -β P H R e |Ψ α (s) (R (s) ), R (s) dR (s) ,
(4.17) where the trace condition imposes periodic boundary conditions, R (P +1) = R 1 and Ψ α (P +1) = Ψ α 1 , on the Trotter discretization parameters s = 1, .., P . It remains to evaluate the following matrix element:

R (s+1) | Ψ α (s+1) (R (s+1) )|e -β P T R e -β P H R e |Ψ α (s) (R (s) ) |R (s) = (4.18) R (s+1) | Ψ α (s+1) (R (s+1) )|e -β P T R |Ψ α (s) (R (s) ) |R (s) e -β P E α (s) (R (s) ) ,
where the effect of the electronic Hamiltonian is expressed by the adiabatic energy, E α (s) (R (s) ), in electronic state Ψ α (s) evaluated at the fixed nuclear configuration R (s) at time slice s and can be pulled out of the integral. At this stage the adiabatic approximation can be applied, which amounts to assuming that the nuclear kinetic energy operator, T R , does not act on the electronic wave functions, Ψ α (s) . This means that the variation of the electronic wave functions as the result of changes in the nuclear coordinates is negligible. In the same spirit, the overlap between the two electronic wave functions at time slice s and s + 1 will yield the orthogonality relation:

Ψ α (s+1) (R (s+1) )|Ψ α (s) (R (s) ) = Ψ α (s+1) (R (s+1) )| |Ψ α (s) (R (s+1) ) (4.19) + I ∇ I |Ψ α (s) (R (s+1) ) (R (s+1) -R (s) ) +O (R (s+1) -R (s) ) 2 ≈ Ψ α (s+1) (R (s+1) )|Ψ α (s) (R (s+1) ) = δ α (s+1) α (s) .
The final expression of the partition function in the adiabatic approximation writes:

Z adia = lim P →∞ P s=1 dR (s) α (s) δ α (s+1) δ α (s) R (s+1) |e -β P T R |R (s) e -β P [E α (s) (R (s) )] = α lim P →∞ P s=1 dR (s) R (s+1) |e -β P T R |R (s) e -β P [Eα(R (s) )] , (4.20)
where as a result of the product of Kronecker deltas, each of them connecting two adiabatic states Ψ (s)

α and Ψ (s+1) α under periodic boundary conditions in imaginary time s = P + 1 = 1, only one sum, α (1) = α , survives. More compactly the eq. 4.20 can be written as the trace,

Z adia = α Tr{e -β( T R +Eα(R)) }. (4.21)
The remaining sum over all electronic eigenstates, Ψ α , means that the nuclear fluctuations are sampled according to Boltzmann's thermal distribution with the weight given by the electronic eigenvalue ∝ exp [-βE α ]. Thus, the contribution of excited states (α > 1) to the partition function is suppressed exponentially as a function of the energy with respect to the ground state, ∝ exp [β (E α -E 0 )], separately for each nuclear configuration R. Therefore, if the electronic energy gap between the excited states and the ground state is large, i.e. E α -E 0 1 β , then the electronic ground state contribution to the partition function will be dominant and the Born-Oppenheimer approximation can be invoked. The sum over α collapses to only a ground electronic state,

Z BO = Tr{e -β( T R +E 0 (R)) }, (4.22)
bringing us to the Born-Oppenheimer approximation of the partition function.

To finalize the formal definition of the path integrals approach, it is worth noticing that the propagator of the nuclear kinetic energy operator in eq. 4.20 can be evaluated analytically in momentum space. By the virtue of momentum's completeness relation and noting that the plane wave can be expressed as R|K = e -iKR . The nuclear kinetic energy operator,

R|e -β P T R |R = dKdK R|K K |e -β P T R |K K|R (4.23) = dKdK e -iK R e -β 2M I P K 2 δ(K -K )e -iKR = M I P 2πβ 3N I 2 e -M I P 2β (R-R ) 2 ,
is represented via the harmonic interactions that couple neighboring imaginary time slices. We will consider, for convenience, a single component system with nuclei of mass M I In the latter we used the fact that the kinetic energy operator in momentum representation has the following form:

K| T R |K = -K 2 2M I δ(K -K ).
The complete Trotter product assembles to,

Z adia = M I P 2πβ 3N I 2 α lim P →∞ P s=1 dR (s) (4.24) × exp - P s=1 M I P 2β (R (s) -R (s+1) ) 2 - β P E α (R (s) ) ,
which yelds the partition function expressed as a path integral in the adiabatic approximation. The result of Trotter discretization is the path of P time slices s = 1, ..., P of "duration" ∆τ = β/P ,

R (s) = R (1) 1 , ..., R (1)
N I ; ...; R The two nuclei {R

(1)

1 , R (2)
1 , ..., R

(P ) 1 } and {R (1) 2 , R (2)
2 , ..., R

(P ) 2 } interact at each time slice (s) via E 0 (R (s) 1 , R (s)
2 )/P obtained from the ground state electronic structure calculations. The interaction within each ring polymer I = 1, 2 is given by harmonic springs

M I P 2β (R (s) I -R (s+1) I ) 2
between nearest neighbors only and is subject to periodic boundary conditions, R

(P +1) I = R (1)
I .

where we have explicitly wrote single nuclei coordinates R (s)

I and has to be closed due to the trace condition R

(P +1) I ≡ R (1)
I , i.e. they are periodic in imaginary time τ . The effective classical partition function eq. 4.24 can be interpreted as the N I polymers each composed of P monomers. Each quantum degree of freedom is represented by a ring polymer with intrapolymer interactions coming from the kinetic energy and consisting of harmonic nearest-neighbor couplings along the closed chain. The interpolymer interaction is given by the scaled potential E (s) α /P , which is the eigenvalue of the electronic hamiltonian ĤR e evaluated at fixed nuclear position R (s) at the corresponding time slice (s). A sketch of such a polymer for two particles is represented on fig. 4.1. Note that this picture is for distinguishable particles (Boltzmannon). If other statistic is important, the picture is different and more complex [START_REF] Ceperley | Path integrals in the theory of condensed helium[END_REF]. It should be reminded that the eq. 4.24 is an exact reformulation of adiabatic partition function, eq. 4.21, in the infinitely fine discritization P → ∞. In order to evaluate the desired expectation values based on eq. 4.24, one can utilize Metropolis Monte Carlo sampling, which was firstly done by Pollock and Ceperley in 1984 [START_REF] Pollock | Simulation of quantum many-body systems by path-integral methods[END_REF] to perform accurate simulations of liquid and solid helium.

In this work a method fully based on quantum Monte Carlo is used to perform PI simulations. The coupled electron-ion Monte Carlo method (CEIMC) is based on PIMC, where the electronic properties are computed by quantum Monte Carlo [START_REF] Pierleoni | Computational methods in coupled electronion Monte Carlo simulations[END_REF]. Details of CEIMC can be found in Appendix D.

Quasiparticle energy gap in a canonical ensemble at finite temperature

In this section, we will discuss the thermal renormalization of the quasiparticle energy gap, defined in section 2.1.3. We will consider the system of N p nuclei and N e = N p + n electrons. The canonical partition function in the Born-Oppenheimer approximation,

Z BO (N e , N p ) = Tr{e -β( T R (Np)+E 0 (R;Ne,Np)) } (4.26) = e -βF (Ne,Np)
now is written in terms of the Helmholtz free energy and the dependency on the number of nuclei and electrons has been indicated explicitly. For convenience, the free energy, F , corresponding to the PES E 0 (R) has been introduced. In the further, we will remove the subscript BO from the partition function and will explicitly state when it is otherwise. We will try now to reformulate the partition function, Z(N e , N p ), of N e = N p + n electrons in terms of N e = N p electrons times some correction. Multiplying and dividing eq. 4.26 by the e -βF (Np,Np) , Z(N e , N p ) = e -βF (Np,Np) e -βF (Np,Np) e -βF (Ne,Np) (

= e -βF (Np,Np) Tr{e -β( T R (Np)+E 0 (R;Np,Np)+∆E 0 (R;n,Np)) } e -βF (Np,Np) = e -βF (Np,Np) e -β∆E 0 (R;n,Np) Np , one can argue that the correction, which is defined via the average over the system with the N e = N p , is given by the increment of the PES, when n electrons are added to the system, ∆E 0 (R; n, N p ) = E 0 (R; N e , N p ) -E 0 (R; N p , N p ). Notice that as long as n << N p , ∆E 0 (R; n, N p ) E 0 (R; N p , N p ), because the energy is extensive property. We assume that the ∆E 0 (R; n, N p ) is normally distributed over the nuclear fluctuations (the verification of the distribution is given on fig. 4.2). The average of the exponent can be expanded in cumulants, for normally distributed variable being, log e -β∆E 0 (R;n,Np)

Np = -β ∆E 0 (R; n, N p ) Np + σ 2 β 2 2 + O(β 3
), (4.28)

σ 2 = ∆E 0 (R; n, N p ) 2 Np -∆E 0 (R; n, N p ) 2

Np

where σ is the variance of normally distributed ∆E α (R; n, N p ). the square root of the variance σ of the electron addition energy, normally distributed over the nuclear configurations, is in fact 10 times smaller than the average. The provided result is for crystalline hydrogen at 200 K. For other materials the variance would decrease as they are heavier elements, which makes this assumption more valid. The average over the exponent can be then translated to the argument, giving,

Z(N e , N p ) e -β(F (Np,Np)+ ∆E 0 (R;n,Np) Np ) . (4.29)
This result implies that the free energy of the system of N e electrons and N p protons,

F (N e , N p ) F (N p , N p ) + ∆E 0 (R; n, N p ) Np , (4.30)
can be to a good approximation represented as the free energy of the system with N p electrons and N p protons corrected by the average of the electronic energy that is required to add n electrons. The electron addition/removal energy in the canonical ensemble can be then computed as the Free energy difference,

∆ Ne±1 c = F (N e ± 1, N p ) -F (N p , N p ) ≈ E 0 (R; N e ± 1, N p ) -E 0 (R; N p , N p ) Np , (4.31)
where the BOA was invoked by considering only the ground state PES and the average, .. Np , is taken using nuclear states sampled on the BO energy surface with N e = N p .

If the electronic problem has been solved within the generalized Kohn-Sham (GKS) approximation and if one assumes that the single particle electron wave function Ψ Ne GKS (r, R) do not change upon addition or removal of electrons, then the Koopmans' theorem [START_REF] Koopmans | Über die zuordnung von wellenfunktionen und eigenwerten zu den einzelnen elektronen eines atoms[END_REF] holds and the electron affinity and ionization can be substituted by the N e + 1 and N e eigenvalue of the GKS system respectively,

∆ Ne+1 c (θ) = (Ne+1)R GKS (θ) Np , (4.32) ∆ Ne-1 c (θ) = (Ne)R GKS (θ) Np ,
where the translation symmetries of the supercell are invoked via the k-point or twist θ, Bloch's vectors of the supercell. Quasiparticle energy gap in the canonical ensemble is then given by,

∆ c (θ, θ ) = ∆ Ne+1 c (θ) -∆ Ne-1 c (θ ). (4.33)
In analogy with section 2.1.3, the fundamental energy gap is then,

∆ c (θ, θ ) = min φ ∆ Ne+1 c (φ) -max φ ∆ Ne-1 c (φ ) , (4.34)
where θ and θ correspond to the Bloch's vectors at which the min and max are realized. If θ = θ the gap is indirect, while when θ = θ the gap can be direct. The connection of the supercell twist angle θ to the crystal momentum k and therefore to the vertical (direct) energy gap will be discussed in section 4.5.

Quasiparticle energy gap in a semigrand canonical ensemble at finite temperature

We can extend the previously defined in section 2.1.6 formalism to compute the gap for quantum crystals at finite temperature in the grand canonical ensemble. Note that only electrons are treated in the grand canonical ensemble, that is why we call it semigrand canonical ensemble. Consider again the system of N p nuclei and N e = N p + n electrons for a given twist angle θ. In the semigrand-ensemble the partition function is

Z(µ, θ) = ∞ Ne=0
e βµNe e -βF (Ne,θ) = e -βΩ(µ,θ) (4. [START_REF] Salpeter | A relativistic equation for bound-state problems[END_REF] where the dependence on temperature, volume, and the number of nuclei has been kept implicit and F and Ω are respectively the Helmholtz free energy and the grand-potential.

For k B T much smaller than the energies of electronic excitations (either at fixed N e or different N e ), electrons can be assumed to be in the ground state and the sum over The Helmholtz free energy includes the average nuclear kinetic energy, the average potential energy over the ground state Born-Oppenheimer (BO) surface of the Ne electron system, and the nuclear entropy. To reduce finite size effects we can average over the twists, as in section 2.1.6, obtaining the free energy and electron densities, ω(µ) = f (n e (µ e ))µ e n e (µ) (4.37)

f (n e (µ)) = 1 M θ V θ F ( N , θ) (4.38
)

n e (µ) = 1 M θ V θ N (µ, θ) (4.39)
As in the ideal crystal case, the fundamental energy gap is

∆ gc = µ + -µ -= df dn e n + p - df dn e n - p . (4.40)
In this work, we analyze nuclear configurations generated during CEIMC calculations performed in the canonical ensemble at N e = N p with twist averaged energies. As we don't have access to free energies, it is tempting to replace them with total energies in eq. 4.40. We assume that both the average nuclear kinetic energy and the nuclear entropy are nearly independent of the specific number of electrons and can be replaced by their values at N e = N p . To justify this we write the partition function in the canonical ensemble as in eq. [START_REF] Kubo | Statistical-mechanical theory of irreversible processes. i. general theory and simple applications to magnetic and conduction problems[END_REF] As before .. Np means that the averages are taken using nuclear states sampled on the BO energy surface with N e = N p . To ensure the convergence of the averages we consider 40 statistically independent nuclear configurations from the CEIMC trajectory. Note that putting the average over nuclear configurations outside the derivatives in eq. 4.45 gives a different value for the fundamental gap. This is the usual procedure to compute electronic properties from nuclear trajectories. Further, we call it the "semiclassical approximation". The application of these results can be found in chapters 5 and 6.

Quasi-momentum of the electronic wave function of quantum crystals

We will focus here on a single adiabatic PES, which is the N th eigenstate Ψ N (r, R) of the electrons for the corresponding ground state wave function χ N 0 (R) of the nuclei and we neglect for a moment any other adiabatic surface entering the wave function expansion in eq. 4.5. Note that it has in general lost any symmetry properties of the perfect crystal. To restore them, the electronic wave function Ψ N (r, R) can be expanded around the eigenstate of the ideal crystal Ψ R 0 kn (r) in the spirit of the Allen-Heine expansion [START_REF] Allen | Theory of the temperature dependence of electronic band structures[END_REF]. Considering U (R) = ĤR e -ĤR 0 e as a small perturbation it can then be Taylor-expanded around the ideal nuclear configuration R 0 ,

U (R) = N I l=1 ∆R l [∇ R l U (R)] R 0l + O(∆R 2 ), (4.46)
where ∆R l = R l -R 0,l and the expansion was truncated after the first order for the reasons that will be apparent later. Perturbation theory can then be used to expand the wave function Ψ N (r, R) up to first order in the nuclear displacements,

Ψ N (r, R) = Ψ R 0 kn (r) + l qm Ψ R 0 qm (r) E R 0 qm -E R 0 kn ∆R l ∇ R l qm|U (R)|kn R 0l (4.47) = Ψ R 0 kn (r) + qm Ψ R 0 qm (r) E R 0 qm -E R 0 kn ∆RA qm kn , with A qm kn,l = ∇ R l qm|U (R)|kn R 0
being the electron-nuclei coupling matrix element, and the sum over the nuclei is inside the scalar product ∆RA qm kn . The electronic eigenvalues, E R 0 kn , and eigenfunctions, Ψ R 0 kn (r), of the ideal crystal are characterized by the crystal momentum k and the Band index n. By integrating the electronic wave function in eq. 4.47 over the nuclei density distribution |χ N 0 (R)| 2 , one can restore the symmetry properties of the perfect crystal,

ΨN (r) ≡ dR|χ N 0 (R)| 2 Ψ N (r, R) (4.48) ≈ dR|χ N 0 (R)| 2 Ψ R 0 kn (r) + qm Ψ R 0 qm (r) E R 0 qm -E R 0 kn (R -R 0 )A qm kn = Ψ R 0 kn (r).
Note that, if the higher order terms had been included in the expansion in eq. 4.46, then the symmetry properties of the perfect crystal would not be restored. Therefore, within the harmonic approximation, the wave function of the perfect crystal is restored when integrating over the nuclear ground state density. In particular, the translational symmetries of the integrated wave function ΨN (r) are the same as of the ideal crystal wave function Ψ R 0 kn (r), so that the crystal wave vector k of Ψ R 0 kn (r) carries over to ΨN (r). However, note that the BO electronic wave function Φ N (r, R) that depends parametrically on R is defined up to an arbitrary phase φ(R) which depends on the nuclear position R, this arbitrariness in the choice of phases is called gauge arbitrariness [START_REF] Grosso | Solid State Physics[END_REF] Φ

N (r, R) Ψ N (r, R)e iφ(R) . (4.49)
Therefore if the phase φ(R) was known one could determine the crystal momentum k by calculating, ΨN (r)

≡ dR|χ N 0 (R)| 2 Φ N (r, R)e -iφ(R) ≈ Ψ R 0 kn (r). (4.50)
A numerical determination of φ(R) could be possible by adiabatic changes of the nuclear positions, R, unfortunately that will be out of the scope of this thesis. Alternatively, the crystal momentum of ΨN (r) can be determined by designing an appropriate combination of matrix elements that will cancel the phase. Consider the overlap between the BO wave function Φ N (r, R) and the perfect crystal wave function Ψ R 0 qm (r),

Ψ R 0 qm |Φ N (R) drΨ R 0 qm Ψ R 0 kn (r) + pl Ψ R 0 pl (r) E R 0 pl -E R 0 kn (R -R 0 )A qm kn e iφ(R) = δ qk δ mn + (1 -δ qk δ mn ) (R -R 0 ) A qm kn E R 0 qm -E R 0 kn e iφ(R) . (4.51)
Then, multiplying the previous equation by the overlap Ψ R 0 qm |Φ N (R) * with the same momentum q and different band index m , one can eliminate the phase φ(R),

Ψ R 0 qm |Φ N (R) Ψ R 0 qm |Φ N (R) * = δ qk δ m n (1 -δ qk δ mn ) (R -R 0 ) A qm kn E R 0 qm -E R 0 kn +δ qk δ mn (1 -δ qk δ m n ) (R -R 0 ) A qm kn * E R 0 qm -E R 0 kn + (1 -δ qk δ mn ) (1 -δ qk δ m n ) (R -R 0 ) 2 (4.52)
× A qm kn E R 0 qm -E R 0 kn A qm kn * E R 0 qm -E R 0 kn .
Further, to get rid of the last term in eq. 4.52, it can be multiplied by (R -R 0 ) and integrated over the nuclei density distribution |χ N 0 (R)| 2 , giving

T N kn (q, m; q, m ) = dR|χ N 0 (R)| 2 (R -R 0 ) Ψ R 0 qm |Φ N (R) Ψ R 0 qm |Φ N (R) * = δ qk δ m n (1 -δ mn ) A qm kn E R 0 qm -E R 0 kn + δ mn (1 -δ m n ) A qm kn * E R 0 qm -E R 0 kn × dR|χ N 0 (R)| 2 (R -R 0 ) 2 . (4.53)
Therefore, a non-vanishing T N kn (q, m; q, m ) means that the Born-Oppenheimer electronic wave function, Φ N (r, R), averaged over the nuclear ground state, has the crystal momentum k. A plot of a square modulus of the quantity defined in eq. 4.53 as a function of q is given in fig. 4.3. It is clear that the maximum overlap is when the crystal momentum of the thermal crystal is equal to the one of the ideal crystal, q = k.

For all the systems considered here we have verified, analyzing the overlap T N kn (q, m; q, m ) within a DFT-HSE wave function, that the information on the translational symmetries of the BO electronic wave function Φ N (r, R) can be obtained from the ideal crystal states

k k + G 1 k + G 2 k + G 3 k + G 4 k + G 5 k + G 6 k + G 7 k + G 8

Band structure at finite temperature

Having the thermally renormalized electronic energies for a different number of electrons E 0,k (R; N e , N p ) Np , where again, .. Np means that the averages are taken using nuclear states sampled on the BO energy surface with N e = N p . By correctly determining the crystal momentum k of extra electrons or holes, it is possible to plot the thermal band structure. Note that as the nuclear simulations are usually performed in large supercells, in order to correctly determine the crystal momentum, it is necessary to "unfold" the electronic energies determined in the reduced Brillouin zone (BZ) of the supercell into the BZ of the primitive cell. More details on the refolding procedure are present in Appendix C. Including the crystal momentum into the eq. 4.31, the band structure at finite temperature is defined as,

Nek ≡ E 0,k (R; N e + 1, N p ) -E 0,k (R; N e , N p ) Np , (4.54)
and in GKS-DFT it is equivalent to eq. 4.32,

nk = R nk Np , (4.55)
and amounts to averaging the electronic eigenvalues over nuclear fluctuations. The reduction of the gap due to the nuclear effects is of the order of 2 eV, we believe that such a strong reduction is mainly caused by the nuclear quantum effects, which are large in hydrogen. More details on the electronic gaps of thermal hydrogen crystals will be given in the next chapter.

X+

Optical properties renormalization

In the section to come, we will discuss how to account for the nuclear motion when computing the optical properties of solids. For practical reasons we only focus on the single particle optical properties, namely computed with the Kubo-Greenwood formula as in eq. 3.38.

Semiclassical averaging

Let us assume the product approximation to the exact electron-nuclear wave function |αn |Ψ R α |χ αn , where {Ψ R α , E R α } is the solution of electronic problem depending parametrically on the nuclear configuration R and |χ αn is the nuclear wave-function on the Born-Oppenheimer energy surface E R α . Assuming the electrons are initially in the ground state, one can write the real part of KG conductivity (and similarly for the imaginary part) as a thermal average over nuclear states,

σ 1 (ω, T ) = 1 Z n e - E 0n k B T σ 1 (ω, n), (4.56)
where E 0n are the eigenvalues of the nuclear motion in the BO ground electronic states, Z = n e -E 0n /k B T is the partition function and k B is Boltzmann constant. In the single electron representation -the Kubo-Greenwood theory, σ 1 (ω, n) takes the form (in SI units),

σ 1 (ω, n) ∝ 1 ω occ. α unocc. β,m χ αn | P R αβ |χ βm × χ βm | P R βα |χ αn δ( αn -βm -ω), (4.57)
where the constant prefactor was omitted for convenience, α indicate Kohn-Sham initial states in the valence band |φ R α , β and m indicate, respectively, final electronic and nuclear states in the conduction band,

P R αβ = φ R α | ∇ |φ R β
is the matrix element of the single electron momentum operator at fixed nuclear configuration R, and αn = χ αn | R α |χ αn are the joint electron-nuclear eigenvalues. The conventional quasiclassical procedure introduced by Williams [START_REF] Williams | An absolute theory of solid-state luminescence[END_REF] and Lax [START_REF] Lax | The franck-condon principle and its application to crystals[END_REF] (WL) substitutes the final nuclear states with a continuum. In practice it replaces the eigenvalues αn in Eq. (4.57) by the eigenvalues evaluated at fixed nuclear configuration R α , a procedure that can be justified as discussed in refs. [START_REF] Lax | The franck-condon principle and its application to crystals[END_REF][START_REF] Patrick | Unified theory of electron-phonon renormalization and phonon-assisted optical absorption[END_REF].

σ W L 1 (ω, n) = 1 ω occ. α unocc. β χ αn | |P R αβ | 2 δ(∆ R α,β -ω) |χ αn (4.58)
Using second order perturbation theory, it can be shown that this expression considers in an effective way the phonon-assisted indirect transitions [START_REF] Patrick | Unified theory of electron-phonon renormalization and phonon-assisted optical absorption[END_REF][START_REF] Zacharias | One-shot calculation of temperaturedependent optical spectra and phonon-induced band-gap renormalization[END_REF]. However, for light nuclei as in the case of hydrogen, replacing the nuclear spectrum by a classical continuum might not be accurate enough.

Quantum averaging

An alternative method is to consider only direct transitions between pairs of electronic states of thermally averaged bands. This procedure will include temperature renormalization of the bands but, assuming momentum conservation does not include indirect transitions. In practice, we replace the eigenvalues in eq. 4.57 by its thermal average.

To justify this approximation, consider, in analogy with Lax [START_REF] Lax | The franck-condon principle and its application to crystals[END_REF], the integral repre-sentation of the delta function in time domain in eq. 4.57,

δ( αn -βm -ω) = 1 +∞ -∞ dt exp [i( αn -βm -ω)t/ ] . (4.59)
Then, the conductivity, σ 1 (ω, n), can be re-written in the form

σ 1 (ω, n) = 1 ω occ. α unocc. β,m +∞ -∞ dt χ αn | P R αβ |χ βm χ βm | P R βα |χ αn × exp [i( αn -βm -ω)t/ ] . (4.60)
Further, considering the low temperature, nuclear states will then occupy only the ground state and the few low lying excitations. The energy spectrum of the nuclei states is less quantum, and therefore the eigenvalue differences between the nuclear states within the one electronic PES are small compared to the energy difference between different electronic PES and can be neglected, resulting in

βm = βn + ( βm -βn ) ≈ βn . (4.61)
Note that, when the nuclei are in their ground state and the nuclear ground state wave function does not change significantly when one electron is excited, the following approximation becomes exact. On the other hand, when the temperature is high, then the thermal occupation of nuclear states grows and the difference ( βmβn ) cannot be neglected anymore and the following approximation breaks down.

The summation over nuclear states m can be replaced with the completeness relation

m |χ βm χ βm | = 1, σ 1 (ω, n) = 1 ω occ. α unocc. β +∞ -∞ dt χ αn | |P R αβ | 2 |χ αn e i( αn-βn -ω)t/ . (4.62)
The thermal average over the initial nuclear states n in eq. 4.56 now applies only to the following part of eq. 4.62,

1 Z n e -βEαn χ αn | |P R αβ | 2 |χ αn e i( αn-βn )t/ = |P R αβ | 2 exp [i∆ αβ t/ ] T (4.63) ≈ |P R αβ | 2 T exp i ∆ αβ T t - σ 2 t 2
σ 1 (ω, T ) = 2πe 2 2 m 2 Ω 1 ω occ. α unocc. β ∞ -∞ dt |P R αβ | 2 T exp i ( ∆ αβ T -ω)t - σ 2 t 2 2 = 2πe 2 2 m 2 Ω 1 ω √ π σ occ. α unocc. β |P R αβ | 2 T exp - ( ∆ αβ T -ω) 2 4σ 2 , (4.64)
where time integral is a Fourier transform of a gaussian, which is itself a Gaussian. We have obtained effectively a Kubo-Greenwood conductivity, where the delta function is represented as a gaussian with the eigenvalues and matrix elements averaged over the nuclear states and the smearing defined by the variance of the eigenvalues. We will further denote this expression as Quantum Averaging. To summarize, the eq. 4.64 is valid only at low temperature and when the transition matrix elements do not correlate with the eigenvalues, which effectively means that the transitions are computed at the thermally averaged electronic eigenvalues. The illustration of the different renormalization procedures of the absorption spectra of solid C2c-24 hydrogen at 200 K is presented in fig. 4.5(a). When looking at the large absorption (inset of fig. 4.5(a)), the difference is negligible, however, looking at absorption onset (below 40000 cm -1 ) we clearly see the difference between the two procedures. Especially at 290 GPa, the semiclassical way predicts the closed gap, i.e. finite absorption at zero frequency, when in reality the gap is ∼ 0.4 eV. Moreover, the quantum procedure brings the absorption onset closer to the experimental values, which is discussed in chapter 5.

Conclusion

In this chapter, we have developed the necessary formalism to treat the excitations and optical properties in the presence of finite temperature nuclei. We have first considered the path integral formalism, defining the partition function in the Born-Oppenheimer and adiabatic approximation. Further, we have discussed the electron addition and removal energies in the canonical and semi-grandcanonical ensemble at finite temperature. We have argued that the correct quasiparticle and fundamental gap can only be obtained if the addition and removal energies are first averaged over the nuclear fluctuations.

To extract the band structure from the supercell calculations at finite temperature, the crystal momentum has to be determined first. We have provided a general formalism of determining the crystal momentum of thermal crystals, which we further use to plot the band structure of solid crystalline hydrogen at 200 K.

Finally, we have introduced an alternative to the semiclassical procedure of renormalizing optical properties obtained with single electron theory. We have argued that the new procedure is more suitable for systems at low temperatures with large zero-point nuclear effects. We call the new procedure quantum averaging. We illustrate the effect of quantum averaging by computing the absorption coefficient of solid crystalline hydrogen at 200 K using single electron Kubo-Greenwood theory.

109

Chapter 5

Metallization of crystalline molecular hydrogen

The theoretical framework developed in chapters 2, 3 and 4 is applied to study optical properties and gap closure with pressure of crystalline molecular hydrogen. We start by discussing the solid hydrogen phase diagram and possible metallization picture. The temperature range considered in this work is between 100 K and 430 K. Depending on the structure, we find that the fundamental indirect gap closes between 380 and 530 GPa for ideal crystals and 330-380 GPa for quantum crystals. Beyond this pressure, the system enters into a bad metal phase where the density of states at the Fermi level increases with pressure up to ∼450-500 GPa when the direct gap closes. To draw a connection between theoretical and experimental work, we compute optical properties using Kubo-Greenwood formalism. Furthermore, a procedure of benchmarking DFT functionals, which is based on mapping DFT and QMC densities of states, is introduced.

For an extensive review of hydrogen, its experimental and theoretical ways of studying, consider the review paper of McMahon [START_REF] Mcmahon | The properties of hydrogen and helium under extreme conditions[END_REF], which is accompanied by the recent update by Goncharov [START_REF] Goncharov | Phase diagram of hydrogen at extreme pressures and temperatures; Updated through 2019 (Review article)[END_REF].

Introduction

The metallization of crystalline hydrogen under pressure has attracted considerable attention over the last century. Predicted to be stable in an atomic bcc lattice around 25GPa, the mechanism for molecular dissociation was first discussed by Wigner and Huntington [START_REF] Wigner | On the Possibility of a Metallic Modification of Hydrogen[END_REF]. The search for its metallization has driven high pressure research until the recent [START_REF] Ranga | Observation of the Wigner-Huntington transition to metallic hydrogen[END_REF], still debated [START_REF] Goncharov | Comment on 'Observation of the Wigner-Huntington transition to metallic hydrogen[END_REF][START_REF] Loubeyre | Comment on: Observation of the Wigner-Huntington transition to metallic hydrogen[END_REF][START_REF] Liu | Comment on " Observation of the Wigner-Huntington transition to metallic hydrogen[END_REF][START_REF] Silvera | Response to Comment on "Observation of the Wigner-Huntington transition to metallic hydrogen[END_REF], observation of reflective samples at 495GPa in a Diamond Anvil Cell (DAC) apparatus. Even though it is the simplest element and H 2 the simplest homonuclear molecule in nature, the study of hydrogen under extreme conditions has uncovered rich and unexpected physics [START_REF] Mcmahon | The properties of hydrogen and helium under extreme conditions[END_REF][START_REF] Mao | Ultrahigh-pressure transitions in solid hydrogen[END_REF][START_REF] Nellis | Dynamic compression of materials: metallization of fluid hydrogen at high pressures[END_REF].

The mechanism by which the insulating crystal transforms into a conducting crystal is still unclear. Experiments have difficulty in determining the crystalline structure and its evolution with pressure because of the low cross section to X-rays [START_REF] Loubeyre | X-ray diffraction and equation of state of hydrogen at megabar pressures[END_REF][START_REF] Ji | Ultrahigh-pressure isostructural electronic transitions in hydrogen[END_REF][START_REF] Dubrovinsky | No evidence of isostructural electronic transitions in compressed hydrogen[END_REF] and the small volume of the samples for neutron scattering. Structural information is obtained indirectly through vibrational spectroscopy while the electronic structure is probed by optical measurements [START_REF] Silvera | Metallic Hydrogen[END_REF]. Direct measurements of static conductivity in the DAC remain inconclusive [START_REF] Eremets | Conductive dense hydrogen[END_REF][START_REF] Nellis | Has Metallic Hydrogen Been Made in a Diamond Anvil Cell?[END_REF][START_REF] Goncharov | To the Editor: comment on Eremets[END_REF][START_REF] Goncharov | Bonding, structures, and band gap closure of hydrogen at high pressures[END_REF][START_REF] Eremets | Low temperature phase diagram of hydrogen at pressures up to 380 GPa. A possible metallic phase at 360 GPa and 200 K. jan[END_REF][START_REF] Eremets | Molecular semimetallic hydrogen[END_REF]. A complex phase diagram comprising up to at least four different molecular phases (from I to IV) with different vibrational spectra has been traced experimentally [START_REF] Mcmahon | The properties of hydrogen and helium under extreme conditions[END_REF] (see figure 5.1). Recent experiments [START_REF] Loubeyre | Synchrotron infrared spectroscopic evidence of the probable transition to metal hydrogen[END_REF][START_REF] Ranga | Observation of the Wigner-Huntington transition to metallic hydrogen[END_REF][START_REF] Eremets | Molecular semimetallic hydrogen[END_REF][START_REF] Dias | Quantum phase transition in solid hydrogen at high pressure[END_REF][START_REF] Ranga P Dias | New insulating low temperature phase in dense hydrogen: The phase diagram to 421 GPa[END_REF] searched for metallization at low temperature (≤ 100K) while raising pressure in phase III. The phase III, also known in the literature as hydrogen-A (H-A) phase occurs above approximately 160 GPa. Spectroscopically determined thermodynamic stability range of this phase extends to pressures beyond 300 GPa and temperatures up to 300 K [START_REF] Sheng Zha | Synchrotron infrared measurements of dense hydrogen to 360 GPa[END_REF]. However, infrared (IR) and Raman spectroscopies have only been able to provide partial information about the molecular orientations in phase III of hydrogen. Recently, an experimental work based on the new high-quality single-crystal x-ray diffraction (XRD) technique was performed [START_REF] Ji | Ultrahigh-pressure isostructural electronic transitions in hydrogen[END_REF]. The authors have predicted that in phase III crystalline hydrogen has a lattice of molecular centers, which is close to hexagonal close packed (hcp). However, the critics show that the data remains very limited and insufficient for the complete structural determination [START_REF] Dubrovinsky | No evidence of isostructural electronic transitions in compressed hydrogen[END_REF]. The fact that the II→III transition pressure is relatively insensitive to the isotope (H 2 and D 2 ) suggests that phase III is mainly determined by the Born-Oppenheimer energy of interacting static molecules (i.e. obeys the classical ordering). This allows the use of DFT as the first approximation to search for the structure of phase III. Indeed, the study of Pickard and Needs [START_REF] Chris | Structure of phase III of solid hydrogen[END_REF] revealed that a monoclinic structure C2/c-24 with 12 molecules per unit cell provides a good match to the experimental vibrational data for phase III and is the lowest-enthalpy phase over the pressure range, where phase III is observed (160-300 GPa). However, by including the nuclear quantum effects a new structure, hexagonal P 6 1 22 closer related to the hcp, has emerged and was shown to be more stable than C2/c-24 below 200 GPa [START_REF] Monserrat | Hexagonal structure of phase III of solid hydrogen[END_REF].

Phase diagram

Considerable attention has also been paid to the higher temperature phase IV since its discovery [148, 154-156, 161, 162]. Predicted above 220 GPa at room temperature [START_REF] Eremets | Conductive dense hydrogen[END_REF][START_REF] Howie | Proton tunneling in phase IV of hydrogen and deuterium[END_REF] with the Raman spectra suggesting the existence of two distinct local environments, phase IV could be a mixture of graphene like layers and unbound hydrogen molecules (see fig. 5.2(c)), which is thought to be a precursor of molecular dissociation. Indeed, repeating the study of Pickard and Needs [START_REF] Chris | Structure of phase III of solid hydrogen[END_REF], but using larger unit cells, several consistent mixed phases were found with space groups Pc [START_REF] Pickard | Density functional theory study of phase IV of solid hydrogen[END_REF][START_REF] Pickard | Erratum: Density functional theory study of phase IV of solid hydrogen[END_REF], with the best prediction being Pc-48 structure. Further increasing pressure it was shown by Dalladay-Simpson et al. [START_REF] Dalladay-Simpson | Evidence for a new phase of dense hydrogen above 325 gigapascals[END_REF], based on the changes in vibrational spectra, that the new phase (V) will appear, which should coexist with phase IV in some P-T range [START_REF] Dias | Quantum phase transition in solid hydrogen at high pressure[END_REF][START_REF] Ranga P Dias | New insulating low temperature phase in dense hydrogen: The phase diagram to 421 GPa[END_REF], because no first order phase transition was observed [START_REF] Eremets | Low temperature phase diagram of hydrogen at pressures up to 380 GPa. A possible metallic phase at 360 GPa and 200 K. jan[END_REF] but a gradual phase change from IV to V.

At higher pressure, the emerging metallization picture is that the transparent insulating molecular phase III transforms into a strongly absorbing (in the visible) molecular phase at ∼ 350-360GPa with different IR frequencies, first named phase V [START_REF] Eremets | Low temperature phase diagram of hydrogen at pressures up to 380 GPa. A possible metallic phase at 360 GPa and 200 K. jan[END_REF] and later H 2 -PRE or phase VI [START_REF] Silvera | Metallic Hydrogen[END_REF][START_REF] Dias | Quantum phase transition in solid hydrogen at high pressure[END_REF], with semiconducting characteristics [START_REF] Eremets | Semimetallic molecular hydrogen at pressure above 350 GPa[END_REF]. Hydrogen finally reaches a metallic phase with the observation of reflective samples at ∼495GPa [START_REF] Ranga | Observation of the Wigner-Huntington transition to metallic hydrogen[END_REF], although disagreement concerning the pressure scale still remains [START_REF] Loubeyre | Comment on: Observation of the Wigner-Huntington transition to metallic hydrogen[END_REF][START_REF] Silvera | Metallic Hydrogen[END_REF][START_REF] Eremets | Comments on the claimed observation of the Wigner-Huntington Transition to Metallic Hydrogen[END_REF]. New synchrotron infrared spectroscopy measurements [START_REF] Loubeyre | Synchrotron infrared spectroscopic evidence of the probable transition to metal hydrogen[END_REF] report a reversible collapse of the IR transmission spectrum at 427GPa, interpreted as a first order transition to the metallic state.

High pressure hydrogen crystal structures

The primary information for theoretical investigations of solids is the crystalline structure. Candidate structures for high pressure phases have been obtained by ab initio Random Structural Search methods [START_REF] Chris | Structure of phase III of solid hydrogen[END_REF][START_REF] Pickard | Density functional theory study of phase IV of solid hydrogen[END_REF][START_REF] Pickard | Erratum: Density functional theory study of phase IV of solid hydrogen[END_REF][START_REF] Monserrat | Structure and metallicity of phase V of hydrogen[END_REF], most of them are formed by bidimensional layers stacked in various ways.

One of the most probable candidate structures for phase III is the layered C2/c-24 structure with 24 atoms in the primitive cell. Indeed, it gives the lowest free energy in ground state QMC calculations assuming harmonic phonons corrections (with DFT-PBE) [START_REF] Mcminis | Molecular to Atomic Phase Transition in Hydrogen under High Pressure[END_REF][START_REF] Azadi | Dissociation of high-pressure solid molecular hydrogen: A quantum monte carlo and anharmonic vibrational study[END_REF]. The four layers are depicted in figure 5.2(a). Stacked in ABCDA fashion, each layer consists of molecules almost parallel to the respective planes. Such arrangement of the molecules creates a non vanishing electric dipole moment, leading to a strong infrared signal [START_REF] Chris | Structure of phase III of solid hydrogen[END_REF], which is qualitatively compatible with experimental results [START_REF] Howie | Proton tunneling in phase IV of hydrogen and deuterium[END_REF].

A structure that is competitive at higher pressures (P > 250 GPa) is the Cmca-12 [START_REF] Mcminis | Molecular to Atomic Phase Transition in Hydrogen under High Pressure[END_REF][START_REF] Azadi | Dissociation of high-pressure solid molecular hydrogen: A quantum monte carlo and anharmonic vibrational study[END_REF] (fig. 5.2(b)). The structure is similar to C2/c-24, with slightly different stacking of the layers in ABAB fashion and with the molecules being completely parallel to the planes.

For Phase IV we consider only Pc-48, since the recently proposed Pca21 structure [START_REF] Monserrat | Structure and metallicity of phase V of hydrogen[END_REF] is found to be rather similar to Pc48 after geometry relaxation with the vdW-DF functional. Depicted in fig. 5.2(c), Pc-48 structure which has 48 atoms in the primitive cell. The structure is layered in ABAB fashion. While the A layers are formed by molecules, the atoms in the B layers form a distorted hexagonal network. The distortion of the hexagonal lattice results in different distances among the first neighbours.

We first consider ideal crystal structures (protons fixed at lattice sites) relaxed at constant pressure with the DFT-vdW-DF functional. Quantum crystals, with protons represented by path integrals at finite temperature, are addressed with CEIMC at constant volume 1 . All systems considered have 96 protons in nearly cubic supercells. Optimized Slater-Jastrow-Backflow trial wave functions have been used for the CEIMC calculations [START_REF] Pierleoni | Liquid-liquid phase transition in hydrogen by coupled electron-ion Monte Carlo simulations[END_REF]; details of the CEIMC simulations are reported in ref. [START_REF] Rillo | Coupled Electron-Ion Monte Carlo simulation of hydrogen molecular crystals[END_REF]. Averages over ionic positions for gaps are obtained using 40 statistically independent configurations from the CEIMC trajectories.

Fundamental energy gap

In this section we investigate the closure of the electronic gap of candidate structures for phase III (Cmca-12 and C2/c-24) and phase IV (Pc48) [START_REF] Chris | Structure of phase III of solid hydrogen[END_REF][START_REF] Pickard | Density functional theory study of phase IV of solid hydrogen[END_REF] within a Quantum Monte Carlo (QMC) framework introduced in section 2.1.6. For ideal structures, the fundamental gap decreases with pressure from ∼ 3-3. [START_REF] Mcminis | Molecular to Atomic Phase Transition in Hydrogen under High Pressure[END_REF]). These structures were optimized with vdW-DF functional. QMC for C2/c-24 optimized with the BLYP from ref. [START_REF] Sam Azadi | Nature of the metallization transition in solid hydrogen[END_REF] (closed blues triangles). GW results from Refs. [START_REF] Liang | Band gaps and effective oscillator models for solid hydrogen and h 2 o ice at high pressure[END_REF][START_REF] Lebegue | Semimetallic dense hydrogen above 260 GPa[END_REF][START_REF] Dvorak | Quasiparticle energies and excitonic effects in dense solid hydrogen near metallization[END_REF] for C2/c-24 (blue) and Cmca-12 (orange) optimized with the PBE functional. Note that pressures from RQMC are 10-15GPa lower than the nominal optimization pressure.

Path Integral approach. Extending the grand canonical method [START_REF] Yang | Electronic band gaps from Quantum Monte Carlo methods[END_REF] to determine the electronic gap of the quantum crystals at finite temperature (see section 4.4), we observe a strong gap reduction of ∼ 2eV due to nuclear quantum effects (NQE) while temperature effects below 300K are minor. At 200K the fundamental indirect gap closes at ∼330GPa for Cmca-12 and ∼380GPa for C2/c-24. Raising the temperature of C2/c-24 to 290K reduces the closure pressure to 340GPa while decreasing it to 100K does not give any noticeable effect. For both structures, the direct gap, as obtained by the unfolding of the supercell bands (see section 4.5), remains open up to ∼470-500GPa. Values for the C2/c-24 structure are in agreement with recent experimental data [START_REF] Loubeyre | Synchrotron infrared spectroscopic evidence of the probable transition to metal hydrogen[END_REF], although we cannot discuss the experimentally observed sudden closure at 427GPa, due to our sparse pressure resolution.

For a given fixed nuclear configuration, the fundamental energy gap is obtained by considering systems with a variable number of electrons n ∈ [-6, 6] where n = N e -N p .

For each system we perform Reptation Quantum Monte Carlo (RQMC) calculations with imaginary-time projection t =2.00 Ha -1 and time step τ =0.01 Ha -1 for up to 6 × 6 × 6 Monkhorst-Pack grid of twists. We check that those values are adequate for converging the band gaps within our resolution. The fundamental gap is obtained, according to section 2.1.6, from grand-canonical twist-averaged boundary conditions (GCTABC) RQMC and corrected for finite size effects in leading and next-to-leading order, as discussed in section 2.1.5.

Extending calculations of the fundamental gaps to quantum crystals, the trace over nuclear degrees of freedom must be taken with care. In the semiclassical approximation (see section 4.3), the fundamental gap is the smallest electronic excitation energy that occurs from quantum or thermal fluctuations of the lattice. Strictly speaking this gap is always closed, since the probability of a proton configuration with a metallic character is never exactly zero. For dense molecular hydrogen phonon energies are ∼ 0.1 -0.5 eV [START_REF] Pickard | Density functional theory study of phase IV of solid hydrogen[END_REF]. ZPM dominates for T ≤ 1000K, so the semi-classical approach is not appropriate. Electronic energies should be averaged over the nuclear configurations according to their thermal distribution. The gap will be given by the minimum of the average excitation energies, always larger than the semiclassical gap. Figure 5.5a illustrates typical results for the integrated density of states as a function of (electronic) chemical potential. The gap of the quantum crystal can be directly read off from the width of the incompressible region. More details are given in section 4.4.

Figure 5.3 shows estimates of the fundamental gap for ideal crystals versus pressure. The gap decreases with pressure in a similar fashion for all structures: Cmca-12 has the smallest gap, followed by C2/c-24 and by Pc48. We find reasonable agreement with the QMC estimates of ref. [START_REF] Sam Azadi | Nature of the metallization transition in solid hydrogen[END_REF] 2 . References [START_REF] Liang | Band gaps and effective oscillator models for solid hydrogen and h 2 o ice at high pressure[END_REF][START_REF] Lebegue | Semimetallic dense hydrogen above 260 GPa[END_REF][START_REF] Dvorak | Quasiparticle energies and excitonic effects in dense solid hydrogen near metallization[END_REF] report smaller values of the gap based on GW. We believe this disagreement is primarily due to the lattice geometry that has been optimized at constant pressure with PBE in refs. [START_REF] Liang | Band gaps and effective oscillator models for solid hydrogen and h 2 o ice at high pressure[END_REF][START_REF] Lebegue | Semimetallic dense hydrogen above 260 GPa[END_REF][START_REF] Dvorak | Quasiparticle energies and excitonic effects in dense solid hydrogen near metallization[END_REF] and with vdW-DF in the present work. It has been previously observed that PBE optimized geometries has longer H 2 bonds and smaller gap values at DFT level [START_REF] Morales | First Principles Methods: A Perspective from Quantum Monte Carlo[END_REF][START_REF] Raymond C Clay | Benchmarking exchange-correlation functionals for hydrogen at high pressures using quantum Monte Carlo[END_REF] (see fig. 3.2b for illustration). This propagates into G 0 W 0 . Indeed, GW results from structures optimized with vdW-DF [START_REF] Mcminis | Molecular to Atomic Phase Transition in Hydrogen under High Pressure[END_REF] are in excellent agreement with our predictions.

Values of the fundamental gap from GCTABC for quantum crystals at various temperatures and pressures are shown in Fig. 5.4a: they are smaller by ∼2eV with respect to the ideal crystal. ZPM is almost entirely responsible for this reduction. Note that the gap hardly changes from 300K to 200K within our estimated errors. Similar to ideal crystals, Cmca-12 gap is smaller than C2/c-24 gap at T=200K, the former closing at ∼340GPa while the latter at higher pressures ∼380GPa. As for the Pc-48 structure at T=430K (phase IV) the gap is slightly below values for C2/c-24 at 200K. Our results show that the electronic gap is fairly independent of the specific crystalline structure of the molecular quantum crystals. We also report gap values for C2/c-24 at T=200K from Path Integral Molecular Dynamics (PIMD) [START_REF] Morales | First Principles Methods: A Perspective from Quantum Monte Carlo[END_REF] with two different DFT functionals, namely HSE [START_REF] Heyd | Energy band gaps and lattice parameters evaluated with the Heyd-Scuseria-Ernzerhof screened hybrid functional[END_REF] and vdW-DF2 [START_REF] Lee | Higher-accuracy van der Waals density functional[END_REF]. As vdW-DF2 underestimates the molecular bond lengths of the ideal crystalline structure [START_REF] Raymond C Clay | Benchmarking exchange-correlation functionals for hydrogen at high pressures using quantum Monte Carlo[END_REF], its PIMD configurations are expected to bias the electronic gap towards larger values. Our results do not agree with predictions of ref. [START_REF] Azadi | Nuclear quantum effects induce metallization of dense solid molecular hydrogen[END_REF] (not Experiments: indirect gap from the Tauc analysis at 100K (phase III), (black squares) [START_REF] Loubeyre | Optical studies of solid hydrogen to 320 GPa and evidence for black hydrogen[END_REF], and at 300K (phase IV), (black triangles) [START_REF] Goncharov | Bonding, structures, and band gap closure of hydrogen at high pressures[END_REF][START_REF] Howie | Mixed molecular and atomic phase of dense hydrogen[END_REF]; direct gap at 100K (black squares) [START_REF] Loubeyre | Synchrotron infrared spectroscopic evidence of the probable transition to metal hydrogen[END_REF][START_REF] Loubeyre | Optical studies of solid hydrogen to 320 GPa and evidence for black hydrogen[END_REF].

shown) yielding a metallic state for C2/c-24 at 300GPa and 300K, and predict substantially larger gap reduction for C2/c-24 quantum crystals than ref. [START_REF] Li | Classical and quantum ordering of protons in cold solid hydrogen under megabar pressures[END_REF]. However, those works are based on less controlled assumptions such as using "scissor corrected" BLYP band structure and an ad hoc procedure for including nuclear motion. For all structures considered the observed fundamental gap is indirect. Estimate of the direct gap can be obtained by unfolding the band structure of the supercell (see section 4.5). Fig. 5.4b shows the direct gap for both C2/c-24 and Cmca-12 structures. While for the indirect gap Cmca-12 is always lower than C2/c-24, the direct gap is systematically larger. The difference between direct and indirect gap is of ∼ 1eV for C2/c-24, and of ∼ 2eV for Cmca-12. Closure of the direct gaps, obtained by linear extrapolation, occurs ∼ 450GP a in C2/c-24 and ∼ 500GP a in Cmca-12. Hence for both structures we observe an intermediate pressure region where the fundamental indirect gap is closed but the direct vertical gap remains open and decreases linearly with pressure. In this region, we expect the density of states around the Fermi level to increases progressively with pressure, as qualitatively reported in ref. [START_REF] Rillo | Coupled Electron-Ion Monte Carlo simulation of hydrogen molecular crystals[END_REF]. This indicates the formation of a bad metal with properties similar to a semi-metal upon closure of the indirect gap, a scenario strongly supporting the recently proposed experimental picture [START_REF] Eremets | Semimetallic molecular hydrogen at pressure above 350 GPa[END_REF](see also refs. [START_REF] Silvera | Metallic Hydrogen[END_REF][START_REF] Dias | Quantum phase transition in solid hydrogen at high pressure[END_REF]). The non-vanishing direct gap naturally explains the reported observation of absorbing (black) hydrogen around 320-360 GPa (depending on the experimental pressure scale) [START_REF] Loubeyre | Optical studies of solid hydrogen to 320 GPa and evidence for black hydrogen[END_REF]. Fig. 5.4b also shows the experimental estimates of both indirect and direct gaps from optical absorption. Measuring indirect gaps is difficult in hydrogen since samples are very thin and the optical signal from phonon-assisted absorption is too low to be detected [START_REF] Goncharov | To the Editor: comment on Eremets[END_REF][START_REF] Eremets | Molecular semimetallic hydrogen[END_REF]. The indirect gap value has been extracted from a Tauc analysis of the absorption profiles at 300K (Phase IV) [START_REF] Goncharov | Bonding, structures, and band gap closure of hydrogen at high pressures[END_REF][START_REF] Howie | Mixed molecular and atomic phase of dense hydrogen[END_REF] and 100K (Phase III) [START_REF] Sheng Zha | Synchrotron infrared measurements of dense hydrogen to 360 GPa[END_REF][START_REF] Loubeyre | Optical studies of solid hydrogen to 320 GPa and evidence for black hydrogen[END_REF] assuming the low energy absorption spectrum can be reliably extrapolated to zero energy. We have re-analyzed the spectra of ref. [START_REF] Loubeyre | Optical studies of solid hydrogen to 320 GPa and evidence for black hydrogen[END_REF] to extract the value of the indirect gap from a Tauc plot [START_REF] Tauc | Optical properties and electronic structure of amorphous germanium[END_REF], as was performed in ref. [START_REF] Goncharov | Bonding, structures, and band gap closure of hydrogen at high pressures[END_REF] for the data from ref. [START_REF] Howie | Mixed molecular and atomic phase of dense hydrogen[END_REF]. Details are further given in section 5.3.2. Conversely the direct gap at 100K (phase III) has been associated with the absorption edge at lower pressure [START_REF] Loubeyre | Optical studies of solid hydrogen to 320 GPa and evidence for black hydrogen[END_REF] or with full absorption at higher pressure [START_REF] Loubeyre | Synchrotron infrared spectroscopic evidence of the probable transition to metal hydrogen[END_REF] and corresponds roughly to the energy where the absorption coefficient equals 30000cm -1 . The direct gap of C2/c-24 structure is in agreement with the experimental data up to 425GPa, where experiments report a collapse of the gap value ascribed to the metallization transition [START_REF] Loubeyre | Synchrotron infrared spectroscopic evidence of the probable transition to metal hydrogen[END_REF]. Our results do not allow to predict this transition, but rule out C2/c-24 and Cmca-12 for this new metallic phase3 . For the indirect gap we predict ∼ 0.3 -0.5eV smaller values than in experiments. However, the Tauc analysis of refs. [START_REF] Goncharov | Bonding, structures, and band gap closure of hydrogen at high pressures[END_REF][START_REF] Howie | Mixed molecular and atomic phase of dense hydrogen[END_REF][START_REF] Loubeyre | Optical studies of solid hydrogen to 320 GPa and evidence for black hydrogen[END_REF] does not consider the energy of the emitted or absorbed phonons, which should be comparable to the observed discrepancy. However, excitonic effects and exciton-phonon coupling, neglected within the present approach, need to be addressed for this level of precision. In agreement with our findings, the experimental indirect gap depends little on both temperature and structure4 .

Optical properties

Next we explore optical properties computed using the Kubo-Greenwood (KG) framework with Kohn-Sham (KS) orbitals. To reduce the bias of the underlying DFT functional, we have benchmarked several XC approximations to reproduce the behavior of the QMC density of states close to the gap. In fig. 5.5a for C2/c-24 at 200K, we compare the electronic excess density, n en p , as a function of electronic chemical potential, µ, from QMC and from DFT-HSE 5 . The observed plateau at n en p = 0 is the signature of the indirect gap. Deviations from the plateau on both sides characterize the density of states of the valence and conduction band close to the band edges. As shown in Fig. 5.5a the HSE approximation provides slightly smaller values of the fundamental gap and reproduces [START_REF] Loubeyre | Optical studies of solid hydrogen to 320 GPa and evidence for black hydrogen[END_REF] (corrected by 20 GPa [START_REF] Loubeyre | Synchrotron infrared spectroscopic evidence of the probable transition to metal hydrogen[END_REF]), 386GPa -magenta filled circles and 406GPa -red filled circles [START_REF] Loubeyre | Synchrotron infrared spectroscopic evidence of the probable transition to metal hydrogen[END_REF])

reasonably well the integrated density of states from GCTABC around the Fermi energy (more details are in the next section 5.3.1). We therefore employed HSE to compute optical properties exploiting the KGEC code [START_REF] Calderín | Kubo-Greenwood electrical conductivity formulation and implementation for projector augmented wave datasets[END_REF] in the QuantumEspresso suite [START_REF] Giannozzi | QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials[END_REF][START_REF] Giannozzi | Advanced capabilities for materials modelling with quantum ESPRESSO[END_REF]. For thermal and quantum crystals considered here, the William-Lax (WL) semiclassical (SC) approximation [START_REF] Williams | An absolute theory of solid-state luminescence[END_REF][START_REF] Lax | The franck-condon principle and its application to crystals[END_REF][START_REF] Patrick | Unified theory of electron-phonon renormalization and phonon-assisted optical absorption[END_REF][START_REF] Zacharias | One-shot calculation of temperaturedependent optical spectra and phonon-induced band-gap renormalization[END_REF][START_REF] Zacharias | Stochastic Approach to Phonon-Assisted Optical Absorption[END_REF] is not appropriate as already discussed in section4.7.2 Instead of a joint density of states based on excitation energies for each nuclear configuration entering the WL expression, we have used the corresponding one based on electronic energies averaged over ionic ZPM, more appropriate for low temperatures. In Fig. 5.5b we compare the absorption profiles for C2/c-24 at T=200K and different pressures 6 to experimental profiles from Refs [START_REF] Loubeyre | Synchrotron infrared spectroscopic evidence of the probable transition to metal hydrogen[END_REF][START_REF] Loubeyre | Optical studies of solid hydrogen to 320 GPa and evidence for black hydrogen[END_REF] at T=100K. We observe a higher absorption than in experiments at comparable pressure, which cannot be explained by the temperature difference. We marked each predicted profile with a red dot at the energy corresponding to the observed direct gap and we report a thick horizontal line at 30000cm -1 the value of the absorption used in the experiments to extract the value of the direct gap. Our results at lower pressures are in reasonable agreement with this criterion. However at the higher pressure absorption at the energy gap is about 2-3 times higher than 30000cm -1 .

QMC benchmark of XC functionals

In this section we report some benchmark of various XC approximations for the band gap. We will focus on the ideal crystal in the C2/c-24 structure at two values of pressure: 300GPa and 500GPa. Figs. (5.6a-5.6b) show the values of the gap using various functionals compared with the gap from RQMC-GCTABC. HSE and vdW-DF2 are of comparable accuracy and perform better than the other approximations when the gap is large (300GPa). At 500GPa all XC approximations provides too small gaps. Again HSE and vdW-DF2 accuracy is comparable which is somehow surprising [START_REF] Morales | First Principles Methods: A Perspective from Quantum Monte Carlo[END_REF][START_REF] Raymond C Clay | Benchmarking exchange-correlation functionals for hydrogen at high pressures using quantum Monte Carlo[END_REF].

In order to get more information about energy bands around the gap edges we compare in Fig. (5.7a) the integrated density of states from vdW-DF2, HSE and RQMC-GCTABC. Besides a slightly larger incompressible region, the deviation of the QMC profile from the plateau is smoother than using either DFT approximations. The number of available states around the band edges is related to the intensity of optical properties. Unfortunately we cannot easily obtain optical properties with QMC methods. Therefore, it is important to try assessing the accuracy of the XC approximations.

As is the common practice, for each DFT approximation we first correct the value of the DFT gap to match the QMC one, and then subtract from the corrected integrated density of states the QMC DOS. This gives an indirect measure of the expected accuracy of optical properties of various XC approximations: the more accurate the approximation the smaller this difference should be. The difference is shown in figure (5.7b) for HSE The difference of the integrated density of states between DFT and QMC for the ideal C2/c-24 crystal at 300GPa. The scissor correction on the horizontal axis from the gap value has been applied to the DFT profiles before subtracting the QMC profile. For each profile, µ F has been assigned at the maximum of the valence band.

and vdW-DF2 functionals. The two functionals perform in similar way but the HSE one has slightly smaller deviations from the QMC profile and it should provide better optical properties. Note, however, that in this discussion we have neglected the excitonic effects, discussed in section 2.2.4. In the following we have used HSE in computing optical spectra.

Tauc analysis of absorption profiles

We have reanalyzed the absorption profiles of ref. [START_REF] Loubeyre | Optical studies of solid hydrogen to 320 GPa and evidence for black hydrogen[END_REF] using the indirect gap Tauc analysis [START_REF] Tauc | Optical properties and electronic structure of amorphous germanium[END_REF][START_REF] Tauc | Optical properties and electronic structure of amorphous Ge and Si[END_REF], similar to the analysis performed in ref. [START_REF] Goncharov | Bonding, structures, and band gap closure of hydrogen at high pressures[END_REF] of the data of ref. [START_REF] Howie | Proton tunneling in phase IV of hydrogen and deuterium[END_REF]. Tauc's relation describes the observed absorption profile α(ω) of the semiconductors: α(ω) ∝ ( ω -E g ) 2 / ω, where E g is the inferred band gap. In figure (5.8) we show the profiles at the four values of pressure reported in the original paper and the linear fit from which have extracted the indirect gap values reported in figure 5.4b of the main text. The pressure values have been adjusted according to the new scale reported in ref. [START_REF] Loubeyre | Synchrotron infrared spectroscopic evidence of the probable transition to metal hydrogen[END_REF]. We did not include the phononic shift, which should appear in the final expression for the onset of indirect phonon-assisted absorption [START_REF] Wooten | Optical properties of Solids[END_REF][START_REF] Grosso | Solid State Physics[END_REF]. This shift will correct the gaps reported from the Tauc analysis, by lowering the estimated gap by typical phonon energies emitted in the indirect transition since at the temperatures of the experiments the protons are in their ground state with few thermally excited phonons. At temperatures low enough, thermal occupation of involved phonons can be neglected, so that the value of the indirect gap extracted from the Tauc analysis should be reduced by the typical phonon energy [START_REF] Loubeyre | Optical studies of solid hydrogen to 320 GPa and evidence for black hydrogen[END_REF] for an indirect gap system. Values of the gap extracted from the intercept of the linear fits and the horizontal axis are reported in the legend and in figure 3 of the main manuscript.

emitted in the transition.

Optical properties: details

We consider 40 statistically independent nuclear configurations from the CEIMC trajectory to perform the thermal averaging. For each configuration, we employ HSE-DFT electronic structure with 8 × 8 × 8 k-grid and 2 × 2 × 2 q-grid to sample the Fock operator. The q-grid convergence was tested on one configuration for C2/c-24 structure at 200K and 400 GPa. Figure 5.9b illustrates the k-grid convergence. We see that going from 6 × 6 × 6 grid to 8 × 8 × 8 does not modify the onset of absorption (i.e. band gap is converged), but for the larger grid, oscillations are smaller, allowing a better comparison to experiment. We note that smearing can be used to improve k-grid convergence, i.e. larger k-grid usually implies that one can use smaller smearing. However, due to dispersion of the eigenvalues, which is of the order of 0.2 eV, we are limited on the resolution of band structure and cannot take a smaller smearing. In figure 5.9a we show the conductivity of the C2/c-24 structure at 343GPa at three different temperatures. We see an increase of conductivity with temperature, hence a semiconducting behaviour, in agreement with the observation of a small but still open gap at this pressure.

In figures (5.10a-5.10b) we show the reflectivity and the optical conductivity of the C2/c-24 structure at various pressures. The peak at ω = 0 comes from the intraband transitions with a gaussian widening of 0.2eV. We see a progressive increase in optical conductivity with pressure. These results are in agreement with ref. [START_REF] Rillo | Coupled Electron-Ion Monte Carlo simulation of hydrogen molecular crystals[END_REF].

Conclusions

Within this chapter the discussion for focused on solid crystalline hydrogen, we have first described the phase diagram of solid hydrogen at pressures between ∼150 GPa and 500 GPa and temperatures from 0 to 430 K. As the structure of crystalline hydrogen is still unknown, we have presented a discussion of the the potential candidate crystalline structures.

Further, we have studied the closure of the fundamental gap with pressure of candidate structures of molecular hydrogen in phase III (C2/c-24 and Cmca-12) and phase IV (Pc48) entirely based on Quantum Monte Carlo. For ideal structures, our gap values are in excellent agreement with GW prediction [START_REF] Mcminis | Molecular to Atomic Phase Transition in Hydrogen under High Pressure[END_REF]. Considering quantum nuclei at finite temperature, we observe a strong reduction of the energy gap with respect to the ideal structures at the same pressure (∼ 2eV) caused by ZPM. At 200K the fundamental (indirect) gap closes at ∼ 370-380GPa for C2/c-24 and at ∼ 340GPa for Cmca-12. We observe a reasonable agreement with experimental determinations of indirect gaps from optical absorption. The direct gap remains open until ∼ 450GPa for C2/c-24 and ∼ 500GP a for Cmca-12. In this range of pressure, the system is a bad metal (or semi-metal) suggesting a scenario that qualitatively supports recent experiments [START_REF] Eremets | Molecular semimetallic hydrogen[END_REF][START_REF] Dias | Quantum phase transition in solid hydrogen at high pressure[END_REF][START_REF] Ranga P Dias | New insulating low temperature phase in dense hydrogen: The phase diagram to 421 GPa[END_REF][START_REF] Eremets | Semimetallic molecular hydrogen at pressure above 350 GPa[END_REF]. In refs. [START_REF] Eremets | Molecular semimetallic hydrogen[END_REF][START_REF] Eremets | Semimetallic molecular hydrogen at pressure above 350 GPa[END_REF] no discontinuities in the Raman vibrational spectrum are reported when entering the semi-metallic phase, while in refs. [START_REF] Dias | Quantum phase transition in solid hydrogen at high pressure[END_REF][START_REF] Ranga P Dias | New insulating low temperature phase in dense hydrogen: The phase diagram to 421 GPa[END_REF] new IR vibron peaks are reported in this pressure range and ascribed to a structural phase transition. They have been tentatively assigned to a transition from the C2/c-24 to the Cmca-12 structure [START_REF] Dias | Quantum phase transition in solid hydrogen at high pressure[END_REF]. Our present results do not disprove this hypothesis. Our predictions for the direct gap are in good agreement with the experimental data at T=100K [START_REF] Loubeyre | Synchrotron infrared spectroscopic evidence of the probable transition to metal hydrogen[END_REF][START_REF] Loubeyre | Optical studies of solid hydrogen to 320 GPa and evidence for black hydrogen[END_REF]. However, our absorption profiles do not agree as well with the experimental behavior. This disagreement remains an open question.

In addition, in this chapter, we have presented conductivity and reflectivity, computed with Kubo-Greenwood formalism, C2c-24 of solid hydrogen at temperature and pressure range considered here.

125

Chapter 6

Metal insulator transition in dense liquid hydrogen

Liquid hydrogen at temperatures between 900 K and 3000 K and pressures from 70 GPa to 220 GPa, the area where the molecular dissociation occurs, is the focus of this chapter. In particular, here we address the fundamental gap closure across the liquid-liquid phase transition in hydrogen. The main finding is that the gap closure coincides with the molecular dissociation transition of liquid hydrogen. The chapter begins with a review of the experimental and theoretical works focusing on studying liquid hydrogen and its metallization. Besides the fundamental gap results, obtained with QMC, we provide a benchmark of some DFT functionals based on the QMC density of states. Finally, using the benchmarked HSE functional, we reanalyze optical absorption results provided in Rillo et al. [START_REF] Rillo | Optical properties of high-pressure fluid hydrogen across molecular dissociation[END_REF].

For a discussion of the physics of liquid hydrogen, I again refer to a review by McMahon [START_REF] Mcmahon | The properties of hydrogen and helium under extreme conditions[END_REF] with an additional discussion of the recent advances by Goncharov [START_REF] Goncharov | Phase diagram of hydrogen at extreme pressures and temperatures; Updated through 2019 (Review article)[END_REF] and Gregoryanz [START_REF] Gregoryanz | Everything you always wanted to know about metallic hydrogen but were afraid to ask[END_REF].

Introduction

Discussion of previous experimental and theoretical works

The insulator-metal (IM) transition in liquid hydrogen has been a longstanding focus in physics. Initially, the first order transition from the insulating molecular to metallic monoatomic fluid, called plasma-phase transition (PPT), was predicted theoretically based on chemical models [START_REF] Ld Landau | On the relation between the liquid and the gaseous states of metals[END_REF][START_REF] Norman | Thermodynamics of a dense plasma[END_REF][START_REF] Ebeling | Plasma phase transition in hydrogen[END_REF][START_REF] Saumon | Fluid hydrogen at high density: The plasma phase transition[END_REF]. Besides, there was proposed that above some critical temperature dielectric liquid should continuously transform to a metal [START_REF] Ld Landau | On the relation between the liquid and the gaseous states of metals[END_REF][START_REF] Ebeling | Plasma phase transition in hydrogen[END_REF]. Experimentally high P-T conditions necessary to observe the PPT can be achieved in two ways: using dynamic or static compression. Dynamically, hydrogen can be compressed Figure 6.1: Phase diagrams of hydrogen and deuterium around the liquid-liquid phase transition (LLPT) line (adapted from [START_REF] Rillo | Optical properties of high-pressure fluid hydrogen across molecular dissociation[END_REF]). Shaded lines (blue for hydrogen and red for deuterium) are the LLPT predicted by CEIMC [START_REF] Pierleoni | Liquid-liquid phase transition in hydrogen by coupled electron-ion Monte Carlo simulations[END_REF]. Filled symbols are estimates of the LLPT from the reflectivity coefficient; open symbols indicate the inception of absorption. Squares correspond to deuterium, circles to hydrogen. Shown are data from sp-DAC (green), Z-machine (orange), NIF (red), and lp-DAC methods (purple). DAC-p, data from sp-DAC corresponding to the temperature plateau from refs. [START_REF] Dzyabura | Evidence of a liquid -liquid phase transition in hot dense hydrogen[END_REF] and [START_REF] Zaghoo | Evidence of a first-order phase transition to metallic hydrogen[END_REF] (T ≤ 1700 K) and from ref. [START_REF] Ohta | Phase boundary of hot dense fluid hydrogen[END_REF] (T ≥ 1700 K); DAC-r, data from sp-DAC at R = 0.3; lp-DAC [START_REF] Jiang | A Spectroscopic Study of the Insulator-Metal Transition in Liquid Hydrogen and Deuterium[END_REF], filled purple points are conducting conditions, and open purple points are nonconducting conditions (for both hydrogen and deuterium); NIF-a, data from NIF when the absorption coefficient > 1µm -1 ; NIF-r, data from NIF at R = 0.3 [START_REF] Peter M Celliers | Insulator-metal transition in dense fluid deuterium[END_REF]; Z-a, data from Z-machine when the sample becomes dark; Z-r, data from Z-machine at the observed discontinuity in reflectivity [START_REF] Knudson | Direct observation of an abrupt insulator-to-metal transition in dense liquid deuterium[END_REF] with shock waves, following the time-varying changes in pressure, the metallic states can be detected via electrical, optical and density measurements [START_REF] Peter M Celliers | Insulator-metal transition in dense fluid deuterium[END_REF][START_REF] Knudson | Direct observation of an abrupt insulator-to-metal transition in dense liquid deuterium[END_REF][START_REF] Nellis | Minimum metallic conductivity of fluid hydrogen at 140 gpa (1.4 mbar)[END_REF][START_REF] Nellis | Electronic energy gap of molecular hydrogen from electrical conductivity measurements at high shock pressures[END_REF][START_REF] Weir | Metallization of fluid molecular hydrogen at 140 GPa (1.4 Mbar)[END_REF][START_REF] Celliers | Shock-induced transformation of liquid deuterium into a metallic fluid[END_REF][START_REF] Loubeyre | Coupling static and dynamic compressions: First measurements in dense hydrogen[END_REF][START_REF] Loubeyre | Extended data set for the equation of state of warm dense hydrogen isotopes[END_REF][START_REF] V E Fortov | Phase transition in a strongly nonideal deuterium plasma generated by quasi-isentropical compression at megabar pressures[END_REF]. Meanwhile, statically metallic liquid hydrogen can be reached in diamond anvil cell (DAC), following the controlled laser heating at constant pressure [START_REF] Dzyabura | Evidence of a liquid -liquid phase transition in hot dense hydrogen[END_REF][START_REF] Zaghoo | Evidence of a first-order phase transition to metallic hydrogen[END_REF][START_REF] Ohta | Phase boundary of hot dense fluid hydrogen[END_REF][START_REF] Jiang | A Spectroscopic Study of the Insulator-Metal Transition in Liquid Hydrogen and Deuterium[END_REF][START_REF] Stewart Mcwilliams | Optical Properties of Fluid Hydrogen at the Transition to a Conducting State[END_REF][START_REF] Zaghoo | Conductivity and dissociation in liquid metallic hydrogen and implications for planetary interiors[END_REF].

The detailed experimental and theoretical phase diagrams of hydrogen and deuterium around the PPT, which is usually called liquid-liquid phase transition (LLPT), is presented on figure 6.1.

The most direct information on the IM transition can be achieved via the conductivity measurements. The first experimental work that directly determined conductivity via the resistance measurements in the liquid hydrogen was carried out using a shock wave compression in the gas-gun experiment [START_REF] Nellis | Minimum metallic conductivity of fluid hydrogen at 140 gpa (1.4 mbar)[END_REF][START_REF] Nellis | Electronic energy gap of molecular hydrogen from electrical conductivity measurements at high shock pressures[END_REF][START_REF] Weir | Metallization of fluid molecular hydrogen at 140 GPa (1.4 Mbar)[END_REF]. To achieve high pressures, the initial shock is split into multiple, relatively weak shocks reverberating in hydrogen between two sapphire anvils. The resistance was measured by putting the electrodes to the hydrogen/anvil interface, the electrodes, in turn, were connected to the oscilloscope through a battery-charged isolated capacitor. When a shock wave transits the liquid hydrogen between the electrodes, the liquid becomes conducting and the capacitors discharge through hydrogen allowing the resistance measurements. Pressure was determined via the measured mass velocity of the initial shock and the Hugoniot equation of state (EOS) of the sapphire anvil. The other thermodynamic parameters such as density and temperature were determined from different EOS of hydrogen [START_REF] Kerley | A Model for the Calculation of Thermodynamic Properties of a Fluid[END_REF][START_REF] Ross | Linear-mixing model for shock-compressed liquid deuterium[END_REF] which result in a large temperature uncertainty. Based on the minimum conductivity of 2000 (Ωcm) -1 the IM transition was placed at 140 GPa and 2600 K. To determine the energy gap the authors fitted the conductivities in the range 93 -120 GPa to the equation for a liquid semiconductor with the thermally activated conductivity that depends on the mobility gap and the limiting value of conductivity.

Shock compression can as well be a laser-driven process [START_REF] Peter M Celliers | Insulator-metal transition in dense fluid deuterium[END_REF][START_REF] Celliers | Shock-induced transformation of liquid deuterium into a metallic fluid[END_REF] (NIF). The setup is almost identical to the one in gas-gun experiment, except that now the shock wave is created by laser irradiation of the pusher (Al or Be or Cu) and then being transmitted to the liquid hydrogen/deuterium. The shock velocity and the reflectance are measured with another lasers. Thermodynamic parameters are inferred from the velocity data based on the known EOSs. Different EOSs result again in the large uncertainty on temperature. The latest results on IM transition in liquid deuterium in this experimental setup predict it to be the first order with the critical temperature in the range of 1100 K< T C <3000 K with the metallization pressure around 200 GPa [START_REF] Peter M Celliers | Insulator-metal transition in dense fluid deuterium[END_REF]. Based on optical measurements, two transition boundaries were identified: first, sample becomes opaque, which corresponds to the onset of absorption, then follows by the increase of the reflectivity to 30%, which is attributed to the IM transition. The reported band gap is based on the empirical relations to the refractive index of semiconductor materials.

In a similar experiment, the shock wave in deuterium was created using the electromagnetic current pulse [START_REF] Knudson | Direct observation of an abrupt insulator-to-metal transition in dense liquid deuterium[END_REF] (Z-machine). The authors see the absorption increase at the same P-T range as Cellier et al. [START_REF] Peter M Celliers | Insulator-metal transition in dense fluid deuterium[END_REF], however their reflectivity increases abruptly at higher pressure, between 280 and 305 GPa. The temperature range (again inferred via the EOS) is between 1000 K and 1800 K. The band gap is not measured directly but based on the energy of absorption onset (∼ 2.3 eV) and qualitatively compared to the reanalyzed data of Weir et al. [START_REF] Weir | Metallization of fluid molecular hydrogen at 140 GPa (1.4 Mbar)[END_REF] and to first-principle density functional theory (DFT) predictions.

Hydrogen is a very diffusive material, therefore it is difficult to achieve high temperatures, required to observe the IM transition, during static compression. However, using short pulsed-laser heating it was possible to reach up to 3000 K in a DAC with compressed liquid hydrogen [179-181, 194, 197] (sp-DAC). By increasing the laser power, a plateau in temperature between 1100-2200 K and 90-160 GPa [START_REF] Zacharias | One-shot calculation of temperaturedependent optical spectra and phonon-induced band-gap renormalization[END_REF][START_REF] Dzyabura | Evidence of a liquid -liquid phase transition in hot dense hydrogen[END_REF][START_REF] Ohta | Phase boundary of hot dense fluid hydrogen[END_REF] accompanied by the increase of reflectivity and decrease in optical transmission [START_REF] Zaghoo | Evidence of a first-order phase transition to metallic hydrogen[END_REF] was interpreted as being due to the latent heat, a signature of the first order phase transition. Although, a finite element analysis (FEA) of the pulsed-laser heated DAC predicts the latent heat necessary to reach the plateau to be rather large (∼ 2eV /atom), in contrast to the theoretical predictions for the latent heat at the PPT (∼ 0.035 eV) [START_REF] Morales | Evidence for a first-order liquid-liquid transition in high-pressure hydrogen from ab initio simulations[END_REF]. The plateaus were alternatively interpreted by other authors as the onset of hydrogen absorption [START_REF] Peter M Celliers | Insulator-metal transition in dense fluid deuterium[END_REF][START_REF] Knudson | Direct observation of an abrupt insulator-to-metal transition in dense liquid deuterium[END_REF][START_REF] Goncharov | Comment on "evidence of a first-order phase transition to metallic hydrogen[END_REF]. Measured reflectivity reaches saturation at higher temperatures than the plateau [START_REF] Zaghoo | Conductivity and dissociation in liquid metallic hydrogen and implications for planetary interiors[END_REF]. Based on the Drude fit, at the saturation hydrogen is predicted to be largely atomic and degenerate, in contrast to the semiconductor model. However, below the saturation, the nature of the liquid is non-free-electron like [START_REF] Zaghoo | Dynamic conductivity and partial ionization in dense fluid hydrogen[END_REF].

Using long pulsed-laser heating, another experimental group observed similar two-stage transition: anomalous temperature behavior and the onset of absorption followed by the rapid increase of the reflectivity [START_REF] Jiang | A Spectroscopic Study of the Insulator-Metal Transition in Liquid Hydrogen and Deuterium[END_REF][START_REF] Stewart Mcwilliams | Optical Properties of Fluid Hydrogen at the Transition to a Conducting State[END_REF] (lp-DAC). However, the P-T conditions ascribed to these transitions are somewhat in disagreement with the previous DAC experiments [179-181, 194, 197]. The authors used Tauc's relation [START_REF] Tauc | Optical properties and electronic structure of amorphous Ge and Si[END_REF] to describe the observed absorption profile of the semiconducting liquid hydrogen: α ∝ ( ω -E g ) 2 / ω, where E g is the inferred band gap.

Overall, most of the experimental works conclude that the metallization of liquid hydrogen occurs in two steps: entering first into the absorbing semiconductor regime which follows by the rapid increase of reflectivity and IM transition. However, it remains uncertain the behavior of the fundamental gap: whether the transition is Mott-like temperature activated and accompanied by the continuous band overlap, or gap closure is discontinuous and coincides with the PPT.

To shine more light on the problem of the IM transition in liquid hydrogen a considerable number of theoretical studies were made in the past. Recent various theoretical investigations [START_REF] Pierleoni | Liquid-liquid phase transition in hydrogen by coupled electron-ion Monte Carlo simulations[END_REF][START_REF] Knudson | Direct observation of an abrupt insulator-to-metal transition in dense liquid deuterium[END_REF][START_REF] Gorelov | Benchmarking vdW-DF firstprinciples predictions against Coupled Electron-Ion Monte Carlo for high-pressure liquid hydrogen[END_REF][START_REF] Morales | Evidence for a first-order liquid-liquid transition in high-pressure hydrogen from ab initio simulations[END_REF][START_REF] Mazzola | Phase Diagram of Hydrogen and a Hydrogen-Helium Mixture at Planetary Conditions by Quantum Monte Carlo Simulations[END_REF][START_REF] Lorenzen | First-order liquid-liquid phase transition in dense hydrogen[END_REF][START_REF] Holst | Electronic transport coefficients from ab initio simulations and application to dense liquid hydrogen[END_REF][START_REF] Morales | Nuclear quantum effects and nonlocal exchange-correlation functionals applied to liquid hydrogen at high pressure[END_REF][START_REF] Norman | Plasma phase transition in warm dense hydrogen[END_REF][START_REF] Lu | Towards the Same Line of Liquid-Liquid Phase Transition of Dense Hydrogen from Various Theoretical Predictions[END_REF][START_REF] Hinz | Fully Consistent Density Functional Theory Determination of the Insulator-Metal Transition Boundary in Warm Dense Hydrogen[END_REF] based on the Born-Oppenheimer molecular dynamics (BOMD) and path integral Monte-Carlo (PIMC) predict the presence of the first order transition between insulating molecular and conductive monoatomic fluid. The location of the transition in P-T space is strongly influenced by the choice of the exchange-correlation (XC) approximation in the DFT driven BOMD or PIMC [START_REF] Morales | Nuclear quantum effects and nonlocal exchange-correlation functionals applied to liquid hydrogen at high pressure[END_REF][START_REF] Lu | Towards the Same Line of Liquid-Liquid Phase Transition of Dense Hydrogen from Various Theoretical Predictions[END_REF][START_REF] Hinz | Fully Consistent Density Functional Theory Determination of the Insulator-Metal Transition Boundary in Warm Dense Hydrogen[END_REF][START_REF] Knudson | Evaluation of exchange-correlation functionals with multiple-shock conductivity measurements in hydrogen and deuterium at the molecular-to-atomic transition[END_REF][START_REF] Ramakrishna | Influence of finite temperature Exchange-Correlation effects in Hydrogen[END_REF][START_REF] Hua | Thermodynamic anomalies and three distinct liquid-liquid transitions in warm dense liquid hydrogen[END_REF]. More reliable QMC-based methods (CEIMC and QMC-based molecular dynamics [START_REF] Pierleoni | Liquid-liquid phase transition in hydrogen by coupled electron-ion Monte Carlo simulations[END_REF][START_REF] Mazzola | Phase Diagram of Hydrogen and a Hydrogen-Helium Mixture at Planetary Conditions by Quantum Monte Carlo Simulations[END_REF]) predict the transition line that is in agreement with the experimental observation of the reflective sample in most of the experiments except by Knudson et al. [START_REF] Knudson | Evaluation of exchange-correlation functionals with multiple-shock conductivity measurements in hydrogen and deuterium at the molecular-to-atomic transition[END_REF].

Electronic properties necessary to identify the IM transition, such as optical conductivity, reflectivity and absorption can be computed within DFT [START_REF] Rillo | Optical properties of high-pressure fluid hydrogen across molecular dissociation[END_REF][START_REF] Pierleoni | Liquid-liquid phase transition in hydrogen by coupled electron-ion Monte Carlo simulations[END_REF][START_REF] Morales | Evidence for a first-order liquid-liquid transition in high-pressure hydrogen from ab initio simulations[END_REF][START_REF] Lorenzen | First-order liquid-liquid phase transition in dense hydrogen[END_REF][START_REF] Holst | Electronic transport coefficients from ab initio simulations and application to dense liquid hydrogen[END_REF][START_REF] Lu | Towards the Same Line of Liquid-Liquid Phase Transition of Dense Hydrogen from Various Theoretical Predictions[END_REF][START_REF] Hinz | Fully Consistent Density Functional Theory Determination of the Insulator-Metal Transition Boundary in Warm Dense Hydrogen[END_REF] by the Kubo-Greenwood formula [START_REF] Kubo | Statistical-mechanical theory of irreversible processes. i. general theory and simple applications to magnetic and conduction problems[END_REF][START_REF] Greenwood | The Boltzmann Equation in the Theory of Electrical Conduction in Metals[END_REF]. Based on the HSE density functional and nuclear trajectory from CEIMC [START_REF] Rillo | Optical properties of high-pressure fluid hydrogen across molecular dissociation[END_REF], the DC conductivity and reflectivity jump coincide with the dissociation transition, which together with the onset of absorption agrees with most experiments [START_REF] Peter M Celliers | Insulator-metal transition in dense fluid deuterium[END_REF][START_REF] Stewart Mcwilliams | Optical Properties of Fluid Hydrogen at the Transition to a Conducting State[END_REF][START_REF] Zaghoo | Striking isotope effect on the metallization phase lines of liquid hydrogen and deuterium[END_REF]. However, changing the XC approximation gives rather different results on optical properties and shifts the IM transition line, which can be mostly explained by the incorrect band gaps. Therefore, considering a correlated many-body theory, such as QMC, can give an accurate prediction of optical properties and might further serve as a benchmark for single electron theories.

In the past, within the QMC and using the many-body Kubo formula [START_REF] Kubo | Statistical-mechanical theory of irreversible processes. i. general theory and simple applications to magnetic and conduction problems[END_REF][START_REF] Lin | Electrical conductivity of high-pressure liquid hydrogen by quantum Monte Carlo methods[END_REF] electrical conductivity, computed for liquid hydrogen at temperatures above the critical, showed a good agreement with the experimental results available at the time [START_REF] Nellis | Minimum metallic conductivity of fluid hydrogen at 140 gpa (1.4 mbar)[END_REF][START_REF] Weir | Metallization of fluid molecular hydrogen at 140 GPa (1.4 Mbar)[END_REF]. However, to address the IM transition it is needed to have the information on the large temperature scale below and above the critical point. With the recently available method to accurately compute energy gaps within the QMC for ideal [START_REF] Yang | Electronic band gaps from Quantum Monte Carlo methods[END_REF] and thermal crystals [START_REF] Gorelov | Energy Gap Closure of Crystalline Molecular Hydrogen with Pressure[END_REF], we perform the next step to the fully consistent characterization of the IM transition in liquid hydrogen within this theory.

Theoretical method

Here we report results of an extensive study of the band gap closure of hydrogen near the LLTP using a recently developed QMC based method (see section 2.1.6 and 4.4) [START_REF] Yang | Electronic band gaps from Quantum Monte Carlo methods[END_REF][START_REF] Gorelov | Energy Gap Closure of Crystalline Molecular Hydrogen with Pressure[END_REF]. We have studied the systems along three isotherms: T = 900, 1500, and 3000 K. Quantum effects are addressed with path integrals at finite temperature using CEIMC at constant volume. All systems considered have 54 protons. Optimized Slater-Jastrow-Backflow trial wave functions have been used for the CEIMC calculations; details of the CEIMC simulations are reported in Ref. [START_REF] Pierleoni | Liquid-liquid phase transition in hydrogen by coupled electron-ion Monte Carlo simulations[END_REF]. Averages over ionic positions for band gaps are obtained using 40 statistically independent configurations from the CEIMC trajectories.

For a given configuration we perform several reptation quantum Monte Carlo (RQMC) calculations with a varying number of electrons N ∈ [-6, 6], where N = N p -N e . We consider the imaginary-time projection t = 2.0 and time step τ = 0.01 Ha -1 with the 6 × 6 × 6 Monkhorst-Pack grid of twists. Electronic size effects on the gap are treated as discussed in [START_REF] Yang | Electronic band gaps from Quantum Monte Carlo methods[END_REF]. The gaps of quantum liquid are computed by first averaging the electronic total energies for different number of electrons according to the thermal distribution over the nuclear configuration and then applying the grand-canonical twist-averaged boundary conditions (GCTABC) [START_REF] Yang | Electronic band gaps from Quantum Monte Carlo methods[END_REF][START_REF] Gorelov | Energy Gap Closure of Crystalline Molecular Hydrogen with Pressure[END_REF]. In the grand-canonical ensemble the fundamental gap is defined as the difference in chemical potentials between adding, µ + , and removing, µ -, the fractional number of electrons. At finite temperature this can be formulated as in eq. 4.45,

∆ gc = µ + -µ - d e(n e (µ)) Np dn e N + p - d e(n e (µ)) Np dn e n - p , (6.1)
where the derivatives are computed at n e = n p = N p /V , e(n e (µ)) Np means that the total energies at different electrons numbers are averaged over the nuclear states, N p . Calculation of optical properties was done within the single electron theory using the linear response Kubo-Greenwood formula [START_REF] Kubo | Statistical-mechanical theory of irreversible processes. i. general theory and simple applications to magnetic and conduction problems[END_REF][START_REF] Greenwood | The Boltzmann Equation in the Theory of Electrical Conduction in Metals[END_REF]. Thermodynamic averages of optical properties, computed with the HSE XC approximations, over 16 uncorrelated configurations from the CEIMC run were performed to ensure the convergence. More details on optical properties are given in the work of Rillo et al. [START_REF] Rillo | Optical properties of high-pressure fluid hydrogen across molecular dissociation[END_REF]. To compare DFT gaps to the QMC ones we rerun some of the HSE-DFT calculations reported in Ref. [START_REF] Rillo | Optical properties of high-pressure fluid hydrogen across molecular dissociation[END_REF] with increased k-point grid (8 × 8 × 8) to ensure the convergence.

To correct the band gap error when computing the optical properties within DFT, one can rigidly shift the unoccupied eigenvalues by the QMC-DFT gap difference, ∆ sc = ∆ QM C -∆ DF T . Alternatively, it is possible to shift the obtained Kubo-Greenwood conductivity directly by the ∆ sc . We verified that the two procedures are, in fact, equivalent. depends on both temperature and density as can be seen on the inset. Below the critical temperature of the LLPT, the gap closure coincides with the beginning of the coexistence region, as indicated by vertical bars of corresponding color. Note that at all temperatures the gap decreases linearly with pressure, with the slope becoming steeper at increasing temperature. The discontinuity at 1500 K is related to the pressure discontinuity (see the inset of fig. 6.3), at 900 K the gap closes earlier and the discontinuity is hidden. In any case, the determination of the transition pressure is a difficult task as the system changes character from molecular to atomic along the simulation. On the same time, the convergence test shows that for the gap, taking 40 configurations or even less is enough to perform the average. We compare our results to the various experimental estimates [START_REF] Peter M Celliers | Insulator-metal transition in dense fluid deuterium[END_REF][START_REF] Nellis | Minimum metallic conductivity of fluid hydrogen at 140 gpa (1.4 mbar)[END_REF][START_REF] Stewart Mcwilliams | Optical Properties of Fluid Hydrogen at the Transition to a Conducting State[END_REF]. Cellier et al. [START_REF] Peter M Celliers | Insulator-metal transition in dense fluid deuterium[END_REF] have extracted the gap based on the empirical relations to the refractive index data The agreement to our results is quite good, although the provided gap is for deuterium, not hydrogen. It is important to note that the gaps of Cellier et al. extrapolate to higher values of the closure pressure.

Results

Fundamental gap

Another estimate of the gap is based on the semiconductor model of thermally activated conductivity, σ, [START_REF] Nellis | Minimum metallic conductivity of fluid hydrogen at 140 gpa (1.4 mbar)[END_REF] σ(ρ, T ) = σ 0 exp(-E g (ρ)/2k B T ), (

where σ 0 is the limiting value of conductivity and E g (ρ) is the energy gap, assumed sampling. Noticeably, below the critical temperature, the gap is closed immediately after the molecules start to dissociate. At higher temperatures, on the contrary, the gap progressively closes with the dissociation. In Ref. [START_REF] Pierleoni | Electron localization properties in high pressure hydrogen at the liquid-liquid phase transition by Coupled Electron-Ion Monte Carlo[END_REF], analyzing the momentum distribution of liquid hydrogen with QMC, authors pointed out to the localized nature of electrons at temperatures above critical for all densities, consistently with our results of open gap at low molecular fraction. On the contrary, below the critical temperature, it was shown that the electronic properties of liquid change its nature, becoming more delocalized and Fermi-like, in accordance with the closed gap.

Our GCTABC method allows us to have an access to the density of states (DOS) near the conduction-valence band edges. Figure 6.3 shows the DOS at three isotherms for densities around the gap closure. Below the critical temperature, at 1500 K and 900 K, we show the DOS at four densities around the LLPT, the equation of state is plotted on the inset as reported in ref. [START_REF] Pierleoni | Local structure in dense hydrogen at the liquid-liquid phase transition by coupled electron-ion Monte Carlo[END_REF]. Remarkably, when the LLPT transition occurs, the DOS changes rapidly from semiconductor shape with just a few states around the Fermi energy to metal like almost straight line shape, which clearly indicates the metal to insulator transition. At high temperature, where the IM transition is continuous we do not see such behaviour, on the contrary, the DOS indicates that when gap closes the deep at Fermi energy will be progressively filled with states. Note that at 3000 K we cannot precisely locate the gap closure density, more calculations will be needed between r s = 1.6, where the gap is 0.8 eV and r s = 1.55, where the gap is closed.

Benchmark of XC approximations

Further, we provide a benchmark for the various XC approximations for the band gap. We will focus only on one isotherm of 1500 K and one density r s = 1.54. Four functionals are considered: non-local and semi-local van-der-Waals density functionals rVV10 [START_REF] Sabatini | Nonlocal van der waals density functional made simple and efficient[END_REF] and vdW-DF2 [START_REF] Lee | Higher-accuracy van der Waals density functional[END_REF], semi-local meta-GGA TPSS functional [START_REF] Tao | Climbing the density functional ladder: Nonempirical meta-generalized gradient approximation designed for molecules and solids[END_REF], and non-local hybrid HSE [START_REF] Heyd | Energy band gaps and lattice parameters evaluated with the Heyd-Scuseria-Ernzerhof screened hybrid functional[END_REF]. Figure 6.4(a) shows the value of the gap using different DFT functionals compared to the thermal gap from RQMC-GCTABS. The discrepancy is on the order of ∼ 1 eV with HSE and vdW-DF2 being the closest to the QMC prediction, scenario similar to solid hydrogen [START_REF] Gorelov | Energy Gap Closure of Crystalline Molecular Hydrogen with Pressure[END_REF]. We try to asses the accuracy of the intensity of optical properties computed with different XC functionals. With QMC we do not have a direct access to the optical properties, however, note that to the large extend they are defined by the density of states. For each DFT approximation we correct the value of the gap to match the QMC one and then we plot the difference in the integrated density of states between QMC and gap-corrected DFT. The results for single density, r s = 1.54, and temperature 1500 K are shown in Figure 6.4(b). Similarly to the gap comparison, HSE and vdW-DF2 perform better than the others with HSE being the best. At the same time note that the vdW-DF2 performs almost as good as HSE at higher and better at lower density at 136 the conduction side. Therefore, considering the computational cost of non-local hybrid functionals, it might be advantageous to use semi-local vdW-DF2 almost without the accuracy loss. Another important conclusion, coming from Fig. 6.4 is that even just correcting the gap error in DFT does not guarantee that the intensities of the spectra are accurate and probably underestimated within the considered XC approximations since the difference between the QMC and DFT DOS is always positive, implying that there are fewer states contributing to the DFT spectrum.

Optical properties

Plotting the HSE and QMC gap versus pressure at three isotherm in Figure 6.6(a), we notice that the shift is uniform for all densities below the critical temperature and of the order of ∼ 1 eV. At high temperature, T = 3000 K, the difference is decreasing with increased pressure. The gap closes at the same values of pressure with DFT and QMC at all temperatures, except the lowest (T = 900 K). Such difference will result in a shift in theoretical absorption profile reported previously [START_REF] Rillo | Optical properties of high-pressure fluid hydrogen across molecular dissociation[END_REF]. Indeed, Figure 6.6(b) shows that when shifting the HSE eigenvalues in order to match the QMC gap, the value of absorption at 2.3 eV plotted as a function of pressure, decreases with the shift being more pronounced at lower densities and lower temperatures.

The reflectivity can be reanalysed in the same manner. Consistently with decreased absorption at lower pressure, reflectivity decreases as well. However, we do not provide the same analyses here for the following reason: we are interested at the IM transition, which characterised by the typical values of reflectivity ∼ 0.3, the gap at this value is either small or already closed. Pressure at which reflectivity is 0.3, according to [START_REF] Rillo | Optical properties of high-pressure fluid hydrogen across molecular dissociation[END_REF] is reported on Fig. 6.6(a) as squares of the corresponding colors. Therefore, correcting the gap for the reflectivity will not give a meaningful result, because of the incorrect determination of the intensity of optical properties within HSE-DFT. As it has been seen in fig. 6.4(b), the DFT optical intensities might be underestimated.

Pressure at which hydrogen turns opaque was attributed by several experimental works to correspond to the absorption of ∼ 1 (µm) -1 [START_REF] Peter M Celliers | Insulator-metal transition in dense fluid deuterium[END_REF][START_REF] Knudson | Direct observation of an abrupt insulator-to-metal transition in dense liquid deuterium[END_REF][START_REF] Stewart Mcwilliams | Optical Properties of Fluid Hydrogen at the Transition to a Conducting State[END_REF]. Based on our QMC shifted HSE absorption, we predict that the onset of absorption line on phase diagram will shift to higher pressure, with respect to the previously reported line [START_REF] Rillo | Optical properties of high-pressure fluid hydrogen across molecular dissociation[END_REF]. This fact again indicates that the absorption intensities might be underestimated within the HSE XC approximation and that there might be an error cancellation between the underestimated band gap and underestimated intensities.

On the other side, closure of the gap, being an indicator of IM transition, perfectly underlines experimentally observed increase of reflectivity and coincides with the previously reported [START_REF] Pierleoni | Liquid-liquid phase transition in hydrogen by coupled electron-ion Monte Carlo simulations[END_REF] structural phase transition. This fact concludes the assumption that the liquid-liquid phase transition is as also an IM transition accompanied as well by the increase of the reflectivity and conductivity.

Conclusion

In this final chapter, we have addressed the fundamental gap closure across the liquidliquid phase transition (LLPT) in hydrogen using the newly developed (see section 4.4) formalism purely within QMC. The main finding was that gap closure is continuous and coincides with the the beginning of the molecular dissociation transition. The gap has been found to depend on temperature. Comparing to the previous work [START_REF] Pierleoni | Electron localization properties in high pressure hydrogen at the liquid-liquid phase transition by Coupled Electron-Ion Monte Carlo[END_REF], we find that below the critical temperature the closing is accompanied by the change of electron localization. We have benchmarked some DFT functional on the basis of the QMC inferred density of states and have found that all considered functionals underestimate the gap and underestimate the optical transition intensities (see fig. 6.4(b)). Finally, we have found that by shifting the HSE gap values to the ones of QMC, the absorption curves also shift, and therefore, a better description of the optical intensities is required.

Chapter 7 Conclusions

It is true that even having access to a "complete set of information" for a theoretical understanding of all real materials, the amount of those understood is finite. The reason essentially is that some materials, for example, hydrogen, require more accurate approximations to the solution of Schrödinger equation 1.1. However, there are methods, like QMC, that are potentially able to reach better accuracy. Nowadays, however, QMC can only reliably predict some ground state properties. The goal of this thesis, therefore, was to extend the QMC to make it as accurate in computing the excited state properties, in particular charged quasiparticle excitations.

The starting point of most of the QMC calculations is the single electron methods. Therefore, before switching to QMC, in chapter 1, the two methods, Hartree-Fock and DFT, were first discussed. Within the QMC we have introduced the variational Monte Carlo and the reptation quantum Monte Carlo. The former provides a way to statistically evaluate most observables for electronic systems using given many-body trial wave function. The latter is an almost exact method for finding the fermionic ground state of a Hamiltonian, with the main approximation being the use of an approximate nodal surface from a trial wave function.

We have extended the ground state QMC to compute the excited states. Firstly, in chapter 2 we have proposed a new way of computing the fundamental gap using QMC only, which relies on grand canonical twist averaging. Most importantly we have shown that for charged systems, finite size supercell calculations are necessarily biased by a finite size error decaying as 1/L, with L being the linear size of the system, where the prefactor is determined by the absolute value of the Madelung constant and the inverse dielectric constant. We have applied this procedure to determine the fundamental gap of molecular hydrogen at high pressure and carbon and silicon in the diamond structure at zero pressure. Secondly, we have discussed the neutral excitations, which are made by promoting an electron from the valence to conduction band and characterized by the coupling between the promoted electron and the hole that was left behind. Within the QMC, we have constructed a wave function to model the excited states, with the generalized variational principle, we can assure that the new excited stated are orthogonal to the lower lying states. Thirdly, in chapter 3 we have made an attempt to compute many-body conductivity with variational Monte Carlo using Kubo [4] formula. Our preliminary results for the C2c-24 ideal solid crystalline hydrogen indicate that the electron-hole interaction energy is about ∼0.4 eV, and concerning the conductivity, further investigations of the convergence of the spectra are necessary.

Bearing in mind the fact that the system of our interest is hydrogen at extreme conditions, we have had to consider the nuclear degrees of freedom and its effects on the electronic structure. In chapter 4 we have introduced the concept of Path Integrals, which is the basis of the coupled electron ion Monte Carlo (CEIMC), the method of preference in this work. We have formally argued the way of treating electron addition and removal energies in canonical and semigrand canonical ensemble at finite temperature. We have provided a general formalism of determining the crystal momentum of thermal crystals, which we further use to plot the band structure of solid crystalline hydrogen having the C2c-24 symmetry at 200 K and 250 GPa. An important extension of the traditional semiclassical way of temperature renormalization of the optical properties was introduced in the end of chapter 4 and named quantum averaging. We illustrate the effect of quantum averaging by computing the absorption coefficient of solid crystalline hydrogen at 200 K using single electron Kubo-Greenwood theory.

Most of the methodological developments done in this work has been driven by the practical goal to study the metallization of liquid and solid hydrogen. In chapter 5, we have studied the closure of the fundamental gap with pressure of candidate structures of solid molecular hydrogen in phase III (C2/c-24 and Cmca-12) and phase IV (Pc48) entirely based on Quantum Monte Carlo. For ideal structures, our gap values are in excellent agreement with GW prediction [START_REF] Mcminis | Molecular to Atomic Phase Transition in Hydrogen under High Pressure[END_REF]. Considering quantum nuclei at finite temperature, we observe a strong reduction of the energy gap with respect to the ideal structures at the same pressure (∼ 2eV) caused by ZPM. At 200K the fundamental (indirect) gap closes at ∼ 370-380GPa for C2/c-24 and at ∼ 340GPa for Cmca-12. We observe a reasonable agreement with experimental determinations of indirect gaps from optical absorption. The direct gap remains open until ∼ 450GPa for C2/c-24 and ∼ 500GP a for Cmca-12. We have found that in this range of pressure, the system is a bad metal (or semi-metal) suggesting a scenario that qualitatively supports recent experiments [START_REF] Eremets | Molecular semimetallic hydrogen[END_REF][START_REF] Dias | Quantum phase transition in solid hydrogen at high pressure[END_REF][START_REF] Ranga P Dias | New insulating low temperature phase in dense hydrogen: The phase diagram to 421 GPa[END_REF][START_REF] Eremets | Semimetallic molecular hydrogen at pressure above 350 GPa[END_REF]. In refs. [START_REF] Eremets | Molecular semimetallic hydrogen[END_REF][START_REF] Eremets | Semimetallic molecular hydrogen at pressure above 350 GPa[END_REF] no discontinuities in the Raman vibrational spectrum are reported when entering the semi-metallic phase, while in refs. [START_REF] Dias | Quantum phase transition in solid hydrogen at high pressure[END_REF][START_REF] Ranga P Dias | New insulating low temperature phase in dense hydrogen: The phase diagram to 421 GPa[END_REF] new IR vibron peaks are reported in this pressure range and ascribed to a structural phase transition. They have been tentatively assigned to a transition from the C2/c-24 to the Cmca-12 structure [START_REF] Dias | Quantum phase transition in solid hydrogen at high pressure[END_REF]. Our present results do not disprove this hypothesis. Our predictions for the direct gap are in good agreement with the experimental data at T=100K [START_REF] Loubeyre | Synchrotron infrared spectroscopic evidence of the probable transition to metal hydrogen[END_REF][START_REF] Loubeyre | Optical studies of solid hydrogen to 320 GPa and evidence for black hydrogen[END_REF]. However, our absorption profiles do not agree as well with the experimental behavior. This disagreement remains an open question.

In the case of liquid hydrogen, in chapter 6, we have found that the gap closure is continuous and coincides with the beginning of the molecular dissociation transition meaning that the metallization and the molecular breaking occur at the same time. Our method of computing the gaps with QMC allows to benchmark some DFT functional basing on the QMC inferred density of states. We have found that for liquid and solid hydrogen almost all considered functionals underestimate the gap and underestimate the optical transition intensities.

Certainly, the goal to have the QMC methods to be as reliable to predict excited state properties as the ground state is a formidable goal and can not be fully reached within a thesis. Therefore, it is important to continue to push the QMC methods for excited states to make it both affordable computationally and accurate at the same time. time consuming, particularly for complex trial functions. The dimension of the parameter space increases rapidly with the complexity of the system and the optimisation can become very cumbersome since it is, in general, a nonlinear optimisation problem.

The density matrix ρ(R, R; β) then becomes,

ρ(R, R; β) = P s=1 exp - M I (R (s) -R (s+1) ) 2 2τ - ij u(R (s) ij , R (s+1) ij
) e -τ P s=1 [E0(R (s) )-V pair (R (s) )]

= P pair (R)e -τ P s=1 [E0(R (s) )-V pair (R (s) )] , (D.5)

where u(R

(s) ij , R (s+1) ij
) is the interpotential action for a single pair of particles. The introduction of the potential V pair is meant to serve two purposes. First, better approximating the density matrix allows us to use a smaller number of beads to converge to the infinite limit. Second, the first product of exponents, indicated as P pair (R) in D.5, is rather cheap to evaluate, leaving the expensive part of evaluating electronic Born Oppenheimer energy to the second exponent, which suggests to use a multilevel Metropolis algorithm to sample the distribution [START_REF] Pierleoni | The Coupled Electron-Ion Monte Carlo Method[END_REF]:

• a move in the P ×3N p configurational space is proposed from point R to R according to a priori transition probability T (R, R )

• the acceptance probability for this move is

A 1 = min 1,
T (R, R )P pair (R ) T (R , R)P pair (R) (D.6)

• if accepted, then a second test is performed with acceptance A 2 = min 1, e -τ P s=1 E 0 (R (s) )-E 0 (R (s) )-V pair (R (s) +V pair (R (s) )

(D.7)

The consequence of separating the acceptance test in two phases is that the first phase (called the "pre-rejection step") is computationally cheap, as it prevents the expensive evaluation of the electronic energy for moves that are not likely to be accepted. Of course, an important role is played by the effective potential V pair : the more it resembles the true energy landscape, the more efficient this process becomes. There is an important issue, however associated with the algorithm, in the CEIMC method, variational Monte Carlo is used to compute the electronic energy differences E 0 (R (s) ) -E 0 (R (s) ), which is evaluated with a statistical error that can bias the outcome of the acceptance test. The simplest solution would be to make this error negligible with respect to the difference E 0 (R (s) ) -E 0 (R (s) ), however this would require sampling a lot of electronic configurations and would be highly inefficient. Instead a different approach can be used called: the penalty method. where u b is the penalty term obtained using an asymptotic expansion that converges as long as χ 2 /n < 1/4. This procedure is known as penalty method [START_REF] Yang | Electronic band gaps from Quantum Monte Carlo methods[END_REF]. For large n the first term is dominating, since χ 2 is a good enough estimator of σ 2 .On the other hand the purpose of the penalty method is to get consistent results even when n is not too large: in this case the additional terms in eq. D.15 provides an extra rejection factor motivated by our ignorance of σ. Comme l'hydrogène est un matériau très léger, les effets quantiques nucléaires sont importants.

Une description précise des effets nucléaires peut être réalisée par la méthode de Monte Carlo à ions et électrons couplés (CEIMC), une méthode de simulation des premiers principes basée sur le QMC. Nous utilisons les résultats de la méthode CEIMC pour discuter des effets quantiques et thermiques des noyaux sur les propriétés électroniques. Nous introduisons une méthode formelle de traitement du gap électronique et de la structure des bandes à température finie dans l'approximation adiabatique et discutons des approximations qui doivent être faites. Nous proposons également une nouvelle méthode pour calculer les propriétés optiques à basse température, qui constituera une amélioration par rapport à l'approximation semi-classique couramment utilisée.

Enfin, nous appliquons l'ensemble du développement méthodologique de cette thèse pour étudier la métallisation de l'hydrogène solide et liquide.

Nous constatons que pour l'hydrogène moléculaire dans sa structure cristalline parfaite, le gap QMC est en accord avec les calculs précédents de GW. Le traitement des effets quantiques nucléaires entraîne une forte réduction du gap ( 2 eV). Selon la structure, le gap indirect fondamental se ferme entre 380 et 530 GPa pour les cristaux idéaux et 330-380 GPa pour les cristaux quantiques, ce qui dépend moins de la symétrie cristalline. Audelà de cette pression, le système entre dans une phase de mauvais métal où la densité des états au niveau de Fermi augmente avec la pression jusqu'à 450-500 GPa lorsque le gap direct se ferme. Notre travail confirme partiellement partie l'interprétation des récentes expériences sur l'hydrogène à haute pression.

Nous explorons également la possibilité d'utiliser une représentation multidéterminante des états excités pour modéliser les excitations neutres et calculer la conductivité via la formule de Kubo. Nous avons appliqué cette méthodologie à l'hydrogène cristallin idéal et limité au niveau de Monte Carlo variationnel de la théorie.

Pour l'hydrogène liquide, la principale conclusion est que la fermeture du gap est continue et coïncide avec la transition de dissociation moléculaire. Nous avons été en mesure d'étalonner les fonctions de la théorie fonctionnelle de la densité (DFT) en nous basant sur la densité QMC des états. En utilisant les valeurs propres des calculs Kohn-Sham corrigé par QMC pour calculer les propriétés optiques dans le cadre de la théorie de Kubo-Greenwood , nous avons constaté que l'absorption optique théorique calculée précédemment s'est déplacée vers des énergies plus faibles. conditions, solid hydrogen Abstract: The hydrogen metallization problem posed almost 80 years ago, was named as the third open question in physics of the XXI century. Indeed, due to its lightness and reactivity, experimental information on high pressure hydrogen is limited and extremely difficult to obtain. Therefore, the development of accurate methods to guide experiments is essential.

In this thesis, we focus on studying the electronic structure, including excited state phenomena, using quantum Monte Carlo (QMC) techniques. In particular, we develop a new method of computing energy gaps accompanied by an accurate treatment of the finite simulation cell error. We formally relate finite size error to the dielectric constant of the material. Before studying hydrogen, the new method is tested on crystalline silicon and carbon diamond, systems for which experimental information on the gap is available. Although finite-size corrected gap values for carbon and silicon are larger than the experimental ones, our results demonstrate that the bias due to the finite size supercell can be corrected for, so precise values in the thermodynamic limit can be obtained for small supercells without need for numerical extrapolation.

As hydrogen is a very light material, the nuclear quantum effects are important. An accurate capturing of nuclear effects can be done within the Coupled Electron Ion Monte Carlo (CEIMC) method, a QMC-based firstprinciples simulation method. We use the results of CEIMC to discuss the thermal renormalization of electronic properties. We introduce a formal way of treating the electronic gap and band structure at a finite temperature within the adiabatic approximation and discuss the approximations that have to be made. We propose as well a novel way of renormalizing the optical properties at low temperature, which will be an improvement upon the commonly used semiclas-sical approximation.

Finally, we apply all the methodological development of this thesis to study the metallization of solid and liquid hydrogen. We find that for ideal crystalline molecular hydrogen the QMC gap is in agreement with previous GW calculations. Treating nuclear zero point effects cause a large reduction in the gap ( 2 eV). Determining the crystalline structure of solid hydrogen is still an open problem. Depending on the structure, the fundamental indirect gap closes between 380 and 530 GPa for ideal crystals and 330-380 GPa for quantum crystals, which depends less on the crystalline symmetry. Beyond this pressure, the system enters into a bad metal phase where the density of states at the Fermi level increases with pressure up to 450-500 GPa when the direct gap closes. Our work partially supports the interpretation of recent experiments in high pressure hydrogen. However, the scenario where solid hydrogen metallization is accompanied by the structural change, for example, a molecular dissociation, can not be disproved.

We also explore the possibility to use a multideterminant representation of excited states to model neutral excitations and compute the conductivity via the Kubo formula. We applied this methodology to ideal crystalline hydrogen and limited to the variational Monte Carlo level of the theory.

For liquid hydrogen, the main finding is that the gap closure is continuous and coincides with the molecular dissociation transition. We were able to benchmark different density functional theory (DFT) functionals based on the QMC electronic density of states. When using the QMC renormalized Kohn-Sham eigenvalues to compute optical properties within the Kubo-Greenwood theory, we found that previously calculated theoretical optical absorption has a shift towards lower energies.

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

6. 1

 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1.1 Discussion of previous experimental and theoretical works . . . . . 6.2 Theoretical method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3.1 Fundamental gap . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3.2 Benchmark of XC approximations . . . . . . . . . . . . . . . . . . . 6.3.3 Optical properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Conclusions A Monte Carlo methods and the Metropolis algorithm B Optimising the wave function C Unfolding the band structure D Coupled electron-ion Monte Carlo D.1 The Penalty method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Bibliography List of abbreviations

  List of publications related with the thesis• V. Gorelov, M. Holzmann, D. M. Ceperley and C. Pierleoni, Energy gap closure of crystalline molecular hydrogen with pressure, Phys. Rev. Lett. 124, 116401, (March 2020) • Y. Yang, V. Gorelov, C. Pierleoni, D. M. Ceperley and M. Holzmann Electronic band gaps from Quantum Monte Carlo methods, Phys.Rev.B 101, 085115, (February 2020) • V. Gorelov, C. Pierleoni and D. M. Ceperley Benchmarking vdW-DF first-principles predictions against Coupled Electron-Ion Monte Carlo for high-pressure liquid hydrogen, Contributions to Plasma Physics, DOI: 10.1002/ctpp.201800185 (February 2019) • V. Gorelov, M. Holzmann, D. M. Ceperley and C. Pierleoni, Electronic energy gap closure and metal-insulator transition in dense liquid hydrogen., accepted to Phys.Rev.B, arXiv:2009.00652, (october 2020) • V. Gorelov, M. Holzmann, D. M. Ceperley and C. Pierleoni, Electronic structure and optical properties of quantum crystals from first principles calculations in the Born-Oppenheimer approximation, accepted to Journal of Chemical Physics, arXiv:2010.01988 (october 2020) 15

Figure 1 :

 1 Figure 1: Hydrogen phase diagram with recent experimental results on hydrogen metallization

Figure 1 . 1 :

 11 Figure 1.1: A schematic illustration of the local energy as a function of Monte Carlo steps for different wave function quality, indicated by parameters a 1 and a 2 The green line illustrates an exact ground state energy.

Figure 1 . 2 :

 12 Figure 1.2: Pictorial representation of the trial moves. In the new configuration (bottom), a new head for the reptile is generated on the right side from the old configuration (top) and the tail is discarded.

Figure 2 . 1 : 2 . 1 . 5 . 1

 212151 Figure 2.1: Left: band structure plot for crystalline C2/c hydrogen at 234 GPa. The dashed lines correspond to the two highest occupied and two lowest unoccupied bands computed with the vdW-DF XC approximation and corrected to match the QMC gap. The blue points are the QMC electron addition/removal energies unfolded into the Brillouin zone -path represented on the right. Right: Brillouin zone (BZ) of the C2/c structure, green line represents the BZ path used for the band structure plot. The BZ path was shifted from the one indicated on the plot by a small constant shift to better match the QMC twist grid.

Figure 2 . 2 :

 22 Figure 2.2:Change in the static structure factor as an electron (upper curves) or a hole (lower curves) is added to the insulating system with N atoms. The lines are fits to the data points. The horizontal lines show the expected k → 0 limit based on the experimental dielectric constants. We have used c = 0.41 for C and c = 0.57 for Si determined from eq. 2.22.

Figure 2 . 3 :

 23 Figure 2.3: Upper bound to the inverse dielectric constant eq. (2.11), where Γ k ≡ 2ωpS Ne (k) k 2

2 1 - 2 1 + 1 .

 121 1 and in the hole channel, |v M | The corrections independent of correspond to the change in the background charge which cancel for the fundamental gap and we obtain eq. (2.23).

Figure 2 . 4 :

 24 Figure 2.4: GCTABC analyses of the C2/c-24 structure of solid hydrogen at r s = 1.38 (234GPa).(a) the electron density n e as a function of the chemical potential µ obtained from HSE functional in comparison to QMC, the inset illustrates the energy density as a function of µ from HSE functional. (b) energy density e 0 as a function of n e using QMC, the inset shows the derivative discontinuity where δn e is the change of the electronic density with respect to the insulating state. Size corrections as discussed in the text are included.

7 ( 1 )

 71 -0.25 (G 0 W 0 )[106] -0.06 (DF T )[START_REF] Monserrat | Comparing electron-phonon coupling strength in diamond, silicon, and silicon carbide: First-principles study[END_REF] 1.17[START_REF] Madelung | Intrinsic Properties of Group IV Elements and III-V, II-VI and I-VII Compounds[END_REF] H 2 r s = 1.38 234GPa

Figure 2 . 6 :

 26 Figure 2.6: Change to static structure factor as an electron (blue filled circle) or a hole (orange open circle) is added to the neutral system. The lines are the fits to obtain asymptotic values at k → 0.

k 1 =Figure 2 . 7 :

 127 Figure 2.7: Inverse dielectric constant upper bound for solid H 2 , where Γ k ≡ 2ωpS Ne (k) k 2

MMFigure 2 . 8 :

 28 Figure 2.8: Left: Comparison of the optical VMC gap and quasiparticle (QP) RQMC gap, computed with QMC for C2/c-24 crystalline hydrogen at approximately 285 GPa at different values of twist angle. Two different basis sizes are considered for the neutral excitations: M = 19 (orange) and 25 (green). The values without size corrections are present for the quasiparticle gap (blue line). (a) Direct/vertical excitations with crystal momentum k f = k i and supercell twist θ. (b) Indirect excitations with crystal momentum k f = k i and supercell twist θ that upon unfolding gives indirect excitations, red line is done by considering only one neutral excitation (M=1) by substituting orbital in Slater determinant by the lowest unoccupied one. Right: excitonic effect energy as defined in eq. 2.50 for different basis sizes.

Figure 3 . 1 :Figure 3 . 2 :

 3132 Figure 3.1: Convergence of the real and imaginary part of Kubo-Greenwood conductivity with respect to the simulation parameters. The spectra are calculated for C2/c-24 hydrogen crystal at 300 GPa with vdW-DF2 XC approximation. Upper panel: k-points with 56 bands and 0.2 eV smearing; middle panel: number of bands with 8 × 8 × 8 k-grid and 0.2 eV smearing; lower panel: smearing factor with 56 bands and 8 × 8 × 8 k-grid.

Figure 3 . 3 :

 33 Figure 3.3: Variational Monte Carlo Kubo electrical conductivity computed for C2c-24 ideal crystalline hydrogen at approximately 285 GPa. Two basis sets consisting of M = 19 and 25 excitations are considered. For comparison we provide Kubo-Greenwood conductivity computed with vdW-DF2 XC approximation.

  where ρ(R, Ψ α ; R, Ψ α ; β) is the density matrix at inverse temperature β = 1/k B T with k B being the Boltzmann constant. Since the nuclear kinetic energy operator does not commute with the rest of the full Hamiltonian, the lowest-order Trotter factorization[START_REF] Trotter | On the product of semi-groups of operators[END_REF],e -β[ T R +H R e

Figure 4 . 1 :

 41 Figure 4.1: Sketch of a diatomic molecule, in the path integral representation ("ring polymer").

Figure 4 . 2 ,Figure 4 . 2 :

 4242 Figure 4.2: Histogram of the distribution of the energies required to add two electrons over the nuclei fluctuations and the normal distribution fit. Analyses is done for the C2/c crystalline hydrogen at 200 K and 300 GPa.

  N e reduces to the values Ne (µ, θ) that minimizes the exponent Ω(µ, θ) = min Ne [F (N e , θ) -µN e ] = F ( Ne , θ)µ Ne . (4.36)

Figure 4 . 3 :

 43 Figure 4.3: Square modulus of the overlap defined in eq. 4.53 as a function of different momentum q. This example is from the DFT-HSE calculations of solid hydrogen in C2/c-24 structure at T = 200K and P = 250GP a. G i are the reciprocal lattice vectors

Figure 4 . 4 :

 44 Figure 4.4: Band structure at finite temperature from QMC-CEIMC calculations (orange points) and for perfect crystal QMC (blue points) of 96 hydrogen atoms compared to the band structure from vdW-DF density functional for a unitary cell of C2/c-24 hydrogen crystal at approximately 250 GPa. Horizontal lines are the corresponding valence band maximum and conduction band minimum. The Brillouin zone path is equivalent to the one in fig. 2.1 and the vdW-DF bands are shifted to match the QMC ideal bands.

Figure 4 .

 4 Figure 4.4 illustrates the QMC highest occupied and lowest unoccupied bands at 200 K and for ideal structure plotted on top of the vdW-DF DFT band structure of the unitary cell of C2/c-24 hydrogen crystal at approximately 250 GPa. The Brillouin zone path and ideal crystal QMC and DFT bands are identical to fig. 2.1.The reduction of the gap due to the nuclear effects is of the order of 2 eV, we believe that such a strong reduction is mainly caused by the nuclear quantum effects, which are large in hydrogen. More details on the electronic gaps of thermal hydrogen crystals will be given in the next chapter.

2 ,

 2 whereσ 2 = ∆ 2 αβ T -∆ αβ

Figure 4 . 5 :

 45 Figure 4.5: (a) Optical absorption for C2/c-24 quantum crystals at 200K using the semiclassical (WL) and quantum (QA) averaging procedures. The dashed vertical lines indicate the values of indirect band-gaps computed using the HSE functional. The inset shows the same plots, but for a larger range, indicating that the difference between two averaged is noticeable only at the absorption edge. (b) Histogram of the distribution of the HSE-DFT LU M O -HOM O over the nuclei fluctuations and the normal distribution fit. Analyses is done for the C2/c crystalline hydrogen at 200 K and 300 GPa.

4. 3

 3 the cumulant expansion was applied, assuming that the excitation energies αβ are normally distributed over the nuclear states. The distribution of the excitation energies LU M O -HOM O over the nuclear states is illustrated for solid hydrogen at 200 K in figure 4.5(b). Inserting the above relation back into eq. 4.56 for the conductivity σ 1 (ω, T ) and recalling the constant one gets,

Figure 5 . 1 :

 51 Figure 5.1: Low temperature experimental phase diagram of H 2 . Phase III to IV Howie et al. [154, 155], Zha et al. [156]. Phase IV to V Dalladay-Simpson et al. [157], Eremets et al.(IV/V) [152]. Phase III to H 2 -PRE Eremetz et al (VI) [152], Dias et al. [158] (red and blue). Phase H 2 -PRE to Metallic H 2 Loubeyre et al. [13]. Phase H 2 -PRE to MH Dias et al. [158] (green)

Figure 5 . 2 :

 52 Figure 5.2: (a) The C2/c-24 structure, (b) the Cmca-12 structure and (c) the Pc-48 structure based on the lattice proposed at P = 300 GPa in the supplementary material of ref. [111]. The number of atoms in the primitive cell is given in the end of the symmetry nomenclature. Right panel: a 3D view, depicting different layers stacked in an (a) ABCD, (b) AB and (c) ABAB fashion. Left panel: (a), (b) is a top view of one layer and (c) is a top view of the whole structure on the right (four layers). (Adapted from ref. [167])

Figure 5 . 3 :

 53 Figure 5.3: Fundamental energy gap for ideal crystals. This work (closed circles): C2/c-24 (blue), Cmca-12 (orange) and Pc48 (green), open GW results for C2/c-24 (open blue circles[START_REF] Mcminis | Molecular to Atomic Phase Transition in Hydrogen under High Pressure[END_REF]). These structures were optimized with vdW-DF functional. QMC for C2/c-24 optimized with the BLYP from ref.[START_REF] Sam Azadi | Nature of the metallization transition in solid hydrogen[END_REF] (closed blues triangles). GW results from Refs.[START_REF] Liang | Band gaps and effective oscillator models for solid hydrogen and h 2 o ice at high pressure[END_REF][START_REF] Lebegue | Semimetallic dense hydrogen above 260 GPa[END_REF][START_REF] Dvorak | Quasiparticle energies and excitonic effects in dense solid hydrogen near metallization[END_REF] for C2/c-24 (blue) and Cmca-12 (orange) optimized with the PBE functional. Note that pressures from RQMC are 10-15GPa lower than the nominal optimization pressure.

Figure 5 . 4 :

 54 Figure 5.4: (a) The fundamental gap of quantum crystals at finite temperature. Closed circles indicate results from this work, for the three structures at various temperature as detailed in the legend. PIMD-DFT results at 200K are obtained with two different XC approximations, namely HSE (downward open triangles) and vdW-DF2(upward open triangles) and the semiclassical averaging are reported for comparison [49]. (b) Direct (closed symbols) and indirect (open symbols) gaps of quantum crystals. GCTABC-RQMC at T=200K: C2/c-24 indirect (blue triangles), direct (blues squares); Cmca-12 indirect (orange triangles), direct (closed squares). Experiments: indirect gap from the Tauc analysis at 100K (phase III), (black squares)[START_REF] Loubeyre | Optical studies of solid hydrogen to 320 GPa and evidence for black hydrogen[END_REF], and at 300K (phase IV), (black triangles)[START_REF] Goncharov | Bonding, structures, and band gap closure of hydrogen at high pressures[END_REF][START_REF] Howie | Mixed molecular and atomic phase of dense hydrogen[END_REF]; direct gap at 100K (black squares)[START_REF] Loubeyre | Synchrotron infrared spectroscopic evidence of the probable transition to metal hydrogen[END_REF][START_REF] Loubeyre | Optical studies of solid hydrogen to 320 GPa and evidence for black hydrogen[END_REF].

1 Figure 5 . 5 :

 155 Figure 5.5: (a) Integrated density of states for C2/c-24 quantum crystals at 200K from GCTABC-RQMC (points) and HSE (smooth lines) at various pressures.(b) Absorption spectra from HSE band structure for C2/c-24 quantum crystals (solid lines) and comparison with the available experimental profiles (opened and filled circles). The spectra from HSE have been shifted in energy by an amount equal to the difference between QMC and HSE direct gap. The reported pressure are as in figure5.5a (see the colors). The red dots indicate the location in energy of the direct gap of figure 5.4b. Experimental pressures are: 296GPa -open orange circles[START_REF] Loubeyre | Optical studies of solid hydrogen to 320 GPa and evidence for black hydrogen[END_REF] (corrected by 20 GPa[START_REF] Loubeyre | Synchrotron infrared spectroscopic evidence of the probable transition to metal hydrogen[END_REF]), 386GPa -magenta filled circles and 406GPa -red filled circles[START_REF] Loubeyre | Synchrotron infrared spectroscopic evidence of the probable transition to metal hydrogen[END_REF])

Figure 5 . 6 :

 56 Figure 5.6: The fundamental gap computed using various XC approximations for the ideal C2/c-24 structure at the nominal pressure of (a) 300GPa and (b) 500 GPa. The orange horizontal bar reports the value of the RQMC-GCTABC gap and its thickness corresponds to its statistical uncertainty.

Figure 5 . 7 :

 57 Figure 5.7: (a) Integrated density of states around the band edges from RQMC-GCTABC (green points with error bars), HSE (orange) and vdW-DF2 (blue) XC approximations for the ideal C2/c-24 crystal at 300GPa. (b)The difference of the integrated density of states between DFT and QMC for the ideal C2/c-24 crystal at 300GPa. The scissor correction on the horizontal axis from the gap value has been applied to the DFT profiles before subtracting the QMC profile. For each profile, µ F has been assigned at the maximum of the valence band.

Figure 5 . 8 :

 58 Figure 5.8: Tauc analysis of the absorption profile reported in ref.[START_REF] Loubeyre | Optical studies of solid hydrogen to 320 GPa and evidence for black hydrogen[END_REF] for an indirect gap

123

 123

Figure 5 . 9 :

 59 Figure 5.9: (a) The optical conductivity of the C2/c-24 quantum crystals at 343GPa and at three temperatures using the QA procedure. (b) K-grid convergence for C2/c-24 structure at 200K and 300 GPa

Figure 5 . 10 :

 510 Figure 5.10: (a) Reflectivity for C2/c-24 quantum crystals at 200K from the QA procedure. (b) The conductivity for C2/c-24 quantum crystals at 200K using the QA procedure. The smearing was 0.2 eV.

  . Two dashed purple lines indicate the inception of absorption (McWilliams-a) and the metallic boundary (McWilliams-m)[START_REF] Stewart Mcwilliams | Optical Properties of Fluid Hydrogen at the Transition to a Conducting State[END_REF]. Brown shaded circles(Weir) show the inception of metallicity from gas gun experiments[START_REF] Nellis | Minimum metallic conductivity of fluid hydrogen at 140 gpa (1.4 mbar)[END_REF]. Blue points are theoretical estimates from Rillo et al.[START_REF] Rillo | Optical properties of high-pressure fluid hydrogen across molecular dissociation[END_REF]: Filled circles show when R = 0.3 for H/vacuum interface; open circles shown when the absorption coefficient equals 1µm -1 . Two slightly different melting lines are reported at low temperature[11] [184].

Figure 6 . 2 :

 62 Figure 6.2: (a) Fundamental energy gap of liquid hydrogen along the isotherms: T = 900 K, 1500 K and 3000 K as a function of pressure. Inset: fundamental energy gap as a function of density. Vertical bar of the corresponding color indicates the coexistence region at the LLPT according to the [213]. Dashed lines are the gaps reported by Cellier et at. [26]. Brown and green squares are the results of Nellis et al. at around 2000-3000 K [28] reanalyzed by [208], green squares correspond to the temperature around 3000 K. Red circle is the gap reported by McWilliams et al. at 2400 K [183]. (b) Molecular fraction, as a function of the fundamental energy gap for three temperatures. Inset: Molecular fraction as a function of density taken from [213].

Figure 6 . 3 :

 63 Figure 6.3: Density of states of liquid hydrogen near the band edge at densities around the gap closure for three isotherms (a) 900 K, (b) 1500 K and (c) 3000 K. Inset shows the equation of states as reported in [213], the dashed and solid red lines indicate the atomic and molecular region. The color code of the DOS matches the colors of points on the inset.

Figure 6 .Figure 6 . 4 :

 664 Figure 6.2a shows the estimates of fundamental gap computed according to eq. 6.1 for different isotherms of liquid hydrogen. The gap gradually decreases with pressure and

Figure 6 . 6 :

 66 Figure 6.6: (a) HSE and QMC band gap at three isotherms. Dashed line are the HSE values and solid line are the QMC results. Squares indicate pressure at which the reflectivity is 0.3 according to ref. [6]. (b) Absorption at ω = 2.3 eV along the T = 1500 K isotherm. Dashed line has the original HSE values reported in Ref. [6] and the solid line is computed considering the QMC corrected ban gap.

  with the new acceptance formula beinga(∆, χ 2 ) = min 1, T (R , R) T (R, R ) e ∆-u b (D.15)

Titre:

  Méthodes de monte carlo quantique pour le calcul des structures électroniques: application à l'hydrogène dans des conditions extrêmes Mots clés: Monte Carlo quantique, structure électronique, méthode CEIMC, conditions ex- trêmes, hydrogène solide Résumé: Le problème de la métallisation de l'hydrogène, posé il y a près de 80 ans, a été désigné comme la troisième question ouverte en physique du XXIe siècle. En effet, en raison de sa légèreté et de sa réactivité, les informations expérimentales sur l'hydrogène à haute pression sont limitées et extrêmement difficiles à obtenir. Il est donc essentiel de mettre au point des méthodes précises pour guider les expériences.Dans cette thèse, nous nous concentrons sur l'étude de la structure électronique, y compris les phénomènes d'état excité, en utilisant les techniques de Monte Carlo quantique (QMC). En particulier, nous développons une nouvelle méthode de calcul pour le gap accompagnée d'un traitement précis de l'erreur induit par la taille finie de la cellule de simulation. Nous établissons un lien formel entre l'erreur de la taille finie et la constante diélectrique du matériau. Avant d'étudier l'hydrogène, la nouvelle méthode est testée sur le silicium cristallin et le carbone de diamant, pour lesquels des informations expérimentales sur le gap sont disponibles. Nos résultats montrent que le biais dû à la supercellule de taille finie peut être corrigé, de sorte que des valeurs précises dans la limite thermodynamique peuvent être obtenues pour les petites supercellules sans avoir besoin d'une extrapolation numérique.

  Following the same reasoning as in eqs. 4.36-4.40 we arrive at our final expression of the fundamental gap∆ gc = µ +µ -=df (n e (µ))

	(4.42) (4.43) with n = n(µ) depending on the chemical potential. The free energy density is now, Using this expression in eq. 4.36 we obtain Ω(µ, θ) F (N p , θ) + min Ne ∆E(n, θ) Np -µN e = F (N p , θ) + ∆E(n, θ) Np -µ(N p + n), f (n e (µ)) 1 θV θ n + p -df (n e (µ)) dn e n -p (4.45) F (N dn e d e(n e (µ)) Np dn e n + p d e(n e (µ)) Np -dn e n -p

.29 Z(N e , θ) e -β[F (Np,θ)+ ∆E(n,θ) Np ] (4.41) which implies F (N e , θ) F (N p , θ) + ∆E(n, θ) Np . p , θ) + ∆E(n(µ), θ) Np (4.44) = f (n p ) + e(n e (µ))e(n p ) Np , where e(n e (µ)) = 1 θV θ E(n(µ) + N p , θ).

k u k ρ k ρ -k /2Ω, and we are only interested in the long-wave length limit, k → 0, of the electron-electron correlation, with wave vectors smaller than the reciprocal lattice

= n

k 2 , 2 = 4πσ/ω = 2nk.

T is the variance. The matrix element, |P R αβ | 2 , and the exponent, exp [i∆ αβ t/ ], were taken to be non-correlated and, by analogy with section

We have checked that the stress tensor in the constant volume CEIMC run remains diagonal with same diagonal elements within our statistical noise.

The observed small difference, in particular at the higher pressure, is probably due to the different XC approximation used for geometry optimization, vdW-DF in our case, BLYP in ref.[START_REF] Sam Azadi | Nature of the metallization transition in solid hydrogen[END_REF] and different size extrapolation.

Our estimates of the direct gap could be biased by ∼ 0.3eV due to the discreteness of our twist grid. Correcting for this bias will place the experimental data in between the C2/c-24 and Cmca-12 predictions.

Note that the pressure values of ref.[START_REF] Loubeyre | Optical studies of solid hydrogen to 320 GPa and evidence for black hydrogen[END_REF] have been recently corrected[START_REF] Loubeyre | Synchrotron infrared spectroscopic evidence of the probable transition to metal hydrogen[END_REF]

This quantity is closely related to the integrated density of states.

To partially correct for HSE inaccuracy, we shifted the energy scale by the difference between the QMC and HSE gap.

Acknowledgements

to depend linearly on the density ρ and independent of the temperature T , k B is the Boltzmann constant. We report results of Nellis et al. [START_REF] Nellis | Minimum metallic conductivity of fluid hydrogen at 140 gpa (1.4 mbar)[END_REF] reanalysed by Knudson et al. [START_REF] Knudson | Evaluation of exchange-correlation functionals with multiple-shock conductivity measurements in hydrogen and deuterium at the molecular-to-atomic transition[END_REF] who used a different equation of state [START_REF] Irwin | Equations of state for hydrogen and deuterium[END_REF] and different σ 0 . They assumed that hydrogen before the transition behaves like a fluid semiconductor, where the conductivity is progressively increased upon the closure of the gap with density. The gap is assumed to weakly depend on temperature. Temperature was not measured in the experiment and, according to the latest equation of state [START_REF] Irwin | Equations of state for hydrogen and deuterium[END_REF], varies between 2000-3000 K, increasing towards the higher pressure [START_REF] Knudson | Evaluation of exchange-correlation functionals with multiple-shock conductivity measurements in hydrogen and deuterium at the molecular-to-atomic transition[END_REF]. Below the critical temperature, our results do not fully support this model, as the QMC density of states increases rapidly at the transition (see fig. 6.3 (a) and (b)) and, in addition, our gap is temperature dependent. Above the critical temperature, we do not have enough data to asses the model, as we would need at least three isotherms. We can conclude that using a semiconductor model can only provide qualitative estimates of gap energies. Lastly, analysing absorption profile with Tauc model [START_REF] Tauc | Optical properties and electronic structure of amorphous Ge and Si[END_REF], McWilliams et al. [START_REF] Stewart Mcwilliams | Optical Properties of Fluid Hydrogen at the Transition to a Conducting State[END_REF] have reported the gap value of 0.9 eV at 2400 K and 140 GPa [START_REF] Stewart Mcwilliams | Optical Properties of Fluid Hydrogen at the Transition to a Conducting State[END_REF]. On figure 6.5 we analysed several theoretical absorption profiles trying to asses the validity of this model. We found that the fitting of the theoretical absorption to the Tauc model slightly overestimates the values of the gaps (by ∼ 0.3 eV), computed at the same level of approximation as optical properties, e.g. DFT -HSE. However, it gives a good agreement, when compared to the QMC gap, indicating the possibility of error cancellation, when calculating the spectra.

To discuss the molecular character of liquid at gap closure, in Figure 6.2a(b) we report it as a function of molecular fraction. The values of molecular fraction are taken from Ref. [START_REF] Pierleoni | Local structure in dense hydrogen at the liquid-liquid phase transition by coupled electron-ion Monte Carlo[END_REF] and represent the probability to find persistent pairs of the atoms during the CEIMC Appendix A

Monte Carlo methods and the Metropolis algorithm

Monte Carlo integration is a general purpose algorithm used to calculate averages of arbitrary functions. Consider the one dimensional integral

where g(x) is any positive definite function such that b a dxg(x) = 1, which can be interpreted as a distribution. The integral can be approximated by an average over a discrete set of points distributed according to g(x). Using the Central Limit Theorem, it is easy to show that the error in the resulting estimate is given by [START_REF] Gaskell | The Collective Treatment of a Fermi Gas: II[END_REF]

where {x i } is distributed according to g(x) and σ 2 f is the variance of f (x) over g(x). Traditional quadrature methods, which would work better in this case, have a very bad scaling with the number of dimensions. Note that because the error in eq. A.2 is independent of the number of dimensions, Monte Carlo integration is able to handle integrals in an arbitrary number of dimensions and is very efficient if the integrand is peaked in a localized region of space.

Extending this simple idea one could calculate the integral in eq. A.1 by sampling any arbitrary distribution, π(x) being any positive definite function that satisfies

This idea is typically called importance sampling. In fact, any definite integral in dx over [a, b] can be simultaneously calculated from a sample of π(x) by "reweighting" the averages according to:

where {x i } is distributed according to π(x). To sample arbitrary distributions, π(x), the Metropolis algorithm [START_REF] Frank | Many-Body Quantum Monte Carlo Study of 2D Materials: Cohesion and Band Gap in Single-Layer Phosphorene[END_REF] is used. It consists on the following three steps:

• Starting from a state x, choose a trial state x with probability T (x → x ), where T is a simple distribution function that can be sampled directly.

• Evaluate the acceptance probability:

• Accept the step x with probability A. For this generate a random number, η, uniformly distributed in [0, 1] and accept the step if η < A.

The algorithm is both general and very simple, the only knowledge needed is of the limiting distribution and the only with requirement for the transition probability, T , to be ergodic. The main drawback of the method is that the steps along the random walk are highly correlated, meaning that it takes a certain number of steps to produce statistically uncorrelated data. This, in turn, implies that longer simulations are needed to reduce error bars in the resulting averages. Nonetheless, the simplicity of the method allows us to apply it to very complicated problems like the Schrödinger equation.

Appendix B Optimising the wave function

The simplest way to find the best trial function is to perform independent VMC runs (possibly in parallel) with different values for the variational parameters. Then one fits the energies as a function of parameters to a simple analytic form, does more calculations around the predicted minimum, and continues until convergence in parameter space is attained. However this is not the most efficient method as close to the minimum the independent statistical errors will mask the variation with respect to the trial function parameters. And what is more important, doing optimisation by hand is limited to a few parameters. There exist several other techniques, for example the reweighting method, that allow efficient wave function optimisation.

Also called correlated sampling method, the reweighting method uses the output of a single random walk to estimate the energy as the function of variational parameters. One has to store the set of configurations {R (1) , R (2) , ..., R (M ) }, sampled from initial distribution Π 0 (R) = |Ψ T (R; a 0 )| 2 . Then the energy for a given set of variational parameters a is written as:

where the weight of the n th configuration is

One then can use a minimiser to find the lowest variational energy or variance (zero variance principle Eq. 1.31) or, as it is usually done, their linear combination as a function of a.

The problem of this method is that if the overlap between current and initial distributions |Ψ T (R (n) ; a)| 2 and |Ψ T (R (n) ; a 0 )| 2 gets too small, the error of ĒV (a) will become large and optimisation becomes unstable. For that reason it is useful to monitor the effective number of configurations N ef f = ( n w n (a)) 2 / n w 2 n (a) and when it gets too small one can reset Π 0 (R) = |Ψ T (R (n) ; a * )| 2 , where a * is the current best estimate of the parameters, regenerate the set of configurations, and proceed.

It is worth mentioning that the optimisation of trial functions for many-body systems is

Appendix C

Unfolding the band structure

If in the simulation we are using the supercell, in order to determine the Bloch crystal momentum ideal structure in the primitive cell we will have to first unfold the bands of the supercell. Single electron wave function can be written using plane waves and a Bloch representation as the following:

where n = 1..N/2 and u n,θ is the periodic function, which has the periodicity of the simulation supercell. In order to unfold the bands of a supercell we can displace the electrons by the primitive cell vectors (R ijk = ia 1 + ja 2 + ka 3 ), the resulting phase will be the crystal momentum, k

In practice to determine the phase one can compute the ratio of the two wave functions respectively displaced by the vector of the unit cell:

then the phase k could be determined from p = |p|e ikR ijk . One should use at least three independent displacements.

147

Appendix D

Coupled electron-ion Monte Carlo

The Coupled Electron-Ion Monte Carlo (CEIMC) [START_REF] Pierleoni | The Coupled Electron-Ion Monte Carlo Method[END_REF] is based entirely on the Monte Carlo method, both for solving the electronic problem, and for sampling the ionic configuration space. In treating finite temperature ions which are coupled with ground state electrons this method relies on the Born-Oppenheimer approximation. The essence of this method is the following: at finite temperature the ionic degrees of freedom (either classical point particles or path integrals) are sampled using a Metropolis Monte Carlo based on the electronic energies computed during independent ground state quantum Monte Carlo calculations.

As discussed in section 4. where β is inverse temperature, R denotes the set of all nuclear positions and ρ(R, R; β) is density matrix. We can then write the density matrix as a chain of P density matrices at τ = β/P with each element being,

In CEIMC, however, to reach the convergence for smaller number of slices, an effective pair potential, V pair (R), called pair action is introduced. The lowest-order Trotter factorisation [START_REF] Trotter | On the product of semi-groups of operators[END_REF] for the difference E 0 (R) -V pair (R) can be then written as,

With the pair potential, the matrix element in eq. D.3 can be evaluated analytically,

D.1 The Penalty method

Consider ∆ and σ 2 ∆ being an average value and a variance of ∆ = -τ P s=1 E 0 (R (s) ) -E 0 (R (s) ). Since ∆ is a random variable with the probability distribution P noise (∆), the acceptance a(∆) is also a random variable, with average

We assume that P noise (∆) is a gaussian

which can be justified by using the central limit theorem if the energy difference is averaged over many samples and has a finite variance. Choosing as the acceptance probability and defining Γ(R, R ) = ln T (R ,R) T (R,R ) , the average acceptance A has the form,

×e -∆(R,R )-Γ(R,R ) , which solves the averaged detailed balance condition, as ∆(R, R ) = -∆(R , R) and Γ(R, R ) = -Γ(R , R). The effect of the statistical noise is increasing the rejection by a factor e -σ 2 2 , as we can see from eq. D.10. However, the variance σ is not known a priori, but is estimated during the run to compute the energy as well. Given n independent samples y i , one has for the estimators of ∆ and σ 2 ∆ = 1 n n i=1 y i (D.12)

(D.13)

The average acceptance in eq. D.8 is then

dχ 2 a(∆, χ 2 )P noise (∆, χ 2 ), (D. [START_REF] Ranga | Observation of the Wigner-Huntington transition to metallic hydrogen[END_REF]