# Spectral and magnetic properties of copper oxides: towards a theoretical description from first principles

Benjamin Bacq-Labreuil

► To cite this version:

Benjamin Bacq-Labreuil. Spectral and magnetic properties of copper oxides: towards a theoretical description from first principles. Strongly Correlated Electrons [cond-mat.str-el]. Institut Polytechnique de Paris, 2022. English. NNT: 2022IPPAX084. tel-04397326

HAL Id: tel-04397326

https://hal.science/tel-04397326

Submitted on 16 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
Spectral and magnetic properties of
copper oxides: towards a theoretical
description from first principles

Thèse de doctorat de l'Institut Polytechnique de Paris
préparée à l'École Polytechnique

École doctorale n°626 de l'Institut Polytechnique de Paris (EDIPP)
Spécialité de doctorat: Physique

Thèse présentée et soutenue à Palaiseau, le 23 septembre 2022, par

BENJAMIN BACQ-LABREUIL

Composition du Jury :

Alexander Lichtenstein
Professeur, Institute for Theoretical Physics, Universität
Hamburg
Président

Michele Casula
Directeur de Recherche, Institut de Minéralogie, de Physique
des Matériaux et de Cosmochimie, Sorbonne Université
Rapporteur

Simon Moser
Group Leader, Physikalisches Institut,
Julius-Maximilians-Universität Würzburg
Rapporteur

Marcello Civelli
Maître de Conférences, Laboratoire de Physique des Solides,
Université Paris-Saclay
Examinateur

Matteo d'Astuto
Chargé de Recherche, Institut Néel
Examinateur

Benjamin Lenz
Maître de Conférences, Institut de Minéralogie, de Physique
des Matériaux et de Cosmochimie, Sorbonne Université
Co-encadrant de thèse

Silke Biermann
Professeure, Centre de Physique Théorique, Ecole
Polytechnique
Directrice de thèse
<figure><img src="image_1.png" /></figure>
Spectral and magnetic properties of copper
oxides: towards a theoretical description from
first principles

Benjamin Bacq-Labreuil

Thèse de doctorat de l'Institut Polytechnique de Paris
préparée au Centre de Physique Théorique
à
L'ÉCOLE POLYTECHNIQUE
<figure><img src="image_2.png" /></figure>
# CONTENTS

<table><tr><td>I</td><td>RÉSUMÉ</td><td>1</td></tr><tr><td>II</td><td>INTRODUCTION</td><td>5</td></tr><tr><td>III</td><td>METHODS</td><td>15</td></tr><tr><td>1</td><td>ELECTRONIC STRUCTURE CALCULATIONS</td><td>17</td></tr><tr><td>1.1</td><td>Density Functional Theory</td><td>17</td></tr><tr><td>1.2</td><td>Mapping to a solvable non-interacting problem: the Kohn-Sham equations</td><td>19</td></tr><tr><td>1.3</td><td>Local Density Approximation</td><td>22</td></tr><tr><td>1.4</td><td>Local Density Approximation + U</td><td>23</td></tr><tr><td>2</td><td>MANY-BODY METHODS</td><td>27</td></tr><tr><td>2.1</td><td>From the band structure to a low-energy effective model</td><td>27</td></tr><tr><td>2.1.1</td><td>The Hubbard model</td><td>27</td></tr><tr><td>2.1.2</td><td>Parameterization</td><td>29</td></tr><tr><td>2.2</td><td>Solving the Hubbard model: Dynamical Mean-Field Theory</td><td>32</td></tr><tr><td>2.2.1</td><td>The limit of infinite dimension</td><td>32</td></tr><tr><td>2.2.2</td><td>The Cavity Method</td><td>36</td></tr><tr><td>2.3</td><td>Beyond local correlations: Cluster Dynamical Mean-Field Theory</td><td>42</td></tr><tr><td>2.4</td><td>Cluster Dynamical Mean Field Theory: Implementation Details</td><td>45</td></tr><tr><td>2.4.1</td><td>Fourier Transformation</td><td>46</td></tr><tr><td>2.4.2</td><td>Initialization and choice of chemical potential</td><td>48</td></tr><tr><td>2.4.3</td><td>Basis transformation</td><td>49</td></tr><tr><td>2.5</td><td>Analytic Continuation</td><td>50</td></tr><tr><td>2.5.1</td><td>Maximum Entropy Method</td><td>52</td></tr><tr><td>2.5.2</td><td>Dealing with matrix-valued functions: element-wise continuation</td><td>55</td></tr><tr><td>2.5.3</td><td>Dealing with matrix-valued functions: Maximum Quantum Entropy Method</td><td>57</td></tr><tr><td>2.6</td><td>An effective spin model for lightly doped cuprates</td><td>58</td></tr><tr><td>IV</td><td>STUDIES OF SELECTED COPPER OXIDES</td><td>65</td></tr><tr><td>3</td><td>STUDY OF THE ANTIFERROMAGNETIC FLUCTUATIONS IN CA<sub>2-x</sub>NA<sub>x</sub>CUO<sub>2</sub>CL<sub>2</sub></td><td>67</td></tr><tr><td>3.1</td><td>Introduction</td><td>67</td></tr></table>
Contents

<table><tbody><tr><td>3.2</td><td>Model and method</td><td>71</td></tr><tr><td>3.2.1</td><td>Band structure of Ca<sub>2</sub>CuO<sub>2</sub>Cl<sub>2</sub></td><td>71</td></tr><tr><td>3.2.2</td><td>Low-energy effective Hubbard model</td><td>72</td></tr><tr><td>3.2.3</td><td>Angle-Resolved Photoemission Spectroscopy measurements</td><td>74</td></tr><tr><td>3.3</td><td>Spectral function</td><td>75</td></tr><tr><td>3.3.1</td><td>Parent compound</td><td>75</td></tr><tr><td>3.3.2</td><td>Hole doped system</td><td>79</td></tr><tr><td>3.4</td><td>Cluster analysis</td><td>83</td></tr><tr><td>3.4.1</td><td>State histogram</td><td>84</td></tr><tr><td>3.4.2</td><td>Origin of the pseudogap</td><td>87</td></tr><tr><td>3.5</td><td>Dynamical Spin Structure Factor: a direct probe of the spin fluctuations</td><td>92</td></tr><tr><td>3.6</td><td>Conclusion</td><td>100</td></tr><tr><td>4</td><td>BULK PROPERTIES OF TETRAGONAL CuO</td><td>103</td></tr><tr><td>4.1</td><td>Introduction</td><td>103</td></tr><tr><td>4.2</td><td>Model and method</td><td>105</td></tr><tr><td>4.2.1</td><td>One-band effective model</td><td>106</td></tr><tr><td>4.2.2</td><td>Methodology</td><td>107</td></tr><tr><td>4.3</td><td>Sublattice decoupling</td><td>108</td></tr><tr><td>4.4</td><td>Spectral function</td><td>111</td></tr><tr><td>4.4.1</td><td>The block construction</td><td>111</td></tr><tr><td>4.4.2</td><td>Results</td><td>112</td></tr><tr><td>4.5</td><td>Temperature dependence</td><td>117</td></tr><tr><td>4.5.1</td><td>Results</td><td>117</td></tr><tr><td>4.5.2</td><td>Effective single site model</td><td>121</td></tr><tr><td>4.6</td><td>Superconductivity</td><td>125</td></tr><tr><td>4.7</td><td>Conclusion</td><td>126</td></tr><tr><td>5</td><td>INTERFACE PROPERTIES OF TETRAGONAL CuO</td><td>129</td></tr><tr><td>5.1</td><td>Puzzling spin and orbital moments</td><td>129</td></tr><tr><td>5.2</td><td>Model and method</td><td>131</td></tr><tr><td>5.2.1</td><td>Model</td><td>131</td></tr><tr><td>5.2.2</td><td>Methodology</td><td>131</td></tr><tr><td>5.3</td><td>Bare interface</td><td>132</td></tr><tr><td>5.4</td><td>Vacancy in the CuO interface layer</td><td>134</td></tr><tr><td>5.5</td><td>Vacancy in the TiO<sub>2</sub> interface layer</td><td>138</td></tr><tr><td>5.6</td><td>Two vacancies in the TiO<sub>2</sub> interface layer</td><td>145</td></tr><tr><td>5.7</td><td>Conclusion</td><td>147</td></tr><tr><td>V</td><td>CONCLUSION</td><td>149</td></tr></tbody></table>
<table><tr><td>VI APPENDIX</td><td>155</td></tr></table>

<table><tr><td>A PROBING THE ELECTRONIC STRUCTURE: ANGLE-RESOLVED PHOTOEMISSION SPECTROSCOPY</td><td>157</td></tr></table>

<table><tr><td>A.1 The photoemission cross-section and the Green's function</td><td>158</td></tr></table>

<table><tr><td>A.2 From angles to momentum</td><td>161</td></tr></table>

<table><tr><td>ACKNOWLEDGMENTS</td><td>163</td></tr></table>

<table><tr><td>ACRONYMS</td><td>165</td></tr></table>

<table><tr><td>BIBLIOGRAPHY</td><td>169</td></tr></table>
<figure><img src="image_3.png" /></figure>
# PART I

## RÉSUMÉ
<figure><img src="image_4.png" /></figure>
Si l'oxyde de cuivre peut être une malédiction pour les anciennes casseroles, il serait plutôt une bénédiction pour le monde scientifique. La découverte des cuprates en 1986 par Bednorz et Müller, une famille de matériaux à haute température critique ayant comme élément central les plans 2D CuO₂, a lancé un champ entier de recherche. D'immenses efforts ont été déployés, menant à la découverte d'un diagramme de phase bien plus riche que le seul dôme supraconducteur. Même s'il est établi que les plans 2D de CuO₂ sont responsables des propriétés électroniques à basses énergies, une compréhension globale théorique fait toujours défaut. La raison principale tient dans le fait qu'il est extrêmement difficile de modéliser les cuprates puisqu'ils font partie de la famille des matériaux à fortes corrélations électroniques. Dans ces systèmes, l'énergie d'intéraction de Coulomb est du même ordre de grandeur que l'énergie cinétique, ce qui rend ineffectives les théories de perturbation, et les approximations de couplage faible/fort. Le travail des théoriciens est donc de développer des modèles plus affinés, mais toujours solubles numériquement ou analytiquement.

Dans cette thèse, nous nous concentrons sur l'étude des propriétés spectrales et magnétiques de deux oxydes de cuivre : Ca₂₋ₓNaₓCuO₂Cl₂ (NaCCOC) et CuO tétragonal (t-CuO). Une grande partie de cette étude est faite en adoptant le modèle de Hubbard à une bande pour capturer l'essentiel de la physique des plans de CuO₂. Malgré son apparente simplicité, ce modèle ne peut être résolu, et nous utilisons donc la théorie du champ dynamique appliquée aux amas (CDMFT). En complément, nous mobilisons la théorie fonctionelle de la densité (DFT). En sus de cet approche théorique, nous avons fait deux campagnes de mesure de spectroscopie par photoémission résolue en angle aux synchrotrons SOLEIL (Saint-Aubain, France) et Photon Factory (Tsukuba, Japon).

Cette combinaison d'outils théoriques et experimentaux nous permet d'étudier en détail les propriétés spectrales de NaCCOC non-dopé, et dopé en trou. Ce matériaux a une structure cristalline simple, et n'est pas supraconducteur à faible dopages ($x < 0.1$, x étant le taux de trous), ce qui le rend intéressant pour l'étude de la phase pseudogap qui survit donc jusqu'à des températures abritrairement basses. Nos calculs et mesures révèlent que les fluctuations antiferromagnétiques de spin jouent un rôle prépondérant. En premier lieu pour le phénomène de "cascade" : nous montrons qu'il s'agit d'un croisement entre une physique de polaron de spin en haut de la bande de valence, et une physique d'électrons localisés à plus haute énergie. Cette explication, mise en évidence par l'utilisation d'un modèle effectif de spin résolu avec l'approximation auto-cohérente de Born pour les polarons de spin et l'approximation Hubbard-I pour les électrons localisés, reste valide à faible dopage en trou. Nous montrons également que les fluctuations de spins sont à l'origine d'interférences destructives supprimant les excitations à basse énergie au point anti-nodal dans la zone de Brillouin. Cela entraîne l'ouverture du pseudogap, que l'on observe également en photoémission. Enfin, nous observons qu'à basse température la fonction de corrélation de spin à temps long ne diminue pas en augmentant le dopage, mais au contraire augmente jusqu'à $x = 0.05$, stagne, puis diminue. Ce phénomène s'étend même aux corrélations de spin entre plus proches voisins, et prochains plus proche voisins.

t-CuO possède une structure cristalline encore plus simple que NaCCOC, puisqu'il est seulement composé de plans 2D CuO qui se superposent dans une symétrie tétragonale sur un substrat de SrTiO₃ (STO). Il a été observé expérimentalement que les plans CuO se décomposaient
en deux sous-réseaux CuO₂ faiblement couplés. En utilisant la CDMFT avec plusieurs résolveurs d'impurité complémentaires, nous expliquons ce phénomène comment étant la conséquence de l'intéraction de Coulomb locale. La suppression des corrélations entre plus proches voisins, i.e. entre sous-réseaux, entraîne de façon effective une renormalisation des termes de saut qui résulte en un fort découplage. La symétrie du paramètre d'ordre est directement impactée : nous prédisons qu'elle est du type $d_{xy}$, au lieu du typique $d_{x^2-y^2}$ pour les cuprates.

Enfin, des expériences ont montré récemment l'émergence de moments pamagnétiques de spin et de moments orbitaux fixes sur les atomes de Cu. Compte tenu des connaissances actuelles sur t-CuO, de tels moments sont surprenants. En utilisant la DFT sur de larges cellules unitaires, nous démontrons qu'il est possible d'expliquer les deux types de moment par la formation de lacunes d'oxygène à l'interface t-CuO/STO. En effet, si les lacunes se forment dans le plan interfacial TiO₂, un gaz d'électrons en 2D est généré dans le plan interfacial CuO. Ce gaz d'électrons serait à l'origine des moments paramagnétiques de spin. De plus, la formation de lacunes engendre de fortes distortions dans le plan CuO à l'interface, au point où la symétrie par rotation locale et globale sont brisée. Cela engendre la fixation du moment orbital dans le plan. Notre scenario étant capable d'expliquer l'observation expérimentale, il reste à le confronter à d'autres expériences supplémentaire pour attester de sa validité : nous proposons donc plusieurs expériences qui permettrait de le vérifier.
# PART II

## INTRODUCTION
<figure><img src="image_5.png" /></figure>
<figure><img src="image_6.png" /><figcaption>Figure 1: (a) Sketch of a typical cuprate's phase diagram as a function of hole doping x and temperature T. The acronym AF refers to an antiferromagnetic charge transfer insulator phase. (b)-(c) Crystal structures of the undoped (x = 0) NaCCOC and t-CuO respectively.</figcaption></figure>

Copper oxide systems belong to the fascinating family of strongly correlated solids. These materials display extremely rich physics, providing a wide range of properties and applications. Some are good metals, other large-gap insulators, they can also be strong ferromagnets, or show antiferromagnetic ordering as for the undoped cuprates. All these different properties originate from a common ground: the electron-electron interactions are strong and can not be approximated as is usually the case for metals. Worse than that, the kinetic energy can neither be neglected since usually the Coulomb interaction is of the order of the bandwidth. In such a case, standard mean-field theories, or perturbative expansions either in kinetic or interaction energy fail to describe properly those materials. The strong competition between the two terms leads to the emergence of unique properties like giant magnetoresistance, ferroelectricity, and high temperature superconductivity to name a few, which are promising and even already used for practical applications.

Cuprates, discovered by Bednorz and Müller in 1986 [1], are best known for their high critical temperature under which superconductivity emerges, the highest at ambient pressure for Hg-based compounds [2], which sparked an immense scientific effort in the hope of finding the underlying mechanism. Indeed, the standard phonon-mediated Bardeen-Cooper-Schrieffer [3] pairing of the electrons is not compatible with the experimental evidences [4, 5]. It has been observed that the superconducting gap has a d-wave symmetry [4, 6], and although phonons can help the pairing caused by another mechanism [7], it can not alone account for this symmetry. Moreover, since in cuprates both the electron-phonon coupling and the phonon frequency responsible for the pairing are weak, the critical temperature derived from this theory could not reach the measured values in cuprates. While the coupling constant could be increased by enhancing the density of states at the Fermi level, this would ultimately lead to a renormalization of the phonons due to screening by the electronic cloud. The predominance of antiferromagnetic fluctuations, as we shall describe hereafter, point to a possible origin of superconductivity from antiferromagnetic correlations, but the mechanism remains mysterious [8–13]. The study is further complicated by the discovery of a complex phase diagram at low-temperature upon hole or electron doping. A typical example for the hole-doped compounds is sketched in Fig. 1(a) [9, 14].
While many different cuprates exist, they all have in common the 2D $CuO_2$ layers responsible for the low-energy properties, eventually sandwiched between the so-called charge reservoir layers (inert electronically), as shown in Fig. 1(b)-(c) where are displayed the crystal structures of $Ca_{2-x}Na_xCuO_2Cl_2$ (NaCCOC) and tetragonal CuO (t-CuO), the two materials on which we focus in this thesis. In the parent compounds, the Cu atoms are in a $3d^9$ configuration, i.e. with one partially filled orbital, while the oxygen atoms are in a $2p^6$ configuration. The half-filled orbital is usually of $d_{x^2-y^2}$ symmetry, but it can also be $d_{xy}$ as is the case for t-CuO. Within a band structure picture, the presence of a half-filled band would mean that the system is metallic. However, due to the large effective Coulomb interaction on these narrow orbitals, it is energetically favorable for the electrons to localize [15, 16]: the parent compounds are charge transfer insulators. The half-filled band is splitted into a lower Hubbard band (LHB) and upper Hubbard band (UHB), and the splitting (Mott gap) is such that the lowest energy excitation is the transfer of an electron from an O site to a Cu site, hence the name charge transfer. The gap is large: usually of the order of 2eV or larger [9]. Even though the on-site Coulomb interaction is large, the hopping of an electron between different Cu sites, mediated by the O atoms, is non-zero, and lead to the coupling of the single electrons on each Cu site via an antiferromagnetic super-exchange term [17, 18]. Undoped cuprates are therefore antiferromagnetically ordered below the Néel temperature $T_N$ [19, 20].

Upon hole doping of the $CuO_2$ layers, i.e. the insertion of holes, the electronic properties of the materials change drastically. Doping is usually achieved by chemical substitution in the charge reservoir layers that separate the 2D $CuO_2$ layers: for instance in NaCCOC, by the substitution of the Ca atoms by Na. At low temperatures, when increasing the hole doping one enters the $d$-wave [4, 6] superconductivity double dome, which vanishes at large dopings where the system becomes a Fermi liquid. Above the second dome, at temperature $T \ge T_c$ is observed the strange metal phase, characterized by a $T$-linear resistivity [14]. Finally, at moderate dopings and $T_c \le T \le T^*$ appears the pseudogap (PG) phase [21–23] characterized in the spectral function by the opening of an anisotropic gap in the Brillouin zone (BZ) [12, 24], as well as charge density wave (CDW) at lower temperatures [25–29].

The complex electronic properties of cuprates, highlighted by their phase diagram in which many different phases are interlaced, is extremely difficult to model theoretically. Indeed, even the undoped antiferromagnetic insulating phase can not be accounted for by the standard density functional theory (DFT), which within the local density approximation (LDA) would predict the system to be metallic since it has a half-filled band. Band structure calculations nevertheless provide a key first step in allowing to identify the partially occupied bands relevant for the description of the low-energy electronic properties [30, 31]. In cuprates, we usually expect a single band close to the Fermi level since only the Cu atoms have a partially filled shell. It is the result of the hybridization of the Cu $d_{x^2-y^2}$ and O $p_{x/y}$ orbitals which are coupled by the hopping amplitude $t_{dp}$. As mentioned above, a simplified view of this problem is that the Cu $d_{x^2-y^2}$ host a hole lying at an energy $E_d$, while the O $p_{x/y}$ is fully occupied and lies at the energy $E_p$, higher than $E_d$. The
lowest energy excitation process is the hopping of the hole from the Cu to O with amplitude $t_{dp}$, but which has an extremely large cost $E_p - E_d$ since the energy of double occupancy is $E_d \sim U_d$, the Coulomb on-site repulsion. The latter being large compared to the hopping amplitude $t_{dp}$, the hole remains localized on the Cu site and the system has a large charge transfer gap [32].

From these considerations, the three-band model, called Emery model, taking into account the Cu $d_{x^2-y^2}$ and O $p_{x/y}$ orbitals seems to be the minimal model possible for cuprates [33, 34]. However, a further simplification can be done by considering a single-band Hubbard model [35–37] at half-filling, with the on-site Coulomb repulsion $U$ being identified to $E_p - E_d$. Similarly to the three-band model, at large $U$ the electrons are localized and a Mott gap opens [15, 16]. In fact, the simplification of the three-band model to a single-band one was first argued by Zhang and Rice, who proposed that upon doping the extra holes form a bound state between the Cu $d_{x^2-y^2}$ and the surrounding O $p_{x/y}$ orbitals, leading to the formation of the isolated Zhang-Rice singlets (ZRS) band at the Fermi level [38]. Hence, the low-energy properties of cuprates can be ultimately represented in a minimal 2D single-band Hubbard model, which we adopt in this thesis for the study of NaCCOC and t-CuO. The one-band Hubbard model captures most of the essential phases of the phase diagram sketched in Fig. 1 when parameterized for cuprates [39]. Still, a debate remains open on whether the one-band model can indeed capture the right physics [40, 41], and if the phenomena observed theoretically are indeed relevant for the real materials.

Although the minimal one-band Hubbard model drastically simplifies the problem, unless for specific cases like the non-interacting and atomic limit, or in 1D, it can not be solved exactly due to the competition between the kinetic and the electron-electron interaction terms. The Hamiltonian reads

$$H = \sum_{\langle i,j \rangle, \sigma} t_{ij} c_{i\sigma}^{\dagger} c_{j\sigma} + U \sum_{i} n_{i\uparrow} n_{i\downarrow}, \qquad (0.1)$$

where $t_{ij}$ are the hopping amplitudes between neighbors $\langle ij \rangle$, and $U$ is the on-site Coulomb interaction. At parameters relevant for the cuprates, the bandwidth is roughly $W \approx 8t$, with $t$ the nearest-neighbor (NN) hopping term, and is of the order of the Coulomb repulsion $U \approx 8 - 10t$ in general [31]. This prevents the use of perturbation theory, as well as weak- or strong-coupling approximations, such that even at half-filling the problem is not tractable analytically. In this context, the development of dynamical mean-field theory (DMFT) [42, 43], sparked by considerations on the properties of the Hubbard model at infinite dimension [44–47], was a huge step forward in the theoretical investigation of the Hubbard model, and hence of cuprates. Based on the mapping of the full interacting lattice problem onto an embedded impurity problem, which can be efficiently solved numerically, DMFT incorporates exactly the local correlations, and is able to capture the Mott transition [43].

Although extremely powerful, DMFT suffers from the approximation of purely local correlations for systems like cuprates, and in fact already at half-filling it overestimates the critical interaction strength $U_c$ at which the Mott transition occurs [48]. A natural extension of DMFT is cluster dynamical mean-field theory (CDMFT), which we will use in this thesis. Within this framework, the single site impurity is replaced by a finite size cluster of sites [49–52], so that the non-local
correlations up to the length scale of the cluster are treated exactly in the impurity problem. The latter is solved by the so-called impurity solvers, a large variety of which are nowadays available. Some examples are exact diagonalization methods [53, 54], renormalization group techniques [55], matrix-product state (MPS) [56, 57]-based methods [58, 59], and continuous-time quantum Monte Carlo (CTQMC) solvers [60–62]. CDMFT belongs to the family of quantum cluster methods [52], which are either used in a self-consistent DMFT scheme like CDMFT and dynamical cluster approximation (DCA) [63–65], or not, like the cluster perturbation theory [66, 67] and the variational cluster approximation (VCA) [68–70]. Similarly to CDMFT, although some conceptual differences, all these techniques aim at incorporating a part of the non-local fluctuations by considering a cluster defined in real or momentum space. Yet another route, not used in this thesis, is to go beyond the purely local correlations by incorporating the non-local fluctuations using diagrammatic expansions around the DMFT solution [71–76].

In addition to the theoretical difficulties in building and solving cuprates' effective models, the comparison to experiments is further complicated by possible effects from the charge reservoir layers, which may shed some ambiguity on the intrinsic properties of the CuO₂ planes. A typical example is the observation in *n*-layer Hg based cuprates of the coexistence of CDW and antiferromagnetic ordering in different CuO₂ layers [77, 78]. These materials are made up of groups of *n* CuO₂ planes, separated by charge reservoir layers. For *n* ≥ 3 the CuO₂ planes can be separated into two groups: the outer planes directly connected to the charge reservoirs which can face large structural distortions, and the inner planes protected from such deformations. While in the former the antiferromagnetic order is quickly suppressed at low doping, the inner planes show a long-range ordering until a relatively large doping *x* = 0.08 [77]. Moreover, it was interpreted recently from quantum oscillations measurements that there could be a coexistence of CDW in the outer planes and antiferromagnetic ordering in inner planes, as signaled by the presence of three different frequencies in the oscillation spectrum [78]. This inner/outer plane differentiation also has strong consequences on the spectral function, as was shown by angle-resolved photoemission spectroscopy (ARPES) measurements on the same family of compounds, where was found a Fermi arc coexisting with a Fermi pocket at the nodal point in the PG phase [79]. Hence, in order to connect the theoretical models to the experiment, it is convenient to study cuprates with simple crystal structures, and even compounds free from charge reservoir layers.

In this thesis, we are particularly interested in NaCCOC and t-CuO, two compounds whose crystal structures are shown in Fig. 1(b)-(c). NaCCOC is composed of single CuO₂ planes separated by Cl and Ca atoms, with the particularity that the apical O atoms are replaced by Cl [80–83]. It sparked some attention when discovered since this was one of the first examples of superconducting samples without oxygen octahedra [80, 81]. Strictly speaking, NaCCOC is not a copper oxide, it is an oxychloride, but the similarity of its structure and its phase diagram [84] to those of usual cuprates makes it a promising candidate to confront the theoretical models to experiments. The parent compound [82] is a charge transfer insulator with a gap 0.5 ≤ Δ ≤ 2.5eV [85–89], and long-range magnetic order [84]. A kink has been observed in the low-energy dispersion [86], reminiscent of the cuprates' *waterfall* feature [24], and which survives at low to intermediate
dopings [90–92]. Upon hole doping, NaCCOC first enters the PG phase at arbitrarily low temperature [26, 90–94], and becomes superconducting for dopings $x \ge 0.1$ [84].

The waterfall feature appears at the nodal point due to the crossover between two different physical regimes at low and high energy. The low-energy part is well understood in terms of spin-polaron, that is to say carriers heavily dressed by the antiferromagnetic spin fluctuations [95–98]. However, the high-energy part has been attributed to different possible mechanisms like spin-charge separation [99, 100], weak-coupling spin-density wave, or to long-range hopping processes [101, 102] within a single magnetic sublattice [98]. The PG is characterized in the spectral function by the opening of an anisotropic gap in the BZ and the appearance of Fermi arcs (or pockets [79]) at the nodal points [12, 24]. It was in fact first discovered in Knight shifts [21] and spin-lattice relaxation time [22] measurements, and slightly later in transport experiments [23]. It has been shown with embedded quantum cluster methods like CDMFT [103–106], and DCA [65, 107–110], that the 2D one-band Hubbard model captures the PG upon hole-doping. Although it remained mysterious for quite some time, a theoretical consensus seems to be reached nowadays that the PG is intimately related to the antiferromagnetic correlations. First signs were obtained by quantum Monte Carlo (QMC) calculations [111], confirmed later on by quantum cluster methods [104, 107, 108] using relatively small clusters, as well as larger clusters [110]. The development of the fluctuation diagnostic reinforced the antiferromagnetic fluctuations scenario, by showing that the dominant channel contributing to the self-energy was the spin channel at $Q = (\pi, \pi)$ [112–114], and was further confirmed by a dual boson expansion around the DMFT solution [75]. The similarities of the theoretical results to the experiments, for instance for the spectral function or the phase diagram, support this interpretation [115].

We will focus on the effects of the antiferromagnetic spin fluctuations on the spectral properties of NaCCOC by simultaneously studying the waterfall feature and the PG. CDMFT [49–52] calculations are performed with 2 × 2 and 8-site clusters, using the hybridization-expansion CT-HYB [116] and the interaction-expansion CT-INT [117] solvers based on the ALPSCore libraries [118], as well as the CT-HYB [119] solver from the TRIQS library [120]. Along with the theoretical calculations, ARPES measurements were carried out at the Photon Factory (Tsukuba, Japan) and SOLEIL (Saint-Aubin, France) synchrotrons to which the author actively participated. These measurements were performed in close collaboration with Matteo d'Astuto, Chafic Fawaz, and Hervé Cercellier from the Institut Néel (Grenoble, France), as well as Yukiko Obata from the Tokyo Tech World Research Hub Initiative (Tokyo, Japan).

We find that at all dopings considered, the one-band Hubbard model solved with a 2 × 2 and 8-site cluster is in astonishing agreement with the experimental measurements. The kink in the dispersion is particularly well reproduced, both in terms of momentum and energy position. Consistently with the existing literature, we identify the low-energy part of the waterfall as stemming from a spin-polaron by showing that the ARPES and CDMFT spectral functions are well accounted for by the propagation of a hole in an antiferromagnetic background, the dispersion of which is calculated within the self-consistent Born approximation (SCBA) [95]. Moreover,
the high-energy part can be unambiguously described by assuming purely local correlations, i.e. within the Hubbard-I approximation, which is consistent with the vision of heavily dressed carriers in a background of localized moments. Furthermore, by studying the state histogram of the simple 2 × 2 cluster, it appears that the strong antiferromagnetic correlations are intimately linked with the opening of the PG, for it promotes as the leading term of the histogram a non-degenerate singlet state. We show, following Ref. 121, that such state leads to the suppression of the low-energy photoemission processes and thus to the opening of the PG at the antinodes, and argue that the picture is not limited to the specific 2 × 2 cluster. Finally, we calculate the dynamical spin structure factor (DSSF) for the 2 × 2 cluster, and find that it is maximum at $Q = (\pi, \pi)$ for all dopings until $x = 0.1$, confirming the role of the antiferromagnetic fluctuations. We also identify a freezing of the inter-site spin-spin correlation function at large times $\tau \to \beta/2$ and at finite dopings, in agreement with recent findings [122, 123], and show that such freezing can be understood as a consequence of the emergence of low-energy spin excitations upon hole-doping. Most interestingly, we argue that such freezing may extend to arbitrarily large inter-site distances for larger systems, given a finite momentum broadening of the $Q = (\pi, \pi)$ peak.

In contrast to NaCCOC, t-CuO has an even simpler crystal structure: it is composed of CuO planes (not CuO$_2$!), which are stacked along the *c* axis with a tetragonal distortion characterized by the ratio of basal and apical distances *c*/*a* ≈ 1.37 [124, 125]. In fact, t-CuO has to be synthesized as a thin film on a TiO$_2$-terminated *SrTiO$_3$* (STO) substrate, and is stable until approximately 10 unit cells of CuO in the *c* direction [124, 126]. Above this critical thickness, CuO returns to its bulk low-symmetry monoclinic phase [127], in which well defined 2D CuO planes are missing. t-CuO has not been experimentally doped yet: the lack of charge reservoir layers prevents the use of the chemical substitution.

Ab initio studies including DFT with hybrid functionals [128–131] and DFT+U [132, 133] gave first insights into the electronic structure of t-CuO and were able to reproduce the experimentally observed tetragonal distortion [129], which is traced back to Jahn-Teller orbital ordering at the Cu d⁹ ions [128, 131]. t-CuO is an insulator with a sizable gap $\Delta > 2.35$ eV [126, 134], accompanied by an antiferromagnetic stripe ordering as found in resonant inelastic x-ray scattering (RIXS) and muon spin resonance (μ-SR) measurements [135, 136], as well as first principle calculations [128, 129, 131]. From heuristic extrapolations based on other binary transition metal oxides (TMO) [124, 137], and DFT calculations [128, 129, 131], the Néel temperature is estimated to be around $T_N \approx 800$K, significantly higher than the critical temperature of its monoclinic counterpart ($T_N \sim 220$K [127]). Most interestingly, it was observed that the ARPES and RIXS spectra are consistent with an effective decoupling of the CuO planes into two weakly-correlated CuO$_2$ sublattices [134, 135]. However, within the simple one-band model derived from RIXS [135], the NN inter-sublattice hopping $t_d$ is evaluated to be about one-fourth of the next-NN (NNN) intra-sublattice hopping *t*. A ratio of one-fourth is clearly too high to explain the sublattice decoupling seen experimentally, given the strong replica features observed in ARPES and the similitude of the magnon spectrum with that of conventional cuprates. Since the model is constructed from a Heisenberg-like Hamiltonian, and therefore lacks the dynamical effects from electron-electron
interactions, this raises the question on the possible role of the dynamic correlations on the sublattice decoupling.

In close collaboration with Max Bramberger, Martin Grundner, Sebastian Paeckel and Ulrich Schollwöck from Ludwig-Maximilians-Universität München (Munich, Germany), we investigate the dynamical influence of the electronic correlations by means of CDMFT [49–52] calculations. Three different solvers were used: the CTQMC [60–62] hybridization-expansion CT-HYB [116] solver for which the author performed the calculations, and two MPS [56, 57]-based impurity solver working at effectively zero temperature and on the imaginary and the real axis [58, 59, 138–144] for which M. Bramberger and M. Grundner did the calculations. The inspection of the different self-energy components reveals a strong suppression of the inter-sublattice correlations, corresponding to the hopping $t_d$. This leads to an effective decrease of the non-interacting $\frac{t_d}{t} \sim \frac{1}{4}$ ratio, and therefore the enhancement of the sublattice decoupling which justifies the view of the CuO layers in t-CuO as two weakly coupled CuO$_2$ layers. Having revealed the underlying mechanism of the sublattice decoupling, we motivate the construction of a super-cluster allowing an improvement of the momentum resolution of our calculations without additional computational cost. This results in a spectral function in excellent agreement with the experimental one [134], whether using the CT-HYB solver or the MPS-based solver directly on the real axis. This *hand-shake* between finite and zero temperature techniques allows to study the influence of the antiferromagnetic spin fluctuations on the properties of t-CuO. Using the CT-HYB solver, we find two different insulating regimes: a paramagnetic insulator at high temperatures, characterized by a momentum selective divergence of the self-energy, and an antiferromagnetically ordered insulator at low temperature, for which the self-energy displays a spin dependent static shift of the order of $U$. This transition is accompanied by a drastic transformation of the spectral function, which agrees with the ARPES measurements and the zero temperature calculations only in the magnetically ordered regime, in which the low-energy dispersion is well accounted for by a SCBA calculation, thus further confirming the role the spin fluctuations. Finally, based on VCA [68–70] calculations performed by Benjamin Lenz, the sublattice decoupling is shown to survive upon hole doping and especially in the superconducting phase: it leads to a $d_{xy}$ superconducting order parameter, which coexists with the antiferromagnetic stripes, although the two phases are in competition [50, 145–148].

Along with the existing literature, our findings would lead to a coherent and complete picture of the electronic properties of t-CuO, if it was not for the recent puzzling observation of paramagnetic spin and pinned orbital moments in x-ray magnetic circular dichroïsm (XMCD) measurements at the Cu L$_{2,3}$ edge [136]. Since t-CuO is an antiferromagnetic insulator with a sizable gap, one would rather expect the XMCD spectrum to be zero in the range of magnetic fields applied (from -6.5 to 6.5T). A minimal model was proposed in Ref. 136, which captures the essential features of the observations, but which is not line with the previous theoretical and experimental works. Indeed, according to their model, t-CuO would be composed of CuO layers with in-plane ferromagnetic ordering, which are antiferromagnetically stacked along the c axis. The last layer would follow the external field orientation and lead to the appearance of the paramagnetic spin
moment, while the one-to-last layer would be pinned and uncompensated, hence giving rise to
the pinned moment. Instead, the presence of a 2D electron gas (2DEG) at the t-CuO/STO in-
terface would be a plausible scenario for the emergence of a paramagnetic spin moment, and
that would especially be compatible with the antiferromagnetic stripe ordering observed both in
experiments [135, 136] and ab initio calculations [128, 129, 131]. However, this possibility was
argued against in Ref. 136 because of the absence of a XMCD Ti³⁺ signal at the Ti L₂,₃ edge,
which was previously seen in LaAlO₃/STO and γ-Al₂O₃/STO heterostructures in the presence
of a 2DEG at the interface [149, 150], therefore leaving the question on the origin of these puz-
zling moments open.

By performing DFT+U [151–153] calculations on large unit cells, we show that a 2DEG emerges
at the t-CuO/STO interface upon the formation of oxygen vacancies. Most importantly, the 2DEG
is invisible at the Ti L<sub>2,3</sub> edge since the valence of the Ti atoms remains unchanged. We there-
fore show that the 2DEG scenario for the paramagnetic spin moments is in fact *not incompatible*
with the XMCD measurements of Ref. 136. Moreover, the formation of oxygen vacancies leads to
strong distortions, mainly located at the interface, which, in combination with the misfit of the
monoxide CuO and the dioxide TiO<sub>2</sub> layers, leads to the breaking of both the local and the global
C<sub>4</sub> rotational symmetries. This induces the pinning of the in-plane component of the orbital
moment, while the out-of plane remains proportional and collinear with the spin moments, in
excellent agreement with the XMCD measurements [136]. Our scenario therefore reconciles the
appearance of the paramagnetic spin and orbital pinned moments with the current knowledge
of t-CuO.

The remainder of the thesis is organized as follows. In Chap. 1, we first introduce DFT, and discuss to which extent it can be suited for the study of the strongly correlated systems' ground states properties. Then, in Chap. 2, we present how a low-energy effective model suitable for cuprates can be constructed starting from a DFT band structure, as well as a possible way of solving it: DMFT and CDMFT. The concept of analytic continuation is introduced, which is a necessary step to transform the Monte Carlo impurity solvers data from the imaginary to the real axis. We also present an effective spin model suited for lightly doped cuprates that can be solved within the SCBA, and which will prove to be useful for the interpretation of our CDMFT calculations. We proceed with the Chap. 3 dedicated to the study of NaCCOC, then Chap. 4 in which are discussed the bulk properties of t-CuO, and finally Chap. 5 focusing on the t-CuO/STO interface properties.
# PART III

## METHODS
<figure><img src="image_7.png" /></figure>
# 1 ELECTRONIC STRUCTURE CALCULATIONS

## 1.1 DENSITY FUNCTIONAL THEORY

Density functional theory (DFT) has imposed itself as one of the most widely used numerical methods to study the ground state properties of solids and molecules. One of its main strengths lies in the fact that almost any solid can in principle be investigated ab initio, i.e. simply from the knowledge of its crystal structure, and with nowadays computational technologies large systems of $10^2-10^3$ atoms can be considered. Another advantage of DFT is that a large range of properties can be determined from the knowledge of the ground state: structural parameters, magnetism, optical properties, etc. In this chapter, we aim at giving an overview of the main equations and concepts of DFT, mainly following the review by Cottenier [154] and Kohn's Nobel lecture [155].

A solid is composed of many electrons moving in a network of positively charged ions. The Hamiltonian of this quantum many-body system reads:

$$ \begin{aligned} \hat{H} &= \hat{T}_{ion} + \hat{T}_{e} + \hat{V}_{e-ion} + \hat{V}_{ion-ion} + \hat{V}_{e-e} \\ &= -\sum_{i} \frac{\nabla_{\mathbf{R}_{i}}^{2}}{2M_{i}} - \sum_{i} \frac{\nabla_{\mathbf{r}_{i}}^{2}}{2m_{e}} - \sum_{i,j} \frac{Z_{i}}{|\mathbf{R}_{i} - \mathbf{r}_{j}|} + \frac{1}{2} \sum_{i,j} \frac{Z_{i}Z_{j}}{|\mathbf{R}_{i} - \mathbf{R}_{j}|} + \frac{1}{2} \sum_{i,j} \frac{1}{|\mathbf{r}_{i} - \mathbf{r}_{j}|}, \end{aligned} \quad (1.1) $$

where $\mathbf{R}_i$, $M_i$ are respectively the ionic position and mass, $\mathbf{r}_i$, $m_e$ the electronic ones, $Z_i$ is the atomic number of the ions, and we used atomic units $e = \hbar = \frac{1}{4\pi\epsilon_0} = 1$. The two first terms of the Hamiltonian ($\hat{T}_{ion}$ and $\hat{T}_e$) are the kinetic energy of the ions and the electrons, they are single particle quantities. The last three terms are two-body operators ($\hat{V}_{e-ion}$, $\hat{V}_{ion-ion}$ and $\hat{V}_{e-e}$) which correspond to the Coulomb interaction between the particles. The Hamiltonian, in its exact form, is not solvable, so methods like DFT have been developed to obtain an approximate solution.

The first step is to notice that the mass of the ions is usually three orders of magnitude larger than that of the electrons $M_i \gg m_e$, such that their dynamic is much slower. The Born-Oppenheimer approximation, in its crudest form, consists in considering the ions as frozen at fixed positions and therefore neglecting their kinetic energy. They act as a positive charge background for the electrons: the ion-electron Coulomb interaction term $\hat{V}_{e-ion}$ can be considered as an external potential acting on the electrons, while $\hat{V}_{ion-ion}$ becomes a constant. Note that this assumption is justified for the compounds under considerations in this thesis which are composed of Cu, O atoms and eventually Ca, Cl and Na atoms, but it can break down for systems like H-based superconductors due to the small mass of the H atom [156]. This approximation results in a simpler Hamiltonian:

$$ \hat{H} = \hat{T}_e + \hat{V}_{e-e} + \hat{V}_{ext}, \quad (1.2) $$
with $\hat{V}_{ext}$ the external potential created by the ionic network. A key aspect of this Hamiltonian is that the two first terms describe a system of interacting electrons independently from the particular solid they belong to, they are universal. The material's specificity is entirely encoded in the external potential $\hat{V}_{ext}$ and the electron filling.

Although the exact Hamiltonian has been drastically simplified, Eq. 1.2 remains not solvable because of the two-body term $\hat{V}_{e-e}$ that couples the electrons all together. Indeed, the eigenfunctions $|\phi(\mathbf{r})\rangle$ would be obtained by Slater determinants [157], i.e. linear combinations of anti-symmetric products of N-electron wave functions. In order to draw a rough estimate of the problem's complexity, the N-electron ground state wave function can be viewed as depending on the total number of parameters M [155]:

$$M \approx p^{3N}, \qquad (1.3)$$

where $p$ is the number of parameters per variable for obtaining a desired accuracy, and at large N the number of variables for the N-body wave function is roughly $3N$ (in 3D space). Let us consider the simple case of $p = 3$. For a usual solid, $N \sim 10^{23}$ and obviously Eq. 1.2 can not be solved. A crystalline solid is periodic, so one could take into account only one unit-cell of the crystal, making the problem much easier. For instance, the parent compound NaCCOC ($x = 0$) that we study later in this thesis can be reduced to a unit cell containing a single atom of Ca, Cu, O and Cl atoms. This amounts to a system of 74 electrons, which would roughly require $M \approx 3^{3 \times 74} \approx 10^{105}$ parameters! This hand wavy argument is of course simplified, but it illustrates the need for tools to approximate the many-body problem.

One method to tackle this issue is DFT, which was formally proposed by Hohenberg and Kohn in 1964 [158]. They formulated two theorems:

1. The ground state density $n_0(\mathbf{r})$ of the many-body electronic system is uniquely determined by the external potential $\hat{V}_{ext}$. It follows that the ground state expectation value of any operator $\hat{O}$ is a unique functional of the ground state density

$$\langle\phi_0|\hat{O}|\phi_0\rangle = O[n_0] \qquad (1.4)$$

2. Taking $\hat{O}$ to be the Hamiltonian $\hat{H}$, the ground state total energy functional reads

$$\begin{aligned} E[n] &= \langle\phi|\hat{T}_e + \hat{V}_{e-e}|\phi\rangle + \langle\phi|\hat{V}_{ext}|\phi\rangle \\ &= F_{HK}[n] + \int n(\mathbf{r})V_{ext}(\mathbf{r})d\mathbf{r}, \end{aligned} \qquad (1.5)$$

where $F_{HK}[n]$ is the universal Hohenberg-Kohn density functional that contains the kinetic and potential energy of the electrons. The ground state energy functional $E[n]$ is minimized at the ground state density $n(\mathbf{r}) = n_0(\mathbf{r})$.

These theorems tell us that the knowledge of the ground state density $n_0(\mathbf{r})$ provides the material's specific external potential $\hat{V}_{ext}$. Since the rest of the Hamiltonian 1.2 is universal, the ground state density determines the full $\hat{H}$, as well as the number of electrons N through its normalization. Therefore, the ground state density contains implicitly the same amount of information
about the energy of the system than the ground state eigenfunction. This is a huge simplification since only a single function $n(\mathbf{r})$ has to be known. It is important to keep in mind that this is only true for the ground state: any other density will determine a different $\hat{V}_{\text{ext}}$ with no physical relevance for the system under investigation. Unfortunately, this step forward does not tell us how to obtain the ground state density $n_0(\mathbf{r})$. A way to its practical determination was given by Kohn and Sham [159] that we shall describe in the following section.

Before proceeding to the Kohn-Sham equations, it is useful to reformulate the second theorem following the framework of Levy [160]. The idea is to start from the variational principle on the wave functions, which gives the ground state energy as:

$$E_0 = \min_{\phi} \langle \phi | \hat{H} | \phi \rangle, \qquad (1.6)$$

where the $|\phi\rangle$ are chosen in a set of normalized wave functions. This minimization procedure can be separated into two steps:

(i) given a density $n(\mathbf{r})$, minimize over all the wave functions that give the density $n(\mathbf{r})$,

(ii) minimize over the density functions $n(\mathbf{r})$.

The energy functional therefore reads:

$$\begin{aligned} E[n] &= \min_{\phi(n)} \langle \phi | \hat{H} | \phi \rangle \\ &= \min_{|\phi\rangle(n)} \langle \phi | \hat{T}_e + \hat{V}_{e-e} | \phi \rangle + \int n(\mathbf{r}) V_{\text{ext}}(\mathbf{r}) d\mathbf{r}, \end{aligned} \qquad (1.7)$$

where the minimization is done over a set of wave functions $|\phi\rangle$ corresponding to the density $n$. Note that the second term involving an integration over the external potential does not need to be minimized with respect to the set of wave functions since they all correspond to the same density $n$. The ground state energy is minimal when $n$ is the ground state density since:

$$E_0 = \min_n \left[ \min_{|\phi\rangle(n)} \langle \phi | \hat{H} | \phi \rangle \right] = \min_n E[n]. \qquad (1.8)$$

The minimization has to be done under the constraint $\int n(\mathbf{r})d\mathbf{r} = N$. Introducing a Lagrange multiplier $\lambda$, we obtain the stationarity condition:

$$\frac{\partial T_e}{\partial n}(\mathbf{r}) + \frac{\partial V_{e-e}}{\partial n}(\mathbf{r}) + V_{\text{ext}}(\mathbf{r}) = \lambda, \qquad (1.9)$$

which tells us that $\hat{V}_{\text{ext}}$ is uniquely determined by the density $n$.

## 1.2 MAPPING TO A SOLVABLE NON-INTERACTING PROBLEM: THE KOHN-SHAM EQUATIONS

The key concept of the Kohn-Sham equations is to map the many-body interacting system onto a fictitious *solvable non-interacting* system which mimics the ground state density $n_0$ [159]. Con-
sider a non-interacting system, in which the electrons are moving inside an effective external field:

$$ \hat{H}^0 = \hat{T}^0 + \hat{V}_{eff}^0, \qquad (1.10) $$

where the superscript 0 denotes the non-interacting quantities. This fictitious system can follow the same minimization procedure as the "true" system, and we can therefore write the following stationarity condition:

$$ \frac{\partial T^0}{\partial n}(\mathbf{r}) + V_{eff}^0(\mathbf{r}) = \lambda', \qquad (1.11) $$

where $\hat{V}_{eff}^0$ is the external potential for the given density $n$. If this stationarity condition was to be equivalent to that of the interacting system 1.9, we then get

$$ \begin{aligned} V_{eff}^0(\mathbf{r}) &= \underbrace{\left[ \frac{\partial T_e}{\partial n}(\mathbf{r}) - \frac{\partial T^0}{\partial n}(\mathbf{r}) \right] + \left[ \frac{\partial V_{e-e}}{\partial n}(\mathbf{r}) - \frac{\partial V_H}{\partial n}(\mathbf{r}) \right] + \frac{\partial V_H}{\partial n}(\mathbf{r})}_{\frac{\partial V_{xc}}{\partial n}(\mathbf{r})} + V_{ext}(\mathbf{r}) \\ &= \frac{\partial V_{xc}}{\partial n}(\mathbf{r}) + \frac{\partial V_H}{\partial n}(\mathbf{r}) + V_{ext}(\mathbf{r}) \end{aligned} \qquad (1.12) $$

where the constant difference between $\lambda$ and $\lambda'$ can be absorbed in the external potential. $\frac{\partial V_{xc}}{\partial n}(\mathbf{r})$ is the so-called exchange-correlation potential, and $V_H$ is the Hartree term defined by:

$$ V_H[n] = \frac{1}{2} \int d\mathbf{r} \int d\mathbf{r}' \frac{n(\mathbf{r})n(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|}. \qquad (1.13) $$

Now, the non-interacting Hamiltonian, which shall be called the Kohn-Sham Hamiltonian hereafter, can be written:

$$ (\hat{H}^0 =) \hat{H}_{KS} = -\sum_i \frac{\nabla_i^2}{2m_e} + \int d\mathbf{r}' \frac{n(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} + \hat{V}_{xc} + \hat{V}_{ext} \qquad (1.14) $$

where the exchange correlation potential is given by:

$$ \hat{V}_{xc} = \frac{\partial V_{xc}}{\partial n}. \qquad (1.15) $$

Therefore, the many-body interacting problem can be mapped onto the non-interacting effective Kohn-Sham Hamiltonian, which consists of non-interacting electrons subject to two potentials: the material's specific $\hat{V}_{ext}$, and the effective exchange-correlation potential $\hat{V}_{xc}$. $\hat{V}_{xc}$ contains all the involved contributions from the electron-electron interaction that are not included in the simple Hartree term. If in principle the mapping is exact, in general $\hat{V}_{xc}$ is unknown and needs to be approximated.

Assuming for a moment that $\hat{V}_{xc}$ is known, the initial variational problem of finding the ground state density of the fully interacting Hamiltonian is reformulated into solving the effective Kohn-Sham Hamiltonian for its ground state. Since it is a non-interacting Hamiltonian, its eigenfunctions can be obtained from:

$$ \hat{H}_{KS}|\phi_i\rangle = \epsilon_i|\phi_i\rangle \qquad (1.16) $$
<figure><img src="image_8.png" /><figcaption>Figure 1.1: Sketch of the DFT self-consistent cycle at the iteration j. When the desired accuracy $\delta_{err}$ is reached, the procedure is stopped and the ground state density $n_0$ can be extracted.</figcaption></figure>

where the $|\phi_i\rangle$ are the Kohn-Sham wave functions. The exact ground state density is then given by:

$$n_0(\mathbf{r}) = \sum_{i}^{N} \phi_i^*(\mathbf{r})\phi_i(\mathbf{r}). \qquad (1.17)$$

This is a significant improvement, since instead of solving the interacting problem we can solve Eq. 1.16. However, one has to note that the Kohn-Sham wave functions are not the eigenfunctions of the original interacting problem. Only the density $n_0(\mathbf{r})$ calculated from these functions using Eq. 1.17 has a physical meaning. Similarly, the Kohn-Sham energies $\epsilon_i$ are not the one of the physical system, although some of them can approximately be identified to the ionization energies [161, 162].

One of the prices to pay for this construction, besides the determination of the exchange correlation that we still assume to be known for now, is that the latter $\hat{V}_{xc}$ and the Hartree potential $\hat{V}_H$ depend on the density $n(\mathbf{r})$. Hence, the Kohn-Sham Hamiltonian that we solve for its density $n(\mathbf{r})$, depends on this very density $n(\mathbf{r})$. We face here a self-consistent problem, that can be tackled by the following procedure (sketched in Fig. 1.1):

(i) Start from an initial guess $n^{j=0}(\mathbf{r})$ (not to be mistaken with the ground state density $n_0(\mathbf{r})$),

(ii) Determine the corresponding Hartree and exchange correlation terms $\hat{V}_H$ and $\hat{V}_{xc}$,

(iii) Solve the Kohn-Sham Hamiltonian 1.16 for its eigenfunctions $|\phi_i\rangle$

(iv) Compute the new density $n^{j+1}(\mathbf{r})$ using Eq. 1.17

(v) if $n^{j+1}(\mathbf{r}) \approx n^j(\mathbf{r})$ up to the desired accuracy, stop the procedure and the ground state density is $n_0(\mathbf{r}) = n^{j+1}(\mathbf{r})$. If not, go back to (ii) and iterate.

Up to now, as mentioned above the procedure is exact (apart from the Born-Oppenheimer approx-
imation). In fact, the interactions were hidden under the carpet since all the unknown terms were put into the exchange correlation potential. For practical calculations however, approximations of the exchange-correlation are needed.

## 1.3 LOCAL DENSITY APPROXIMATION

One of the most popular approximation to the exchange-correlation functional is the so-called local density approximation (LDA), which was already proposed by Kohn and Sham [159]. It is based on the observation that if the density of a system is sufficiently slowly varying, then the exchange correlation function can be written:

$$V_{xc}^{LDA}[n] = \int n(\mathbf{r})\epsilon_{xc}(n(\mathbf{r}))d\mathbf{r}, \qquad (1.18)$$

where $\epsilon_{xc}(n(\mathbf{r}))$ is the exchange-correlation energy of the homogeneous electron gas at a particular density $n(\mathbf{r})$, which is known numerically. The main idea behind is the following: if the density $n(\mathbf{r})$ is slowly varying, then it can be divided into infinitesimal volumes where the density is locally constant. In each of these small volumes, at position $\mathbf{r}$, the electron gas can be considered as homogeneous, and therefore the exchange-correlation energy in each individual volume can be approximated by the well known exchange-correlation energy of an homogeneous gas of density equal to $n(\mathbf{r})$ (the local density of the true system). The total exchange-correlation functional is then obtained by integrating over all these infinitesimal volumes, yielding Eq. 1.18.

Within LDA, $V_{xc}[n]$ is purely local since the contribution to the functional at each point $\mathbf{r}$ is given by $\epsilon_{xc}$ evaluated at $n(\mathbf{r})$, and which is independent from the surrounding density. We expect this approximation to behave well for solids with delocalized electrons, like metals, for the wave functions in these systems extend over large distances thus leading to smoothly varying densities. In fact, it has proven to be efficient for many other systems beyond metals. It is for this reason that sometimes the Kohn-Sham energy bands are directly compared to ARPES measurements for instance, as will be done later in this thesis: for weakly correlated orbitals, the correction from the exchange-correlation potential is small and the electronic properties are rather lead by the crystalline potential which is accurately taken into account within DFT.

LDA can be extended to the spin polarized case with the local spin density approximation (LSDA), which can be written:

$$V_{xc}^{LSDA}[n] = \int n(\mathbf{r})\epsilon_{xc}(n^{\uparrow}(\mathbf{r}), n^{\downarrow}(\mathbf{r}))d\mathbf{r}, \qquad (1.19)$$

where $n^{\uparrow/\downarrow}$ is the density for down/up spins, and $\epsilon_{xc}$ is now the exchange-correlation energy of the homogeneous spin polarized gas.

However, for materials like TMOs, such as NaCCOC or t-CuO that we study in this thesis, for which the valence electrons tend to be rather localized on atomic-like orbitals, the assumption of slowly varying density breaks down. While non-magnetic LDA predicts NaCCOC and t-CuO to be metallic (see Fig. 3.1(c) and Fig. 4.2), experiments show that they are insulators with a quite
sizable gap [85, 124]. The electronic correlations are responsible for the localization of the electrons in such system: a more refined treatment is needed.

Many improvements to the local approximation have been formulated, a natural one being the generalised gradient approximation (GGA) in which the dependence to the gradient of the density is incorporated [163, 164], and which remains computationally cheap. Another possibility, more computationally demanding, is to include explicitly a part of the exact Hartree-Fock exchange

$$V_{HF} = -\frac{1}{2} \sum_{i,j} f_i f_j \int d\mathbf{r} \int d\mathbf{r}' \frac{\phi_j^*(\mathbf{r})\phi_i(\mathbf{r})\phi_i^*(\mathbf{r}')\phi_j(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|} \quad (1.20)$$

which depends on the Kohn-Sham orbitals, i.e. is indirectly linked to the density. This class of approximation are called hybrid functionals [165], and they differ by the amount of mixing between the exact exchange and the LDA/GGA exchange correlation energy. We do not aim here at listing all the different approximations available which have been more or less successfully applied to real materials studies. We will instead focus on one of these many extensions of LDA in the next section, called LDA+U [151–153], that we have used in this thesis to study the interface properties of t-CuO.

## 1.4 LOCAL DENSITY APPROXIMATION + U

TMOs and rare-earth compounds showing the Mott-insulator transition usually contain partially filled d or f orbitals. The Coulomb interaction is effectively much stronger in these narrow orbitals as compared to the other shells since the electron-electron distance is shorter, and the inter-site overlaps smaller. Such materials can be insulating although they have a partially filled shell, in contrast to the usual band insulators. They are called Mott insulators [15, 16]. Not only LSDA usually predicts these materials to be metallic, but it also fails to capture the magnetic ground state. For instance for t-CuO, or NaCCOC ($x = 0$), the magnetic order vanishes in LSDA while there is experimental evidence for long-range magnetic order [136, 166]. The failing of LSDA for these system can be traced back to two main reasons [151, 167–170]:

1. the LSDA potential is orbital independent since it is built on the homogeneous gas exchange-correlation potential.

2. In LSDA, the splitting between two different spin states is governed by the homogeneous gas exchange parameter J, while in the Mott insulators it should be lead by the Coulomb interaction $U \gg J$ as will be shown in the following [151].

In 1991, Anisimov et al. proposed in Ref. 151 to overcome the LSDA shortcomings by adding an effective potential based on the mean-field solution of the Hubbard model [35–37], which is a simple model known to well account for strongly correlated systems [171]. Its simplest version consists of electrons moving on an effective single-orbital lattice:

$$H = \sum_{\langle i,j \rangle, \sigma} t_{ij} c_{i\sigma}^{\dagger} c_{j\sigma} + U \sum_{i} n_{i\uparrow} n_{i\downarrow}, \quad (1.21)$$
where $c_{i\sigma}^{\dagger}$, $c_{i\sigma}$ are the creation and annihilation operators of an electron on site $i$ with spin $\sigma$, $n_{i\sigma}$ is the occupation operator, $t_{ij}$ is the hopping amplitude from site $j$ to $i$ which is summed over the pairs of neighboring sites $\langle i, j \rangle$, and $U$ is the effective on-site Coulomb interaction.

If $t_{ij} \gg U$, the electrons delocalize and the ground state energy dispersion is given by the Fourier transformation of the hoppings $t_{ij}$. In such a case, LSDA is a good approximation since the electronic wave functions are delocalized and the limit of the free electron gas is recovered. If $U \gg t_{ij}$ however (atomic limit), at half-filling, the electrons localize such that each orbital is singly occupied. The hopping processes are effectively suppressed because of the strong energy barrier $U$ which makes the double occupation energetically too costly. This situation is similar to that of TMOs: the electronic wave functions are localized, they show strong local variations leading to the breakdown of the LSDA approximation.

The LDA+U approach [151–153] consists in treating the exchange-correlation potential with a mean-field solution of the Hubbard model in the atomic limit for the subset of $d$ and/or $f$ shells, while the remaining delocalized $s$, $p$ electrons are still treated within LSDA. We represent the $d$ or $f$ shell by a $N$ orbital-degenerate system, for which the Hamiltonian can be written down for one site only in the atomic limit $U \gg t_{ij}$:

$$H = \frac{U}{2} \sum_{\substack{\alpha, \alpha' \\ \alpha \neq \alpha'}} n_{\alpha} n_{\alpha'} \qquad (1.22)$$

where $\alpha = \{m, \sigma\}$ combines the spin ($\sigma$) and orbital ($m = 1, \dots, N$) indices. Following Ref. 151, we show in the following that this simple Hamiltonian can overcome the two limitations of LSDA mentioned earlier in this section.

In the single orbital $N = 1$ case at half-filling, if the system is paramagnetic $n_{\uparrow} = n_{\downarrow} = 0.5$ then its energy at the mean-field level is $E_{PM} \approx \frac{U}{4}$. This is not optimal compared to the spin-polarized case with $n_{\uparrow} = 1$ (or $n_{\downarrow} = 1$) which yields $E_{SP} \approx 0$. The system favors the magnetic polarized solution and the energy splitting between the two spin states is $\approx U$, i.e. is governed by the Coulomb interaction, in contrast to the LSDA approximation where it is lead by the homogeneous gas' exchange parameter. This exchange parameter is to be identified with the Hund's exchange, which is usually order order 1eV, much smaller than the effective Coulomb interaction $U \approx 10$eV [151].

Now consider a $N$ orbital degenerate system with $n$ electrons. In the paramagnetic (both in terms of spin and orbital) limit, each of the $2N$ energy levels are filled by $\frac{n}{2N}$ electrons. Thus, the mean-field estimate for the energy is:

$$E_{PM} \approx \frac{U}{2} 2N(2N-1) \left(\frac{n}{2N}\right)^2 = \frac{U}{2} \left(1 - \frac{1}{2N}\right) n^2 \qquad (1.23)$$
If spin polarization only is allowed, then the majority states will be on average occupied by $\frac{n}{N}$ (1) electrons if $n < N$ ($n > N$), and the minority states by $\frac{(n-N)}{N}$ if $N < n \le 2N$, and 0 otherwise. This leads to:

$$ 
\begin{aligned}
E_{SP} &\approx \frac{U}{2}N(N-1)\left(\frac{n}{N}\right)^2 = \frac{U}{2}\left(1-\frac{1}{N}\right)n^2 & \text{if } 0 < n \le N \\
E_{SP} &\approx \frac{U}{2}\left(1-\frac{1}{N}\right)n^2 + U(n-N) & \text{if } N < n \le 2N.
\end{aligned} 
\quad (1.24) $$

The spin polarized solution improves the energy compared to the paramagnetic case, but it is still not optimal. Most importantly, because of the orbital mixing the Slater gap is only between states of different spins. Therefore, there will be a gap only at half-filling, since for lower (higher) occupations states with the same spin flavor as the last occupied ones are available, providing gapless excitation channels. This underlines the necessity of including orbital polarization, since if we consider that the n electrons preferably occupy a certain set of orbitals, the ground state energy is approximately:

$$ E_{SOP} \approx \frac{U}{2}n(n-1). \quad (1.25) $$

Orbital and spin polarization therefore leads to the most optimal ground state energy, and the gap between the higher occupied state and the lower unoccupied is $\approx U$, in contrast with standard LSDA in which the splitting is governed by Hund's exchange $J \ll U$ [151, 172]. Similarly to the spins, the orbital degree of freedom also polarizes due to the interactions, what can not be captured by the LSDA orbital-independent potential. To this effect can be added the crystal field splitting acting on the orbitals, which is usually weaker than the Coulomb interaction driven polarization.

LDA+U incorporates these effects in practice by defining a rotationally invariant potential for the correlated orbitals [153] (written for a single site for simplicity):

$$ V^{U}[n] = \frac{1}{2} \sum_{m_1, m_2, \sigma} U_{m_1 m_2} n_{m_1 - \sigma} n_{m_2 \sigma} + \sum_{m_1 > m_2, \sigma} (U_{m_1 m_2} - J_{m_1 m_2}) n_{m_1 \sigma} n_{m_2 \sigma} \quad (1.26) $$

where $U_{m_1 m_2} = \langle m_1 m_2 | \hat{V}_{Coulomb} | m_1 m_2 \rangle$ and $J_{m_1 m_2} = \langle m_1 m_2 | \hat{V}_{Coulomb} | m_2 m_1 \rangle$ are the screened Coulomb interaction matrix elements. The LSDA potential is corrected via:

$$ V[n] = V^{LSDA}[n] + V^{U}[n] - V^{dc}[n] \quad (1.27) $$

where $V^{dc}[n]$ ensures that we recover the LSDA potential in absence of orbital polarization and is defined by:

$$ V^{dc}[n] = \frac{1}{2}Un(n-1) - \frac{1}{2}J(n^{\uparrow}(n^{\uparrow}-1) + (n^{\downarrow}(n^{\downarrow}-1)) \quad (1.28) $$

with U and J the screened Coulomb and exchange parameters.

One of the main advantage of LDA+U is that it is computationally cheap. It will prove extremely useful in the study of the interface properties of the t-CuO/STO junction in Chap. 5 of this thesis.
For these kind of systems, **DFT** is the only computationally accessible technique due to the system size to be considered and **LSDA** alone would fail to give a correct description. Indeed **t-CuO** would be predicted to be a non-magnetic metal, preventing us to draw scenarios on the appearance of magnetic moments at the **t-CuO/STO** interface.

A clear limitation of **LDA+U** is that it requires spin and orbital polarization to open an electronic gap. This is related to the fact that this method amounts to add a constant and purely real correction term (the self-energy) to the non-interacting Hamiltonian. Such correction is non zero only if the system displays an orbital and/or spin polarization. The self-energy, that we discuss in more details in the next chapter, incorporates the modification of the non-interacting Hamiltonian from the presence of electronic interactions. Since it is fully real within **LDA+U**, the Kohn-Sham energies are purely real, which means that they have infinite lifetime. We remain in a band structure picture. **LDA+U** would not be able to capture a correlation-driven gap in the half-filled paramagnetic case as we will observe later in Chap. 4. In this situation, the metal-insulator transition is related to the divergence of the imaginary part of the self-energy at the Fermi level which suppresses the quasiparticle excitations [171] and decreases drastically their lifetime. All these effects showed to be essential to capture the spectral properties that can be probed for instance in photoemission measurements, and for which the knowledge of the excited states is essential but can not be obtained from a **DFT** calculation.

For all these reasons, one has to go beyond **DFT(+U)** and treat the interaction effects in the correlated subsystem: the *d* orbitals for the systems under investigation in this thesis. The Hubbard model has proven to be useful at the simple mean-field level within **LDA+U**, but a more involved treatment is needed to obtain the dynamic electronic properties. One powerful method is the so-called **DMFT** and its cluster extension **CDMFT** that we describe in the next chapter.
# 2 MANY-BODY METHODS

## 2.1 FROM THE BAND STRUCTURE TO A LOW-ENERGY EFFECTIVE MODEL

Despite the developments of new functionals, band-theory remains strongly limited for the description of correlated systems like the family of the TMOs [5], to which the copper oxides belong. The Coulomb interaction $U$ being of the same order of magnitude than the bandwidth $W$, a mean-field treatment of the interaction fails and perturbation theory can not be applied neither in $U$, nor in $W$. Moreover, in order to interpret experimental probes like ARPES, the inclusion of the dynamic effects is essential (see App. A). Indeed, in a metal the quasi-free electrons have a long life-time which leads to sharp experimental bands close to the theoretical band picture. In contrast, strongly correlated systems can yield dramatic lifetime suppression, spectral weight redistributions, resulting in a much different spectral function as we shall see later in this thesis. Effective models, which focus on a restricted set of degrees of freedom, are needed to describe the relevant mechanisms in those materials.

In this section, we first motivate the use of the Hubbard model as an effective system that aims at describing the low-energy physics of copper oxides. We then present a way to parameterize the model such that it is suited to capture the essential low-energy physics of a specific material.

### 2.1.1 THE HUBBARD MODEL

The rich phase diagram of copper oxides at low temperature is related to the low-energy properties, i.e. to the electronic states close to the Fermi level. Since we have seen that the manifold of low-energy bands for such materials is rather related to localized orbitals, we can model it by a simple lattice of atomic-like orbitals. The overlaps between these orbitals allow the electrons to hop between neighboring sites, which can be represented by a tight-binding (TB) Hamiltonian. In the case of copper oxides, the $d$ orbitals of the Cu atoms are in a $d^9$ configuration due to the hybridization with the surrounding O atoms: all $d$ orbitals are filled but one, usually the $d_{x^2-y^2}$ which is half-filled. It has been proposed by Zhang and Rice that the low-energy properties of such system can be approximated by a single-orbital model due to the formation of ZRS between the Cu $d_{x^2-y^2}$ and O $p_{x/y}$ orbitals centered on the Cu sites [38], which then leads to a single-band TB Hamiltonian:

$$H_0 = \sum_{\langle i,j \rangle, \sigma} t_{ij} c_{i\sigma}^{\dagger} c_{j\sigma}, \qquad (2.1)$$

where the indices $i$ refer to the lattice sites, and $t_{ij}$ is the hopping amplitude between two neighbors $\langle i, j \rangle$. $H_0$ gives rise to energy bands similarly to DFT. The ZRS form a 2D lattice of single
orbitals, and for the sake of simplicity we assume that the hopping terms are limited to the NN term $t$. Using a Fourier transformation, one easily finds the following energy dispersion:

$$ \epsilon(\mathbf{k}) = -2t[\cos(k_x) + \cos(k_y)], \qquad (2.2) $$

where we set the lattice constant $a = 1$. This dispersion (Fig. 2.1(c)) can, for instance, be compared with the band structure of NaCCOC (Fig. 2.1(a)) calculated using the Wien2k package [173, 174] along the the $\Gamma - X - M - \Gamma$ k-path. As expected from the ZRS approximation, the low-energy part of NaCCOC is composed of an almost isolated single band crossing the Fermi energy. LDA predicts the system to be metallic since the ZRS band is half-filled, as discussed in the previous chapter. The simple TB results looks quite alike the DFT band that crosses the Fermi energy (see Fig. 2.1(c)), and in fact by including longer-range hoppings and choosing the right hopping parameters, one can adjust the TB dispersion so to fit as best a possible the DFT bands (see Fig. 2.1(c)). We show later in this section that this parameterization can be done rigorously by a transformation of the delocalized Bloch-like wavefunctions used in DFT (Fig. 2.1(d)) into a set of atomic-like localized Wannier functions (Fig. 2.1(e)). By providing the non-interacting part $H_0$ of the effective Hamiltonian, the band structure calculations with the LDA functional is a good starting point.

To our TB hamiltonian we can now add the Coulomb interaction, which in principle is a long-ranged interaction. The problem is simplified by considering that, due to electronic screening of the Coulomb potential, the effective on-site interaction $U$ is much larger than the inter-site interaction $V$, which is then neglected in our modeling. This is verified in cuprates as shown by Hirayama et al. for different compounds: $U$ is typically almost an order of magnitude larger than $V$ [31]. Moreover, as is explained in Chap. II, 3 and 4, we will restrict our effective model to a single effective ZRS orbital on each site. With the addition of the on-site interaction $U$, we get to the well-known Hubbard model [35–37]:

$$ H = \sum_{\langle i,j \rangle, \sigma} t_{ij} c_{i\sigma}^{\dagger} c_{j\sigma} + U \sum_{i} n_{i\uparrow} n_{i\downarrow}, \qquad (2.3) $$

which we already mentioned in the context of LDA+U. This model can not be solved exactly because of the interplay between the one-body kinetic term and the two-body interaction one. To solve it beyond the mean-field solution, we use the powerful method called DMFT [42, 43] that is described later in this chapter. The combination of LDA to obtain the band structure and DMFT to treat the correlated subsystem is usually called single-shot LDA+DMFT, and was first apply shortly after DMFT was invented under the name of LDA++ [49, 152].

One might be puzzled here since in principle a part of the electronic interactions is taken into account at the DFT level, and which seems to be ignored in the TB Hamiltonian, such that if we add interactions to $H_0$ we shall face a double-counting problem. The interactions within DFT lead to a constant self-energy for each band, i.e. static band-dependent energy shifts which have to be corrected before adding correlations to the model [152]. In this thesis, the problem of double-counting is avoided since we consider only the single-band Hubbard model for which the
<figure><img src="image_9.png" /><figcaption>Figure 2.1: (a) Crystal structure of NaCCOC ($x = 0$) with a sketch of the **k** path in the first BZ. (b) Band structure obtained with a LDA calculation, along with the Wannier band obtained with the MLWF method (blue points). (c) Wannier band (blue points), and its approximated TB band obtained by neglecting the hoppings terms beyond the NNNN. The orange dashed-dotted line shows the 2D dispersion of Eq. 2.2 calculated with only $t$. (d) Sketch of a delocalized Bloch wavefunction $\phi_{nk}(\mathbf{r})$: the black lines represent the periodic function $u_{nk}$ weighted by the $e^{ikr}$ term (green). (e) Sketch of localized Wannier functions $w_n(\mathbf{r})$ centered around the atomic positions $\{n\}$.</figcaption></figure>

constant shift can be absorbed in the chemical potential.

Now that we have motivated the use of the Hubbard model as an effective system for the correlated manifold, we are left with describing how to parameterize the hopping and interaction parameters to obtain a realistic description of the low-energy properties of the material of interest (the subject of the remainder of this section), and how efficiently solve this model (the remainder of this chapter).

### 2.1.2 PARAMETERIZATION

The DFT bands can be expressed with a set of delocalized Bloch wavefunctions (illustrated in Fig. 2.1(d))

$$ \phi_{nk}(\mathbf{r}) = u_{nk}(\mathbf{r})e^{ikr}, \qquad (2.4) $$

where $u_{nk}(\mathbf{r})$ is a periodic function of $\mathbf{r}$, i.e. $u_{nk}(\mathbf{r} + \mathbf{R}) = u_{nk}(\mathbf{r})$ with $\mathbf{R}$ a lattice vector. To parameterize the non-interacting part of Eq. 2.3, we need to switch to a localized basis set $\{|w_i\rangle\}$ and compute the hopping terms as the overlaps between neighboring sites:

$$ t_{ij} = \int d\mathbf{r} \langle w_i | H_{KS} | w_j \rangle. \qquad (2.5) $$
There exist different methods to perform the transformation to a localized basis: one can use linear muffin-tin orbitals [49, 152, 175], N-th order muffin-tin orbitals [176], Wannier functions constructed by projection [177, 178], or projection of the full-potential linear augmented plane-wave onto a local basis [179, 180], just to name a few. In this thesis, we have used another method called the maximally localized Wannier functions (MLWF) [181, 182] that we present in the following.

From an isolated manifold of N occupied Bloch bands $\phi_{nk}(\mathbf{r})$, a set of N Wannier functions $w_{nR}(\mathbf{r} - \mathbf{R})$, where $\mathbf{R}$ is a lattice vector, can be constructed:

$$|w_{n\mathbf{R}}\rangle = \frac{V}{(2\pi)^3} \int_{BZ} \left( \sum_{m=1}^{N} U_{mn}^{\mathbf{k}} |\phi_{m\mathbf{k}}\rangle \right) e^{-i\mathbf{k}\mathbf{R}} d\mathbf{k}. \qquad (2.6)$$

The integration is done over the BZ, and $U^k$ is a unitary matrix mixing the bands at each k point. The Wannier functions are not unique since:

(i) the Bloch wave functions are defined up to a phase factor $e^{i\theta(k)}$,

(ii) At each k point the set of occupied bands may be unitary transformed by some $U^k$.

Instead of being an issue, this non-uniqueness of the Wannier functions can be exploited to construct MLWF as was first proposed by Marzari et al. [181]. The strategy consists in performing a variational search of the unitary matrix $U^k$ that minimizes the spread function of the Wannier orbitals around their center:

$$\begin{align}
\Omega &= \sum_{n=1}^{N} \langle (\mathbf{r} - \langle \mathbf{r} \rangle_n)^2 \rangle_n = \sum_{n=1}^{N} \left( \langle \mathbf{r}^2 \rangle_n - |\langle \mathbf{r} \rangle_n|^2 \right) \tag{2.7} \\
\langle \hat{A} \rangle_n &= \langle w_{n0} | \hat{A} | w_{n0} \rangle. \nonumber
\end{align}$$

It is easily shown that the spread function $\Omega$ can be separated into two parts:

$$\begin{align}
\Omega &= \sum_{n=1}^{N} \left( \langle \mathbf{r}^2 \rangle_n - \sum_{m\mathbf{R}} |\langle w_{m\mathbf{R}} | \mathbf{r} | w_{n0} \rangle|^2 \right) + \sum_n \sum_{m\mathbf{R} \neq n0} |\langle w_{m\mathbf{R}} | \mathbf{r} | w_{n0} \rangle|^2 \tag{2.8} \\
&= \Omega_I + \tilde{\Omega}. \nonumber
\end{align}$$

It turns out that $\Omega_I$ is gauge invariant [181], that is to say that it does not depend on the choice of the unitary matrix $U^k$, such that the minimization procedure is only done on $\tilde{\Omega}$ by evolving $U^k$ using a standard minimization algorithm. Note that the set of matrices $U^k$ can be constrained if one expects a MLWF of a certain symmetry, for more details the interested reader can refer to the wannier90 package documentation [183, 184].

The main drawback of this method is that it is constructed on the assumption that there is an isolated set of N Bloch wavefunctions at our disposal. However, this is usually not the case, and especially when constructing simpler effective models like the single band Hubbard model we use in this thesis. Indeed, by inspecting Fig. 2.1(b), we notice that the band of interest is entangled with the manifold of bands lying around -1.5eV, preventing us from using the above method directly since it is designed to extract N Wannier functions from a set of N bands.
To overcome this issue, a disentanglement procedure was proposed by Souza et al. in Ref. 182, which amounts to perform an additional minimization of the gauge invariant function $\Omega_I$. The first step is to select an energy window containing $N_k \le N$ bands. At each $k$ point, there exist multiple $N$-dimensional subspaces $S^{(k)}$ which can be spanned by unitary transformation of the $N_k$ states:

$$|u_{nk}^{opt}\rangle = \sum_{m \in N_k} U_{mn}^{dis(k)} |u_{mk}\rangle, \qquad (2.9)$$

where $U^{dis(k)}$ is a rectangular $N_k \times N$ matrix, said to be unitary since $U^{dis(k)}(U^{dis(k)})^{\dagger} = \mathbb{1}$, and $|u_{nk}\rangle$ refers to the periodic part of the Bloch functions. Since $\Omega_I$ is invariant within a given subspace $S^{(k)}$, it can be considered as a functional of the subspace $\Omega_I[S^{(k)}]$. For each $k$ point, the latter is minimized by choosing the optimal subspace $S_{opt}^{(k)}$, and then $\tilde{\Omega}$ is optimized within each $S_{opt}^{(k)}$.

For a physical interpretation of this additional minimization, it is insightful to re-write $\Omega_I$ as follows [182]:

$$\begin{align}
\Omega_I &\propto \sum_{k,k'} w_{k'} \sum_{m=1}^{N} \left( 1 - \sum_{n=1}^{N} |M_{mn}^{(k,k')}|^2 \right), \tag{2.10} \\
M_{mn}^{(k,k')} &= \langle u_{mk} | u_{mk+k'} \rangle, \nonumber
\end{align}$$

where $w_{k'}$ is a weight arising from the discretization linked with finite differences [181]. Intuitively, we notice that $\Omega_I$ is minimum when the overlaps $M_{mn}^{(k,k')}$ are maximized, i.e. when $S^{(k)}$ is the subspace with minimum mismatch as $k$ is varied. For instance, in the attempt to apply such procedure to the case of NaCCOC, $S^{(k)}$ is optimal when it keeps the single band crossing the Fermi energy and excludes the bands from the $-1.5$eV manifold.

MLWF can be routinely obtained from a band structure calculation with the wannier90 package [183, 184]. Most importantly, hopping parameters can be extracted and an effective TB Hamiltonian is obtained. In Fig. 2.1(b)-(c), we show an example of an effective MLWF band obtained for NaCCOC by fitting the ZRS band that crosses the Fermi energy, and its approximated TB band for which the hopping terms were kept until the next-NNN (NNNN) neighbor. The Kohn-Sham band is extremely well reproduced by the Wannier band, as well as the simpler TB band: we now dispose of a low-energy non-interacting model with a set of hopping parameters $t_{ij}$ obtained ab initio.

Having determined the non-interacting part of our effective model, we are left with the Coulomb interaction $U$. In principle, the Coulomb matrix (in the case of a multi-orbital system) can be determined ab initio [185] using method like the constrained random phase approximation [186, 187], or constrained DFT [188–190], but in this thesis we either determined it via a fit to an experimental measurement, or empirically. We therefore do not describe the procedure here, and rather go on with the description of a way to solve our realistic Hubbard model: the dynamical mean-field theory (DMFT).
## 2.2 SOLVING THE HUBBARD MODEL: DYNAMICAL MEAN-FIELD THEORY

There is no exact solution to the Hubbard model, except for limiting cases like the atomic limit ($t_{ij} = 0$), or the non-interacting limit ($U = 0$), or in 1D [191]. As we have seen in the previous chapter, static mean-field approaches can lead to constant and purely real corrections to the (almost) non-interacting solution, but it is insufficient for instance to capture the physics of the Mott insulating transition. In contrast, DMFT [42, 43] in its construction retains the full dynamics of the local quantum fluctuations. It is based on a mapping from the full lattice interacting problem onto an impurity model in which a site is dynamically embedded into a self-consistent electron bath. The impurity model can be seen as an Anderson's impurity model (AIM) [18] in which all the local fluctuations are treated exactly (while non-local correlations are discarded), and which can be efficiently solved numerically. In contrast to the well-known Weiss molecular-field theory for the Ising Model [192], in which the spin interactions are replaced by the action of a static mean-field, in DMFT the occupation of the effective impurity site is not kept constant since electrons can hop back and forth from the electron bath, i.e. it is a time-dependent problem. In principle, DMFT makes no further approximation and is valid for all ranges of temperatures, interaction strength, and similarly to the Weiss theory it becomes exact in the limit of infinite dimension $d \to \infty$, which amounts to $z \to \infty$ with z the coordination number.

In fact, DMFT was developed after the pioneering works of Metzner and Vollhardt [44–47] who studied the properties of the Hubbard model in infinite dimensions, that we present first in this section. Following these findings, Ohkawa [193] and Georges and Kotliar [42] showed that in the $d \to \infty$ limit there is an exact mapping between the Hubbard model and the AIM with an effective non-interacting electronic bath. We detail one of the possible derivations, the so-called cavity method [42, 43], and finally present the set of self-consistent equations at the heart of DMFT.

### 2.2.1 THE LIMIT OF INFINITE DIMENSION

Similarly to the Ising model, the Hubbard model in infinite dimension yields peculiar properties [44–46] that will lead to drastic simplifications of the problem. We first recall the Hamiltonian for the single-orbital case:

$$H = \sum_{\langle i,j \rangle, \sigma} t_{ij} c_{i\sigma}^{\dagger} c_{j\sigma} + U \sum_{i} n_{i\uparrow} n_{i\downarrow} - \mu \sum_{i} n_{i}, \qquad (2.11)$$

where we keep, for the sake of simplicity, only the nearest-neighbor hopping term $t$, and we added the chemical potential ($\mu$) term. Since the Hubbard $U$ interaction is purely local, only the kinetic part of the Hamiltonian depends on the dimension $d$ of the problem. Assuming that the system is a $d$-dimensional hypercubic lattice, i.e. the coordination number is $z = 2d$, the non-interacting dispersion can be written:

$$\epsilon_k = -2t \sum_{i=1}^{d} \cos(k_i), \qquad (2.12)$$
where we set the lattice constant $a = 1$. The density of states (DOS) associated to $\epsilon_k$ is given by

$$D(\omega) = \sum_{\mathbf{k}} \delta(\omega - \epsilon_{\mathbf{k}}), \qquad (2.13)$$

where $\delta$ is the delta Dirac function. The DOS can be interpreted as the probability density for $\omega = \epsilon_k$ given a random choice of $\mathbf{k}$: the dispersion (Eq. 2.12) is simply a sum over random variables $-2t \cos(k_i)$ with $\mathbf{k} = (k_1, ..., k_d)$. The summation being over the dimension $d$, in the limit of $d \to \infty$ one can apply the central limit theorem to the DOS $D(\omega)$ giving

$$D(\omega) \approx \frac{1}{\sqrt{2\pi t \sqrt{2d}}} e^{-\frac{1}{2} \left(\frac{\omega}{t\sqrt{2d}}\right)^2} = \frac{1}{\sqrt{2\pi t \sqrt{z}}} e^{-\frac{1}{2} \left(\frac{\omega}{t\sqrt{z}}\right)^2}. \qquad (2.14)$$

The variance $\sigma = t\sqrt{z}$ becomes infinite in the $d \to \infty$ limit if the hopping $t$ is not rescaled. The DOS would therefore become flat and featureless, without any physical relevance. Rescaling the hopping $t$ is necessary for having a well-defined DOS in infinite dimension, one of the possible choices is

$$t \to \frac{t}{\sqrt{z}}, \qquad (2.15)$$

which results in [44]

$$D(\omega) = \frac{1}{\sqrt{2\pi t}} e^{-\frac{1}{2} \left(\frac{\omega}{t}\right)^2}. \qquad (2.16)$$

When $d \to \infty$, with this choice of rescaling, the DOS becomes a Gaussian. Interestingly, one observes that the DOS converges quickly to this limit as $d$ is increased [47]: already at $d = 3$, i.e. $z = 6$, the shape is very similar to a Gaussian. Note that this choice of rescaling is not the only one possible [47], but it has the virtue of keeping the problem non-trivial at infinite dimension since the kinetic and Coulomb parts of the Hamiltonian remain in competition.

This scaling has important consequences on the many-body observables of the system like the Green's function or the self-energy. As a brief reminder, the Matsubara Green's function is defined as follows

$$G_{ij,\sigma}(\tau - \tau') = -\langle T_{\tau} c_{i\sigma}(\tau) c_{j\sigma}^{\dagger}(\tau') \rangle \qquad (2.17)$$

where $i, j$ denotes the site indices, $\tau$ is the imaginary time, $T_\tau$ is the time ordering operator, and $\langle \cdot \rangle$ denotes the thermodynamic average. It measures the propagation of an electron through the correlated system from site $j$ at time $\tau'$ to site $i$ at time $\tau$. In the case of a non-interacting system, say the Hubbard model restricted to its kinetic part, the non-interacting Green's function can be written as a function of Matsubara frequency and momentum $\mathbf{k}$ as:

$$G_{\mathbf{k},\sigma}^{0}(i\omega_{n}) = \frac{1}{i\omega_{n} + \mu - \epsilon_{\mathbf{k}}} \qquad (2.18)$$

where $\omega_n = \frac{(2n+1)\pi}{\beta}$ are the fermionic Matsubara frequencies, $\beta$ is the inverse temperature. When including interactions, the Hamiltonian can not be diagonalized anymore and the correlation ef-
fects can be conveniently encoded in the Self-energy $\Sigma_{\mathbf{k},\sigma}(i\omega_n)$, which is related to the interacting Green's function by:

$$ 
\begin{align} 
G_{\mathbf{k},\sigma}(i\omega_n) &= \frac{1}{i\omega_n + \mu - \epsilon_{\mathbf{k}} - \Sigma_{\mathbf{k},\sigma}(i\omega_n)} \tag{2.19} \\ 
\Sigma_{\mathbf{k},\sigma}(i\omega_n) &= [G_{\mathbf{k},\sigma}^0(i\omega_n)]^{-1} - [G_{\mathbf{k},\sigma}(i\omega_n)]^{-1}. \nonumber 
\end{align} 
$$

The second equation, that relates the Self-energy to the interacting and non-interacting Green's function, is the so-called Dyson equation.

To see the effects of the infinite dimension limit on the many-body observables, it is instructive to start from the following sum rule for a non-interacting system [45]:

$$ \sum_i (\langle c_{i\sigma}^{\dagger} c_{j\sigma} \rangle^0)^2 = n_{j,\sigma}^0 \qquad (2.20) $$

where the superscript 0 denotes the observables evaluated for the non-interacting system. Now, the sum may be re-written by grouping the $i$ indices with respect to their distance from the site $j$, that is to say we group them by groups of NN, NNN, etc. Taking care of the one term $i = j$, Eq. 2.20 may be written:

$$ 
\begin{align} 
\sum_i (\langle c_{i\sigma}^{\dagger} c_{j\sigma} \rangle^0)^2 &= (n_{j,\sigma}^0)^2 + \sum_{l=1}^{\infty} z_l (\langle c_{l\sigma}^{\dagger} c_{j\sigma} \rangle^0)^2 \tag{2.21} \\ 
\sum_{l=1}^{\infty} z_l (\langle c_{l\sigma}^{\dagger} c_{j\sigma} \rangle^0)^2 &= n_{j,\sigma}^0 - (n_{j,\sigma}^0)^2 \nonumber 
\end{align} 
$$

where $l$ denotes the neighbor rank (NN, NNN, etc), and $z_l$ the number of neighbors of rank $l$. Therefore, we can obtain an upper bound for the non-interacting Green's function:

$$ |\langle c_{i\sigma}^{\dagger} c_{j\sigma} \rangle^0| = |G_{ij,\sigma}^0(\tau \to 0)| \le \left[ \frac{n_{j,\sigma}^0 (1 - n_{j,\sigma}^0)}{z_l} \right]^{1/2}, \quad \text{for } \langle i, j \rangle \text{ of rank } l. \qquad (2.22) $$

In an hypercubic lattice, each site has $z$ NN, each of them having $z-1$ neighbors (excluding the first site), and etc. This means that in the $d \to \infty$ limit, $z_l \propto z^l$, such that one may write

$$ G_{ij,\sigma}^0(\tau \to 0) \propto \frac{1}{Z^{\|\mathbf{R}_i - \mathbf{R}_j\|/2}} \qquad (2.23) $$

where $\|\mathbf{R}_i - \mathbf{R}_j\|$ denotes the "Manhattan" distance between site $i$ and $j$. Although this was demonstrated only for $\tau \to 0$, the Green's function $G_{ij,\sigma}^0$ must follow the same scaling at all times since it does not depend on the time evolution. We therefore have, for a general couple of sites $(i, j)$:

$$ G_{ij,\sigma}^0 \propto \frac{1}{Z^{\|\mathbf{R}_i - \mathbf{R}_j\|/2}} \qquad (2.24) $$
<figure><img src="image_10.png" /><figcaption>Figure 2.2: Sketch of a second order diagram in real space contributing to the fully irreducible self-energy $\Sigma_{ij}$ (considering $\langle ij \rangle$ nearest-neighbor for simplicity). In the $d \to \infty$ limit all the diagrams become purely local due to the scaling discussed in the text.</figcaption></figure>

Despite this scaling of the Green's function, one must not discard the inter-site elements since these terms still contribute to the kinetic energy:

$$E_{kin}^0 = - \sum_{\langle i,j \rangle, \sigma} t_{ij} \langle G_{ij,\sigma}^0 (\tau \to 0^-) \rangle^0. \qquad (2.25)$$

For a couple of neighbors $\langle ij \rangle$, the Green's function scales as $z^{-\|\mathbf{R}_i - \mathbf{R}_j\|/2}$, and the hopping as $t_{ij} \propto z^{-\|\mathbf{R}_i - \mathbf{R}_j\|/2}$, but the sum brings a factor $z^{\|\mathbf{R}_i - \mathbf{R}_j\|}$, the number of equivalent ways to consider two neighbors $\langle ij \rangle$ with the same distance $\|\mathbf{R}_i - \mathbf{R}_j\|$.

From the Green's function scaling one can show the most important simplification at $d \to \infty$: the self-energy becomes purely local [44, 45]. One suggestive way to show that is to consider the diagrammatic expansion of the self-energy. We can reduce the discussion to diagrams of order $n > 2$ since the Hartree term is just a constant and local correction, while the Fock term vanishes in the case of a Hubbard model with only the on-site U interaction. A typical second order diagram that contributes to the fully irreducible self-energy $\Sigma_{ij}$ is shown in Fig. 2.2. Each line of the diagram, representing a non-interacting propagator between site i and j, will contribute to $z^{-\|\mathbf{R}_i - \mathbf{R}_j\|/2}$ in the $d \to \infty$ limit. In total, this amounts to a scaling of $z^{-3\|\mathbf{R}_i - \mathbf{R}_j\|/2}$, but the sum over all neighbors $\langle ij \rangle$ brings only a factor $z^{\|\mathbf{R}_i - \mathbf{R}_j\|}$, i.e. the overall contribution to $\Sigma_{ij}$ scales as $z^{-\|\mathbf{R}_i - \mathbf{R}_j\|/2}$. The same can be shown for higher order diagrams, such that $\Sigma_{ij} \propto z^{-\|\mathbf{R}_i - \mathbf{R}_j\|/2}$ in the $d \to \infty$ limit. In other words, only the local component of the self-energy survives at infinite dimension, it is therefore purely local:

$$\lim_{d \to \infty} \Sigma_{ij,\sigma}(i\omega_n) = \delta_{ij} \Sigma_{ij,\sigma}(i\omega_n). \qquad (2.26)$$

It follows that the self-energy becomes momentum independent:

$$\Sigma_{\mathbf{k},\sigma}(i\omega_n) = \Sigma(i\omega_n) \qquad (2.27)$$

and the local Green's function is calculated via:

$$G_{\sigma}(i\omega_n) = \sum_{\mathbf{k}} \frac{1}{i\omega_n + \mu - \epsilon_{\mathbf{k}} - \Sigma_{\sigma}(i\omega_n)} \qquad (2.28)$$
illustrated in Fig. 2.3, leads to an effective representation of the many-body interacting problem by an isolated site embedded into a self-consistent non-interacting electron bath. By identifying the hybridization function of the AIM to the Weiss field of the cavity method, one observes that the solution of the Hubbard model in $d \to \infty$ can be obtained by solving an impurity problem. Similarly to DFT with the density, here the effective hybridization function is not known a priori, but can be found with a self-consistent procedure.

The cavity construction relies on the action formulation of the Hubbard model. We start from the partition function $Z$ of the lattice system which is described by the action $S$:

$$Z = \int \prod_{i,\sigma} Dc_{i\sigma}^{\dagger} Dc_{i\sigma} e^{-S}, \qquad (2.29)$$

where $c_{i\sigma}^{\dagger}, c_{i\sigma}$ are the Grassman variables [200], and the action $S$ is given by:

$$S = \int_{0}^{\beta} d\tau \left[ \sum_{i,\sigma} c_{i\sigma}^{\dagger}(\tau) \left( \frac{\partial}{\partial\tau} - \mu \right) c_{i\sigma}(\tau) - \sum_{i,j,\sigma} t_{ij} c_{i\sigma}^{\dagger}(\tau) c_{j\sigma}(\tau) + U \sum_{i} c_{i\uparrow}^{\dagger}(\tau) c_{i\uparrow}(\tau) c_{i\downarrow}^{\dagger}(\tau) c_{i\downarrow}(\tau) \right]. \quad (2.30)$$

We proceed by separating the action into three terms: $S_0$ which contains only the site 0, $S^{(0)}$ for the rest of the lattice, and $S_{hyb}$ containing the hopping terms between the lattice and the site 0.

$$S = S_0 + S_{hyb} + S^{(0)}$$

$$S_0 = \int_0^\beta d\tau \left[ \sum_\sigma c_{0\sigma}^\dagger(\tau) \left( \frac{\partial}{\partial\tau} - \mu \right) c_{0\sigma}(\tau) + U c_{0\uparrow}^\dagger(\tau) c_{0\uparrow}(\tau) c_{0\downarrow}^\dagger(\tau) c_{0\downarrow}(\tau) \right]$$

$$S_{hyb} = - \int_0^\beta d\tau \sum_{i \neq 0, \sigma} \left[ t_{i0} c_{i\sigma}^\dagger(\tau) c_{0\sigma}(\tau) + t_{0i} c_{0\sigma}^\dagger(\tau) c_{i\sigma}(\tau) \right] = - \int_0^\beta d\tau S_{hyb}(\tau)$$

$$S^{(0)} = \int_{0}^{\beta} d\tau \left[ \sum_{i \neq 0, \sigma} c_{i\sigma}^{\dagger}(\tau) \left( \frac{\partial}{\partial\tau} - \mu \right) c_{i\sigma}(\tau) - \sum_{i,j \neq 0, \sigma} t_{ij} c_{i\sigma}^{\dagger}(\tau) c_{j\sigma}(\tau) + U \sum_{i \neq 0} c_{i\uparrow}^{\dagger}(\tau) c_{i\uparrow}(\tau) c_{i\downarrow}^{\dagger}(\tau) c_{i\downarrow}(\tau) \right]. \quad (2.31)$$

We can now insert Eq. 2.31 into Eq. 2.29, and explicitly separate the integration over the site 0 and the rest of the lattice:

$$\begin{aligned} Z &= \int \prod_{\sigma} Dc_{0\sigma}^{\dagger} Dc_{0\sigma} e^{-S_0} \int \prod_{i \neq 0, \sigma} Dc_{i\sigma}^{\dagger} Dc_{i\sigma} e^{-S_{hyb}-S^{(0)}} \\ &= Z_{latt} \int \prod_{\sigma} Dc_{0\sigma}^{\dagger} Dc_{0\sigma} e^{-S_0} \langle e^{-S_{hyb}} \rangle_{latt} \end{aligned} \qquad (2.32)$$

where the second line is obtained by integrating out the rest of the lattice, yielding the thermodynamic average $\langle \cdot \rangle_{latt}$ and the partial partition function $Z_{latt}$. The exponential $\langle e^{-S_{hyb}} \rangle_{latt}$ can be expanded in power series of $S_{hyb}$, and because it is averaged over the lattice without the 0 site,
only the even terms in the expansion are non-zero (since odd terms do not conserve the number of particle in the lattice with site 0 excluded). The partition function is then written:

$$Z = \int \prod_{\sigma} Dc_{0\sigma}^{\dagger} Dc_{0\sigma} e^{-S_0} \left[ 1 + \frac{1}{2!} \int_{0}^{\beta} d\tau_1 \int_{0}^{\beta} d\tau_2 \langle S_{hyb}(\tau_1)S_{hyb}(\tau_2) \rangle_{\text{latt}} \right. \\ \left. + \frac{1}{4!} \int_{0}^{\beta} d\tau_1 \int_{0}^{\beta} d\tau_2 \int_{0}^{\beta} d\tau_3 \int_{0}^{\beta} d\tau_4 \langle S_{hyb}(\tau_1)S_{hyb}(\tau_2)S_{hyb}(\tau_3)S_{hyb}(\tau_4) \rangle_{\text{latt}} + \dots \right]. \quad (2.33)$$

In the limit $d \to \infty$ the expansion of the exponential is drastically simplified thanks to the scaling of the propagators and the hopping terms. This is clearly seen by re-writing for instance the second order term of the expansion:

$$\frac{1}{2!} \int_0^\beta d\tau_1 \int_0^\beta d\tau_2 \langle S_{hyb}(\tau_1)S_{hyb}(\tau_2) \rangle_{\text{latt}} = -\frac{1}{2!} \sum_{ij \neq 0, \sigma} t_{i0} t_{0j} \int_0^\beta d\tau_1 \int_0^\beta d\tau_2 G_{ij,\sigma}^{(0)}(\tau_1 - \tau_2) c_{0\sigma}^\dagger(\tau_1) c_{0\sigma}(\tau_2), \quad (2.34)$$

where $G_{ij,\sigma}^{(0)}$ is the Green's function of the lattice without the site 0, which is defined by:

$$G_{ij,\sigma}^{(0)}(\tau_1 - \tau_2) = -\langle T_{\tau} c_{i\sigma}(\tau_1) c_{j\sigma}^{\dagger}(\tau_2) \rangle_{\text{latt}}. \quad (2.35)$$

Similarly to Eq. 2.34, the $n^{th}$ (n even) term involves n hopping terms and a n-particle Green's function $G_{i_1, \dots, i_n, \sigma}^{(0)}$. In the second order term, the two hopping terms contribute to the scaling as $\propto \frac{1}{d}$, the sums $\propto d^2$, while the Green's function scales as $\propto \frac{1}{d^{\|R_i - R_j\|/2}}$. Since $\|R_i - R_j\| \ge 2$ in the second order term, the scaling is overall $O(1)$. This is different for the fourth order term, which involves four hopping terms $\propto \frac{1}{d^2}$, four sums $\propto d^4$, and a four-particle Green's function $G_{ijkl,\sigma}^{(0)}$ which scales as $\propto \frac{1}{d^{\|R_i - R_j\|/2}} \frac{1}{d^{\|R_j - R_k\|/2}} \frac{1}{d^{\|R_k - R_l\|/2}}$. Hence the scaling is at least $\propto \frac{1}{d}$, such that the fourth order and also all higher order terms are suppressed in infinite dimension. We end up with the effective local DMFT action:

$$S_{\text{eff}} = \int_0^\beta d\tau \sum_\sigma c_{0\sigma}^\dagger(\tau) \left(\frac{\partial}{\partial\tau} - \mu\right) c_{0\sigma}(\tau) + U \int_0^\beta d\tau c_{0\uparrow}^\dagger(\tau) c_{0\uparrow}(\tau) c_{0\downarrow}^\dagger(\tau) c_{0\downarrow}(\tau) \\ + \sum_{ij \neq 0, \sigma} t_{i0} t_{0j} \int_0^\beta d\tau_1 \int_0^\beta d\tau_2 G_{ij,\sigma}^{(0)}(\tau_1 - \tau_2) c_{0\sigma}^\dagger(\tau_1) c_{0\sigma}(\tau_2), \quad (2.36)$$

where all the hopping terms are implicitly rescaled using Eq. 2.15 since we are in the $d \to \infty$ limit. We now introduce the hybridization function $\Delta_\sigma$ and the Weiss field $G_\sigma$ as:

$$\Delta_{\sigma}(\tau_1 - \tau_2) = \sum_{ij \neq 0} t_{i0} t_{0j} G_{ij,\sigma}^{(0)}(\tau_1 - \tau_2) \\ G_{\sigma}^{-1}(\tau_1 - \tau_2) = -\delta(\tau_1 - \tau_2) \left(\frac{\partial}{\partial \tau_1} - \mu\right) - \Delta_{\sigma}(\tau_1 - \tau_2). \quad (2.37)$$

$$G_{\sigma}^{-1}(i\omega_n) = i\omega_n + \mu - \Delta_{\sigma}(i\omega_n)$$

The hybridization function $\Delta_\sigma(\tau_1 - \tau_2)$ describes the coupling of the impurity electrons to the effective bath. It encodes the propagation of an electron that hops from the impurity to the bath at time $\tau_2$, then propagates through the lattice without the impurity site via the lattice Green's function
<figure><img src="image_11.png" /><figcaption>Figure 2.4: Each line of a diagram represents a hopping term, and each vertex a cumulant (which does not scale with dimension). Examples of diagrams that contribute to the Green's function $G_{ij,\sigma}$ at $d \to \infty$ (a), and that vanishes (b). (c) In infinite dimension, the diagrams contributing to $G_{ij,\sigma}$ and that pass by the specific site 0 can be separated into contributions $G_{i0,\sigma}$ and $G_{0j,\sigma}$, such that $G_{ij,\sigma} \propto G_{i0,\sigma}G_{0j,\sigma}$. One observes that such decomposition induces a double counting of the decorations (in green) on site 0, evaluated to be $G_{00,\sigma}$, which is therefore corrected by writing $G_{ij,\sigma} = \frac{G_{i0,\sigma}G_{0j,\sigma}}{G_{00,\sigma}}$.</figcaption></figure>

$G_{ij,\sigma}^{(0)}(\tau_1 - \tau_2)$, and hops back to the impurity at time $\tau_1$. It therefore allows a time-dependent fluctuation of the electronic density on the impurity. When connected to a bath, the non-interacting Green's function of the impurity electrons, given by Eq. 2.18, is then corrected by the hybridization and yields the Weiss field $G_{\sigma}^{(0)}i\omega_n$, also called bath Green's function.

Using these definitions, the effective action may be written:

$$S_{eff} = - \int_{0}^{\beta} d\tau_1 \int_{0}^{\beta} d\tau_2 \sum_{\sigma} c_{\sigma}^{\dagger}(\tau_1) G_{\sigma}^{-1}(\tau_1 - \tau_2) c_{\sigma}(\tau_2) + U \int_{0}^{\beta} d\tau c_{\uparrow}^{\dagger}(\tau) c_{\uparrow}(\tau) c_{\downarrow}^{\dagger}(\tau) c_{\downarrow}(\tau), \quad (2.38)$$

where the subscript 0 was removed for readability. Therefore, as advertised above, in $d \to \infty$ the action of the Hubbard model can be exactly mapped onto that of an interacting impurity embedded into an effective bath of electrons. Eq. 2.37 provides the link between the Weiss field $G_{\sigma}$ and the interacting Green's function $G_{ij,\sigma}^{(0)}$ of the full lattice with one site removed. The last missing element is the relation between $G_{ij,\sigma}^{(0)}$ and $G_{ij,\sigma}$, the Green's function of the original problem.

For this last step, it is convenient to follow the ideas of Metzner in Ref. 46: the Green's function is expanded in perturbation theory over the hopping terms, starting from the atomic limit $U \gg t$. The terms of the expansion can be written as diagrams, for which each line corresponds to a hopping term, and each vertex to a cumulant which does not scale with the dimension. At $d \to \infty$, the only diagrams that contribute to $G_{ij,\sigma}$ are fully two-particle reducible: they connect site i and j via a single line and may have local polygon decorations, as shown in Fig. 2.4(a).
Indeed, we will show that the diagrams which are not fully two-particle reducible, like the one shown in Fig. 2.4(b), give a vanishing contribution to the kinetic energy [46]

$$E_{kin} \propto - \sum_{\langle i,j \rangle, \sigma} t_{ij} G_{ij,\sigma}. \qquad (2.39)$$

For simplicity, we consider only the nearest-neighbor hopping (the argument is easily generalized to longer-range hopping terms [46]). Because there are five lines, each corresponding to a hopping term, this diagram scales like $G_{ij,\sigma} \propto d^{-5/2}$. Since i and j are at a distance of 3 sites, $t_{ij} \propto d^{-3/2}$ and there are $d^3$ terms in the sum. In total, its contribution to the kinetic energy scales as $\propto d^{-1}$ and therefore vanishes at $d \to \infty$. Similarly, one can easily show that any polygon decoration as those shown in Fig. 2.4(a)-(c) contribute as $O(1)$.

One immediately sees that for each contribution to $G_{ij,\sigma}$ which passes by the specific site 0, the diagram can be separated into two subdiagrams going from i to 0, and from 0 to j, as shown in Fig. 2.4(c). Since the Green's function $G_{ij,\sigma}$ of the full lattice can be obtained by the sum of all the paths avoiding the site 0 ($G_{ij,\sigma}^{(0)}$), and all the contributions explicitly passing by 0 ($\propto G_{i0,\sigma}G_{0j,\sigma}$), then we have [46]:

$$G_{ij,\sigma} = G_{ij,\sigma}^{(0)} + \frac{G_{i0,\sigma}G_{0j,\sigma}}{G_{00,\sigma}} \qquad (2.40)$$

where we have divided by $G_{00,\sigma}$ to avoid the double counting of the diagram decorations on site 0 (see Fig. 2.4(c)). Such rescaling by $G_{00,\sigma}$ is necessary since in the decomposition into two subdiagrams the site 0 appears twice, hence all the processes starting and ending at site 0 ($G_{00,\sigma}$) are counted twice. Interestingly, this formula was already used by Hubbard [201]. $G_{00,\sigma}$ is nothing else but the local Green's function of the full lattice $G_{\sigma}$. Combining Eq. 2.40 with 2.37, one can show the relation between the local Green's function of the original problem and the Weiss field [43]:

$$G_{\sigma}^{-1}(i\omega_n) = G_{\sigma}^{-1}(i\omega_n) - \Sigma_{\sigma}(i\omega_n), \qquad (2.41)$$

where $\Sigma_{\sigma}$ is the purely local ($d \to \infty$) self-energy of the full lattice. Using the momentum-dependent Green's function definition, the Weiss field and the self-energy can be directly related:

$$G_{\sigma}^{-1}(i\omega_n) = \Sigma_{\sigma}(i\omega_n) + \sum_{\mathbf{k}} \frac{1}{i\omega_n + \mu - \epsilon_{\mathbf{k}} - \Sigma_{\sigma}(i\omega_n)}. \qquad (2.42)$$

Therefore, with Eqs. 2.38, 2.41 and 2.42, we have a set of self-consistent equations where the only unknown observable is the self-energy $\Sigma_{\sigma}$. The DMFT self-consistency cycle, illustrated in Fig. 2.5, is as follows:

(i) Start from an initial guess for the self-energy $\Sigma_{\sigma}$,

(ii) From Eq. 2.42 compute the corresponding Weiss field,

(iii) Construct the effective action with Eq. 2.38,

(iv) Solve the impurity problem for its Green's function $G_{imp,\sigma}$, which is identified to the local Green's function $G_{\sigma}(i\omega_n)$, and obtain a new self-energy from the Dyson Eq. 2.41,
<figure><img src="image_12.png" /><figcaption>Figure 2.5: Sketch of the DMFT self-consistent equations.</figcaption></figure>

(v) Check for self-consistency, and if needed, iterate again from (ii).

Therefore, until here we have shown how to parameterize the Hubbard model for the simple case of the single-band model. We derived the DMFT self-consistent equations, which provides a way to solve the Hubbard model at infinite dimension while keeping all the dynamic fluctuations. It is important to note that despite the local nature of the self-energy, the lattice Green's function remains momentum dependent via the dispersion relation $\epsilon_k$. In practice, DMFT relies on the approximation that $1/z$ is small, which for 3D systems can be justified but becomes questionable for 2D systems like the CuO planes of cuprates. In fact, for the 2D Hubbard model the non-local fluctuations are essential: even though DMFT can resolve the Mott insulator transition at half-filling since the self-energy is frequency dependent, it overestimates the critical interaction strength $U_c$ at which the transition occurs [48] compared to other methods partially including non-local fluctuations. In order to obtain a correct theoretical description of such systems, generalizations of DMFT to include a part of non-local interactions are therefore needed. One way is to start from the DMFT solution, which accounts well for the local correlations and perform a diagrammatic expansion around the DMFT result [71–73, 76]. In this thesis we rather rely on quantum cluster methods, which are based on generalizing the single-site impurity problem to a cluster impurity either containing multiple sites, or multiple patches in momentum space, and therefore allow to consider non-local self-energies on the length scale of the cluster [52].
## 2.3 BEYOND LOCAL CORRELATIONS: CLUSTER DYNAMICAL MEAN-FIELD THEORY

Remaining in the mindset of the cavity construction, a natural way to include non-local fluctuations would be to derive the effective action $S_{eff}$ not only for a single-site impurity, but for a cluster of sites. A first step towards this successful method, called cluster dynamical mean-field theory (CDMFT) [49–52], was first proposed by Lichtenstein and Katsnelson [49] who performed a single-site two-orbital LDA+DMFT calculation. Although not being strictly speaking a finite size cluster, the two-orbital system can be thought of as a dimer in real-space in terms of bonding/anti-bonding states. Most importantly, considering a two-orbital system requires to formulate the DMFT equations with matrix-valued functions, just as the "true" CDMFT method whose first implementation was realized shortly afterwards by Lichtenstein and Katsnelson in the study of the interplay between superconductivity and antiferromagnetic fluctuations in the hole-doped 2D Hubbard model using a 2 × 2 plaquette cluster [50]. Soon after, Kotliar et al. formally expressed CDMFT in terms of the cavity construction [51].

Translational symmetry requires to construct a superlattice made of the individual clusters, as shown in Fig. 2.6(a). To be consistent, the cluster geometry and tiling has to be chosen such that it reproduces the original lattice. Following the notations of the lecture notes by Sénéchal [202], the position of a site may be written as:

$$ \mathbf{r} = \mathbf{R} + \tilde{\mathbf{r}}, \qquad (2.43) $$

where $\tilde{\mathbf{r}}$ is the cluster position, and $\mathbf{R}$ the site position inside the cluster (see Fig. 2.6(a)-(b)). The local (matrix-valued) Green's function is obtained by:

$$ \mathbf{G}_{RR',\sigma}^{-1}(i\omega_n) = \sum_{\tilde{\mathbf{k}}} \left[ (i\omega_n + \mu - \epsilon_{\tilde{\mathbf{k}}} - \Sigma_{\sigma}(i\omega_n))^{-1} \right]_{RR'} \qquad (2.44) $$

where $G_{\sigma}$, $\Sigma_{\sigma}$, $\epsilon$ are matrix-valued quantities, $\tilde{\mathbf{k}}$ is defined in the reduced BZ of the superlattice (see Fig. 2.6(b)) and corresponds to the Fourier transform with respect to $\tilde{\mathbf{r}}$. The self-energy is local on the length-scale of the cluster, i.e. some of the diagonal terms $\Sigma_{ij,\sigma}$ are kept and explicitly calculated within the Anderson impurity model. All the DMFT self-consistent equations remain valid, but should be expressed in terms of matrix-valued functions. Therefore the same self-consistent scheme as the one illustrated in Fig. 2.5 can be used.

However, the division of the original lattice into finite-size clusters leads to a breaking of the original lattice translational symmetry, which can be understood easily by thinking of a 1D chain. Assuming a decomposition of the 1D chain in terms of dimers, say made out of the couples $\langle 2i, 2i + 1 \rangle$, then one immediately sees that $\Sigma_{2i,2i+1} \neq 0$ while $\Sigma_{2i+1,2i} = 0$ since the latter is not included in the cluster, as shown in Fig. 2.6(c). In other words, one bond is renormalized by the interactions, while the other is left bare (i.e. at the non-interacting level). One could thus spark a Peierls gap not necessarily relevant for the system under investigation. For the 2D Hubbard model for instance, it has been shown that it can induce biased CDW [203]. There exist
<figure><img src="image_13.png" /><figcaption>Figure 2.6: (a) Definition of the superlattice for the case of 2 × 2 plaquette clusters, the definition of the position vectors **r**, **R**, and ˜**r** (see text) are illustrated. (b) Sketch of the reduced BZ (red squares) on top of the original BZ (black square), and illustration of the momentum vectors **k**, **K**, and ˜**k**. (c) Illustration of the translation symmetry breaking inherent from the clusterization of the lattice.</figcaption></figure>

another successful extension of DMFT which was developed at the same time as CDMFT, called DCA [63–65], which we do not present here. In contrast to CDMFT, it preserves the translational symmetry by considering cluster defined directly in momentum space. However, this comes to the price of averaging the self-energy over patches in momentum space and is therefore a good approximation when the self-energy has a smooth momentum dependence. For phases like the PG in cuprates, which is studied in detail in Chap. 3, this can be a problem if too few patches are chosen since the self-energy has a strong momentum dependence. Moreover, in the study of the correlation effects on the sublattice decoupling in t-CuO, the real-space formulation of CDMFT allows to easily access to the inter- and intra-sublattice components of the self-energy. We therefore adopt the CDMFT framework in this thesis.

The direct consequence of the translational symmetry breaking is that the original lattice Green's function is no longer diagonal in momentum, i.e. $G \rightarrow G(\mathbf{k}, \mathbf{k}')$. Since by construction the superlattice translation invariance is preserved, then the only off-diagonal elements allowed are the ones sharing the same momentum $\tilde{\mathbf{k}}$ in the superlattice BZ but which differ by the reciprocal superlattice momentum $\mathbf{K}$:

$$G(\mathbf{k}, \mathbf{k}', i\omega_n) = G(\tilde{\mathbf{k}} + \mathbf{K}, \tilde{\mathbf{k}} + \mathbf{K}', i\omega_n), \quad (2.45)$$

where $G$ is the (scalar) Green's function of the original lattice. This can be directly seen by making the link between the latter $G(\tilde{\mathbf{k}} + \mathbf{K}, \tilde{\mathbf{k}} + \mathbf{K}', i\omega_n)$ and the matrix-valued Green's function on the superlattice $G_{RR'}(\tilde{\mathbf{k}}, i\omega_n)$. To this aim, we can use the relation [202]:

$$c(\tilde{\mathbf{k}} + \mathbf{K}) = \frac{1}{L} \sum_{\mathbf{K}', \mathbf{R}} e^{-i\mathbf{R}(\tilde{\mathbf{k}} + \mathbf{K} - \mathbf{K}')}(\mathbf{c}_{\mathbf{K}'}(\tilde{\mathbf{k}})), \quad (2.46)$$
where L is the number of cluster sites, $c(\tilde{\mathbf{k}} + \mathbf{K})$ is the annihilation operator for an electron with momentum $\tilde{\mathbf{k}} + \mathbf{K}$, and $c_{\mathbf{K}'}(\tilde{\mathbf{k}})$ the matrix element $\mathbf{K}'$ of the cluster annihilation operator $c(\tilde{\mathbf{k}})$. It follows that the original lattice Green's function can be re-written:

$$ 
\begin{aligned}
G(\tilde{\mathbf{k}} + \mathbf{K}, \tilde{\mathbf{k}} + \mathbf{K}', i\omega_n) &= \frac{1}{L^2} \sum_{\mathbf{K}_1\mathbf{K}_2, \mathbf{RR}'} e^{-i\mathbf{R}(\tilde{\mathbf{k}}+\mathbf{K}-\mathbf{K}_1)} e^{-i\mathbf{R}'(\tilde{\mathbf{k}}+\mathbf{K}'-\mathbf{K}_2)} G_{\mathbf{K}_1\mathbf{K}_2}(\tilde{\mathbf{k}}, i\omega_n). \\
&= \frac{1}{L} \sum_{\mathbf{RR}'} e^{-i\mathbf{R}(\tilde{\mathbf{k}}+\mathbf{K})} e^{-i\mathbf{R}'(\tilde{\mathbf{k}}+\mathbf{K}')} G_{\mathbf{RR}'}(\tilde{\mathbf{k}}, i\omega_n).
\end{aligned} 
\quad (2.47) $$

Without translational symmetry the relation can not be further simplified. However, since the DOS and the spectral function that we compare to ARPES involves only the diagonal part of the Green's function (see App. A), the so-called periodization procedure [67, 202], in which only the diagonal elements of the lattice Green's function are kept, is justified although still being an approximation. In Eq. 2.47, this amounts to pick the terms with $\mathbf{K} = \mathbf{K}'$:

$$ G(\mathbf{k}, i\omega_n) = G(\tilde{\mathbf{k}} + \mathbf{K}, i\omega_n) = \frac{1}{L} \sum_{\mathbf{RR}'} e^{-i(\mathbf{R}-\mathbf{R}')\tilde{\mathbf{k}}} G_{\mathbf{RR}'}(\tilde{\mathbf{k}}, i\omega_n). \quad (2.48) $$

This equation allows, at the price of neglecting the off-diagonal elements, to relate the translationally invariant lattice Green's function to the super-lattice matrix-valued one.

For a system without translation symmetry breaking this is an efficient way to overcome the CDMFT problem, however it can be questionable when there is, for instance, a long range magnetic ordering, or a CDW. Indeed, Eq. 2.48 restores the translational symmetry in the original lattice and would ignore the possible existence of, say, an antiferromagnetic order. Later in this thesis, we encounter such a situation when studying t-CuO, for which spin polarized CDMFT calculations are performed, and Eq. 2.48 is used to obtain the spectral function that is compared to ARPES. Although t-CuO has been shown to display a long-range antiferromagnetic stripe order [135, 136], the theoretical spectral function obtained with the above periodization procedure is in excellent agreement with the photoemission measurements [134]. The main reason behind is that in photoemission, if the incoming light spot size is larger than the characteristic magnetic correlation length, the measurement averages over the different magnetic configurations possible, similarly to the periodization procedure. This is true for the spectral function, but might be a limitation for other observables, and in such a case a refinement of the procedure to include the eventual symmetry breaking may be needed.

Finally, we presented here the periodization of the Green's function, but it is not the only possibility to restore the translational symmetry. Indeed, since the origin of the problem is that the inter-site self-energy is taken into account for the intra-cluster bonds only, one way would be to
periodize the self-energy instead [51]. Yet another proposition is to periodize the cumulant [204]. The three different procedures can be summarized as follows [203]:

$$ 
\begin{aligned}
G_{\sigma}(\mathbf{k}, i\omega_n) &= \frac{1}{L} \sum_{RR'} e^{-i(\mathbf{R}-\mathbf{R}')\mathbf{k}} [(i\omega_n + \mu - \epsilon(\mathbf{k}) - \Sigma_{\sigma}(i\omega_n))^{-1}]_{RR'}, \\
G_{\sigma}^{\Sigma}(\mathbf{k}, i\omega_n) &= \left( i\omega_n + \mu - \epsilon(\mathbf{k}) - \frac{1}{L} \sum_{RR'} e^{-i(\mathbf{R}-\mathbf{R}')\mathbf{k}} [\Sigma_{\sigma}(i\omega_n)]_{RR'} \right)^{-1}, \quad (2.49) \\
G_{\sigma}^{M}(\mathbf{k}, i\omega_n) &= \left( \left[ \frac{1}{L} \sum_{RR'} e^{-i(\mathbf{R}-\mathbf{R}')\mathbf{k}} [(i\omega_n + \mu - \Sigma_{\sigma}(i\omega_n))^{-1}]_{RR'} \right]^{-1} - \epsilon(\mathbf{k}) \right)^{-1},
\end{aligned}
$$

where $G^{\Sigma}$ and $G^{M}$ denote the reperiodized Green's function obtained respectively with the self-energy and the cumulant periodization. The periodized Green's function, self-energy, or cumulant preserve some features of the symmetry breaking in that they are all obtained by neglecting the off-diagonal momentum elements. The main difference comes from their connection to the real Green's function of the lattice $G(\tilde{\mathbf{k}}+\mathbf{K}, \tilde{\mathbf{k}}+\mathbf{K}', i\omega_n)$. When periodizing the Green's function, there is a direct link since the periodized Green's function at $\mathbf{k}$ is simply the diagonal elements of the lattice Green's function at $(\mathbf{k}, \mathbf{k})$. However, if instead the self-energy, or the cumulant are periodized, then the off-diagonal elements are neglected before performing the matrix inversion, leaving the relation between the periodized Green's function and the lattice Green's function unclear. Since in this thesis we mostly use the periodization procedure to study the momentum resolved spectral function, i.e. the imaginary part of the lattice Green's function, we will rather adopt the Green's function periodization (Eq. 2.48).

So far we have presented the theoretical framework of DMFT and CDMFT, as well as the main set of self-consistent equations. Although the self-consistent loop is simple and imply seemingly straightforward relations between the self-energy, Green's function, and Weiss field, in practice the implementation of CDMFT requires a few tricks necessary not to ruin the final result. In the next section, we briefly describe some tips that may be useful for readers interested in the implementation of the CDMFT self-consistent scheme.

## 2.4 CLUSTER DYNAMICAL MEAN FIELD THEORY: IMPLEMENTATION DETAILS

We review here some of the important details one faces when implementing a CDMFT code. Note that we do not consider here the impurity solver part, but only the set of self-consistent equation, given a solver that already works. We first present the problem of Fourier transforming functions from the Matsubara frequency axis to imaginary time, then we discuss the choice of the starting point for the CDMFT self-consistent calculations, as well as how to chose the chemical potential. Finally we present a simple way to switch to an advantageous basis, which is necessary for the hybridization-expansion [116, 119] CTQMC [60–62] solvers in order to reduce the sign problem.
<figure><img src="image_14.png" /><figcaption>Figure 2.7: Fourier transform of the Green’s function defined in Eq. 2.52 using different number of Matsubara frequencies, with and without the high-frequency expansion correction, calculated at $\beta = 40\text{eV}^{-1}$. The number of $\tau$ points is always $N_{\tau} = 3 \times N_{\omega}$.</figcaption></figure>

### 2.4.1 FOURIER TRANSFORMATION

If the impurity solver used requires the hybridization function as a function of imaginary time, one faces the problem of Fourier transforming a function from the Matsubara to the imaginary time axis. More generally, the filling can also be obtained from the Green’s function at the $\tau \to 0^-$ limit:

$$ 
\begin{aligned}
G_{ii}(\tau \to 0^{-}) &= \lim_{\tau \to 0^{-}} -\langle T_{\tau} c_{i}(\tau) c_{i}^{\dagger}(0) \rangle \\
&= \langle c_{i}^{\dagger}(0) c_{i}(0) \rangle \\
&= n_{i},
\end{aligned} 
\quad (2.50) $$

where $n$ denotes the filling, which is necessarily checked during the self-consistent cycle.

The Fourier transformation reads

$$ G_{ij}(\tau) = \frac{1}{\beta} \sum_{n=-\infty}^{+\infty} e^{-i\tau\omega_n} G_{ij}(i\omega_n), \quad (2.51) $$

and implies a sum over an infinite number of Matsubara frequencies. The problem is that we have to choose a certain cutoff $N$ for the number of Matsubara frequencies, but the diagonal components of the Green’s function decay slowly $G_{ii}(i\omega_n) \propto \frac{1}{i\omega_n}$. A large number of frequency $N$ would then be required to improve the precision, but at the same time one wants to limit as much as possible the number of frequency to keep the impurity solver computational complexity reasonable. In fact, as shown in Fig. 2.7 for the simple Green’s function in the Hubbard-I approximation [35]

$$ G(i\omega_n) = \frac{1}{i\omega_n - \frac{U^2}{4i\omega_n}}, \quad (2.52) $$

at $\beta = 40\text{eV}^{-1}$, even with a relatively large number of frequencies and $\tau$ points, a simple truncation of Eq. 2.51 leads to substantial errors especially around $\tau \approx 0$, which is an important region in $\tau$ space since this is where the filling can be extracted for instance.

A relatively simple solution is to use the high-frequency tail expansion of the Green’s function,
that will allow to approximate the Green's function beyond the truncated frequencies. The expansion of a component of the Green's function is defined by

$$G_{ij}(i\omega_n) \approx M_0 + \frac{M_1}{i\omega_n} + \frac{M_2}{(i\omega_n)^2} + \frac{M_3}{(i\omega_n)^3} + \frac{M_4}{(i\omega_n)^4}, \quad (2.53)$$

which we truncated at the fourth order and is valid at large $\omega_n$. The coefficients $M_0, \dots, M_4$ are easily found using regular fitting methods on the last 20 to 30% Matsubara frequencies. The first coefficient $M_0$ has to be 0, and for a diagonal component the coefficient $M_1$ should be equal to 1. The later condition can be used to check the accuracy of the high-frequency expansion fit.

The Fourier transform may now be written

$$\begin{aligned} G_{ij}(\tau) = & \frac{1}{\beta} \sum_{n=-N}^{N} e^{-i\tau\omega_n} G_{ij}(i\omega_n) + \sum_{n=[-\infty, -N[}^{[N, +\infty]}} e^{-i\tau\omega_n} G_{ij}(i\omega_n) \\ \approx & \frac{1}{\beta} \sum_{n=-N}^{N} e^{-i\tau\omega_n} \left( G_{ij}(i\omega_n) - \frac{M_1}{i\omega_n} - \frac{M_2}{(i\omega_n)^2} - \frac{M_3}{(i\omega_n)^3} - \frac{M_4}{(i\omega_n)^4} \right) \quad (2.54) \\ & + \frac{1}{\beta} \sum_{n=-\infty}^{\infty} e^{-i\tau\omega_n} \left( \frac{M_1}{i\omega_n} + \frac{M_2}{(i\omega_n)^2} + \frac{M_3}{(i\omega_n)^3} + \frac{M_4}{(i\omega_n)^4} \right). \end{aligned}$$

The great advantage of this reformulation is that the second term, which implies a sum over an infinite number of Matsubara frequencies, can be solved analytically (the detailed derivation can be found in the excellent Ref. [205]):

$$\begin{aligned} & \frac{1}{\beta} \sum_{n=-\infty}^{\infty} e^{-i\tau\omega_n} \frac{M_1}{i\omega_n} = -\frac{M_1}{2}\mathrm{sign}(\tau) \\ & \frac{1}{\beta} \sum_{n=-\infty}^{\infty} e^{-i\tau\omega_n} \frac{M_2}{(i\omega_n)^2} = \frac{M_2}{4}(2\tau - \beta) \\ & \frac{1}{\beta} \sum_{n=-\infty}^{\infty} e^{-i\tau\omega_n} \frac{M_3}{(i\omega_n)^3} = \frac{M_3}{4}\tau(\beta - \tau)\mathrm{sign}(\tau) \\ & \frac{1}{\beta} \sum_{n=-\infty}^{\infty} e^{-i\tau\omega_n} \frac{M_4}{(i\omega_n)^4} = \frac{M_4}{48}(2\tau - \beta)(2\tau^2 - 2\tau\beta - \beta^2). \end{aligned} \quad (2.55)$$

By inserting these relations in Eq. 2.54, and using the symmetry $G_{ii}(-i\omega_n) = G_{ii}^*(i\omega_n)$ which holds for the diagonal components of the Green's function, one obtains:

$$\begin{aligned} G_{ii}(\tau) = & \frac{2}{\beta} \sum_{n=0}^{N} \left( \mathrm{Re}[G_{ii}(i\omega_n)] \cos(\omega_n \tau) + \mathrm{Im}[G_{ii}(i\omega_n)] \sin(\omega_n \tau) \right. \\ & \left. + M_1 \frac{\sin(\omega_n \tau)}{\omega_n} + M_2 \frac{\cos(\omega_n \tau)}{\omega_n^2} - M_3 \frac{\sin(\omega_n \tau)}{\omega_n^3} - M_4 \frac{\cos(\omega_n \tau)}{\omega_n^4} \right) \\ & - \frac{M_1}{2}\mathrm{sign}(\tau) + \frac{M_2}{4}(2\tau - \beta) + \frac{M_3}{4}\tau(\beta - \tau)\mathrm{sign}(\tau) + \frac{M_4}{48}(2\tau - \beta)(2\tau^2 - 2\tau\beta - \beta^2). \end{aligned} \quad (2.56)$$

This correction drastically improves the Fourier transform, as is shown in Fig. 2.7: already at $N_\omega = 100$ the result is smooth close to $\tau \approx 0$ and not noisy as compared to the direct transformation.
mation. The high-frequency expansion is therefore a precious trick which allows to significantly reduce the number of Matsubara frequencies. Most importantly, if we showed here an example for the diagonal components of the Green's function, this formalism applies also to the hybridization function, and in fact to other Matsubara functions as long as the first coefficient $M_0$ is zero. Note that in general the off-diagonal components decay faster at high frequency than the diagonal ones, such that in practice for these elements we use the direct transformation (Eq. 2.51).

Similarly, the electron filling at some site/orbital can be determined from the corresponding diagonal component of the Green's function. Indeed, at $\tau \to 0^-$ Eq. 2.56 reduces to:

$$G_{ii}(\tau = 0^{-}) = n_i = \frac{M_1}{2} - \frac{M_2}{4}\beta + \frac{M_4}{48}\beta^3 + \frac{2}{\beta}\sum_{n=0}^{N}\left(\mathrm{Re}[G_{ii}(i\omega_n)] + \frac{M_2}{\omega_n^2} - \frac{M_4}{\omega_n^4}\right). \quad (2.57)$$

This relation is extremely useful in practice, especially when the chemical potential has to be adapted during the self-consistency to ensure that the system remains at the desired filling, as we discuss in the following.

### 2.4.2 INITIALIZATION AND CHOICE OF CHEMICAL POTENTIAL

At half-filling, the U-T phase diagram of the 2D one-band Hubbard model is characterized by a first order metal to insulator transition at low temperature, where a coexistence region is found in which both the metallic and insulating phases are possible solutions [206]. Therefore, in general one must be careful with the initialization of the CDMFT self-consistent loop. An example of possible consequences on the spectral function is discussed later in this thesis, in the study of NaCCOC in Chap. 3.

In our implementation, we find it convenient to initialize the calculation with the self-energy. There is an infinite number of possible choices, but we present here the three we used within this thesis. To start from a metallic solution, the most simple choice is to put all components of the self-energy to zero, such that the initial lattice Green's function is the non-interacting one. For an insulating starting point, the diagonal components for paramagnetic calculations are defined from the Hubbard-I approximation [35]

$$\Sigma_{ii,\sigma} = \frac{U^2}{4i\omega_n}, \quad (2.58)$$

and for spin polarized calculations a spin dependent constant self-energy is used

$$\Sigma_{ii,\sigma} = U/2 \qquad (2.59)$$
$$\Sigma_{ii,-\sigma} = -U/2.$$

In both cases, the off-diagonal components are initialized to zero.

Since the chemical potential $\mu$ determines the overall filling, to which the properties of the Hubbard model are extremely sensitive, it is crucial to precisely keep track of it during the calculation.
tions. In a metallic system, there is no ambiguity on the value of $\mu$, and the latter has to be adapted during the self-consistent cycle to make sure that the system remains at the desired filling. To do so, it is convenient to use the impurity self-energy and calculate the local Green's function at a given $\mu$ via

$$G_{\sigma}(i\omega_n) = \sum_{\mathbf{k}} \frac{1}{i\omega_n + \mu - \epsilon_{\mathbf{k}} - \Sigma_{\sigma}(i\omega_n)} \quad (2.60)$$

and compute the filling using Eq. 2.57. The right $\mu$ can be found using a regula-falsi method: it is determined iteratively using

$$\mu^{i+1} = \frac{\mu^{i-1}(n^i - n) - \mu^i(n^{i-1} - n)}{n^i - n^{i+1}}, \quad (2.61)$$

until $n^i$ corresponds to the desired filling $n$.

If the result is expected to be insulating, this method might lead to spurious jumps of $\mu$ since it can take any value inside the gap without changing the filling. It is therefore convenient to fix the chemical potential somewhere in the gap, for instance in the middle. As a practical example, when using the insulating spin-polarized initialization (Eq. 2.59), a natural choice for the chemical potential would be $\mu = 0$.

### 2.4.3 BASIS TRANSFORMATION

Having described general concepts of the implementation, we now turn to a solver specific aspect: the Monte-Carlo sign problem can be greatly reduced in the hybridization-expansion CTQMC solvers when transforming the cluster basis into an advantageous one in which the hybridization function is as diagonal as possible. For specific clusters like the dimer, or the 2 × 2 plaquette, a diagonal basis can be easily found, but for more general cluster such a basis can be difficult to find. Hence it is more convenient to find a suitable basis on the fly.

In order to determine it, we diagonalize a constant quantity, instead of the hybridization function at a given frequency. We chose the on-site energy matrix, which can be obtained from the first coefficient $M_0$ of the high-frequency expansion for each of the component of the inverse of the Weiss field $G_{\sigma}^{-1}$, and may be understood in terms of orbital-dependent chemical potential. This quantity is then diagonalized and we obtain the transfer matrix $\mathbf{P}$, from which the hybridization function can be "diagonalized" at each frequency

$$\tilde{\Delta}(i\omega_n) = \mathbf{P}^{-1}\Delta(i\omega_n)\mathbf{P}. \quad (2.62)$$

The interaction tensor $\mathbf{V}_{i,j,k,l}$ must be rotated accordingly. The transformation to the new basis set may be written:

$$\tilde{\mathbf{V}}_{a,b,c,d} = \sum_{i,j,k,l} \mathbf{P}_{a,i}^{-1} \mathbf{P}_{b,j}^{-1} \mathbf{V}_{i,j,k,l} \mathbf{P}_{k,c} \mathbf{P}_{l,d} \quad (2.63)$$
where a, b, c, d are the indices defined in the new diagonal basis. Defining the super-indices X = (a, b), Y = (c, d), Z = (i, j) and Z' = (k, l), the problem reduces to a simple matrix product:

$$
\begin{align}
\tilde{\mathbf{V}}_{a,b,c,d} &= \sum_{i,j,k,l} \mathbf{P}_{a,i}^{-1} \mathbf{P}_{b,j}^{-1} \mathbf{V}_{i,j,k,l} \mathbf{P}_{k,c} \mathbf{P}_{l,d} \nonumber \\
\tilde{\mathbf{V}}_{(a,b)(c,d)} &= \sum_{(i,j)(k,l)} \mathbf{T}_{(a,b)(i,j)}^{-1} \mathbf{V}_{(i,j)(k,l)} \mathbf{T}_{(k,l)(c,d)} \tag{2.64} \\
\tilde{\mathbf{V}}_{XY} &= \sum_{Z,Z'} \mathbf{T}_{XZ}^{-1} \mathbf{V}_{ZZ'} \mathbf{T}_{Z'Y} \nonumber \\
\tilde{\mathbf{V}}_{XY} &= (\mathbf{T}^{-1} \mathbf{V} \mathbf{T})_{XY} \nonumber
\end{align}
$$

Note that this derivation is general and does not imply any assumption on the form of the intrac-
tion tensor.

This concludes the presentation of DMFT and its cluster extension CDMFT. The latter allows
to include the non-local fluctuations on short scales, and will prove extremely useful in the study
of the undoped and doped phases of NaCCOC and t-CuO. All the formalism has been written
in Matsubara imaginary time/frequency, but in principle all the CDMFT equations are valid on
the real frequency axis. We chose the Matsubara formulation since the impurity solvers used
during this thesis were based on CTQMC [60–62] techniques, which work at finite temperature
on the imaginary frequency axis. We use the hybridization-expansion CT-HYB [116] solver, and
the interaction-expansion CT-INT [117], based on the ALPSCore libraries [118]. We do not detail
here how these solvers work, and instead refer the reader to the review by Gull et al. [62], except
for the fact that they require the Matsubara imaginary time/frequency formalism. Therefore, in
practice we face the problem of obtaining the real-frequency data (that can be compared with the
experiments) only from the knowledge on the imaginary frequency axis. This is the subject of
the next section.

2.5 ANALYTIC CONTINUATION

The purpose of the analytic continuation is to obtain the value of a function on the real axis given
its value on the imaginary axis (frequency, or time), which amounts to:

$$
G(i\omega_n) \rightarrow G(\omega \pm i\delta) \tag{2.65}
$$

where ω denotes real frequencies, and δ is a small parameter used in practice to prevent the dis-
continuities of the function G on the real axis. As we first briefly show, such a transformation
could be done directly by using the Kramers-Kroning relations [207, 208], but is extremely sen-
sitive to noise and is in-fact ill defined. To illustrate that, for the sake of simplicity we consider
only the diagonal component of the Green’s function for which the spectral function is always
positive:

$$
A(\omega) = -\frac{1}{\pi} \mathrm{Im}[G(i\omega_n)] \ge 0. \qquad (2.66)
$$
<figure><img src="image_15.png" /><figcaption>Figure 2.8: Analytic continuation of the Green’s function $G(i\omega_n) = \frac{1}{i\omega_n}$ using (a) the direct inversion without noise, (b) the direct inversion with the addition of a random noise of $10^{-5}$ amplitude, (c) the MEM with and without noise. The inset in (c) is a zoom around $\omega = 0$.</figcaption></figure>

In principle, one can formally relate the Green’s function on the Matsubara axis to its real axis counterpart using Kramers-Kroning relations [207, 208]:

$$ 
\begin{aligned}
G(i\omega_n) &= \frac{1}{\pi} \int_{-\infty}^{+\infty} \frac{\mathrm{Im}[G_{\mathrm{real}}(\omega)]}{\omega - i\omega_n} d\omega, \\
&= \frac{1}{\pi i} \int_{-\infty}^{+\infty} \frac{\mathrm{Re}[G_{\mathrm{real}}(\omega)]}{\omega - i\omega_n} d\omega,
\end{aligned} 
\quad (2.67) 
$$

where $G_{\mathrm{real}}$ denotes the Green’s function on the real axis, and $G(i\omega_n)$ is restricted to $i\omega_n$ in the upper complex plane. However, usually we need to obtain $G_{\mathrm{real}}(\omega)$ from $G(i\omega_n)$: Eq. 2.67 must be inverted. To this aim, the integral may be discretized on a frequency mesh $[\omega_1, \dots, \omega_N]$, with N the number of real frequencies:

$$ 
\begin{aligned}
G(i\omega_n) &= \lim_{N \to \infty} \sum_{l=1}^{N} \frac{1}{\pi} \frac{\omega_{l+1} - \omega_l}{\omega_l - i\omega_n} \mathrm{Im}[G_{\mathrm{real}}(\omega_l)] \\
[\mathbf{G}]_n &= [\mathbf{M} \cdot \mathbf{G}_{\mathrm{real}}]_n,
\end{aligned} 
\quad (2.68) 
$$

where the Green’s functions are considered as vectors and the matrix $\mathbf{M}$ is defined by $\mathbf{M}_{ln} = \frac{1}{\pi} \frac{\omega_{l+1} - \omega_l}{\omega_l - i\omega_n}$. Therefore, by inverting this matrix equality, one directly obtains the Green’s function on the real axis with the pseudo-inverse of the matrix $\mathbf{M}_{ln}$.

In Fig. 2.8(a) is shown the result obtained from the direct inversion for a simple function $G(i\omega_n) = \frac{1}{i\omega_n}$, which should transform to a delta function on the real frequency axis. We notice many spurious features, although the main delta peak is reproduced. This is linked to the fact that the matrix elements $\mathbf{M}_{ln} = \frac{1}{\pi} \frac{\omega_{l+1}-\omega_l}{\omega_l-i\omega_n}$ vanish for n large, such that the inverse of the matrix $\mathbf{M}^{-1}$ is almost irregular and ill-conditioned. When inverted, it becomes extremely sensitive to small variations in the input $G(i\omega_n)$, which can induce large numerical errors as seen in Fig. 2.8(a). Moreover, in practice we need to perform an analytic continuation of Monte Carlo data, which is noisy and incomplete. The order of magnitude of the noise can be estimated from $\frac{1}{\sqrt{N_{sample}}}$. For a Green’s function obtained with $N_{sample} = 10^7$, a quite high accuracy for impurity solvers, the sampling noise can be estimated to $\sim 3 \cdot 10^{-4}$, and it would lead to an extremely noisy result using the direct inversion, as shown in Fig. 2.8(b) where a small noise of amplitude $10^{-5}$ was added to
$$G(i\omega_n) = \frac{1}{i\omega_n}.$$

Because of the vanishing elements $M_{ln} = \frac{1}{\pi} \frac{\omega_{l+1}-\omega_l}{\omega_l-i\omega_n}$, the direct inversion is an ill-defined problem. Several methods have been proposed to overcome this issue, such as the Padé approximant [209, 210], the stochastic method [211], or the Nevanlina approach [212]. In this thesis we used methods based on the maximum entropy method (MEM) [213, 214], which we present in the following. This method was first developed for only the diagonal components of the Green's function since it is based on the positivity of the spectral function $A(\omega) = -\frac{1}{\pi}\text{Im}[G_{\text{real}}(\omega)]$, such that it is not fully suited for the matrix-valued Green's function. Indeed, the off-diagonal elements are in general not fully positive, and they require a special treatment. Two different strategies to deal with the off-diagonal elements are presented. The first is the "Poor Man's" matrix procedure, in which the off-diagonal elements are treated apart from the diagonal ones [215, 216]. The second, called maximum quantum entropy method (MQEM), consists in generalizing the MEM to treat directly the matrix-valued function as a whole [217].

### 2.5.1 MAXIMUM ENTROPY METHOD

A widely used method to perform analytic continuation is the maximum entropy method (MEM) [213, 214], which attempts to overcome the problem of fitting noisy functions on the imaginary axis by regularizing the trial spectral functions with a smooth default model. It is first convenient to re-write Eq. 2.67 in imaginary time:

$$ 
\begin{align}
G(\tau) &= \frac{1}{\beta} \sum_{-\infty}^{+\infty} e^{-i\omega_n \tau} G(i\omega_n) \nonumber \\
&= \frac{1}{\pi} \int_{-\infty}^{+\infty} \text{Im}[G_{\text{real}}(\omega)] \left[ \frac{1}{\beta} \sum_{-\infty}^{+\infty} \frac{e^{-i\omega_n \tau}}{\omega - i\omega_n} \right] d\omega \nonumber \\
&= \frac{1}{\pi} \int_{-\infty}^{+\infty} \text{Im}[G_{\text{real}}(\omega)] \frac{e^{-\tau\omega}}{e^{-\beta\omega} + 1} d\omega \tag{2.69} \\
&= - \int_{-\infty}^{+\infty} A(\omega) \frac{e^{-\tau\omega}}{e^{-\beta\omega} + 1} d\omega. \nonumber
\end{align}
$$

We restrict the discussion here to the diagonal components of the Green's function, for which $A(\omega)$ is purely real and positive. From the knowledge of a spectral function, the Green's function on the imaginary time axis can therefore be calculated by integrating the spectral function weighted by the kernel $K(\tau, \omega) = -\frac{e^{-\tau\omega}}{e^{-\beta\omega}+1}$. This is called the back-continuation, and is not ill-conditioned contrary to the direct inversion presented previously.

Assuming that we have a measured Green's function $\bar{G}(\tau)$ from a Monte-Carlo calculation, and
a spectral function $A(\omega)$ supposed to reproduce the measured data, the goodness of the fit can be evaluated with the likelihood function

$$ \begin{aligned} \chi^2 &= \int_0^\beta \frac{1}{\sigma^2(\tau)} \left| \bar{G}(\tau) - \int_{-\infty}^{+\infty} A(\omega)K(\tau, \omega)d\omega \right|^2 d\tau \\ &= \int_0^\beta \frac{1}{\sigma^2(\tau)} \left| \bar{G}(\tau) - G(\tau)d\omega \right|^2 d\tau, \end{aligned} \quad (2.70) $$

where $\sigma(\tau)$ is a measure of the noise, and $G(\tau)$ is the back-continued Green's function from $A(\omega)$. Using the likelihood function, an optimal spectral function can be searched via a least-square method, but the noise in the measured data will lead to nonphysical results and over-fitting.

In order to avoid that, the proposed spectral function can be simultaneously optimized and regularized by exploiting the Bayesian interference. As mentioned earlier, the spectral function is always positive and is normalized to 1, so that it can be thought of as a probability density. From Baye's theorem we have that for two events a, b, the probability of a given b $P[a|b]$ is:

$$ P[a|b] = \frac{P[b|a]P[a]}{P[b]}. \quad (2.71) $$

For the analytic continuation, the event a corresponds to the spectral function $A(\omega)$, and b to the measured Green's function $\bar{G}(\tau)$. $P[A|\bar{G}]$ can therefore be interpreted as the probability that $A(\omega)$ is the true spectral function given the measured $\bar{G}(\tau)$. Similarly, $P[\bar{G}|A]$ is the probability that $\bar{G}(\tau)$ is the right measured function given the spectral function, i.e. it is related to the likelihood function $\chi^2$. Taking $P[\bar{G}]$ as a constant that can not be improved within the continuation process, we have:

$$ P[A|\bar{G}] \propto P[\bar{G}|A]P[A]. \quad (2.72) $$

The MEM consists in maximizing the probability $P[A|\bar{G}]$ by searching for a spectral function that maximized the product $P[\bar{G}|A]P[A]$.

$P[\bar{G}|A]$ is connected to the likelihood function, and it can be shown that [213, 214]

$$ P[\bar{G}|A] \propto e^{-\frac{\chi^2}{2}}, \quad (2.73) $$

so that $P[\bar{G}|A]$ is maximized by minimizing the likelihood function. As stated above, focusing only on this procedure would lead to nonphysical results due to the noise. In order to restrict the possible spectral functions, a smooth default model $D(\omega)$ can be defined and used as a reference for the spectral function by defining the following entropy:

$$ S = - \int A(\omega) \ln\left(\frac{A(\omega)}{D(\omega)}\right) d\omega. \quad (2.74) $$
This is a useful quantity in the sense that it measures how "far" is the spectral function from the noise-free default model which is constructed from the prior knowledge of the expected characteristics of the spectral function on the real axis. $P[A]$ is related to this entropy via [213, 214]

$$P[A] \propto e^{\alpha S}, \qquad (2.75)$$

where $\alpha$ is a constant introduced to weight $P[A]$ with respect to $P[\bar{G}|A]$, i.e. $\alpha$ controls how much the final spectral function should resemble the default model since the entropy is maximum when $A(\omega) = D(\omega)$. Finally, we have

$$P[A|\bar{G}] \propto e^{-\frac{\chi^2}{2} + \alpha S}, \qquad (2.76)$$

such that the search of the best $A(\omega)$ given the measured $\bar{G}(\tau)$ consists of a combined minimization of the likelihood function $\chi^2$ and maximization of the entropy $S$.

It is clear from Eq. 2.76 that the final result of the probability maximization will depend on the value of $\alpha$, with the two limiting cases being:

(i) $\alpha = 0$, the procedure only relies on the least-square fitting procedure and is noise sensitive,

(ii) $\alpha \to \infty$, the spectral function is equal to the default model.

Different methods have been proposed to choose the optimal $\alpha$: the classic MEM selects the $\alpha$ which maximizes the probability $P[\alpha|\bar{G}]$ [214], or the Bryan's method which consists of a weighted average over different spectra calculated for a range of $\alpha$ [218]. Recently, another approach has been proposed in which the likelihood function $\chi^2[\alpha]$ is calculated for a range of $\alpha$, and the optimal choice is determined at the maximum curvature of $\log \chi^2[\log \alpha]$ [219].

In Fig. 2.8(c) is shown the resulting spectral function for the Matsubara Green's function $G(i\omega_n) = \frac{1}{i\omega_n}$ obtained with MEM using the MaxEnt code [220] with the Bryan's method [218]. We use a default flat model, and the noise $\sigma(\tau)$ was set to $10^{-5}$. A clear improvement is obtained compared to the direct inversion method, especially when inserting a random numerical noise in the "measured" Matsubara Green's function. This is, of course, an extremely simple function and for more involved functions the result of the MEM can be (severely) parameter dependent: the choice of the $\alpha$ range, of how to choose the optimal $\alpha$, which default model is used, etc. Nevertheless, the analytic continuation remains mandatory when using Quantum Monte Carlo methods, and in this context MEM brought a valuable improvement.

Note that the MEM can also be applied to the self-energy, which is convenient to directly obtain the momentum-resolved Green's function, or the periodized Green's function on the real frequency axis. To this aim, the self-energy has to be transformed into a Green's function-like quantity, and again here we restrict the discussion to the diagonal components. A simple approach can be to simply correct the high-frequency tail of the self-energy so that it has the correct $\frac{1}{i\omega_n}$ behavior:

$$\Sigma^{\text{aux}}(i\omega_n) = \frac{\Sigma(i\omega_n) - \Sigma(\infty)}{M_1}, \qquad (2.77)$$
where $\Sigma(\infty)$ is the asymptotic value of the self-energy at $\omega_n \to \infty$, and $M_1$ is the first coefficient of the high frequency expansion:

$$ \Sigma(i\omega_n) \approx \Sigma(\infty) + \frac{M_1}{i\omega_n} + \frac{M_2}{(i\omega_n)^2} + \dots \qquad (2.78) $$

$\Sigma^{aux}(i\omega_n)$ behaves like a Green's function and can therefore be continued using the MEM. The true self-energy is recovered by transforming back directly on the real frequency axis with the coefficients $M_1$ and $\Sigma(\infty)$. Another method that we have used for the spectral function of t-CuO is to construct an auxiliary Green's function from the self-energy:

$$ G^{aux} = \frac{1}{i\omega_n - \Sigma(i\omega_n) + C}, \qquad (2.79) $$

where $C$ is a constant shift. The analytic continuation is performed for $G^{aux}$ and the self-energy is recovered by inverting the previous equation directly on the real axis.

### 2.5.2 DEALING WITH MATRIX-VALUED FUNCTIONS: ELEMENT-WISE CONTINUATION

When using CDMFT we face the problem that not only the diagonal elements, but also the off-diagonal ones have to be analytically continued to the real axis. The MEM procedure described above can not straightforwardly be applied to the off-diagonal components since, contrary to their diagonal counterparts, they are not fully positive, preventing the use of the entropy in Eq. 2.74. One possible way is to transform the Green's function (or self-energy) to a favorable basis set in which it is almost diagonal and in which the off-diagonal elements can be neglected. The diagonal elements are continued separately, and the function can be rotated back. However, because of the systematic error of the MEM, the result becomes at least basis-dependent and it can even lead to causality problems when rotated back to the original basis (the rotation matrix is no longer suited for the continued quantity). Another strategy is based on the mixing of the off-diagonal elements with the diagonal ones already continued, so to create new fictitious positive spectral functions [221]. The limitation comes from the fact that already continued spectral functions enter the calculation for the off-diagonal components. In this thesis, we have used two different ways to go beyond these two procedures, a first in which each element of the matrix are treated separately with a modified entropy and specific default model for the off-diagonal elements, and a second one in which the entropy is redefined to perform the continuation on the whole matrix.

The "Poor Man's" matrix procedure [215, 216] is implemented within the TRIQS package [120]. The diagonal components are treated separately with the usual MEM described above. In order to treat the off-diagonal elements, the spectral function for each element $A$ can be divided into two positive functions $A^+$ and $A^-$:

$$ A(\omega) = A^{+}(\omega) - A^{-}(\omega). \qquad (2.80) $$
The default model is accordingly separated into two parts $D^+(\omega)$, $D^-(\omega)$ and the modified entropy reads [215]:

$$S^{\pm} = \int \left[ A^{+}(\omega) - D^{+}(\omega) - A^{+}(\omega) \ln \frac{A^{+}(\omega)}{D^{+}(\omega)} + A^{-}(\omega) - D^{-}(\omega) - A^{-}(\omega) \ln \frac{A^{-}(\omega)}{D^{-}(\omega)} \right] d\omega, \quad (2.81)$$

which is equivalent to Eq. 2.74. The MEM consists in maximizing the function

$$-\frac{\chi^2(A = A^+ - A^-)}{2} + \alpha S(A^+, A^-), \quad (2.82)$$

in which we note that the likelihood function depends only on the difference between $A^+(\omega)$ and $A^-(\omega)$. The maximization must therefore be done under the constraint $A(\omega) = A^+(\omega) - A^-(\omega)$, which is reached when [215]:

$$A^{+}(\omega) = \frac{\sqrt{A(\omega)^2 + 4D^{+}(\omega)D^{-}(\omega)} + A(\omega)}{2} \quad (2.83) \\ A^{-}(\omega) = \frac{\sqrt{A(\omega)^2 + 4D^{+}(\omega)D^{-}(\omega)} - A(\omega)}{2}.$$

The entropy $S^{\pm}$ can then be transformed to [215]:

$$S^{\pm} = \int \left[ \sqrt{A(\omega)^2 + 4D^+(\omega)D^-(\omega)} - D^+(\omega) - D^-(\omega) \right. \\ \left. - A(\omega) \log \frac{\sqrt{A(\omega)^2 + 4D^+(\omega)D^-(\omega)} + A(\omega)}{2D^+(\omega)} \right] d\omega. \quad (2.84)$$

In practice, the diagonal elements are first continued so to obtain the set of {$A_{ii}$}, which are then used to construct the default models for the off-diagonal elements

$$\mathbf{D}_{ij}(\omega) = \sqrt{\mathbf{A}_{ii}(\omega)\mathbf{A}_{jj}(\omega)} + \delta, \quad (2.85)$$

where $\delta$ is a small non-zero constant to prevent the default model from being zero. This choice of default model, instead of a flat one, has shown to improve the results [215].

An example is shown in Fig. 2.9(a)-(b), where the analytic continuation of the self-energy of an isolated dimer cluster is performed with the Poor Man's method [215, 216], and compared to the exact result [222]. The two sites are connected with the hopping term $t = 0.4$eV, and the on-site Coulomb term is $U = 3.5$eV. The self-energy was transformed to an auxiliary Green's function following Eq. 2.79, and then continued element-wise. Note that the results depend on the choice of the constant $C$, which can sometimes drastically determine the final result. For this specific example, the two poles of the exact self-energy are well reproduced, although a spurious shoulder can be noticed in the occupied region in Fig. 2.9. Most interestingly, the off-diagonal part of the self-energy is well captured by the Poor Man's method, so that the resulting local Green's function is in fair agreement with the exact result.
<figure><img src="image_16.png" /><figcaption>Figure 2.9: (a)-(b) Analytic continuation of the diagonal and the off-diagonal elements of the isolated dimer’s self-energy, compared with the exact result [222]. The analytic continuation is performed with the MQEM and with the Poor Man’s methods. For the latter the self-energy was transformed to an auxiliary Green’s function following Eq. 2.79. (c) Corresponding local Green’s function obtained from the self-energy on the real axis.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 2.10: (a)-(b) Analytic continuation of the diagonal and the off-diagonal elements of a dimer cluster calculated within CDMFT and parameterized for t-CuO. The analytic continuation is performed with the MQEM and with the Poor Man’s methods. For the latter the self-energy was transformed to an auxiliary Green’s function following Eq. 2.79. (c) Corresponding spectral function at $\mathbf{k} = (\pi, 0)$ obtained by periodizing the Green’s function directly in real frequency using the continued self-energy. Note the non-causality of the Poor’s Man method in the positive frequency region.</figcaption></figure>

In a more complicated situation however, treating the off-diagonal elements of the self-energy separately can lead to non-causality issues when applying the Dyson equation to obtain the Green's function, due to the systematic error of the MEM method. This is illustrated in Fig. 2.10, where is shown the analytic continuation for a dimer self-energy obtained with CDMFT for a single-band Hubbard model parameterized for t-CuO. Although there is no sign of non-causality in the local self-energy obtained with the Poor's Man method in Fig. 2.10(a), the spectral function at $\mathbf{k} = (\pi, 0)$ shows a negative weight in the occupied part, as shown in Fig. 2.10. In contrast, the MQEM method that we describe hereafter ensures the causality by a simultaneous treatment of the self-energy (Green's function) matrix.

### 2.5.3 DEALING WITH MATRIX-VALUED FUNCTIONS: MAXIMUM QUANTUM ENTROPY METHOD

The principle of maximum quantum entropy method (MQEM) [217] is still based on a probability maximization procedure, but a new definition for the entropy is introduced which is defined for matrix-valued functions. MQEM is therefore free of basis-dependence, as well as non-causality
issues. The new entropy is defined by generalizing Eq. 2.74 to matrix-valued spectral function and default models [217]:

$$S = \int \mathrm{Tr} \mathbf{A}(\omega)[\ln \mathbf{A}(\omega) - \ln \mathbf{D}(\omega)]. \qquad (2.86)$$

The likelihood function can also be extended to matrix-valued functions, and the rest of the procedure is similar to the MEM method.

The default model is matrix-valued, and is constructed from the coefficients of the high-frequency tail expansion of each element of the input Green's function:

$$\mathbf{G}(i\omega_n) = \frac{\mathbf{M}_0}{i\omega_n} + \frac{\mathbf{M}_1}{(i\omega_n)^2} + \dots \qquad (2.87)$$

where $\mathbf{M}_i = \int \omega^i \mathbf{A}(\omega) d\omega$ is the ith moment of the spectral function. In the case of a self-energy, one can easily transform it into a Green's function respecting such high-frequency expansion by using Eq. 2.78. Sim et al. proposed to find the default model which maximizes

$$S_D = \int \mathbf{D}(\omega) \ln \mathbf{D}(\omega) d\omega, \qquad (2.88)$$

with the constraint

$$\mathbf{M}_i = \int \omega^i \mathbf{D}(\omega) d\omega. \qquad (2.89)$$

Restricting the constraint to the first three moments leads to the stationarity condition:

$$\mathbf{D}(\omega) = e^{\sum_{i=0}^{2} \lambda_i \omega^i} \qquad (2.90)$$

where $\lambda_i$ is a Lagrange multiplier introduced to enforce the constraint.

In Fig. 2.9 is shown the result for the isolated dimer self-energy which is good agreement with the Poor Man's method and the exact result. Both the local and NN components of the self-energy are well reproduced, and the local Green's function exhibits the expected four peak structure at the right frequencies. Note that the two methods seem to underestimate the difference in peak height between the outer and inner peaks.

In Fig. 2.10 is shown the analytic continuation using the MQEM method for the t-CuO dimer cluster. The results are in good agreement with the Poor's Man method, although some features seem washed out. Most importantly, the spectral function's causality is ensured in the whole frequency range by the simultaneous treatment of the matrix-valued functions, as can be seen in Fig. 2.10(c).

## 2.6 AN EFFECTIVE SPIN MODEL FOR LIGHTLY DOPED CUPRATES

DMFT and CDMFT are powerful techniques to treat exactly a part of the electron-electron interactions, either fully at the local level, or at some finite length scale. If they have been successful for studying various properties of the Hubbard model, like the Mott-insulator transition [43] or the
d-wave superconductivity [50, 52], it remains extremely difficult to extract the underlying mechanisms. For instance, while d-wave superconductivity has been shown to exist in the 2D Hubbard model [50, 52, 145–148, 223], its origin remains elusive. Different strategies exist, which aim at interpreting the results by identifying the key ingredient(s) responsible for the studied feature. An example is the recently developed fluctuation-diagnostic, which allows to analyze the self-energy by a decomposition into different channels by using two-particle quantities [112, 113]. Although powerful and versatile, this method requires the computationally challenging calculation of two-particle quantities. Another strategy is to construct simpler effective models designed to capture only the relevant physics, and to be easier to solve as well as to interpret. We present here such an effective spin model, which is expected to capture at least partially the physics of undoped and lightly-doped cuprates.

Typical values of on-site interaction for cuprates are relatively large, of the order $U \sim 8–10t$. Although $U$ is of the same order than the bandwidth $W \simeq 8t$ (for a 2D square lattice), the ground state of the half-filled single-band Hubbard model can be approximated by a state consisting of singly-occupied sites. It is well described by an effective t-J model in which states with double-occupations are projected out [224–226]. The kinetic term is restricted to the hopping processes that do not create a double occupancy, and the Coulomb on-site interaction term is transformed into an Heisenberg exchange interaction between neighboring electrons mediated by the coupling $J \simeq 4t^2/U$. Considering for now the case with only NN hopping, the Hamiltonian reads

$$H_{t-J} = -t \sum_{\langle ij \rangle, \sigma} \tilde{c}_{i\sigma}^{\dagger} \tilde{c}_{j\sigma} + J \sum_{\langle ij \rangle} (\mathbf{S}_i \cdot \mathbf{S}_j - \frac{1}{4} n_i n_j), \quad (2.91)$$

where $\tilde{c}_{i\sigma}^{\dagger}$, $\tilde{c}_{i\sigma}$ are the restricted creation/annihilation operators to prevent from double occupancy $\tilde{c}_{i\sigma}^{\dagger} = c_{i\sigma}^{\dagger}(1 - n_{i-\sigma})$, and $\mathbf{S}_i$ is the spin operator on site i [224–226]. For simplicity, the spin exchange is assumed isotropic, but this does not restrict the generality of the method.

At half-filling, we can notice that the kinetic term vanishes since it is impossible for an electron to hop without generating a double occupancy. Therefore the ground state corresponds to the one of the Heisenberg model. The determination of the spectral function is however nontrivial: initially a hole is created, and will propagate through the lattice since it can hop to neighboring sites without generating double occupancies. The hole is dressed by the antiferromagnetic spin fluctuations, and forms a spin polaron [95, 227–229], similarly to an electron/hole coupled to phonons [230, 231]. Interestingly, if there is a single spin polaron, the spectral function can be approximated by an analytical formula which allows to compute it directly on the real axis and at the thermodynamic limit [95]. To show that, the kinetic part of the Hamiltonian is rewritten as a coupling term between spinless holes and spin operators, and the spin fluctuations are treated as spin waves. This provides a Hamiltonian describing a hole coupled to magnons, whose self-energy can be approximated using the self-consistent Born approximation (SCBA) [232, 233]. Since we consider a single hole (we neglect the hole-hole coupling), the method is expected to hold for low doping values, and should quickly break down as more holes are introduced. A clear advantage of the SCBA is that it is computationally cheap, and allows to easily treat systems at the
thermodynamic limit, thus free from finite size effects. Note that SCBA is not the only possibility to calculate the hole propagation in an antiferromagnetic background. For instance, a successful recent development is the parton theory [227–229] in which the spin polaron is regarded as the bound state of a spinon and a holon. We however restrict our work to the simple framework of SCBA which will show to be sufficient to interpret the theoretical and experimental spectral functions, and which can be thought of as an approximation to the more involved treatments as the parton theory.

The first step is to re-write the spin operators as bosonic operators using the Holstein-Primakoff transformation [234] (for half-integer spins):

$$ 
\begin{aligned}
S_i^+ &= \sqrt{1 - a_i^\dagger a_i} a_i \sim a_i, \\
S_i^- &= a_i^\dagger \sqrt{1 - a_i^\dagger a_i} \sim a_i^\dagger, \\
S_i^z &= \frac{1}{2} - a_i^\dagger a_i
\end{aligned} 
\qquad (2.92) 
$$

where $a_i^\dagger, a_i$ are Bose operators, and the linear approximation is used. Then, the fermionic operators $c_{i\uparrow}, c_{i\downarrow}$ are decomposed to define the spinless hole operator:

$$ 
\begin{aligned}
c_{i\uparrow} &= h_i^\dagger, \\
c_{i\downarrow} &= h_i^\dagger S_i^+.
\end{aligned} 
\qquad (2.93) 
$$

As defined, there is a discrepancy between the local fermion Hilbert space, composed of 3 states $|0\rangle|\uparrow\rangle, |\downarrow\rangle$ and the local spin-hole Hilbert space which is composed of the four product states $|0, \uparrow\rangle, |1, \uparrow\rangle, |0, \downarrow\rangle, |1, \downarrow\rangle$. This can be cured by noting that there can not be on the same site a hole and a down spin state from the definition of $c_{i\downarrow} = h_i^\dagger S_i^+$, and may be directly taken into account in the SCBA. The constraint on the double occupation is automatically fulfilled:

$$ \tilde{c}_{i\uparrow}^{\dagger} |\downarrow\rangle = (1 - n_{i\downarrow}) c_{i\uparrow}^{\dagger} |\downarrow\rangle = 0 \qquad (2.94) $$

since $h_i |0, \downarrow\rangle = 0$.

From this definition, Eq. 2.91 may be written in terms of hole and bose operators:

$$ H_{t-J} = -t \sum_{\langle ij \rangle} (h_i h_j^\dagger a_j + h.c.) + \frac{J}{4} \sum_{\langle ij \rangle} h_i h_i^\dagger (a_i a_j + a_i^\dagger a_j^\dagger + a_i^\dagger a_i + a_j^\dagger a_j) h_j h_j^\dagger + \text{const.} \qquad (2.95) $$

The first term describes the hopping of the hole which simultaneously couples with the spin degree of freedom. In the second term, the factors $h_i h_i^\dagger$ and $h_j h_j^\dagger$ ensure the loss of magnetic energy by the presence of holes inside the lattice. As it is, the second term is non-trivial and prevents from solving the model directly. It can approximated in a mean-field fashion: the hole operators are replaced by their mean-field average, thus transforming to a factor $(1 - \delta)^2$ where $\delta$ is the concentration of holes. This amounts to take into account the hole doping as a simple damping of the magnon dispersion by a factor $(1 - \delta)^2$. An additional difficulty arising from the
second term, even within the mean-field approximation, is that it does not conserve the number of bosonic paricles. This may be solved by applying a Bogoliubov transformation, after having performed a Fourier transform [95]. The Hamiltonian transforms into:

$$H_{t-J} = \frac{zt}{\sqrt{N}} \sum_{\mathbf{k},\mathbf{q}} \left[ h_{\mathbf{k}} h_{\mathbf{k}-\mathbf{q}}^{\dagger} \alpha_{\mathbf{q}} (u_{\mathbf{q}} \gamma_{\mathbf{k}-\mathbf{q}} + v_{\mathbf{q}} \gamma_{\mathbf{k}}) + \text{h.c.} \right] + \sum_{\mathbf{q}} \omega_{\mathbf{q}} \alpha_{\mathbf{q}} \alpha_{\mathbf{q}}^{\dagger}, \quad (2.96)$$

where $\gamma_{\mathbf{k}} = \frac{1}{2}(\cos(k_x) + \cos(k_y))$ is the usual dispersion of a square lattice, z is the coordination number, $\omega_{\mathbf{q}} = SzJ(1 - \delta)^2\nu_{\mathbf{q}}$ is the magnon dispersion with $\nu_{\mathbf{q}} = \sqrt{1 - (\gamma_{\mathbf{q}})^2}$. $\alpha_{\mathbf{q}}$, $\alpha_{\mathbf{q}}^{\dagger}$ are defined via the Bogoliubov transformation as

$$\begin{align}
\alpha_{\mathbf{q}} &= u_{\mathbf{q}} a_{\mathbf{q}} - v_{\mathbf{q}} a_{-\mathbf{q}}^{\dagger} \nonumber \\
u_{\mathbf{q}} &= \left( \frac{1 + \nu_{\mathbf{q}}}{2\nu_{\mathbf{q}}} \right)^{\frac{1}{2}} \tag{2.97} \\
v_{\mathbf{q}} &= -\text{sign}(\gamma_{\mathbf{q}}) \left( \frac{1 - \nu_{\mathbf{q}}}{2\nu_{\mathbf{q}}} \right)^{\frac{1}{2}}. \nonumber
\end{align}$$

Eq. 2.96 now has the simple form of a propagating hole which couples to the spin degree of free-dom with the coupling $M_{\mathbf{k},\mathbf{q}} = (u_{\mathbf{q}}\gamma_{\mathbf{k}-\mathbf{q}} + v_{\mathbf{q}}\gamma_{\mathbf{k}})$. An important difference with the usual electron-phonon coupling has to be noted at this stage: in the case of a spin polaron, there is no bare kinetic term for the hole, its motion is entirely controlled by the coupling to the spins with a bandwidth on the scale of J. Therefore the usual picture of mass enhancement due to the coupling of "free" particles to phonons, i.e. in a weak-coupling regime, does not apply here. We face a strong-coupling problem since the hopping t is larger than the exchange coupling J, such that the dominating term in Eq. 2.96 is the coupling between the hole and the spin excitations.

The absolute value of the coupling $M_{\mathbf{k},\mathbf{q}}$ is shown in Fig. 2.11: at all k points it vanishes at $\mathbf{q} = (0,0)$ and $\mathbf{q} = (\pi, \pi)$, so it is the short-wavelength spin fluctuations that are important. While at $\mathbf{k} = (0,0)$ and $\mathbf{k} = (\pi, \pi)$ the coupling is relatively small at all $\mathbf{q}$, it shows a strong momentum dependence for $\mathbf{k} = (\pi, 0)$ and $\mathbf{k} = (\frac{\pi}{2}, \frac{\pi}{2})$, i.e. in the region showing the largest spectral weight in the hole's spectral function which is discussed later in this section. Note that here $\mathbf{q}$ and $\mathbf{k}$ are both defined in the lattice BZ.

An analytic solution to this Hamiltonian can be obtained using the SCBA [95, 232, 233], from which the following self-energy can be derived:

$$\Sigma(\mathbf{k}, \omega) = \frac{z^2 t^2}{N} \sum_{\mathbf{q}} \frac{|M_{\mathbf{k},\mathbf{q}}|^2}{\omega - \omega_{\mathbf{q}} - \Sigma(\mathbf{k}-\mathbf{q}, \omega - \omega_{\mathbf{q}}) + i\eta}, \quad (2.98)$$

where $\eta$ is a small broadening number. The hole's Green's function is then given by:

$$G^h(\mathbf{k}, \omega) = \frac{1}{\omega - \Sigma(\mathbf{k}, \omega) + i\eta}. \quad (2.99)$$
<figure><img src="image_18.png" /><figcaption>Figure 2.11: Absolute value of the coupling function $|M_{k,q}|$ calculated as a function of $q$ at different $k$ points represented on the BZ sketch.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 2.12: (a) Illustration of the different exchange couplings on a 2D square lattice. (b) Sketch of the BZ and the k-path along which is calculated the hole's spectral function shown in (c). The green wedge depicts the irreducible part of the BZ.</figcaption></figure>

In practice the latter is obtained by first converging Eq. 2.98 in a self-consistent loop. Moreover, since in our Hubbard one-band effective models for NaCCOC and t-CuO we include the longer-ranged hopping amplitudes $t'$ and $t''$, we use a modified magnon dispersion which incorporates the corrections from $t'$ and $t''$ as well as higher order processes in $t$, following Ref. [235]. This leads to correct the NN exchange $J$ as:

$$J = \frac{4t^2}{U} - \frac{24t^4}{U^3} \qquad (2.100)$$

and to include the NNN $J' = \frac{4t'^2}{U} + \frac{4t^4}{U^3}$, NNNN $J'' = \frac{4t''^2}{U} + \frac{4t^4}{U^3}$, as well as the cyclic $J_c = \frac{80t^4}{U^3}$ exchange couplings. The latter is essential to capture the difference in energy in the magnon dispersion between $\mathbf{k} = (\pi, 0)$ and $\mathbf{k} = (\pi/2, \pi/2)$ which is seen in experiment [166]. These couplings are illustrated in Fig. 2.12(a) and correspond to the following Hamiltonian:

$$H_{J-J',-J'',-J_c} = J \sum_{\langle ij \rangle} \mathbf{S}_i \cdot \mathbf{S}_j + J' \sum_{\langle\langle ij \rangle\rangle} \mathbf{S}_i \cdot \mathbf{S}_j + J'' \sum_{\langle\langle\langle ij \rangle\rangle\rangle} \mathbf{S}_i \cdot \mathbf{S}_j \\ + J_c \sum_{\langleijkl\rangle} [(\mathbf{S}_i \cdot \mathbf{S}_j)(\mathbf{S}_k \cdot \mathbf{S}_l) + (\mathbf{S}_i \cdot \mathbf{S}_l)(\mathbf{S}_k \cdot \mathbf{S}_j) - (\mathbf{S}_i \cdot \mathbf{S}_k)(\mathbf{S}_j \cdot \mathbf{S}_l)], \qquad (2.101)$$

where $\langleijkl\rangle$ refer to a square plaquette of four neighboring sites $i, j, k, l$. With these corrected exchange parameters, the magnon dispersion is given by [236]:

$$\begin{aligned} \omega_q &= Z_c(\mathbf{q})(1 - \delta)^2 \sqrt{A_q^2 - B_q^2} \\ A_q &= 4JS + 4J'S(\cos(q_x) \cos(q_y) - 1) + 2J''(\cos(2q_x) + \cos(2q_y) - 4J_cS^3(\cos(q_x) \cos(q_y) + 1)) \\ B_q &= 2JS(\cos(q_x) + \cos(q_y)) - 4J_cS^3(\cos(q_x) + \cos(q_y)) \end{aligned} \qquad (2.102)$$

where the factor $Z_c(\mathbf{q})$ incorporates the effects of the quantum fluctuations [235, 236]. Moreover, since the $t'$ and $t''$ hopping terms do not disturb the antiferromagnetic order, they simply result in
a dispersion of the hole without coupling to the spin degree of freedom, leading to a modification of Eq. 2.98 as follows:

$$ \Sigma(\mathbf{k}, \omega) = \frac{z^2 t^2}{N} \sum_{\mathbf{q}} \frac{|M_{\mathbf{k},\mathbf{q}}|^2}{\omega - \omega_{\mathbf{q}} - \epsilon_{\mathbf{k}-\mathbf{q}} - \Sigma(\mathbf{k}-\mathbf{q}, \omega - \omega_{\mathbf{q}}) + i\eta}, \quad (2.103) $$

where $\epsilon_k = 4t' \cos(k_x) \cos(k_y) + 2t''(\cos(2k_x) + \cos(2k_y))$ (note the positive sign since it is a hole dispersion).

In the practical implementation, the parameter $\eta$ is first chosen relatively large and reduced as the self-consistent accuracy is reached. Moreover, one can make use of the high symmetry of the magnon dispersion of cuprates: the BZ can be reduced to the irreducible part as shown in green in Fig. 2.12(b). A typical SCBA spectrum obtained with parameters relevant for NaCCOC ($x=0$) is shown in Fig. 2.12(c): a coherent quasi-particle peak can be identified, with replica of the dispersion at lower energy which are increasingly incoherent. The system size is 30 × 30, way larger than what can be treated within CDMFT or other techniques treating explicitly the Hubbard interaction U, and which allows to eliminate finite size effects.
# PART IV

## STUDIES OF SELECTED COPPER OXIDES
<figure><img src="image_20.png" /></figure>
# 3 STUDY OF THE ANTIFERROMAGNETIC FLUCTUATIONS IN Ca$_{2-x}$Na$_x$CuO$_2$Cl$_2$

## 3.1 INTRODUCTION

Although being a cuprate, strictly speaking Ca$_{2-x}$Na$_x$CuO$_2$Cl$_2$(NaCCOC) is not a copper oxide, but an oxychloride. It has a simple tetragonal crystal structure [80–83] illustrated in Fig. 3.1(a)). The main difference is the replacement of the apical O by Cl atoms, which confers a strong 2D character to the CuO$_2$ layers, making it an interesting candidate to connect to the ZRS picture [38]. In fact NaCCOC attracted the scientific attention since it was one of the first cuprates without oxygen octahedra to show superconductivity upon hole-doping [80, 81]. The introduction of hole carriers in the parent compound [82] can be performed by replacing the Ca by Na atoms [80, 81, 83], or by vacancies [237]<fnref n="1" />.

NaCCOC is antiferromagnetically ordered [84] (up to $x = 0.02$ doping) and shows a Mott gap $0.5 \le \Delta \le 2.5$eV [85–89]. Photoemission measurements identified a *kink* feature in the energy dispersion [86] at low energy, similar to the well-known *waterfall* feature of cuprates [24]. Upon hole-doping, NaCCOC becomes superconducting at $x \sim 0.10$ until $x = 0.20$ [80, 81, 84] with a superconducting gap of *d*-wave symmetry [240], and interestingly there is a region of the phase diagram around $x \sim 0.05$ where both the superconductivity and the long-range antiferromagnetic order are suppressed. In this low-doping region, a spin-glass like behavior has been identified [84], as well as a CDW [25–27] of period of 4 lattice constants, but without long-range order [241, 242]. The famous PG phase, present in all cuprate compounds [12, 14, 24], characterized in the spectral function by the opening of an anistropic gap in the BZ, is also observed at low to intermediate dopings [26, 90–94]. The maximum hole concentration experimentally available is $x = 0.20$, which is still in the superconducting phase, and hence prevents to infer on the existence of the Fermi liquid and the strange metal phases [14] in the overdoped regime. Though, the NaCCOC phase diagram until the optimal doping is in excellent agreement with the one of cuprates, with the peculiarity that the PG exists down to $T = 0$ around $x = 0.05$ since both the long-range antiferromagnetic order and the superconductivity are suppressed. Therefore, along with the fact that NaCCOC has a quite simple structure, it appears as a promising candidate to confront the cuprates’ low-energy models to experiments.

In this chapter, we are particularly interested in the physical origin of the waterfall feature and the PG, which are both seen in ARPES measurements. The waterfall manifests itself as a kink

<footnote n="1">¹Note that two other chemical substitutions are possible: replacing the Ca by K [238], and Cl by Br [239], but very few studies are dedicated to these options.</footnote>
in the energy dispersion at low-energy [24, 86], and is induced by the cross-over between low- and high-energy sectors yielding different physics. It is observed in the parent NaCCOC compound ($x = 0$) and survives at low and intermediate dopings [90–92]. A similar kink is observed in the one-band Hubbard model which is often used as an effective model for cuprates, and the low-energy part has been theoretically identified as stemming from a spin-polaron [95–98], i.e. a hole propagating in an antiferromagnetic background. A consensus is however not reached for the high-energy part, which has been attributed to spin-charge separation [99, 100], to weak-coupling spin-density wave, or to long-range hopping processes [101, 102] within a single magnetic sublattice [98].

The PG is a universal feature of cuprates, and was first experimentally observed in the Knight shifts [21], spin-lattice relaxation time [22], transport measurements [23], photoemission [243–245], and others [12]. Its most well-known feature in hole-doped cuprates is the opening of a gap at the anti-nodes ($\pi, 0$) while at the node ($\pi/2, \pi/2$) a quasiparticle peak survives, leading to a discontinuous Fermi surface. Since the experimental discovery, it has been shown with embedded quantum cluster methods like CDMFT [103–106], and DCA [65, 107–110] that the 2D one-band Hubbard model display a PG upon hole-doping. It is linked to the Mott insulating transition since the PG occurs at finite doping for $U$ values above the critical $U_c$ for the insulating transition at half-filling. Though, it is a distinct phase: the value of the PG does not extrapolate to the Mott gap as the doping $x \to 0$ [106]. Early on, from QMC calculations it was proposed that the PG originates from short range antiferromagnetic fluctuations [111], which was corroborated later with the finding that quantum cluster methods with small cluster sizes (2-site and $2 \times 2$ plaquette) were also able to capture the right behavior [104, 107, 108]. Gull et al. showed that it was not an artifact from the finite size of the embedded cluster, and that considering larger systems does not significantly change the results [110]. More recently, a complementary analysis called fluctuation diagnostic, based on the decomposition of the contributions to the self-energy into different channels (spin, charge, and particle-particle), pointed out the dominance of the spin channel at the wavevector $Q = (\pi, \pi)$ in the PG phase from the analysis of the self-energy obtained with DCA [112]. The same technique was also applied to self-energies obtained by diagrammatic Monte Carlo, whose recent developments allowed to perform calculation in a specific parameter regime relevant for the PG physics [113, 114]. In agreement to previous works, the anti-nodal self-energy is dominated by commensurate wavevector $Q = (\pi, \pi)$ spin fluctuations while at the node the spin channel is also leading, but at an incommensurate wavevector $Q = (\pi + \delta, \pi + \delta)$. Similar results were also obtained with the dual boson technique which includes long-range correlation effects by expanding around the DMFT solution [75]. Moreover, the correlation length was measured to be finite from DCA calculations [246]. Therefore, by the combination of many different techniques, a consensus seems to be reached that the PG in the hole-doped one-band Hubbard model, and by extension in cuprates, originates from short-range spin fluctuations [75, 104, 107, 111–114, 246].

Both the waterfall and the PG stem from the antiferromagnetic fluctuations, with the difference that the waterfall occurs already in the parent compound which hosts long-range antiferro-
magnetic order, whereas the PG is expected to emerge from short range fluctuations upon hole-doping. In this chapter, we study the emergence of both phenomena in NaCCOC using CDMFT on a 2 × 2 plaquette and 8-site cluster, confronting the results to the existent experimental literature and ARPES measurements performed in synchrotron facilities for $x \approx 0.05$ and $x \approx 0.10$ doped samples, to which we actively participated in collaboration with Matteo d'Astuto, Chafic Fawaz and Hervé Cercellier (Institut Néel, Grenoble France). We also compare our numerical results to RIXS measurements performed by Matteo d'Astuto, Blair Lebert et al. on undoped and doped NaCCOC [166].

Calculations for the parent compound provide a spectral function in excellent agreement with ARPES measurements [85–88], and show a clear waterfall feature. By performing SCBA calculations for an effective model as described in Sec. 2.6, we show that the low-energy part of the spectrum is well accounted for by a spin-polaron, in agreement with previous observations [95–98]. We find that the high-energy part is well captured at the atomic limit, i.e. by a self-energy calculated with the Hubbard-I [35] approximation, consistent with the picture of a hole propagating into a lattice of localized electrons (spins), which are antiferromagnetically ordered. Our analysis therefore provides a simple and unambiguous explanation for the high-energy part of the waterfall feature.

Upon hole-doping $x > 0$, we identify an emerging PG characterized by the opening of a gap at the anti-node ($\pi, 0$) in the calculations, as well as the measured spectra, both being in good qualitative agreement. Similarly to observations of Shen et al. [92] at low doping, our ARPES measurements on $x \approx 0.05$ and $x \approx 0.1$ show a suppression of the coherent quasiparticle peak at the nodal point ($\pi/2, \pi/2$), in sharp contrast to our calculations in which the peak is always present, challenging the simple picture of a perfect 2D lattice model for the underdoped region. However, we identify a double peak structure at the nodal point which is in excellent agreement, especially for the 8-site cluster, with previous experimental observations [91, 92].

We analyze the PG phase by dissecting the state histogram of the embedded 2 × 2 plaquette using a modified cluster basis in contrast to previous works [106, 121, 247]. We show that due to the presence of two degenerate orbitals in the 2 × 2 plaquette, which correspond to the cluster momenta $\mathbf{K} = (\pi, 0)$ and $\mathbf{K} = (0, \pi)$ and are weakly coupled to the self-consistent bath, the most probable state until $x = 0.12$ is a non-degenerate singlet state $S = 0$, in agreement with previous observations [106, 121, 247]. This singlet nature of the ground state leads to a cancellation of the zero-energy excitations for the two degenerate orbitals, that we explicitely show following Ref. 121 that we adapt to our new choice of basis, such that the lowest excitation channel carries a spin $S = 1$, momentum $\mathbf{Q} = (\pi, \pi)$ and a non-zero energy. This leads to the opening of the PG at the anti-nodes. The 2 × 2 plaquette is the minimum cluster containing this mechanism, and by comparison to the 8-site cluster we show that the picture remains valid for larger cluster.

Moreover, this observation is linked to the observation of a dominant $\mathbf{Q} = (\pi, \pi)$ peak in the DSSF, that we directly compare to RIXS measurements for both the undoped and doped systems.
<figure><img src="image_21.png" /><figcaption>Figure 3.1: (a) Tetragonal crystal structure of the parent compound NaCCOC (x=0). (b) Sketch of the BZ and the k-path along which is plotted the band structure obtained from a LDA calculation shown in (c), along with the orbital projected DOS. (d) The LDA orbital projected band structure overlayed on the ARPES spectrum taken at SOLEIL synchrotron with a 10% doped sample. The LDA bands were all shifted by the same constant energy δε to match the measurement.</figcaption></figure>
Consistently with the experimental findings for NaCCOC [166] and other cuprates [248–250], we show that in the doped sample clear magnon excitations survive although having a shorter life-time than in the parent compound. By inspecting the spin-spin correlation function, we find a similar behavior as Werner et al. [61, 123]: upon hole doping, at long imaginary times $\tau \to \beta/2$ the correlation function becomes larger than for the parent compound and does not decay with distance (on the scale of the 2 × 2 cluster). We show that this behavior, for which no clear explanation was found so far, is directly linked to the existence of low-energy spin excitations and that the "freezing" of the correlation function can extend to long-range distances as the temperature is lowered, what could be an interesting subject of further studies.

Overall, this chapter aims at depicting a unified description of the spectral and magnetic properties of NaCCOC by simultaneously focusing on the waterfall feature and PG phase with a combined theoretical and experimental approach. We emphasize the robustness of the 2×2 plaquette results by a systematic comparison to the 8-site cluster, and find a peculiar behavior of the long-range spin-spin correlation function at long times $\tau \to \beta/2$.

## 3.2 MODEL AND METHOD

Similarly to cuprates, NaCCOC is composed of CuO₂ planes, of 2D character [251], intercalated between charge reservoir layers made of Ca and Cl atoms, as shown in Fig. 3.1(a). The low energy window, relevant for the waterfall feature and the PG, is ruled by the physics of these 2D CuO₂ layers. In the following, we first describe the effective Hubbard model that we use to represent NaCCOC. We then present the clusters used for solving the model with CDMFT, and we give details about the ARPES measurements carried out at the Photon Factory (Tsukuba, Japan) and at SOLEIL (Saint-Aubin, France) synchrotrons.

### 3.2.1 BAND STRUCTURE OF Ca₂CuO₂Cl₂

Although DFT should not be used to study the spectral properties of this correlated material, it provides an interesting starting point for building an effective model of NaCCOC. Indeed, from electronic structure calculations with the LDA potential one can extract a TB Hamiltonian as explained in Sec. 2.1. In Fig. 3.1(c) is shown the band structure along the **k**-path shown in Fig. 3.1(b) as obtained from a LDA calculation with Vienna ab initio simulation package (VASP) [252–255] for the parent compound ($x = 0$). As expected from a LDA calculation (see Sec. 1.4), a single band crosses the Fermi energy, consistent with the valence of the Cu atoms which are in $d^9$ electronic configuration, i.e. one of the *d* orbitals is half-filled. Note that actually two of these bands can be seen in the band structure, which is simply due to the fact that the unit cell contains two equivalent Cu atoms. A striking feature is that the single band is well isolated from the rest of the Cu *d*/ O *p*/ Cl *p* manifold which lies below (−1.5eV) and is in rather good agreement with the ARPES measurement done on a *x* = 0.1 hole-doped sample (see Fig. 3.1(d)). The isolated band, also seen in the experimental spectrum, can be traced back to be of Cu *d*<sub>*x*<sup>2</sup>−*y*<sup>2</sup></sub> and O *p*<sub>*x*</sub>/*y* character as shown by the projected DOS in Fig. 3.1(c), consistent with a ZRS state [38].
<figure><img src="image_22.png" /><figcaption>Figure 3.2: (a) The CuO$_2$ planes are modeled with an effective Hubbard model. The hopping parameters were obtained from the fit of the single band dispersion shown in Fig. 3.1(c), and the on-site U interaction is extracted from a fit of the experimental magnon dispersion (left panel of (b)) obtained from a RIXS measurement carried on an undoped sample by Lebert et al. [166]. The Heisenberg spin model used for the fit is shown in the right panel of (b). (c) The two different clusters considered for the CDMFT calculations: the 2 × 2 plaquette and the 8-site cluster.</figcaption></figure>

The effective TB model can be restricted to a square lattice in two dimensions, for the single band shows a weak dispersion along the Γ-Γ$_Z$ path, in agreement with transport measurements [251]. By performing a Wannier fit using the wannier90 [183, 184] package, we obtain the hopping pa- rameters t = 0.425eV, t' = −0.076eV ~ −0.2t and t'' = 0.05eV ~ 0.1t, which are close to the typical cuprates’ hopping term t'/t ≃ −0.2/0.3. By restricting the 2D TB model to these three hopping terms only, we obtain the Wannier band shown in Fig. 3.1(c) which represents well the single DFT band. Note that a key feature of the dispersion is included in the 2D effective TB model: the dispersion becomes flat at the X point, i.e. at the antinode (π, 0), yielding a van Hove singularity which is of importance in the PG physics [246, 256].

3.2.2 LOW-ENERGY EFFECTIVE HUBBARD MODEL

Based on these observations, we model NaCCOC by a single-band Hubbard model [35–37] which aims at capturing the low-energy physics of the CuO$_2$ planes. The model is illustrated in Fig. 3.2(a), the Hamiltonian reads:

$$
H = U \sum_i n_{i\uparrow} n_{i\downarrow} - \mu \sum_i n_i - t \sum_{\langle ij \rangle, \sigma} c_{i\sigma}^{\dagger} c_{j\sigma} - t' \sum_{\langle\langle ij \rangle\rangle, \sigma} c_{i\sigma}^{\dagger} c_{j\sigma} - t'' \sum_{\langle\langle\langle ij \rangle\rangle, \sigma} c_{i\sigma}^{\dagger} c_{j\sigma}, \quad (3.1)
$$

where *U* is the Hubbard on-site interaction term, *μ* the chemical potential, 〈·〉, 〈·,·〉, 〈〈·,·〉〉 denote respectively the NN, NNN and NNNN associated with the hopping terms *t* = 0.425eV, *t'* = −0.076eV ~ −0.2*t* and *t''* = 0.05eV ~ 0.1*t*. As described above, the latter are obtained by a fit of the LDA single band using MLWF [181, 182], and we are left with determining the *U*
value. Instead of calculating it *ab initio*, we rather rely on the magnon dispersion obtained by a RIXS measurement on the parent compound which was performed by our collaborators Lebert, d'Astuto et al. [166]. Since undoped NaCCOC is an antiferromagnetic insulator with long-range magnetic order [84], the magnon dispersion can approximately be represented by a spin-wave calculation based on the Heisenberg model shown in the right panel of Fig. 3.2(b). The spin model Hamiltonian is given by

$$ H_{\text{spin}} = J \sum_{\langle ij \rangle} \mathbf{S}_i \mathbf{S}_j + J' \sum_{\langle\langle ij \rangle\rangle} \mathbf{S}_i \mathbf{S}_j + J'' \sum_{\langle\langle\langle ij \rangle\rangle\rangle} \mathbf{S}_i \mathbf{S}_j \\ + J_c \sum_{\langleijkl\rangle} \left[ (\mathbf{S}_i \mathbf{S}_j)(\mathbf{S}_k \mathbf{S}_l) + (\mathbf{S}_i \mathbf{S}_l)(\mathbf{S}_k \mathbf{S}_j) - (\mathbf{S}_i \mathbf{S}_k)(\mathbf{S}_j \mathbf{S}_l) \right] \quad (3.2) $$

where $\langleijkl\rangle$ refer to a square plaquette of four neighboring sites i, j, k, l (see Fig. 3.2(b)), and the exchange couplings are: $J = \frac{4t^2}{U} - \frac{24t^4}{U^3}$, $J' = \frac{4t'^2}{U} + \frac{4t^4}{U^3}$, $J'' = \frac{4t''^2}{U} + \frac{4t^4}{U^3}$ and $J_c = \frac{80t^4}{U^3}$ [235, 236]. Keeping the hopping parameters fixed to the value obtained *ab initio*, the on-site interaction $U$ is treated as a fitting parameter. We stress the importance of the corrections of the exchange parameters from the $t^4$ processes, and especially the cyclic exchange $J_c$ which allows to capture the difference in the magnon energy between $(\pi, 0)$ and $(\pi/2, \pi/2)$ [166]. This procedure results in the fitted magnon dispersion shown in Fig. 3.2(b), which is in excellent agreement with the experimental data, especially by considering only $U$ as fitting parameter. We obtain $U \approx 4.335\text{eV} = 10.2t$ which will be our default choice for the theoretical calculations (the only exception is for the state histogram discussed later in this chapter which were obtained with $U \approx 4.25\text{eV} = 10t$).

The on-site Hubbard interaction $U$ is of the same order of magnitude than the bandwidth $W \approx 3.5\text{eV}$, preventing from the use of perturbation theory in either the kinetic or potential energy. To solve this many-body strongly correlated system, we perform CDMFT calculations with two different clusters presented in Fig. 3.2(c). Most of the calculations and results presented in this chapter are obtained with the 2×2 plaquette which contains the minimum ingredients to capture the PG in the hole-doped single band Hubbard model [121, 247]. The effective impurity model is solved using CTQMC [60–62]-based solvers working at finite temperature: the hybridization-CT-HYB [116] and the interaction-expansion CT-INT [117] solvers based on the ALPSCore libraries [118]. Additionally, we also used the hybridization-expansion solver of the TRIQS toolkit [119, 120] in order to obtain the state histograms. Thanks to the limited size of the plaquette cluster, it is possible to perform calculations at low temperature, deep in the PG phase. In order to check that the results we obtain are not an artifact from the specific properties of the plaquette cluster, we also perform calculations for the 8-site cluster shown in Fig. 3.2(c).

We ensured that all the theoretical CDMFT results presented in this chapter were well converged, and we performed extra self-consistent loops with a high Monte-Carlo sampling: the number of counts was always $\ge 2 \cdot 10^6$ for all temperatures and cluster sizes. To obtain the momentum-resolved spectral function, we first perform an analytic continuation of the cluster self-energy using the MQEM [217] presented in Sec. 2.5. From the self-energy, we compute the periodized Green's function using the Green's function periodization procedure (Eq. 2.48) [67]. For the
<figure><img src="image_23.png" /><figcaption>Figure 3.3: In-situ photography taken at SOLEIL synchrotron of the sample in the experiment chamber under UHV (left). A ceramic tip (highlighted by the red dotted frame) is glued on the sample surface to allow for cleaving under vacuum, resulting in a clean 2D surface as shown in the right panel. The sides of the samples are covered with a graphite spray (darker zone around the sample) to ensure the electrical contact with the copper sample holder. On the upper left side of the latter, one can notice the gold sample used for the Fermi energy reference.</figcaption></figure>

DSSF, we use a different procedure: from the two-time two-particle imaginary time Green's function, the spin-structure factor is first calculated at the different Q points of interest, and then an analytic continuation with a bosonic kernel is performed using the MaxEnt package [220]. More details on the link between the DSSF and the two-time two-particle Green's function are given in the dedicated Sec. 3.5.

### 3.2.3 ANGLE-RESOLVED PHOTOEMISSION SPECTROSCOPY MEASUREMENTS

In addition to the theoretical calculations, we performed two ARPES measurements at synchrotron facilities: the Photon Factory (Tsukuba, Japan) in close collaboration with Matteo d'Astuto and Yukiko Obata, and SOLEIL (Saint-Aubin, France) along with Matteo d'Astuto, Chafic Fawaz and Hervé Cercellier. The samples measured were synthesized by Matteo d'Astuto, David Santos-Cottin, Hajime Yamamoto and Ikuya Yamada. The characterization and preparation of the samples for the experiments were done by Matteo d'Astuto, Yukiko Obata, and Chafic Fawaz. We do not detail the synthesis here, the interested reader can obtain more details in Ref. 83. The samples' characteristic size is about 2 × 2 × 0.2mm (width×height×thickness). The doping is estimated by measuring the superconducting transition temperature $T_c$, a standard procedure for this material [90] but not that precise, and complemented by a measure of the lattice parameters. For the $x = 0.1$ samples, the $T_c$ measurement leads to an estimation of the doping close to the expected value, but the error on the Na concentration is of the order of 1–2%, i.e. $x = 0.1 \pm 0.02$. This is even worse for the $x = 0.05$ sample since the sample is not superconducting and one can only rely on the lattice parameters measurements, from which it is estimated that the real doping is somewhere around $x \sim 0.05$.

The crystalline planes' orientation was characterized using x-ray diffraction, and the samples were then glued accordingly on the copper sample holder with a silver epoxy to ensure the electrical contact (see Fig. 3.3), always manipulating the material in a dry environment to prevent