from surface pollution and oxidation. To further make sure that there is a good contact between the sample and the holder, to avoid surface charging during the measurement, graphite has been sprayed over. Finally, a ceramic tip is glued over the sample (see Fig. 3.3) and is only removed under **ultra high vacuum (UHV)** once in the experiment chamber at the synchrotron. The **UHV** was of the order of $10^{-11}$mbar to avoid surface pollution as well as scattering of the emitted photoelectrons. The samples are cleaved in this environment, resulting in a clean 2D surface as shown in Fig. 3.3. At both synchrotrons the measurements were carried out at low temperature, $15\text{K} \le T \le 30\text{K}$. Unfortunately, at the Photon Factory the setup did not allow to control the temperature precisely, and at SOLEIL the temperature controller broke down during the beam time, so that a precise temperature dependent measurement could not be carried out.

Since the maximum momentum resolution is obtained along the slit of the analyzer, we first corrected the orientation of the sample with respect to the analyzer. In practice, the sample is rotated so that either the slit is parallel to the nodal direction (1, 1, 0) or to the anti-nodal direction. Finally, the Fermi energy is determined using the gold sample attached to the sample holder. The author carried the data analysis of the experimental data and created all the experimental plots shown in the following chapter. We used the pyARPES package [257] which the author adapted to the beamlines where the measurements were performed. For more details about the ARPES formalism, its link to the Green's function calculated in CDMFT, and how the raw data is transformed back to momentum space, we refer the reader to App. A.

## 3.3 SPECTRAL FUNCTION

Both the **PG** and the waterfall feature can be observed in the momentum-resolved spectral function, which is directly calculated in **CDMFT** and obtained in **ARPES** (albeit the additional matrix-element effects [258]). In the following, we start by studying the parent compound **NaCCOC** ($x=0$) and show that the waterfall feature seen by Ronning *et al.* [86] is well captured by **CDMFT** calculations. Using two different simpler effective models, we interpret the waterfall feature as a cross-over between low-energy spin polarons and higher-energy localized electrons. Then, we proceed with the hole-doped system for which the waterfall feature is still observed and a clear **PG** is captured by **CDMFT**, independently of the cluster size.

### 3.3.1 PARENT COMPOUND

In Fig. 3.4 is shown the spectral function calculated with the 2×2 plaquette cluster (see Fig. 3.2(c)) at $\beta = 60\text{eV}^{-1}$ with two different initializations, obtained by performing an analytic continuation of the self-energy using the MQEM method [217] presented in Sec. 2.5. All the one-particle spectra shown in this chapter were obtained using the same analytic continuation procedure. In both cases, either a metallic or insulating starting point, the final system is a Mott insulator with a gap $\Delta \simeq 2.5\text{eV}$, in good agreement with the scanning tunneling microscopy (STM) estimate [89]. In the occupied region, apart from the global shift of the chemical potential, the resulting spectra show drastic differences: in Fig. 3.4(a1)-(a3) we observe a low-energy dispersion around −0.5eV which is well separated from the LHB which lies around −2eV. In contrast, we see in Fig. 3.4(b1)-
<figure><img src="image_24.png" /><figcaption>Figure 3.4: Spectral functions calculated with the plaquette cluster at $\beta = 60\text{eV}^{-1}$ with different initializations: starting (a1)-(a3) from a metallic state, (b1)-(b3) from an insulating state. Note that our set of parameters is close to the coexistence region where both the metallic and insulating phases are possible solutions. The final chemical potential is different, hence the static shift in energy of the two Hubbard bands depending on the starting point. The solid and dotted lines were extracted from the Hubbard-I and SCBA spectra and were shifted by a constant energy so to match the CDMFT dispersion. The orange circles depict the experimental dispersion measured by Ronning et al. on undoped NaCCOC ($x = 0$) and are extracted from Ref. 86. The circles have been shifted in energy for the panels (b3) to compensate for the difference in chemical potential.</figcaption></figure>
(b3) that the insulating starting point leads to a rather continuous LHB dispersing from −2.5eV to −1.0eV, except at the nodal point between Γ and M where a low-energy dispersion distinct from the LHB appears.

Despite these differences, both calculations capture the waterfall feature at the nodal point which is signaled by a kink in the LHB dispersion: around −0.5eV and −1.25eV respectively for the metallic and insulating solution. As shown in Fig. 3.4(a2),(b2), it can be understood as a cross-over between two different energy regimes. At high-energy, the LHB is well captured by a Hubbard-I approximation [35] in which the self-energy is approximated by the atomic limit result, given for a general filling and chemical potential by [259]:

$$ \Sigma_{\sigma} = \frac{1}{2}nU + \frac{\frac{1}{2}n(1 - \frac{1}{2}n)U^2}{\omega + \mu - (1 - \frac{1}{2}n)U + i\delta} \quad (3.3) $$

where n is the electron filling. Within this approximation of purely local self-energy, the LHB is well reproduced around Γ where the energy is the lowest. Since this approximation is valid in the atomic limit, we observe in Fig. 3.4 that it works best for the insulating starting point.

At lower energy however this approximation of the self-energy is not valid: the effective mass of the electrons is increased leading to a kink in the dispersion. The purely local self-energy of the Hubbard-I approximation can not capture such a behavior since it is related to non-local antiferromagnetic fluctuations [95–98]. To illustrate that, we perform SCBA calculations following the formalism described in Sec. 2.6 using the specific parameters of NaCCOC. The low-energy dispersion of the SCBA spectrum is extracted and compared to the CDMFT calculations in Fig. 3.4(a2)-(b2). The spin-polaron dispersion is strongly renormalized due to the interaction with the surrounding spins and the dispersion becomes rather flat in comparison to the Hubbard-I approximation. The bandwidth of the spin-polaron feature is of the order of ~ 0.4eV. The low-energy part of the spectrum around the nodal point (π/2, π/2) is well captured by the spin-polaron physics, thus providing a clear explanation for the waterfall feature in terms of a cross-over between a high-energy region where the correlations are purely local and a low-energy part governed by antiferromagnetic spin fluctuations. Note that while the SCBA calculation relies on long-range antiferromagnetic order, the CDMFT calculations here only incorporate short-range fluctuations, hence showing that the kink feature can survive in the hole-doped region in which the long-range spin correlations are suppressed. Note also that while in the SCBA spectrum as shown in Fig. 2.12 the spin-polaron has a large lifetime at the nodal point, in the CDMFT spectrum the feature is less coherent: the system is not a perfect antiferromagnet. Although we did not perform this analysis here, the comparison of the lifetime of the spin-polaron in theory and experiment could provide interesting information on the system.

Most interestingly, the CDMFT spectrum is in excellent agreement with the experimental measurements [85, 86, 91]. Both starting points very well reproduce the ARPES data of Ronning et al. [86] as shown in Fig. 3.4(a3)-(b3). While the high-energy part of the Γ-M dispersion seems to be better captured by the metallic starting point, the spin-polaron part is astonishingly agreeing for the in-
<figure><img src="image_25.png" /><figcaption>Figure 3.5: Spectral functions calculated with the plaquette cluster at $\beta = 20\text{eV}^{-1}$ with different initializations: starting (a) from a metallic state, (b) from an insulating state. The final chemical potential is different, hence the static shift in energy of the two Hubbard bands depends on the starting point.</figcaption></figure>

sulating initialization. Both in the theoretical and experimental spectra we observe a downturn of the low energy dispersion at the nodal point $(\pi/2, \pi/2)$, i.e. at the crossing with the antiferromagnetic BZ boundary, thus being another sign of the role of these fluctuations in the low-energy physics of undoped NaCCOC. The bandwidth of the spin-polaron is another argument for the identification of the experimental waterfall: The experimental and theoretical (both for CDMFT and SCBA) bandwidths ($\sim 0.4\text{eV}$) are in excellent agreement. One of the main differences between the two theoretical spectra is seen around the anti-nodal point, where CDMFT predicts either a low-energy flat dispersion at the same-energy as the nodal point (metallic starting point), or a high-energy dispersion stemming from the atomic limit (insulating starting point). The second option seems to better account for the experimental data of Ref. [85] (most easily visualized in Fig. 2 of Ref. [91]) since a clear difference in energy of roughly $-0.5\text{eV}$ is observed between the nodal and anti-nodal point. Finally, one can note that the waterfall feature is similar to the spin-charge signature observed with ARPES in 1D spin chains [260]. Despite the fact that these two phenomena are related to the cross-over between two different energy regimes (atomic limit vs spin polaron, and spinon vs holon), the physics at play is different: in NaCCOC the motion of the hole is lead by the coupling to the spin degrees of freedom, while in the 1D spin chains the charge of the hole become entirely decoupled from the spin excitation.

Fig. 3.5 shows the spectral functions calculated with the two different initialization but for a higher temperature ($\beta = 20$) at which the spin fluctuations start to be suppressed. With the metallic initialization (Fig. 3.5(a)), we still observe a low-energy feature but there is a strong spectral weight transfer to the LHB as compared to the $\beta = 60\text{eV}^{-1}$ calculation. It is even more striking for the insulating starting point (Fig. 3.5(b)) where the kink in the dispersion at the nodal point has completely disappeared. This is an additional sign that the waterfall feature observed in the parent compound originates from the spin fluctuations, and shows the utmost importance of performing the calculations at sufficiently low temperature to find an agreement with the low-temperature ARPES measurements.
<figure><img src="image_26.png" /><figcaption>Figure 3.6: Experimental and theoretical momentum resolved spectral functions. (a) ARPES measure-
ment on the x = 0.05 and calculated spectra for x = 0.05 and x = 0.1 at β = 60. The Hubbard-I
approximation is no longer shifted in energy, but the SCBA dispersion still is. (b) Zoom at the
nodal point to compare the experimental and theoretical results for both dopings. The pur-
ple circles denote the maximum of the theoretical spectral function at each k-point, and the
relative spectral weight is depicted by the circle size.</figcaption></figure>

3.3.2 HOLE DOPED SYSTEM

At finite hole doping, the chemical potential (and the spectral function) is no longer dependent
on the starting point. Indeed, the coexistence region vanishes with doping: only the metallic
solution survives. All the theoretical results presented in the following for the hole-doped NaC-
COC are obtained with the insulating (at half-filling) starting point. The experimental data has
been obtained with a x = 0.05 sample at the Photon Factory synchrotron, and a x = 0.1 sample
measured at SOLEIL synchrotron. For the x = 0.05 sample, the experimental setup lacked of a
gold sample for the Fermi energy reference, such that the later was obtained by a fit of the occu-
pied bands of NaCCOC. To be consistent, the Fermi energy for the x = 0.1 sample was corrected
according to the method used for the x = 0.05.

The momentum-resolved spectral functions along Γ − X − M − Γ are shown in Fig. 3.6(a) for
the x = 0.05 sample and calculations at β = 60eV⁻¹ for x = 0.05, 0.1 with the plaquette cluster.
We observe a clear coherent dispersion at the nodal point for the x = 0.05 sample, along with a
gapped weaker spectral weight around the anti-nodal point. This is in good agreement with the
calculated spectrum for x = 0.05: a coherent quasiparticle peak is seen at the nodal point and
two split flat dispersions appear at the anti-nodal point below and above the Fermi level. This is
a clear sign of the opening of a PG, which starts to close at higher doping concentrations x = 0.1
in the calculations. Interestingly, the waterfall feature survives in the Γ − M dispersion since the
spectrum is still well described by a cross-over between a Hubbard-I approximation and a SCBA
dispersion. Note that the Hubbard-I dispersion needs no longer to be shifted to adapt to the right
<figure><img src="image_27.png" /><figcaption>Figure 3.7: (a) Theoretical and experimental Fermi surface for the $x = 0.05$ and $x = 0.1$ samples, the ARPES measurement is shown in the top right corner of the BZ. The white rotated square illustrates the antiferromagnetic BZ. (b) EDC around the Fermi arc (the definition of the angle $\alpha$ is shown in (a)) for the $x = 0.05$ sample. (c) EDC around the nodal point on the path $\Gamma - M$. The orange (blue) arrows highlight the high- (low-) energy peaks.</figcaption></figure>

chemical potential, hence showing the robustness of this qualitative picture. At the anti-nodal point where we observe a PG, neither the Hubbard-I nor the SCBA give a good agreement with the CDMFT spectrum since yet another physical mechanism is at play and which is the subject of the next section.

To further compare the experimental and theoretical spectra, we show in Fig. 3.6(b) a zoom on the $\Gamma - M$ dispersion. Similarly to the undoped case, we observe an excellent agreement in the position of the dispersion and the kink which appears around 0.3eV, showing that the antiferromagnetic fluctuations survive at those values of doping. The increase of spectral weight as $k$ gets closer to the nodal point ($\pi/2, \pi/2$) is seen in all the spectra. Note that these data show an excellent agreement with the previous ARPES measurements [26, 90–93]. However, one can notice that in the experimental spectrum, the coherent and intense quasiparticle at $\omega = 0$ captured by CDMFT is not seen, which was also observed by Shen et al. for dopings below $x = 0.10$ [26, 92]. This is best seen by inspecting the energy distribution curve (EDC) around the nodal point, as shown in Fig. 3.7.

We show in panel (a) the Fermi surfaces for the two dopings for which we observe the famous Fermi arc, hallmark of the PG in photoemission data [24], both in the experimental and the theoretical data. The opening of the PG can be directly observed in the EDC around the nodal point, which are presented in Fig. 3.7(b) for the $x = 0.05$ sample. It is clear that in contrast to the CDMFT calculation, no coherent peak is observed for the $x = 0.05$ sample at the nodal point, although a PG is obtained due to a shift and broadening of the spectral function at low energy as going towards the anti-nodal point. Similarly, the quasiparticle peak in the ARPES measurements of the $x = 0.1$ sample is also absent (not shown here), which is a sign that the doping could be overestimated, based on the observations of Shen et al. in Ref. 92. The missing quasi-
<figure><img src="image_28.png" /><figcaption>Figure 3.8: (a) Spectral function at the anti-nodal point (π, 0) for different dopings (0.0 ≤ x ≤ 0.1) at β = 60eV⁻¹. (b) Spectral function at the nodal (π/2, π/2) and anti-nodal point (π, 0) (multiplied by a factor 10) at two temperatures below and above the PG transition for x = 0.05.</figcaption></figure>

particle was attributed to a coupling to phonons which are not included within our theoretical
framework. We do not infer on the validity of this proposition since this would require a better
characterization of many samples and a comparison between a larger number of them, as well
as the modification of our simple effective model to allow for the coupling of the electrons to
additional degrees of freedom.

When a quasiparticle peak is present in the photoemission data, a double peak structure at the nodal point is observed [91, 92]. As shown in Fig. 3.7(c), the double peak structure is reproduced, although the peak separation is much sharper than in experiment. The quasiparticle peak forms at low-energy as the dispersion goes from Γ to M, like in the measured data, and it is backfolded when crossing the nodal point. The latter is not clearly seen in experiment, which may come from the quasi-particle suppression, or an overestimation of the strength of the antiferromagnetic spin fluctuations in the 2 × 2 plaquette.

Though, overall the CDMFT calculations yield an excellent agreement with the ARPES mea-
surements, and capture the essential ingredients of the hole-doped NaCCOC spectral function:
the opening of a PG at low doping and a waterfall feature at the nodal point at low-energy. It is
furthermore insightful to study the doping and temperature dependence of the PG and the quasi-
particle peak at the nodal point. In Fig. 3.8(a) is shown the doping dependence of the spectral
function at the anti-nodal point (π, 0) at β = 60eV⁻¹. The undoped system shows a broad peak
around −0.5eV as we have seen previously. As the doping is increased, and even for the lowest
doping (x = 0.01), we see an asymmetric double peak structure appearing with two peaks located
away from the Fermi level, roughly at ω = ±0.15eV. The difference between the occupied and
<figure><img src="image_29.png" /><figcaption>Figure 3.9: (a) calculated spectral function for $x = 0.05$ at $\beta = 40\text{eV}^{-1}$. The Hubbard-I approximation is slightly shifted for the sake of comparison in energy, as well as the SCBA dispersion. (b) Comparison between the experimental and theoretical spectrum for the $x = 0.05$ sample around the nodal point. The purple circles denote the maximum of the theoretical spectral function at each k-point, and the relative spectral weight is depicted by the circle size.</figcaption></figure>

empty regions come from the particle-hole asymmetry, and is therefore reduced as the doping is increased. Upon increasing the hole concentration, the gap fills up until an almost closing of the PG, but the two peaks below and above the Fermi level do not collapse, consistently with previous observations [261]. Note that the peaks position may vary between the different dopings due to systematic errors from the analytic continuation.

The temperature dependence of the spectral function provides useful information about the possible origin of the PG. In Fig. 3.8(b) is shown the spectral function both at the nodal and anti-nodal point for $x = 0.05$ at two different temperatures: $\beta = 20, 60\text{eV}^{-1}$. Similarly to the undoped case, we observe that for high temperature the anti-nodal PG disappears, as well as the double peak structure at the nodal point. This may indicate an influence of the spin fluctuations on the PG [75, 104, 107, 111–114, 246] since the double peak structure at the nodal point is related to the spin-polaron physics as shown in Fig 3.6(a) where we observe that the coherent quasiparticle peak dispersion corresponds well to that of a spin-polaron. It is the subject of the next section, in which we dissect the CDMFT result to search for the physical origin of the PG in the 2D Hubbard model, and thus, by extension, in cuprates.

Before proceeding to such an analysis, it is important to note that the waterfall feature and the PG are not a peculiar consequence of the 2 × 2 cluster. In Fig. 3.9 is shown the CDMFT spectral function obtained with an 8-site cluster (see Fig. 3.2(c)) at $\beta = 40\text{eV}^{-1}$ for $x = 0.05$, and its comparison to the experimental spectrum. Because the 8-site cluster is breaking the $x/y$ symmetry since it is elongated along one or the other direction (see Fig. 3.2(c)), we restore the symmetry when performing the Green's function reperiodization by averaging between the x- and y-oriented clusters. We first notice that both the waterfall feature and the PG are captured by this larger cluster, in agreement with the findings by Gull *et al.* that the cluster size (though in DCA) does not influence much the main characteristics of the spectrum [110]. The Hubbard-I approximation accounts well for the high-energy part of the dispersion, although a small constant shift of 0.15eV has been applied. Similarly, the spin-polaron captures rather well the low-energy part
<figure><img src="image_30.png" /><figcaption>Figure 3.10: (a) EDC along the Fermi arc following the convention of Fig. 3.7(b) obtained with the 8-site cluster at $\beta = 40\text{eV}^{-1}$ for $x = 0.05$. (b) EDC around the nodal point along the path $\Gamma - M$. The orange (blue) arrows highlight the high- (low-) energy peaks.</figcaption></figure>

of the spectrum, although it seems to overestimate slightly the effective mass of the electrons, which is to be expected since in the doped system the antiferromagnetic fluctuations are weaker than at half-filling. The $\Gamma - M$ dispersion is again in excellent agreement with the experimental measurements with a clear kink in the dispersion around 0.3eV.

In Fig. 3.10 we show the EDC calculated with the 8-site cluster at $\beta = 40\text{eV}^{-1}$ for $x = 0.05$. Interestingly the EDC along the Fermi arc show a different qualitative behavior than that of the plaquette cluster: for the latter a secondary peak appears at the nodal point around $-0.1\text{eV}$ while the low-energy peak is already formed, whereas for the 8-site cluster the low-energy peak arises at the nodal point from the higher-energy feature which disperses between $-0.2, -0.1\text{eV}$. We come back on this difference later in the chapter. The 8-site cluster is in better agreement with the measured spectrum of Ref. 91,92 along the path $\Gamma - M$, as shown in Fig. 3.10(b). Indeed, there is no longer a clear separation between the two peaks at the nodal point, and the absence in experiment of backfolding of the low-energy peak when crossing the nodal point can be explained by the fact that it happens in the occupied part of the spectrum, as shown in Fig. 3.10(b). This is a clear improvement of the 8-site cluster as compared to the plaquette, since with the latter the backfolded quasiparticle peak always lies below the Fermi level.

## 3.4 CLUSTER ANALYSIS

The analysis of the waterfall feature in terms of a crossover between a localized limit at high energy and a spin-polaron at low-energy, as well as the correspondence between the closing of the PG and the disappearance of the waterfall feature at high temperature point to the importance of
the short-range spin fluctuations in describing the low-energy physics of the 2D Hubbard model, and by extension to that of the cuprates. These observations are in line with the existing literature finding hallmarks of dominating short-range spin fluctuations [75, 104, 107, 111–114, 246]. In the following, we aim at deepening further the analysis of the CDMFT calculations by investigating the cluster state histogram which can be extracted from the impurity solver. The histogram provides the probability for the plaquette connected to a self-consistent electron bath to be in specific states. We dissect the results for the 2 × 2 plaquette which is simple enough to provide insightful information at a low cost, and we show, following the argument of Merino and Gunnarsson [121], that its simple state histogram can provide a scenario for the opening of a PG and that it is not limited to the peculiar 2 × 2 plaquette cluster. We emphasize the importance of this analysis: it is essential to understand the mechanisms leading to the opening of the PG in the CDMFT calculations since it will allow to extract precious information about the PG phase in the 2D Hubbard model, and therefore in cuprates.

### 3.4.1 STATE HISTOGRAM

In order to obtain the state histogram, we performed calculations using the CT-HYB solver from TRIQS [119, 120]. Note that the state histograms were obtained with $U = 10t$, in contrast to the usual $U = 10.2t$ in this chapter. This small difference in $U$ has no consequence on the validity of the analysis. All the other quantities discussed in the following were obtained for $U = 10.2t$ as found from the fit of the RIXS data.

Instead of using the site basis, it is more convenient to switch to the cluster momenta basis in which the observables are diagonal in the 2 × 2 cluster. It is defined as follow:

$$ 
\begin{aligned}
c_{(0,0),\sigma}^{\dagger} &= \frac{1}{2}(c_{1,\sigma}^{\dagger} + c_{2,\sigma}^{\dagger} + c_{3,\sigma}^{\dagger} + c_{4,\sigma}^{\dagger}) \\
c_{(\pi,\pi),\sigma}^{\dagger} &= \frac{1}{2}(c_{1,\sigma}^{\dagger} - c_{2,\sigma}^{\dagger} + c_{3,\sigma}^{\dagger} - c_{4,\sigma}^{\dagger}) \\
c_{(\pi,0),\sigma}^{\dagger} &= \frac{1}{2}(c_{1,\sigma}^{\dagger} + c_{2,\sigma}^{\dagger} - c_{3,\sigma}^{\dagger} - c_{4,\sigma}^{\dagger}) \\
c_{(0,\pi),\sigma}^{\dagger} &= \frac{1}{2}(c_{1,\sigma}^{\dagger} - c_{2,\sigma}^{\dagger} - c_{3,\sigma}^{\dagger} + c_{4,\sigma}^{\dagger}),
\end{aligned} 
\quad (3.4) $$

where 1, 2, 3, 4 denotes the plaquette sites labeled clockwise, see Fig. 3.11(a). Although this basis is convenient, in this chapter we prefer to replace the usual two $c_{(\pi,0),\sigma}^{\dagger}$, $c_{(0,\pi),\sigma}^{\dagger}$ by their bonding/anti-bonding combination which are defined by:

$$ 
\begin{aligned}
c_{D13,\sigma}^{\dagger} &= \frac{1}{\sqrt{2}}(c_{(\pi,0),\sigma}^{\dagger} + c_{(0,\pi),\sigma}^{\dagger}) = \frac{1}{\sqrt{2}}(c_{1,\sigma}^{\dagger} - c_{3,\sigma}^{\dagger}) \\
c_{D42,\sigma}^{\dagger} &= \frac{1}{\sqrt{2}}(-c_{(\pi,0),\sigma}^{\dagger} + c_{(0,\pi),\sigma}^{\dagger}) = \frac{1}{\sqrt{2}}(c_{4,\sigma}^{\dagger} - c_{2,\sigma}^{\dagger})
\end{aligned} 
\quad (3.5) $$

where the $D13$, $D42$ orbitals are diagonal bonding orbitals between sites 1, 3 and 4, 2. We therefore analyze the plaquette state histogram with the basis set $\{(0, 0), (\pi, \pi), D13, D42\}$.

For the isolated cluster, these four orbitals have different on-site energies illustrated in Fig. 3.11(a):
<figure><img src="image_31.png" /><figcaption>Figure 3.11: (a) Sketch of the on-site energies of the isolated plaquette cluster. (b) Doping and temperature dependent state histogram calculated with $U = 10t$. (c) Sketch of the three leading groups of states.</figcaption></figure>

at half-filling the (0, 0) orbital is expected to be doubly occupied, while the orbitals D13 and D42 will share the two electrons left since they are degenerate, and the orbital $(\pi, \pi)$ is empty. In agreement with this simple picture, and with previous works [48, 247, 262], we find that the state histogram is dominated by only a few states which can be grouped into three categories, as shown in Fig. 3.11(b)-(c). These states, reduced to their most important components, can be approximately written

$$ 
\begin{aligned}
|N = 4, S = 0\rangle &\approx |\uparrow\downarrow, 0\rangle \otimes (|\uparrow, \downarrow\rangle - |\downarrow, \uparrow\rangle) \\
|N = 4, S = 1\rangle &\approx |\uparrow\downarrow, 0\rangle \otimes |\uparrow, \uparrow\rangle \\
|N = 3, S = 1/2\rangle &\approx |\uparrow\downarrow, 0\rangle \otimes |\uparrow, 0\rangle
\end{aligned} 
\qquad (3.6) 
$$

where we adopt the convention $|\phi\rangle = |(0, 0), (\pi, \pi)\rangle \otimes |D13, D42\rangle$, i.e. we separate $(0, 0), (\pi, \pi)$ and D13, D42 orbitals for convenience, and we show only one of the two (four) possible configurations for the states $|N = 4, S = 1\rangle$ ($|N = 3, S = 1/2\rangle$) (see Fig. 3.11(c)).

The most likely state is $|N = 4, S = 0\rangle$, which corresponds to a superposition of resonating valence bond (RVB) singlets [263] between NN [121, 247]. Interestingly, within our new choice of basis, this state involves a singlet between the D13 and D42 orbitals. As expected from the isolated plaquette limit, this state dominates the histogram at low doping, and its contribution increases with decreasing temperature, showing that the short range antiferromagnetic fluctuations play a crucial role in the physics of the PG. At high temperature ($\beta = 30eV^{-1}$), the second group of states is $|N = 4, S = 1\rangle$ which corresponds to a spin flip excitation: going from $|N = 4, S = 0\rangle$ to $|N = 4, S = 1\rangle$ involves for instance the spin flip term $c_{D13,\sigma}^{\dagger}c_{D42,-\sigma}$, which carries a spin $S = 1$ and a momentum $Q = (\pi, \pi)$. As the doping is increased, i.e. as electrons are removed from the
lattice, this group's contribution shows an increase until $x \sim 0.1$ where it saturates and slowly drops as the doping is further increased. Therefore, at intermediate dopings the $N = 4$ configurations are increasingly populated by these $S = 1$ excitations, showing that the singlet state breaks into $S = 1$ triplet states. At large dopings, the dominating group of states is $|N = 3, S = 1/2\rangle$, which simply corresponds to the removal of an electron in the D13, D42 orbitals.

Since the D13, D42 are directly related to the $(\pi, 0), (0, \pi)$ orbitals, and since they are half-filled, we can expect that the low-energy physics close to the anti-nodal point is related to the physics of these orbitals which forms a singlet in the state $|N = 4, S = 0\rangle$. Similarly, from this simple picture we expect the $(0, 0)$ orbital to account for the LHB dispersion around the $\Gamma$ point, and from the $(\pi, \pi)$ orbital to capture the unoccupied part of the spectrum close to $M$. In fact, this can be directly visualized in the spectral function by decomposing it into different contributions from the $(0, 0), (\pi, \pi)$, D13, D42 orbitals. To this aim, the site-basis component of the Green's function $\tilde{G}_{\sigma}$ can be expressed with the diagonal components of the Green's function $\tilde{G}_{\sigma}$ in the rotated basis:

$$ G_{\sigma} = P \begin{pmatrix} \tilde{G}_{(0,0),\sigma} & 0 & 0 & 0 \\ 0 & \tilde{G}_{D13,\sigma} & 0 & 0 \\ 0 & 0 & \tilde{G}_{D42,\sigma} & 0 \\ 0 & 0 & 0 & \tilde{G}_{(\pi,\pi),\sigma} \end{pmatrix} P^{-1} \quad P = \frac{1}{2} \begin{pmatrix} 1 & \sqrt{2} & 0 & -1 \\ 1 & 0 & -\sqrt{2} & 1 \\ 1 & -\sqrt{2} & 0 & -1 \\ 1 & 0 & \sqrt{2} & 1 \end{pmatrix}. \quad (3.7) $$

The Green's function reperiodization (Eq. 2.48) can then be re-written in terms of the diagonal $\tilde{G}_{(0,0),\sigma}, \tilde{G}_{D13,\sigma}, \tilde{G}_{D42,\sigma}, \tilde{G}_{(\pi,\pi),\sigma}$ components as:

$$ G_{\sigma}^{\text{per}}(\mathbf{k}, i\omega_n) = \sum_{m} a_{m}(\mathbf{k}) \tilde{G}_{m,\sigma}, \quad (3.8) $$

where $m$ denote the orbitals $\{(0, 0), (\pi, \pi), D13, D42\}$ and $a_m(\mathbf{k})$ are momentum-dependent weights for each of the components.

The result of such decomposition is shown in Fig. 3.12 for the $x = 0.05$ and $x = 0.10$ that we have already discussed previously, i.e. calculated at $\beta = 60\text{eV}^{-1}$ with the plaquette cluster. As expected, the $(\pi, \pi)$ contribution is almost limited to the occupied region around the $M$ point, some incoherent spectral weight can be noticed in the occupied part due to the small but non-zero filling of this orbital. The $(0, 0)$ contribution is substantial around $\Gamma$ and is therefore responsible for the atomic limit behavior. Interestingly, it also accounts for most of the spectral weight at the nodal point $(\pi/2, \pi/2)$ and a clear waterfall feature is seen, i.e. the $(0, 0)$ also hosts the coherent spin-polaron at low-energy. Most importantly, its contribution at the anti-nodal point $(\pi, 0)$ vanishes, leaving the D13 and D42 as the only non-zero terms in this region. This provides a hint for why the SCBA dispersion is in good agreement at the nodal point, but always seems to be shifted in energy at the anti-nodal point: these two regions are in fact determined by electrons living on different orbitals. The D13, D42 orbitals are responsible for the PG physics: they display two rather flat dispersions around the anti-nodal point. Note that they also have a non-zero contribution at the nodal point, such that the double peak feature observed in the $2 \times 2$ plaquette results is a consequence of the periodization procedure, which at the nodal point superposes the
<figure><img src="image_32.png" /><figcaption>Figure 3.12: Momentum-resolved spectral function and its decomposition into the different diagonal $\tilde{G}_{(0,0),\sigma}$, $\tilde{G}_{D13,\sigma}$, $\tilde{G}_{D42,\sigma}$, $\tilde{G}_{(\pi,\pi),\sigma}$ components for (a) $x = 0.05$ and (b) $x = 0.1$, both calculated at $\beta = 60\text{eV}^{-1}$ with the plaquette cluster.</figcaption></figure>

contribution from (0, 0) and the *D13*, *D42* orbitals. We shall get back to this point later with the 8-site cluster.

This simple analysis allows us to make a difference between the key features of the spectral func-
tion. Most importantly, it shows that the **PG** is related to the behavior of the *D13*, *D42* orbitals
which have the tendency to form a singlet state at low doping. If the cluster is isolated, this singlet
state naturally leads to an opening of a gap at the anti-nodal point (one can think of it as a dimer).
However, when it is connected to a bath of non-interacting electrons, it is not clear why the gap
would remain open: this is the case of the (0, 0) orbital for instance, which in the isolated cluster
is supposed to be doubly occupied but shows a large quasiparticle peak at the Fermi energy in the
**CDMFT** spectral function! In the following, based on the derivation of Merino *et al.* [121] that we
adapt to our specific choice of basis, we show that the gap remains open due to the very nature of
the non-degenerate singlet state using a simple effective model. Moreover, we also discuss why
this mechanism should hold for the 8-site cluster, and therefore for larger systems.

3.4.2 ORIGIN OF THE PSEUDOGAP

Despite the fact that the isolated *D13*, *D42* orbitals would show a gap since they form a singlet state, it is not straightforward that the gap remains open when coupling to a continuous bath of electrons. Indeed, the latter allows photoemission processes from a state with *N* electrons on the cluster to the same state with *N* electrons on the cluster and one missing bath electrons, thus giving a zero energy contribution. Following the derivation of Ref. 121 that we specialize for the singlet state |*N* = 4, *S* = 0⟩, we show in this section that in the limit of large *U*, the non-degenerate singlet ground state leads to a destructive interference of such processes, so that the gap remains non-zero even if the system is coupled to a continuous bath of electrons.

To derive a simple effective model, we will approximate the *D13*, *D42* orbitals as disconnected
from the (0, 0), (π, π) orbitals such that they form a two-level system connected to a continuous
bath of non-interacting electrons. This is a reasonable approximation given that the Green's func-
<figure><img src="image_33.png" /><figcaption>Figure 3.13: Non-interacting and interacting (obtained from the last CDMFT calculation) imaginary part of the hybridization function on the Matsubara axis $\Delta(i\omega_n)$ and DOS on the real axis, for both the plaquette and 8-site cluster.</figcaption></figure>

tion in this basis is purely diagonal. The details of the Hamiltonian governing the isolated D13, D42 orbitals are not important, as long as they lead to a ground state of the form (for the isolated orbitals):

$$|GS\rangle = \frac{1}{\sqrt{2}}(|\uparrow, \downarrow\rangle - |\downarrow, \uparrow\rangle), \qquad (3.9)$$

consistently with the state $|N = 4, S = 0\rangle$. In the $2 \times 2$ cluster, it is the short range antiferromagnetic fluctuations that causes the emergence of the RVB-like singlet state.

The coupling to the continuous bath of non-interacting electrons can be written:

$$H_{hyb} = \sum_{m\sigma,\epsilon} V_{\epsilon}(d_{m\sigma,\epsilon}^{\dagger}c_{m\sigma} + \text{h.c.}) \qquad (3.10)$$

where $d_{m\sigma,\epsilon}^{\dagger}$, $d_{m\sigma,\epsilon}$ are the creation and annihilation operators of a bath electron of orbital $m$, spin $\sigma$ and energy $\epsilon$, $c_{m\sigma}^{\dagger}$, $c_{m\sigma}$ the ones of the cluster electrons, and $V_{\epsilon}$ the coupling between the bath and the cluster, taken to be orbital independent since they are equivalent. The orbital index $m$ refers either to D13 or D42. Since we are in the limit $U \gg t$, we can also consider that $U \gg V_{\epsilon}$, such that the hybridization to the bath may be treated in perturbation theory. This is justified by inspecting the hybridization function of the different orbitals of the plaquette both for the non-interacting case, and for the interacting case obtained after the CDMFT self-consistent cycle (i.e. with correlated effects already incorporated), as shown in Fig. 3.13. The hybridization function has a maximum of -0.5 for the D13, D42 orbitals in the non-interacting case, which is further reduced after the CDMFT self-consistent calculations, and is in both cases much lower than $U$. While the non-interacting DOS shows a van Hove singularity for the D13, D42 orbital [247], after the CDMFT cycle it shows a PG, while the (0, 0) orbital shows a coherent peak at the Fermi energy.

In the ground state (Eq. 3.9), the two orbitals are half-filled and there are two cluster electrons. When connecting the cluster to the bath, the number of electrons on the cluster is not fixed, but
that of the full system is conserved. At first order, the correction to the ground state originates from states formed by the hopping of a cluster electron to the bath, or the reverse, i.e. states of the form:

$$ \begin{aligned} d_{m\sigma,\epsilon}^{\dagger} c_{m\sigma} |GS\rangle \\ c_{m\sigma}^{\dagger} d_{m\sigma,\epsilon} |GS\rangle, \end{aligned} \qquad (3.11) $$

which are states with $N-1$ and $N+1$ electrons on the cluster, and where we take the convention that the $|GS\rangle$ is a product between the cluster and the bath states (not explicitly written for the sake of simplicity). The corrections to the ground state will involve matrix elements of the type:

$$ \begin{aligned} \langle GS | d_{m\sigma,\epsilon}^{\dagger} c_{m\sigma} V_{\epsilon} c_{m'\sigma',\epsilon}^{\dagger} d_{m'\sigma',\epsilon} | GS \rangle &= \delta_{mm',\sigma\sigma',\epsilon\epsilon'} \langle GS | d_{m\sigma,\epsilon}^{\dagger} c_{m\sigma} V_{\epsilon} c_{m'\sigma',\epsilon}^{\dagger} d_{m'\sigma',\epsilon} | GS \rangle \\ &= V_{\epsilon} \langle GS | n_{m\sigma,\epsilon}^{d} (1 - n_{m\sigma}^{c}) | GS \rangle, \end{aligned} \qquad (3.12) $$

for the $N-1$ contribution, and for the $N+1$:

$$ V_{\epsilon} \langle GS | n_{m\sigma}^{c} (1 - n_{m\sigma,\epsilon}^{d}) | GS \rangle. \qquad (3.13) $$

We consider in the following:

$$ V_{\epsilon} \langle GS | n_{m\sigma,\epsilon}^{d} (1 - n_{m\sigma}^{c}) | GS \rangle \approx V_{\epsilon} \langle GS | n_{m\sigma}^{c} (1 - n_{m\sigma,\epsilon}^{d}) | GS \rangle = \tilde{V}_{\epsilon}, \qquad (3.14) $$

such that both $N+1$ and $N-1$ correction terms have the same amplitude, which is justified since the self-consistent electron bath is supposed to be at the same relative filling as the cluster.

We can now write the corrected ground state as:

$$ \begin{aligned} |\psi\rangle &= |GS\rangle + \sum_{m\sigma,\epsilon} \frac{\tilde{V}_{\epsilon}}{\Delta_{+} - \epsilon} c_{m\sigma}^{\dagger} d_{m\sigma,\epsilon} |GS\rangle + \sum_{m\sigma,\epsilon} \frac{\tilde{V}_{\epsilon}}{\Delta_{-} + \epsilon} d_{m\sigma,\epsilon}^{\dagger} c_{m\sigma} |GS\rangle \\ &= |GS\rangle - \sum_{m\sigma,\epsilon} \frac{\tilde{V}_{\epsilon}}{\Delta_{+} - \epsilon} d_{m\sigma,\epsilon} c_{m\sigma}^{\dagger} |GS\rangle + \sum_{m\sigma,\epsilon} \frac{\tilde{V}_{\epsilon}}{\Delta_{-} + \epsilon} d_{m\sigma,\epsilon}^{\dagger} c_{m\sigma} |GS\rangle, \end{aligned} \qquad (3.15) $$

where $\Delta_{\pm}$ corresponds to the difference in cluster energy between the ground state and the $N+1/N-1$ states, and we ordered the operators in the second term such that the cluster creation/an-nihilation operators always appear to the right with respect to the bath operators. In principle, as mentioned above, the coupling to the bath authorizes photoemission processes between the ground state $|GS\rangle$, and a final state like $\sim d_{m\sigma,\epsilon}|GS\rangle$ which, if it is a bath electron at the Fermi level, leads to a non-zero spectral weight at zero energy since $|GS\rangle$ and $\sim d_{m\sigma,\epsilon_F}|GS\rangle$ would have the same energy. To show that such process is in fact canceled when the ground state is of the nature of Eq. 3.9, we take the example (not restrictive) of a photoemission process that involves the removal of an electron $c_{D_{13}\uparrow}$ on the cluster. We need to calculate a matrix element of the type:

$$ \langle \phi_f | c_{D_{13}\uparrow} | \psi \rangle, \qquad (3.16) $$
where we consider the final state $|\phi_f\rangle$ composed of the cluster in its ground state and a removed bath electron at the Fermi energy:

$$|\phi_f\rangle = d_{m\sigma, \epsilon_F} |GS\rangle - \sum_{m\sigma, \epsilon} \frac{\tilde{V}_\epsilon}{\Delta_+ - \epsilon} d_{m\sigma, \epsilon} c_{m\sigma, \epsilon_F}^\dagger d_{m\sigma, \epsilon_F} |GS\rangle + \sum_{m\sigma, \epsilon} \frac{\tilde{V}_\epsilon}{\Delta_- + \epsilon} d_{m\sigma, \epsilon} c_{m\sigma} d_{m\sigma, \epsilon_F} |GS\rangle. \quad (3.17)$$

In order to calculate this matrix element, it is useful to re-write the corrected ground state using the explicit expression of $|GS\rangle$:

$$|\psi\rangle = \frac{1}{\sqrt{2}}(|\uparrow, \downarrow\rangle - |\downarrow, \uparrow\rangle) \\ + \frac{1}{\sqrt{2}} \sum_{\epsilon} \frac{\tilde{V}_{\epsilon}}{\Delta_{+} - \epsilon} [d_{D13\uparrow,\epsilon}|\uparrow\downarrow, \uparrow\rangle + d_{D13\downarrow,\epsilon}|\uparrow\downarrow, \downarrow\rangle + d_{D42\uparrow,\epsilon}|\uparrow, \uparrow\downarrow\rangle + d_{D42\downarrow,\epsilon}|\downarrow, \uparrow\downarrow\rangle] \\ + \frac{1}{\sqrt{2}} \sum_{\epsilon} \frac{\tilde{V}_{\epsilon}}{\Delta_{-} + \epsilon} [d_{D13\uparrow,\epsilon}^\dagger|0, \downarrow\rangle - d_{D13\downarrow,\epsilon}^\dagger|0, \uparrow\rangle - d_{D42\uparrow,\epsilon}^\dagger|\downarrow, 0\rangle + d_{D42\downarrow,\epsilon}^\dagger|\uparrow, 0\rangle]. \quad (3.18)$$

The action of the annihilation operator $c_{D13\uparrow}$ on the corrected ground state gives

$$c_{D13\uparrow}|\psi\rangle = \frac{1}{\sqrt{2}}|0, \downarrow\rangle + \frac{1}{\sqrt{2}}\sum_{\epsilon} \frac{\tilde{V}_{\epsilon}}{\Delta_{+} - \epsilon} [-d_{D13\uparrow,\epsilon}|\downarrow, \uparrow\rangle + \dots] + \dots \quad (3.19)$$

where we kept only the terms that will contribute to the matrix element. Similarly, keeping only the relevant elements, the final state can be written

$$|\phi_f\rangle = d_{D13\uparrow,\epsilon_F} \frac{1}{\sqrt{2}}(|\uparrow, \downarrow\rangle - |\downarrow, \uparrow\rangle) - \sum_{\epsilon} \frac{\tilde{V}_{\epsilon}}{\Delta_{-} + \epsilon} [d_{D13\uparrow,\epsilon}^\dagger d_{D13\uparrow,\epsilon} |0, \downarrow\rangle + \dots] + \dots \quad (3.20)$$

Combining Eq. 3.19 and 3.20, we obtain the following approximated expression for the matrix element [121]

$$\langle\phi_f|c_{D13\uparrow}|\psi\rangle = \frac{1}{2}\sum_{\epsilon} \frac{\tilde{V}_{\epsilon}}{\Delta_{+} - \epsilon} \langle\downarrow, \uparrow|d_{D13\uparrow,\epsilon}^{\dagger}d_{D13\uparrow,\epsilon}|\downarrow, \uparrow\rangle - \frac{1}{2}\sum_{\epsilon} \frac{\tilde{V}_{\epsilon}}{\Delta_{-} + \epsilon} \langle 0, \downarrow|d_{D13\uparrow,\epsilon}^{\dagger}d_{D13\uparrow,\epsilon}|0, \downarrow\rangle. \quad (3.21)$$

Since $\langle\downarrow, \uparrow|d_{D13\uparrow,\epsilon}^\dagger d_{D13\uparrow,\epsilon}|\downarrow, \uparrow\rangle = \langle 0, \downarrow|d_{D13\uparrow,\epsilon}^\dagger d_{D13\uparrow,\epsilon}|0, \downarrow\rangle = \delta_{\epsilon, \epsilon_F}$, in the symmetric case ($\Delta_+ = \Delta_-$) the two terms cancel each other. The very same derivation can be done for $c_{D13\downarrow}, c_{D42\uparrow}, c_{D42\downarrow}$, hence the zero-energy spectral weight is entirely suppressed. This mechanism is at the origin of the PG observed at the anti-nodal point for the 2 × 2 plaquette.

The first non-zero excitation channel will be the one relating the singlet ground state $|N = 4, S = 0\rangle$ to the lowest excited state, which is the triplet $|N = 4, S = 1\rangle$. The PG is therefore related to the breaking of the RVB-like singlets by excitations carrying a spin $S = 1$ and momentum $Q = (\pi, \pi)$. This is in agreement with the observations made via the fluctuation diagnostic technique [112–114], as well as dual boson calculations [75]. At relatively high doping, this picture breaks down since the state histogram shows substantial contributions from the degenerate set of states $|N = 3, S = 1/2\rangle$ which do not lead to a cancellation of the zero-energy excitation [121].
An important point is whether this mechanism is an artifact of the plaquette cluster, or if this
holds for larger systems. In particular, for the 8-site cluster there are not only the $(\pi, 0)/(0, \pi)$
(i.e. $D13/D42$) orbitals which are degenerate and related by $Q = (\pi, \pi)$, but also the $(\pm\pi/2, \pm\pi/2)$
so that one may expect that the same mechanism occurs also for these orbitals thus inducing a
gap at the nodal point. As shown previously, we observe a clear quasi-particle peak at the nodal
point, so that the scenario described above for the $D13$, $D42$ orbitals seems not to apply to the
$(\pm\pi/2, \pm\pi/2)$ orbitals. The reason behind this is illustrated in Fig. 3.13: the hybridization of the
$(\pm\pi/2, \pm\pi/2)$ orbitals to the bath is twice larger than the one of the $D13$, $D42$ orbitals. The iso-
lated ground state on the cluster is no longer a good approximation for these orbitals. Instead,
the $(\pm\pi/2, \pm\pi/2)$ orbitals hybridize with the bath leading to a clear quasi-particle peak in the
interacting DOS obtained after the CDMFT cycle. Interestingly, the difference between the $D13$,
$D42$ and $(\pm\pi/2, \pm\pi/2)$ orbitals is already present at the non-interacting level (see Fig. 3.13), due
to the difference in the energy dispersion [121].

From these considerations we can understand the qualitative differences in the EDC shown in
Fig. 3.7 and Fig. 3.10, respectively for the plaquette and the 8-site cluster. At the nodal point, it
is shown in Fig. 3.12 that the spectral function calculated with the plaquette is the superposition
of the contributions from the (0, 0) and $D13/D42$ orbitals. These two sets of orbitals are almost
independent from each other, and their superposition at the nodal point is the product of the pe-
riodization procedure, such that the two-peak structure observed in Fig. 3.7(c) is "artificial" and
overestimates the peak splitting. It also has a large influence on the dispersion along the Fermi
arc: in Fig. 3.7(b), the low-energy peak at the Fermi level is already present at the anti-nodal
point (although shifted in energy) and smoothly evolves into the sharp quasiparticle peak at the
nodal point, along with the appearance of a sub feature around -0.1eV. This can be linked with
a smooth evolution from a $D13/D42$ to a (0, 0) contribution as getting closer to the nodal point,
along with the appearance of a well separated satellite at the nodal point due to the $D13/D42$
orbitals, as can be seen in Fig. 3.12.

For the 8-site cluster the scenario is different: at the nodal point the periodization procedure
ensures that the cluster momenta orbitals (±π/2, ±π/2) are the only contributors. Towards the
nodal point, the spectral function smoothly evolves from a (0,0) to a (±π/2, ±π/2) contribution
along the path Γ−M, and from $D13/D42$ to (±π/2, ±π/2) along the Fermi arc (by varying α). Since
the (±π/2, ±π/2) orbitals are also subject to the $Q = (π, π)$ fluctuations and are more sensitive
than the (0,0) orbital (the hybridization to the bath is weaker, see Fig. 3.13), the resulting quasi-
particle peak at the nodal point is broader and the double peak are less distinct, in agreement
with the experiment. Especially along the Fermi arc, the behavior of the EDC as a function of α
is different from the plaquette cluster: the low-energy peak emerges close to the nodal point but
does not come from a smooth evolution of the PG peak at the anti-nodal point. This is linked to
the explicit existence of the (±π/2, ±π/2) orbitals which becomes the main contribution around
the nodal point, in contrast with the plaquette cluster where it is a combination of $D13/D42$ and
(0, 0) orbitals. The PG peak at the anti-nodal point seems to be smoothly absorbed into a shoulder of the quasiparticle peak at the nodal point.

## 3.5 DYNAMICAL SPIN STRUCTURE FACTOR: A DIRECT PROBE OF THE SPIN FLUCTUATIONS

The results we have shown and analyzed so far point to an understanding of the waterfall and PG physics in a unified framework of antiferromagnetic spin fluctuations. It is therefore interesting to probe directly the collective spin fluctuations: in experiment this can be done with RIXS measurements which have proven useful to extract information about the magnon excitations in cuprates [166, 248–250]. Contrary to ARPES, calculating a RIXS spectrum is extremely challenging [264], and it is easier to use indirect comparison between theory and experiments. RIXS is sensitive to the DSSF [264], which is at the origin of the low-energy peaks usually identified as (para)-magnons.

The DSSF can be determined within CDMFT via the two-time two-particle Green's function which is accessible through some impurity solvers. In this thesis, we used the CT-HYB [116] solver to obtain the two-time two-particle Green's function. After a the self-consistent CDMFT loop is converged using the CT-INT [117] solver, an extra loop is done with CT-HYB. In order to improve the sampling of the two-particle quantity, the eigen states of the isolated cluster having an energy larger than 5eV with respect to the ground state were projected out [116]. We ensured that the projection did not change significantly the results, but simply improved the statistics of the two-particle quantity sampling. Unfortunately, the CT-HYB was not able to manage the 8-site cluster for it was impossible to perform the calculation in a reasonable amount of time, and since the two-time two-particle Green's function can not be obtained with the CT-INT (not implemented yet), we are limited to the 2 × 2 plaquette.

In order to relate the two-time two-particle Green's function to the DSSF, we start from the definition of the latter in imaginary time (since the impurity solver works in Matsubara space):

$$S(\mathbf{k}, \tau) = \frac{1}{L} \sum_{\mathbf{R}_i \mathbf{R}_j} e^{i\mathbf{k}(\mathbf{R}_j - \mathbf{R}_i)} \langle \mathbf{S}_i(\tau) \mathbf{S}_j(0) \rangle, \quad (3.22)$$
where L is the number of cluster sites, $\mathbf{R}_i$ denotes the position of the site $i$, and $\mathbf{S}_i(\tau)$ the spin operator acting on site $i$ at time $\tau$. The correlation functions in the sum can be expanded in terms of fermionic creation/annihilation operators.

$$ \begin{aligned} \langle \mathbf{S}_i(\tau) \mathbf{S}_j(0) \rangle &= \langle \mathbf{S}_i^z(\tau) \mathbf{S}_j^z(0) \rangle + \frac{1}{2} [ \langle \mathbf{S}_i^\dagger(\tau) \mathbf{S}_j^-(0) \rangle + \langle \mathbf{S}_i^-(\tau) \mathbf{S}_j^+(0) \rangle ] \\ &= \frac{1}{4} \langle (n_{i\uparrow}(\tau) - n_{i\downarrow}(\tau))(n_{j\uparrow}(0) - n_{j\downarrow}(0)) \rangle \\ &\quad + \frac{1}{2} [ \langle c_{i\uparrow}^\dagger(\tau) c_{i\downarrow}(\tau) c_{j\downarrow}^\dagger(0) c_{j\uparrow}(0) \rangle + \langle c_{i\downarrow}^\dagger(\tau) c_{i\uparrow}(\tau) c_{j\uparrow}^\dagger(0) c_{j\downarrow}(0) \rangle ] \quad (3.23) \\ &= \frac{1}{4} ( G_{i\uparrow,i\uparrow,j\uparrow,j\uparrow}(\tau) + G_{i\downarrow,i\downarrow,j\downarrow,j\downarrow}(\tau) - G_{i\uparrow,i\uparrow,j\downarrow,j\downarrow}(\tau) - G_{i\downarrow,i\downarrow,j\uparrow,j\uparrow}(\tau) ) \\ &\quad + \frac{1}{2} ( G_{i\uparrow,i\downarrow,j\downarrow,j\uparrow}(\tau) + G_{i\downarrow,i\uparrow,j\uparrow,j\downarrow}(\tau) ), \end{aligned} $$

where $G_{i\sigma_i, j\sigma_j, k\sigma_k, l\sigma_l} = \langle c_{i\sigma_i}^\dagger(\tau) c_{j\sigma_j}(\tau) c_{k\sigma_k}^\dagger(0) c_{l\sigma_l}(0) \rangle$ is the two-time two-particle Green's function which can be directly obtained from the CT-HYB solver [116]. Because of the limited size of the 2 × 2 cluster, we have only access to $\langle \mathbf{S}_i(\tau) \mathbf{S}_j(0) \rangle$ up to the NNN, which is a quite drastic approximation to the Eq. 3.22. In fact, this equation corresponds to the periodization procedure applied to the momentum independent spin-spin correlation functions. For the cluster momenta $(0,0), (\pi,0), (0,\pi), (\pi,\pi)$, it amounts to a basis transformation from the site basis. To compare to the RIXS measurements, we again need to perform an analytic continuation of the DSSF. In practice, due to the noise in the data even with a large sampling number (for $\beta = 60\,\text{eV}^{-1}$, the lowest number of counts obtained for $x = 0.1$ was ~ 50 · 10⁶), it is not possible to first perform an analytic continuation of each $\langle \mathbf{S}_i(\tau) \mathbf{S}_j(0) \rangle$ individually, and then apply Eq. 3.22. Instead, the DSSF is calculated at each desired k point in imaginary time, and then we apply the MEM method with a bosonic kernel using the MaxEnt package [220].

The results are shown in Fig. 3.14(a) at a few specific k points as a function of doping, calculated at $\beta = 60\,\text{eV}^{-1}$. In agreement with our analysis above and with the literature [75, 111–113], the DSSF is maximum at $Q = (\pi, \pi)$, hence directly showing the importance of the antiferromagnetic spin fluctuations. Interestingly, the $Q = (\pi, \pi)$ sector is dominant for all the doping considered (up to $x = 0.1$), and it is stronger in the undoped system as expected. The DSSF's peak at $Q = (\pi/2, \pi/2)$ is located at lower energy than at $Q = (\pi, 0)/(0, \pi)$, consistently with the experimental observation of Ref. 166 for the parent NaCCOC compound ($x = 0$). However, the difference in energy is overestimated: it is roughly about 100meV for $x = 0$, while in experiment it is measured to be of the order of 40meV [166]. Given that the experimental resolution was 130meV, and that the $Q = (\pi/2, \pi/2)$ point is not explicitly a cluster momentum, the overall qualitative agreement is satisfying. Most importantly, the CDMFT calculations are able to account for the asymmetric line shape of the peaks which is even more visible at larger dopings as the spectrum broadens, in agreement with the RIXS measurements. As we will show later with the DSSF calculated for a spin-only model, this asymmetric line shape originates from the coupling of the spin excitations to the charge degrees of freedom.

A problem can nevertheless be identified in the doping dependence of the DSSF. Indeed, for
<figure><img src="image_34.png" /><figcaption>Figure 3.14: (a) DSSF at different **k** points as a function of doping (x = 0, x = 0.05, x = 0.1), all calculated at **β** = 60eV⁻¹. The same error estimate (10⁻²) was used for all spectra to perform the analytic continuation. (b) Local, NN, and NNN cluster spin-spin correlation functions as a function of imaginary time, for different dopings and calculated at **β** = 60eV⁻¹. Note that the NN component has been multiplied by -1.</figcaption></figure>

all the **Q** showed here but **Q** = (0, 0), we observe a softening of the DSSF peak as the doping is increased. It is roughly of the order of 50meV, whereas it was observed experimentally and theoretically with different methods that the peak position should almost not move as a function of doping, or even increase [75, 248–250]. Given that some of these results were obtained after an analytic continuation procedure, the culprit is most probably the finite size of our plaquette cluster. Indeed, the **Q** = (π, π) excitation can only connect the D13, D42 orbitals, or the (0, 0), (π, π) orbitals, and therefore lacks the important (±π/2, ±π/2) channel, which would be included in the 8-site cluster. This is of importance, since the D13, D42 orbitals are gapped, while the (±π/2, ±π/2) orbitals show a coherent quasiparticle peak at the Fermi energy (see Fig. 3.13). Also, one can notice in Fig. 3.14(b) that as the doping is increased, the spin-spin correlation functions saturate at τ = β/2 and even exceed the value obtained for the undoped compound. The NNN correlation function at τ = β/2 is even slightly larger than the NN component. This phenomenon has already been observed previously [122, 123], and shows that the longer ranged spin-spin correlation functions should be incorporated since they do not show a fast decay with distance at long times τ ~ β/2.

The non-zero value of the spin-spin correlation functions at large time $\tau \to \beta/2$ is directly linked to the existence of low-energy spectral weight in the DSSF, i.e. to low-energy excitation modes. Indeed, similarly to the Eq. 2.69 linking the fermionic imaginary time Green's function to its real frequency counterpart, we can relate the imaginary time spin-spin correlation function to its spec-
tral function on the real axis. For the sake of simplicity, we write the imaginary time correlation function as $S_{ij}(\tau) = \langle S_i(\tau)S_j(0) \rangle$, and its relation to the real frequency axis $S_{ij}(\omega)$ reads

$$ 
\begin{aligned}
S_{ij}(\tau) &= \frac{1}{2} \int_{-\infty}^{+\infty} S_{ij}(\omega) \frac{e^{-\omega(\beta-\tau)} + e^{-\omega\tau}}{1 - e^{-\omega\beta}} d\omega \\
&= \int_{0}^{+\infty} \frac{S_{ij}(\omega)}{\omega} \omega \frac{e^{-\omega(\beta-\tau)} + e^{-\omega\tau}}{1 - e^{-\omega\beta}} d\omega,
\end{aligned} 
\quad (3.24) $$

where we use the asymmetry of the function $S_{ij}(\omega)$ to obtain the second line, in which we artificially divided and multiplied by $\omega$ for convenience. At $\tau = \beta/2$, the above expression can be drastically simplified:

$$ 
\begin{aligned}
S_{ij}(\beta/2) &= \int_{0}^{+\infty} \frac{S_{ij}(\omega)}{\omega} \omega \frac{e^{-\omega(\beta/2)} + e^{-\omega\beta/2}}{1 - e^{-\omega\beta}} d\omega \\
&= \frac{1}{2} \int_{0}^{+\infty} \frac{S_{ij}(\omega)}{\omega} \frac{\omega}{\sinh(\omega\beta/2)} d\omega \\
&= \frac{2}{\beta^2} \int_{0}^{+\infty} \left( \frac{S_{ij}(2\omega/\beta)}{2\omega/\beta} \right) \frac{\omega}{\sinh(\omega)} d\omega.
\end{aligned} 
\quad (3.25) $$

Therefore, at $\tau = \beta/2$ there is a weighting function in the integral of the form $\frac{\omega}{\sinh(\omega)}$, which very quickly goes to zero as $\omega$ is increased, whereas its low frequency limit is given by:

$$ \lim_{\omega \to 0} \frac{\omega}{\sinh(\omega)} = \frac{1}{1 + \frac{\omega^2}{3!}}. \quad (3.26) $$

This means that $S_{ij}(\beta/2)$ is finite only if $\frac{S_{ij}(\omega)}{\omega}$ is non-zero at $\omega \to 0$, i.e. if there are low-energy spin excitations.

The plateau observed at long imaginary times in Fig. 3.14 for the doped system is easily explained by these considerations. In the undoped case, the system has a large gap which prevents the low-energy spin excitations, such that at all cluster momenta the DSSF is either suppressed or gapped. However, for both $x = 0.05$ and $x = 0.1$ hole concentration, at all cluster momenta the DSSF is non-zero at low energies, therefore leading to an increase of the spin-spin correlation function at $\tau = \beta/2$. At $\beta = 20$, we can notice that at $Q = (\pi, \pi)$ and $Q = (\pi/2, \pi/2)$ the DSSF is not gapped even at $x = 0$, which leads to a much stronger value at $\tau = \beta/2$ for the undoped system, resulting in the absence of the spin freezing phenomena. Note that despite this kind of changes in the DSSF, the temperature also has direct implications on the value of the spin-spin correlation function at $\tau = \beta/2$: as can be seen in Eq. 3.25, a lower temperature results in a lower $\frac{1}{\beta}$ constant, as well as a narrowing of the energy window which effectively contributes to the integral since the frequency dependence of $S_{ij}$ is scaled by a $\frac{1}{\beta}$ factor. Therefore, in general we can expect the value at $\tau = \beta/2$ to decrease with the temperature. It is indeed what is observed if we compare for instance the spin-spin correlation function calculated at $\beta = 40\text{eV}^{-1}$ and $\beta = 60\text{eV}^{-1}$, two inverse temperatures at which we observe the PG. The results are shown in Fig. 3.15: as expected, the value of the correlation function at long times decreases with decreasing temperature. In-
<figure><img src="image_35.png" /><figcaption>Figure 3.15: Local, NN, and NNN cluster spin-spin correlation functions as a function of imaginary time, for different dopings and calculated at $\beta = 60\text{eV}^{-1}$ and $\beta = 40\text{eV}^{-1}$. Note that the NN component has been multiplied by -1.</figcaption></figure>

terestingly, one can notice that for the doping $x = 0.1$, although its absolute value is lower at $\beta = 60\text{eV}^{-1}$ than at $\beta = 40\text{eV}^{-1}$, it is closer to that of $x = 0$ and $x = 0.05$ at lower temperature, hence being a sign that the freezing propagates to larger dopings with decreasing temperature, which is in agreement with the results of Ref. 123.

If we can explain the increase of the correlation function at $\tau = \beta/2$ from the undoped to the doped system, it is still not clear why it freezes to almost the same value independently of the distance between the sites i and j. It is of utmost importance to understand this point, since it is a clear sign that at long times the correlation function is distance-independent, i.e. is long-ranged. A hint can be found by reversing the Fourier transform of Eq. 3.22:

$$S_{ij}(\omega) = \frac{1}{L} \sum_{\mathbf{k}} S(\mathbf{k}, \omega). \qquad (3.27)$$

With the $\mathbf{k}$ sum restricted to the cluster momenta, we can obtain the local, NN and NNN correlation function for the 2 × 2 plaquette:

$$S_{\text{loc}}(\omega) = \frac{1}{4}(S((0,0),\omega) + S((\pi,0),\omega) + S((0,\pi),\omega) + S((\pi,\pi),\omega)) \\ S_{\text{NN}}(\omega) = \frac{1}{4}(S((0,0),\omega) - S((\pi,0),\omega) + S((0,\pi),\omega) - S((\pi,\pi),\omega)) \qquad (3.28)$$

$$S_{\text{NNN}}(\omega) = \frac{1}{4}(S((0,0),\omega) - S((\pi,0),\omega) - S((0,\pi),\omega) + S((\pi,\pi),\omega)).$$

The key point is to notice that the $Q = (\pi, \pi)$ term clearly dominates all the others at low-frequency, for small $\omega$ all the spin-spin correlation functions in absolute value $|S_{\text{loc}}|$, $|S_{\text{NN}}|$, $|S_{\text{NNN}}|$ are roughly of the same order: $\frac{1}{4}|S((\pi, \pi), \omega)|$. Combining this observation with the previous ar-
gument, it is clear that the spin-spin correlation functions at $\tau = \beta/2$ are equal because they all share the same contribution from $S((\pi, \pi), \omega)$, such that we have identified the origin of the spin freezing observed by Werner et al. [122, 123]. Upon further increasing the hole doping, the $Q = (\pi, \pi)$ becomes less dominant because the spin-spin fluctuations strongly weaken so that the freezing mechanism disappears, coherent with the temperature dependence observed in Fig. 3.15 and the findings of Werner et al. [123]. This result is even more general than the ferromagnetic fluctuations between NNN [123]: all longer ranged spin-spin fluctuations will have the same contribution from $Q = (\pi, \pi)$ such that the freezing extends to long-range correlation functions. We emphasize that such a mechanism can occur only if there are low-energy spin excitations, and it is even enhanced when there are zero-frequency spin excitations. In this regard, the freezing mechanism can be expected to be even larger in the 8-site cluster! Interestingly, it is at least confirmed in the 8-site cluster for the NN and NNN components in Ref. 123 using DCA.

One may wonder at this stage if this freezing might be an effect of the finite size of the cluster. Indeed, the $S((\pi, \pi), \omega)$ is in fact weighted by $1/L$, with $L$ the number of cluster sites, such that its contribution at low energy is suppressed for a larger number of sites, hence canceling the freezing mechanism. This is true only if the large spin excitations have no momentum broadening, i.e. if only the $Q = (\pi, \pi)$ response is strong, yielding a scaling factor $S((\pi, \pi), \omega)/L$ going to zero for large systems. However, in reality a relatively large momentum broadening is observed, as can be seen for instance in the diagrammatic Monte Carlo calculations of Wu et al. [113]. Therefore, in a certain length scale $\delta_k$ the DSSF can be considered as large as at $Q = (\pi, \pi)$. For a square cluster made of $\sqrt{L} \times \sqrt{L}$ sites, this implies that each spin-spin correlation function $S_{ij}(\omega)$ has $(\delta_k)^2 L$ contributions of the order of $S((\pi, \pi), \omega)/L$, such that the overall contribution is

$$S_{ij}(\omega) \propto (\delta_k)^2 S((\pi, \pi), \omega). \qquad (3.29)$$

If the broadening is sufficiently large, and if the $Q = (\pi, \pi)$ strongly dominates, there is a freezing of the spin-spin correlation function at large distances in the doped system. This is an intriguing result that asks for further investigations. It is important to point here that we do not claim here that this observation is linked with the spin-glass-like phase observed experimentally [84] which occurs at lower temperature ($T \sim 6K$). Calculations at lower temperature, or even zero-temperature, would be necessary to further explore a possible link.

Another sign that this spin-freezing is related to the PG can be seen in Fig. 3.16, where we show essentially the same plots as in Fig. 3.14 but calculated at $\beta = 20eV^{-1}$ where we do not observe a PG. The DSSF spectra are quite similar to those calculated at $\beta = 60eV^{-1}$, although the peaks at higher dopings are shifted to higher energies. The softening of the DSSF peak maximum seems to be slightly reduced. Most importantly, the spin-spin correlation functions show a drastically different behavior: they continuously decrease with increasing doping, without any sign of freezing, in sharp contrast to the $\beta = 60eV^{-1}$ results. This can be understood from the fact that in the undoped system the DSSF at $Q = (\pi, \pi)$ is already non-gapped thus leading to a different behavior than at lower temperature.
<figure><img src="image_36.png" /><figcaption>Figure 3.16: (a) DSSF at different **k** points as a function of doping ($x = 0, x = 0.05, x = 0.1$), all calculated at $\beta = 20\text{eV}^{-1}$. The same error estimate ($10^{-2}$) was used for all spectra to perform the analytic continuation. (b) Local, NN, and NNN cluster spin-spin correlation functions as a function of imaginary time, for different dopings and calculated at $\beta = 60\text{eV}^{-1}$. Note that the NN component has been multiplied by -1.</figcaption></figure>

The usual antiferromagnetic magnon dispersion goes to zero at $Q = (\pi, \pi)$, while in our calculations the peak rather appears around $\omega \sim 0.1\text{eV}$ even for the undoped system. Since the finite size of our cluster induces strong limitations as we have discussed above, it is legitimate to wonder if the non-zero frequency of the peak at $Q = (\pi, \pi)$ is another spurious finite size effect. In order to check that, we performed exact diagonalization (ED) calculations for a simple Heisenberg spin model with only NN exchange interactions, at half filling:

$$H = J \sum_{\langle ij \rangle} \mathbf{S}_i \cdot \mathbf{S}_j + h \sum_i (-1)^{x_i+y_i} S_i^z \quad (3.30)$$

with $J = \frac{4t^2}{U} - \frac{24t^4}{U^3}$ [235], and where $h$ is a staggered field giving the possibility to lift the degeneracy between the two antiferromagnetic configuration possible. We solve the model for a $2 \times 2$ and a $4 \times 4$ plaquette with ED using the QuSpin package [265, 266] with periodic boundary conditions. Once the system solved for its eigen states and energies, the DSSF at zero temperature can be calculated by a modified Lehmann representation:

$$S(\mathbf{k}, \omega) = \sum_{\mathbf{R}_i \mathbf{R}_j} e^{i\mathbf{k} \cdot (\mathbf{R}_j - \mathbf{R}_i)} \sum_m \frac{\langle GS | \mathbf{S}_i | m \rangle \langle m | \mathbf{S}_j | GS \rangle}{\omega + i\delta + E_{GS} - E_m}, \quad (3.31)$$

where $m$ runs over the excited states. In practice, the vector product of the two spin operators is expanded before introducing the closure relation.

The results are presented in Fig. 3.17(a) for the $2 \times 2$ and $4 \times 4$ system with and without staggered
<figure><img src="image_37.png" /></figure>

Figure 3.17: (a) DSSF of the spin-only model calculated for a 2 × 2 and 4 × 4 system with eventually a non-zero staggered field. The vertical dotted lines in the right panel highlight the position of the peak at **k** = (0, 0) (right panel) for a non-zero staggered field. (b) Momentum resolved DSSF obtained from an antiferromagnetic (left) and paramagnetic (right) CDMFT calculation at $\beta$ = 20eV⁻¹. The red line is a guide to the eye following the peak maximum at each momentum, we emphasize that it should be thought as a dispersion.

field. Without the staggered field, the degeneracy of the two possible antiferromagnetic configurations lead to a suppression of the peak at (0, 0). At $(\pi, \pi)$, we observe that the peak is shifted to lower energy as the system size is increased, in agreement with the expectations mentioned earlier. When turning on the staggered field, we observe a non-zero peak at (0, 0) since the two antiferromagnetic configurations are not equivalent anymore. As highlighted by the dotted vertical lines in the right panel of Fig. 3.17(a), the peak at (0, 0) appears at the exact same energy as the peak at $(\pi, \pi)$. Again, although a small shift of the peak energy due to the non-zero staggered field, increasing the cluster size also leads to a decrease in the energy at which the maxima appears when the antiferromagnetic configuration degeneracy is lifted.

The fact that the peaks at $(\pi, \pi)$ and (0, 0) in the spin-only model appear at the exact same energy, and that the peak line shape is perfectly symmetric in the spin-only model, are direct evidences of the role of the charge degrees of freedom. To show that, we performed an antiferromagnetic CDMFT calculation at $\beta$ = 20eV⁻¹, that we compare to the usual paramagnetic calculations, as shown in Fig. 3.17(b). As for the spin-only model, in the antiferromagnetic calculations a sizable peak at (0, 0) at finite frequency appears, in sharp contrast to the paramagnetic case which shows a suppression of the DSSF at (0, 0), in agreement with the spin-only model. However, the peak at $(\pi, \pi)$ is significantly shifted as compared to the one at (0, 0). Hence the peak position at $(\pi, \pi)$ is not only influenced by finite size effects: it is also governed by the interplay between charge and spin fluctuations, even though the latter dominate the former [112, 113].
Therefore, despite the fact that the **PG** is well described by limited short-range interactions, the **DSSF** seems to require the incorporation of longer range correlation functions, as well as the ex-citation channel between the nodal points **Q** = (±π/2, ±π/2), to be in good agreement with the experiment. If the line shape and the increasing broadening as a function of doping is consistent with the experimental and theoretical findings [75, 166, 248–250], the softening of the **DSSF** peak is not well captured within our method. A future step is to use the 8-site cluster, which we expect to significantly change the results.

3.6 CONCLUSION

By combining state-of-the-art theoretical and experimental techniques, as well as simpler effec-
tive models that capture the essential mechanisms, we show in this chapter that the spectral prop-
erties of NaCCOC from zero to intermediate doping, and by extension that of the cuprates, can be
well understood in a unified framework in which the short range antiferromagnetic fluctuations
play a key role. The waterfall feature, seen both in the undoped and the doped system can be de-
scribed as a crossover between a spin-polaron physics at low-energy [95–98], and an atomic-like
behavior with purely local correlations at higher energy. We obtain an excellent agreement be-
tween the CDMFT calculations and the ARPES measurements, which further confirms that the
kink observed in photoemission is the consequence of the dressing of electrons by the spin fluc-
tuations. Interestingly, at low-doping the waterfall feature is also well described by the SCBA,
showing that strong antiferromagnetic fluctuations are still present. The CDMFT calculations
performed on the small plaquette cluster display a kink position and spectral weight dispersion
that matches the experimental ones. This demonstrates that the waterfall feature at all dopings
can be caused by short range antiferromagnetic fluctuations. We moreover argue that this is not
a consequence of the specific shape and size of the plaquette cluster by performing calculations
for a 8-site cluster which confirm the results.

Our analysis, in agreement with previous findings [75, 104, 107, 111–114, 246], show that the short-range antiferromagnetic fluctuations are also responsible for the opening of the **PG**. By exploiting the simple state histogram of the plaquette cluster, it is possible to draw a scenario for the opening of the **PG**: the non-degenerate singlet nature of the most probable state $|N = 4, S = 0\rangle$ leads to a suppression of the zero-energy excitation channel for the *D*13, *D*42 orbitals which are solely responsible for the spectral weight at the anti-node. The first excitation channel available is the transition to $|N = 4, S = 1\rangle$ state in which the *D*13, *D*42 are in a triplet state. The **PG** is therefore linked to the destruction of the *RVB*-like singlets into triplets for the electrons with $(\pi, 0), (0, \pi)$ momentum. This scenario holds for the 8-site cluster: even though the $(\pm\pi/2, \pm\pi/2)$ are susceptible to form similar singlet states as the *D*13, *D*42 orbitals, they rather couple to the non-interacting electron bath since they show a larger hybridization function either the interactions are included or not.

Finally, to directly probe the spin fluctuations we calculate the **DSSF** for the 2 × 2 plaquette.
The doping dependence of the peak position is not in line with the experimental and theoretical
previous findings [75, 248–250], and we argue that this is probably due to a finite size effect. As
the doping is increased we observe in the **PG** phase a freezing of the spin-spin correlation func-
tions at $\tau \to \beta/2$ which prevents from a fast decay of the correlators with distance. We found
an explanation for the spin freezing phenomenon, which is induced by the dominance of the
**Q = (\pi, \pi)** sector at low energy in the **DSSF**. The mechanism can be generalized to large cluster
sizes given a certain momentum broadening of the strong **DSSF** peak around **Q = (\pi, \pi)**, such
that at long imaginary times $\tau \to \beta/2$ the spin freezing extends to long distances. We believe that
this may have important consequences on the electronic properties and that it would be interest-
ing to deepen the analysis in further studies. It will be also crucial to calculate the **DSSF** for the
8-site (and larger) cluster which contains explicitly the non-gapped $(\pm\pi/2, \pm\pi/2)$ orbitals, and
which could be a first test case with **CDMFT** for the generalization of the spin-freezing mecha-
nism.
<figure><img src="image_38.png" /></figure>
# 4 BULK PROPERTIES OF TETRAGONAL CuO

The following chapter is based on Ref. 267, which was carried out in close collaboration with Max Bramberger, Martin Grundner, Sebastian Paeckel and Ulrich Schollwöck (Arnold Sommerfeld Center of Theoretical Physics, Department of Physics, University of Munich). The figures extracted from the publication are signaled explicitly in the caption. The results shown here that were not calculated by the present author are also clearly indicated.

## 4.1 INTRODUCTION

As we have seen for NaCCOC, the two-dimensional CuO₂ layers are mainly responsible for the cuprates' low-energy physics. Sometimes, however, an emergent phenomenon appearing in the CuO₂ planes can be the consequence of distortions or disorder effects due to the surrounding ligand atoms. For instance, it was shown in Hg-based *n*-layer cuprates, in which *n* CuO₂ layers are sandwiched between charge reservoir layers, that long-range antiferromagnetic order coexists with CDW order in the different CuO₂ layers [77–79]. The CDW is attributed to the CuO₂ outer planes which are closer to the charge reservoir layers and which undergo larger distortions due to chemical substitution, while the antiferromagnetic order exists in the cleaner inner planes which are protected from the structural distortions [77, 78]. This can result in significantly different outcomes in the spectral function, as was observed in Ref. 79 where Fermi arcs and pockets were observed and attributed again to inner and outer CuO₂ planes, as well as for the Fermi surface topology [78]. There is therefore some ambiguity on the properties of the pure CuO₂ layers, which motivates for the search of an *ideal* cuprate without ligand atoms to better connect the theoretical effective models with experiments.

In this context, polymorphs of pure CuO planes were considered [268], but it turns out that bulk CuO has a low-symmetry monoclinic crystal structure [269], far from the ideal 2D planes. Fortunately, when thin films of CuO are grown on a STO substrate, they stack along the *c* axis in a tetragonal crystal structure [124–126], as shown in Fig. 4.1(a). Above a critical thickness, of approximately 10 unit cells [124, 126], the crystal returns in its monoclinic phase. In this distorted rocksalt structure the Cu-O distances for basal and apical oxygens differ by a factor 1.37 [124, 125], which amounts to well separated 2D layers. First principles studies including DFT with hybrid functionals [128–131] and DFT+U [132, 133] gave first insights into the electronic structure of tetragonal CuO (t-CuO) and were able to reproduce the experimentally observed tetragonal distortion [129], which could be traced back to Jahn-Teller orbital ordering at the Cu d⁹ ions [128, 131].

Ab initio calculations also proposed an antiferromagnetic stripe order [128, 129, 131], which is
in agreement with the experimental findings from RIXS [135]. Extrapolation from other binary
TMOs [124, 137] and estimates from first principles calculations [128, 129, 131] place the Néel
temperature around ~ 800K, which is much higher than the critical temperature of its mono-
clinic bulk phase (T<sub>N</sub> ~ 220K [127]). Such a high estimate of the Néel temperature is promis-
ing regarding superconductivity since it signals a high exchange coupling J, which was recently
claimed to be correlated with the superconducting order parameter in cuprates [223, 270–272].
However, hole- or electron-doping of t-CuO has not been experimentally achieved yet, and some
questions remain open for the undoped compound.

t-CuO is an insulator with quite sizeable gap, $\Delta > 2.35$ eV, of which the electronic structure was measured via ARPES [134] and used to construct effective three- and one-band $t-J$ models [134, 273, 274]. Although t-CuO does not escape from the long-standing debate [40, 41] on whether or not ZRS [38] can describe the low-energy spectral features of cuprates, the one-band model derived from RIXS in Ref. 135 is in qualitative agreement with the ZRS one [274]. We do not aim at contributing to this discussion, but we rather study the signature of sublattice decoupling which was observed in ARPES and RIXS [134, 135], and illustrated in Fig. 4.1(b). Strong replica features were found outside of the single sublattice BZ, corroborating with the fact that the t-CuO RIXS spectrum is reminiscent of $Sr_2CuCl_2O_2$ (SCOC), a standard cuprate made of CuO$_2$ layers. These observations suggest that each of the CuO planes constituting t-CuO can be seen as two weakly coupled CuO$_2$ sublattices, but the question about the physical origin remains open.

From 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 NNN intra-sublattice hopping $t$, a surprisingly high value considering the strong replica features observed in ARPES and the similitude of the magnon spectrum with conventional cuprates. Such model is constructed from a Heisenberg-like Hamiltonian [135], and therefore lacks the dynamical effects from electron-electron interactions. In this chapter, we investigate the dynamical influence of the inter-sublattice hopping $t_d$ and the electronic correlations by the means of CDMFT [49–52], which has been successful in studying the Hubbard model and especially cuprates [50, 103–106]. Our key finding is that the inter-sublattice correlations are heavily suppressed as compared to local and short-range intra-sublattice correlations, which formally justifies to regard t-CuO as weakly-coupled interlaced CuO$_2$ lattices. While $t$ hopping is renormalized by the intra-sublattice correlations, $t_d$ is let bare since the inter-sublattice correlations are vanishing, such that the effective $t_d/t$ is decreased by a factor $\approx 2$. With the sublattices being weakly connected to each other, we can motivate an efficient block-construction scheme for our cluster calculations, allowing to increase for free the momentum resolution of the calculated spectral functions. In Sec. 4.4, using a MPS [56, 57]-based impurity solver working directly on the real axis [58, 141–144] and at effectively zero temperature, the equal energy maps and momentum resolved spectral functions are reproduced, in remarkable agreement with ARPES measurements without the need for analytic continuation. These calculations, performed by M. Bramberger and M. Grundner, are also in excellent agreement with results obtained at finite temperature by the author using the hybridization-expansion CTQMC [60, 62] solver CT-HYB [116] based on the ALPSCore libraries [118]. Then in Sec. 4.5.1, the magnetic
<figure><img src="image_39.png" /><figcaption>Figure 4.1: *Figure taken from Ref. 267.* (a) Rock salt crystal structure of tetragonal CuO. (b) slab of CuO within the a-b plane. Bright (dark) red atoms indicate the sublattice A and B of our model. (c) Two identical Cu-sublattices and indication of the hoppings *t*<sub>d</sub>, *t*, *t*′ and *t*″ included in the model. The arrows sketch the stripe order considered throughout the chapter. Highlighted in blue and green are the magnetic sublattices that correspond to the stripe order. (d-g) Clusters including different hopping terms as discussed in the text.</figcaption></figure>

ordering in **t-CuO** is analyzed as a function of temperature and two driving mechanisms for the
insulating phase are identified. We also find an increase of the local correlations close to the crit-
ical temperature *T*<sub>c</sub> at which the clusters antiferromagnetically order, and show with a minimal
model derived by the author that it is a pure effect of thermal and spin fluctuations. Finally in
Sec. 4.6, despite that **t-CuO** has not been experimentally doped, i.e. that superconductivity has
not been observed yet, the emergence and symmetry of the superconducting order parameter
upon hole-doping is studied using **VCA** [68], a complementary cluster technique, for which the
calculations were carried out by B. Lenz. The sublattice decoupling has direct consequences even
in the hole-doped regime: a coexistence of magnetic stripe order and superconductivity of *d*<sub>xy</sub>
symmetry is found, whereas the usual cuprate *d*<sub>x<sup>2</sup>-y<sup>2</sup></sub> symmetry is strongly suppressed.

4.2 MODEL AND METHOD

We start by presenting the band structure of *t*-CuO, from which we can motivate the construction
of an effective one-band Hubbard model, similarly to NaCCOC. We then detail how the model is
solved using CDMFT and VCA calculations.
4.2.1 ONE-BAND EFFECTIVE MODEL

To construct a minimal model for t-CuO, it is instructive to analyze its band structure obtained from a simple DFT calculation done within LDA using the Wien2k package [173, 174], shown in Fig. 4.2. As expected, the system is predicted to be metallic since the Cu atoms are in a d⁹ configuration (and no symmetry is broken). The d<sub>xy</sub> orbital is half-filled, in contrast to usual cuprates like NaCCOC for which it would be the d<sub>x²-y²</sub>. In t-CuO, the d<sub>x²-y²</sub> orbital is favored against d<sub>xy</sub> by the crystal field due to the O atoms' position. Regardless of these differences, there are strong similarities between t-CuO and the previously studied NaCCOC:

(i) we find a single band crossing the Fermi energy (actually two here because there are two Cu atoms in the unit cell),

(ii) this band is only of Cu d<sub>xy</sub> (equivalently d<sub>x²-y²</sub> for NaCCOC) and O p<sub>x/y</sub> character,

(iii) it has a weak dispersion along k<sub>z</sub>.

We observe another set of bands close to X and (π/2, π/2), which are reminiscent of the spectral weight seen in ARPES and labeled β by Moser et al. [134]. They are of Cu d and O p<sub>z</sub> character, and they have no contribution from the d<sub>xy</sub> nor the p<sub>x/y</sub> orbitals.

Close to the Fermi level, t-CuO's electronic structure is composed of a quasi 2D single half-filled band of mixed Cu d<sub>xy</sub>/O p<sub>x/y</sub> character, which we identify to the single band observed in ARPES. It is moreover the only partially filled band, i.e. the one contributing mostly to the magnon part of the RIXS spectrum. One can notice signs of the sublattice decoupling at the DFT level: the single band dispersion between M and Γ is almost symmetric with respect to the red zone boundary which represents the single sublattice BZ, and the band energy at those two points only differs slightly. However, this observation can not be put on equal footing to the experimental observations, since ARPES is sensitive to spectral properties and RIXS to collective fluctuations, whereas DFT provides the ground state properties. Therefore, an explicit treatment of the electronic correlations is needed: both to obtain the expected insulating behavior, as well as to incorporate the collective dynamic correlations.

To this aim, we consider a single CuO layer as shown in Fig. 4.1(a)-(b), and model it by a single-band Hubbard model [35–37]:

$$
H = U \sum_i n_{i\uparrow} n_{i\downarrow} + \sum_{\substack{i,j,\sigma \\ |i-j|=a}} t_d c_{i\sigma}^{\dagger} c_{j\sigma} + \sum_{\substack{i,j,\sigma \\ |i-j|=\sqrt{2}a}} t c_{i\sigma}^{\dagger} c_{j\sigma} + \sum_{\substack{i,j,\sigma \\ |i-j|=2a}} t' c_{i\sigma}^{\dagger} c_{j\sigma} + \sum_{\substack{i,j,\sigma \\ |i-j|=2\sqrt{2}a}} t'' c_{i\sigma}^{\dagger} c_{j\sigma}. \quad (4.1)
$$

Since we are not interested in the fine tuning of the model, we chose to take the hopping param-
eters t<sub>d</sub> = -0.1eV, t = 0.44eV, t' = -0.2eV, t'' = 0.075eV (see Fig. 4.1(c)) from the fit of the
RIXS magnon spectrum of Ref. 135. As mentioned above, the O atoms favor the NNN hopping
t rather than the NN t<sub>d</sub>, in contrast to the usual CuO<sub>2</sub> layers. We set the on-site Hubbard inter-
action U = 7eV, significantly higher than the one obtained from the magnon fitting of Ref. 135,
but necessary for obtaining a gap that is larger than the ARPES lower bound [134], but still lower
than the STM upper bound [126], and consistent with hybrid DFT [129] (2.7eV).
<figure><img src="image_40.png" /><figcaption>Figure 4.2: t-CuO band structure on a k-point illustrated on the BZ (**left**) obtained with a non spin polarized DFT calculation, along with the projected DOS. The red shaded area highlights the BZ of a single CuO₂ sublattice.</figcaption></figure>

### 4.2.2 METHODOLOGY

To solve this effective interacting many-body problem, we use CDMFT [50–52], which was introduced in Sec. 2.2. We perform antiferromagnetic calculations and chose the cluster tilings such that they are in line with the antiferromagnetic stripes. The calculations are initialized using a spin-polarized (constant) self-energy as presented in Sec. 2.4, and the chemical potential is kept fixed in the middle of the gap. In Sec. 4.3, we first study the sublattice decoupling using the *t*-dimer and *t* − *t*<sub>d</sub> clusters of Fig. 4.1(d)-(e). The *t*-dimer is composed of two sites belonging to the same sublattice, i.e. connected via the NNN hopping *t*. The *t* − *t*<sub>d</sub> cluster is composed of two of those dimers explicitly connected through the *t*<sub>d</sub> hopping, thus allowing to directly study the effect of including *t*<sub>d</sub> on the correlations. Following our results of Sec. 4.3, we then motivate the use of the block construction (Fig. 4.1(g)) in Sec. 4.4 to compute the spectral function of *t*-CuO with an increased momentum accuracy. Similarly to the *t* − *t*<sub>d</sub> cluster, the super-cluster is composed of two 2 × 2 plaquette clusters on each sublattice, but this time the inter-sublattice hopping *t*<sub>d</sub> is treated at the non-interacting level and is not included in the effective impurity. In other words, based on our observations in Sec. 4.3, the inter-sublattice correlations are neglected such that we solve the CDMFT impurity problem for only one of the two plaquette clusters. All the details of this construction are given in Sec. 4.4.

Having defined our impurity clusters, we are left with choosing the solver. As mentioned in the introduction of this chapter, three different complementary techniques were used: two MPS-based solvers working both on the imaginary [59, 138–140] and real [58, 141–144] frequency axis for which M. Bramberger and M. Grundner performed the calculations, and the hybridization-expansion CTQMC [60–62] solver CT-HYB [116] based on the ALPSCore libraries [118] for which the author performed the calculations. The combination of these techniques allows us to combine zero and finite temperature analysis on the imaginary axis, as well as to obtain a direct information on the real axis (without the need for analytic continuation) that can be directly compared to experiments.
Finally, in Sec. 4.6, VCA [68–70] calculations, performed by B. Lenz, are used to investigate the possible superconducting solutions. This variational quantum cluster technique is well-suited to check for different symmetry-breaking orders of the interacting lattice system [145, 275]. Here, the self-energy functional which is to be minimized has several variational parameters. It includes the chemical potential of the cluster [276], and the one of the lattice [277], which are first determined. Then, the different symmetry-breaking Weiss fields are sequentially added with an adiabatic switching to allow for long-range order [145, 275, 276]. The antiferromagnetic stripe order with ordering wavevector $Q = (0, \pi)$ is triggered by a suitable field on the cluster:

$$ \mathcal{H}_{\text{AFS}}^{\text{Weiss}} = M \sum_{\mathbf{R}} (-1)^{R_y \cdot \pi} (n_{\mathbf{R}\uparrow} - n_{\mathbf{R}\downarrow}), $$

where we denote the cluster sites with $\mathbf{R}$, $M$ is the field strength determined via the variational principle. Likewise, superconducting pairing fields are added via

$$ \mathcal{H}_{\text{SC}}^{\text{Weiss}} = D \sum_{i,j} \Delta_{i,j} (c_{i\uparrow} c_{j\downarrow} + \text{h.c.}), $$

where $D$ denotes the variational parameter and $\Delta_{i,j}$ is chosen such that it amounts to pairing with $d_{xy}$ or $d_{x^2-y^2}$ symmetry. The cluster self-energy is calculated using exact diagonalization with a band Lanczos algorithm [278]. Moreover, a Nambu transformation is employed to include the superconductivity fields [145].

## 4.3 SUBLATTICE DECOUPLING

At the single particle level it is hard to argue for the decoupling of the two sublattices since the nearest-neighbor hopping $t_d$ is of the same order of magnitude as the next-nearest neighbor hopping ($t_d \sim -\frac{t}{4}$). Therefore it is important to also take into account the self-energy which captures the modification of the non-interacting Hamiltonian due to the presence of electronic interactions in the correlated material.

As we have seen in Chap. 2, within the framework of CDMFT the local interactions, hopping terms on the given cluster and dynamical fluctuations to an electronic reservoir are taken into account exactly, while longer-ranged exchange with the rest of the lattice is included on the single-particle level and enters via the self-consistency loop [279]. The cluster self-energy $\Sigma(\omega)$ is a matrix-valued quantity in terms of combined cluster-spin indices. It links the non-interacting and interacting cluster Green's functions, $G_0(\mathbf{k}, \omega)$ and $G(\mathbf{k}, \omega)$, via the Dyson equation

$$ \Sigma(\omega) = G_0(\mathbf{k}, \omega)^{-1} - G(\mathbf{k}, \omega)^{-1}. $$

Besides the local component, $\Sigma_{\text{loc}}(\omega)$, non-local self-energies within the cluster are accessible, which we denote with respect to the hopping term connecting the corresponding sites, e.g. $\Sigma_t(\omega)$, $\Sigma_{t_d}(\omega)$. By choosing different impurity cluster geometries, namely the t-dimer (Fig. 4.1 (d)) and the t-t_d (Fig. 4.1 (e)) clusters, we are able to probe the influence of the inter-sublattice com-
<figure><img src="image_41.png" /></figure>

<figure><img src="image_42.png" /></figure>

<figure><img src="image_43.png" /></figure>

<figure><img src="image_44.png" /></figure>

<figure><img src="image_45.png" /></figure>

<figure><img src="image_46.png" /></figure>

Figure 4.3: Panels (d)-(f) were extracted from Ref. 267. Comparison between selected elements of the self-energy computed on two different clusters using (a)-(c) the CT-HYB solver at $\beta = 50\text{eV}^{-1}$, and (d)-(f) the MPS-based solver on the imaginary axis. Note the difference in scales between the panels (c),(f) and (a)-(b),(d)-(e). The components shown belong to the up-spin block of the self-energy.
ponent $\Sigma_{t_d}(\omega)$ on the other components $\Sigma_{\text{loc}}(\omega)$ and $\Sigma_t(\omega)$, i.e. on the coupling between the two sublattices. The following results have been obtained by a MPS-based impurity solver [59, 138, 140] working on the imaginary axis at effectively zero temperature ($T = 0$ K), and by the CT-HYB [116, 118] solver at low temperature $\beta = 50\text{eV}^{-1}$.

We show in Fig. 4.3(a)-(b),(d)-(e) the elements of the self-energy already included in the dimer cluster both calculated with the MPS-based and the CT-HYB solver, which do essentially not change by incorporating explicitly the $t_d$ hopping in the impurity cluster. Both the real and imaginary part of $\Sigma_{\text{loc}}(\omega)$ and $\Sigma_t(\omega)$ computed with the dimer are equal to the ones computed with the $t-t_d$ cluster. This is related to the inter-sublattice component $\Sigma_{t_d}(\omega)$ being strongly suppressed, as shown in Fig. 4.3(c),(f). It is two to three orders of magnitude smaller than the intra-sublattice element $\Sigma_t(\omega)$ (Fig. 4.3(b),(e)) depending on which solver is used. Since the inter-sublattice hopping $t_d$ is roughly about one fourth of the leading order hopping ($|t_d| \approx |t_4|$), this strong suppression is far from trivial and indicates that electronic correlations effects favor intra-sublattice hopping $t$. Indeed, from a hand-wavy view of the self-energy as a static correction to the Hamiltonian, $\text{Re}[\Sigma_t(\omega = 0)] \approx -0.5\text{eV}$ would roughly double $t$, while $\text{Re}[\Sigma_{t_d}(\omega = 0)] \approx 0$ lets $t_d$ bare, so that the $t_d/t$ ratio is effectively decreased by the correlations. Note that the CT-HYB $\Sigma_{t_d}(\omega)$ is noisy since its value is close to the sampling error, estimated to be roughly about $0.1 \cdot 10^{-3}$. Apart from this statistical noise in $\Sigma_{t_d}(\omega)$, the result of the self-energy components are in excellent agreement with the results at zero temperature and confirms that the sublattice decoupling survives at finite temperature.

The driving mechanism behind the formation of sublattices is linked to the fact that the hopping elements $t_d, t, t', t''$ are not decreasing monotonically with distance. Independent of whether it is the NN, NNN or any longer ranged hopping, it is the leading order one that is favored by the electronic correlations. Indeed, from a simple bonding viewpoint, it is energetically favorable for the electrons to form bonding states between sites connected via the larger hopping term. Moreover, if the NN hopping would be the largest, then all the sites of the lattice would be connected via processes including only the favored hopping (see Fig. 4.4(a)), leading to a smoothly decaying self-energy with distance. However, in t-CuO the NNN is the leading order hopping due to the position of the oxygen atoms which are located between NNN Cu sites and hence favor $t$ against $t_d$. Since the $t$ hopping term connects sites belonging to the same sublattice, as shown in Fig. 4.4(b), the NN sites are not connected by higher-order hopping processes $t^n$, leading to a strong suppression of the inter-sublattice self-energy $\Sigma_{t_d}(\omega)$.

These results are of importance, since they prove that thinking of t-CuO as two weakly coupled sublattices is correct, and it reveals the physical origin of this behavior. Since hallmarks of the decoupling had been seen in the ARPES spectra of Ref. 134, we calculate in the following the momentum-resolved spectral function of t-CuO. As we have pointed out in Sec. 2.2, CDMFT provides a Green's function of limited momentum resolution due to the finite cluster size, and the only way to improve it is to consider a larger cluster. In order to enlarge our cluster while keeping the impurity problem numerically tractable, we show in the following that the sublattice
<figure><img src="image_47.png" /><figcaption>Figure 4.4: Sketch of higher order hopping processes. In (a), the leading order hopping can connect to every site, thus leading to a smoothly decaying self-energy, contrary to (b) where t is the largest hopping, hence leading to a suppression of $\Sigma_{t_d}$ since it only connects sites of the same sublattice.</figcaption></figure>

decoupling can be used for constructing an 8-site super-cluster (Fig. 4.1(g)) that requires solving a 4-site impurity problem only. This technique, that we call hereafter the *block construction*, allows us to increase *for free* the momentum resolution of our calculations.

## 4.4 SPECTRAL FUNCTION

In this section, we first describe the block-construction framework, and then present the resulting spectral functions.

### 4.4.1 THE BLOCK CONSTRUCTION

An illustration of the block construction principle is given in Fig. 4.5. As mentioned in the previous section, our goal is to solve the CDMFT equations for the 8-site cluster of Fig. 4.5(a) for its self-energy $\Sigma_{8site}$, which would be computationally very expensive. Grouping the indices that belong to the same sublattices, $\Sigma_{8site}$ can be written as a block matrix:

$$ \Sigma_{8site} = \begin{pmatrix} \Sigma_A & \Sigma_{AB} \\ \Sigma_{AB} & \Sigma_B \end{pmatrix}, \qquad (4.2) $$

where the diagonal blocks $\Sigma_A$ ($\Sigma_B$) correspond to the self-energy of the 4-site plaquette on the sublattice A (B), and the off-diagonal $\Sigma_{AB}$ is the inter-sublattice self-energy. We have seen that the inter-sublattice self-energy between NN sites is suppressed, so that we can expect any longer ranged inter-sublattice self-energy to also vanish, i.e. $\Sigma_{AB} \sim 0$. The latter will therefore be neglected, and $\Sigma_{8site}$ is approximated by the block construction self-energy:

$$ \Sigma_{8site} \approx \Sigma_{\text{block-constr.}} = \begin{pmatrix} \Sigma_A & 0 \\ 0 & \Sigma_B \end{pmatrix}. \qquad (4.3) $$
<figure><img src="image_48.png" /><figcaption>Figure 4.5: Sketch of the block construction: we start from the full 8-site cluster in (a), made of two 4-site plaquettes on each sublattice A and B. Based on the observations of Sec. 4.3, the inter-sublattice self-energy $\Sigma_{AB} \approx 0.0$ is neglected and accordingly the cluster impurity problem in (b) is solved only for one of the two 4-site plaquettes.</figcaption></figure>

We are left with only solving the impurity problem for one of the two sublattices so to obtain $\Sigma_{A/B}$, and use Eq. 4.3 to get the 8-site cluster impurity self-energy. After closing the self-consistency loop, we project down onto one of those blocks and obtain the hybridization function for one of the two sublattices only.

Within this approximation, we treat the two largest hopping terms $t = 0.440\text{eV}$ and $t' = -0.2\text{eV}$ explicitly in the correlated effective impurity, while the two lowest ones, $t_d = 0.1\text{eV}$ and $t'' = 0.075\text{eV}$, are treated at the non-interacting level. Similar block construction techniques have been already used successfully within CDMFT [280, 281].

### 4.4.2 RESULTS

Having introduced the block construction scheme, we now compare the calculated spectral functions using the block-construction to the ARPES data of Ref. 134. They were obtained using the MPS-based impurity solver on the real axis [58, 141–144] at zero temperature and the CT-HYB [116, 118] solver at $\beta = 40\text{eV}^{-1}$ for which an additional analytic continuation was performed on the cluster self-energy using the Poor’s Man method [215, 216].

In Fig. 4.6(a),(c) are shown equal energy cuts on the top of the valence band ($E = -2.2\text{eV}$) which are in good agreement with the experimental energy map (Ref. 134, Fig. 1(a)): the strong maxima are recovered in the middle of the BZ (which are rotated by 90° with respect to the experimental data). The replica features outside the single-sublattice BZ (dashed black line), which experimentally justified the assumption of sublattice decoupling, are also reproduced by our calculations. It is interesting to emphasize here that in the case of entirely decoupled sublattices ($t_d = 0$), the spectral features inside the first BZ of a single sublattice would be exactly replicated outside this
<figure><img src="image_49.png" /><figcaption>Figure 4.6: Panels (c)-(d) were extracted from Ref. 267 Spectral function $A(\mathbf{k}, \omega)$ obtained with the (a)-(b) CT-HYB and the (c)-(d) MPS-based solver. (a),(c) Equal energy maps at $E = -2.2$ eV where the dashed black line depicts the first BZ of a single sublattice. (b),(d) $A(\mathbf{k}, \omega)$ along high-symmetry $\mathbf{k}$-path as computed with the block-construction scheme and compared to the experimentally measured dispersion (purple circles in inset) extracted from Ref. 134 and shifted by 0.4 eV in order to align the chemical potentials. All heat maps are normalized to the maximal value displayed and averaged over the possible orientations in the block-construction (see Fig. 4.7).</figcaption></figure>
BZ (i.e. outside the dashed black line in Fig. 4.6(a)-(c)). Here, since $t_d \neq 0$ the outside replica features are slightly different from the ones inside. However, we find a close resemblance since the self-energy decreases the $t_d/t$ ratio by a factor ~2 such that $t$ is effectively 10 times stronger than $t_d$. Unlike in ARPES there are no matrix-element effects [258] in our calculations, i.e. no additional intensity modulation, contrary to the ones seen in the experimental energy map of Ref. 134 across the different BZ. Still, one can notice an asymmetry between x and y axis, which is not due to the asymmetry of the cluster itself since we average over the cluster orientations (Fig. 4.7), but rather to the antiferromagnetic stripes' orientation (Fig. 4.8). We come back in more details on this point later in this section. As in the previous section, we find a striking resemblance of the CT-HYB and the MPS-based solver, even though an additional analytic continuation step is needed for the former.

The agreement between the two solvers is confirmed in the momentum resolved spectral function of the valence band using the block construction, as shown in Fig. 4.6(b),(d). Comparing our results to ARPES (cf. Fig. 2(a) in Ref. 134) we find an overall good agreement. In particular, we observe a Zhang-Rice-like band separated from the lower Hubbard band at higher binding energy, which coincides with the experimental low-energy spectrum, as shown in the inset of Fig. 4.6(b),(d). A clear separation is seen for the MPS solver which works directly on the real axis, but is blurred out by the analytic continuation for the finite temperature results. As we shall justify in the next section with finite temperature calculations (Fig. 4.10), we identify this band to stem from a spin-polaron, i.e. a hole propagating in an antiferromagnetic background. We find that around the *M* and *Γ* high-symmetry points, there is an incoherent and dispersive spectrum without well-defined structures captured by both solvers, consistently with ARPES. Moreover we reproduce the experimentally observed missing spectral weight at the *X* point, a feature which was not obtained within a SCBA calculation based on a ZRS spin-model [38, 274]. Hence, the overall agreement between our model and the experiment is striking.

An obvious feature that the calculations presented in this work can not reproduce are the contributions from a lower lying band marked with $\beta$ in the experimental data [134], which is not included in our low-energy one-band model as we have seen in Sec. 4.2. Although in the undoped system it is justified to neglect the $\beta$ band since they do not directly participate to the magnetic order, nor the sublattice decoupling, one may wonder if this is still a good approximation when hole-doping the material. It is argued by Moser. et al. that the $\beta$ band originates from a Cu *d* orbital manifold of different symmetry than the low-energy Zhang-Rice-like band of $d_{xy}$ character, which is corroborated by the LDA band structure as discussed previously. Upon hole-doping in cuprates, it is experimentally observed that the hole is located on the O sites [282–284], and from a theoretical standpoint the doped holes are expected to form ZRS states between the Cu $d_{xy}$ and O $p_{x/y}$ orbitals. Therefore, it is unlikely that these holes will populate manifolds of different symmetry such as the one at the origin of the $\beta$ band, so that the one-band model should remain a good approximation upon hole-doping.

As mentioned above, the equal energy maps and momentum-resolved spectral functions shown
<figure><img src="image_50.png" /></figure>

Figure 4.7: Figure taken from Ref. 267, the results were obtained with the MPS-based solver on the real frequency axis. Sketch of the two possible orientations within a given magnetic stripe order for (a) the dimer and (b) the block construction. In panels (c) and (d) are shown the corresponding equal energy maps obtained at $E = -2.2$eV using these cluster orientations as well as their average. The dashed black line indicates the BZ of a single sublattice.
<figure><img src="image_51.png" /><figcaption>Figure 4.8: *Figure taken from Ref. 267, the results were obtained with the MPS-based solver on the real frequency axis. Comparison of equal energy maps and momentum resolved spectral functions for the two different directions of the stripe order: (a) along x, (b) along y.*</figcaption></figure>

in Fig. 4.6 were obtained by averaging over the cluster orientations. Similarly, all the results shown in this chapter which were obtained using a dimer cluster were always averaged over the different cluster orientations, following the Oriented-Cluster DMFT scheme [285]. Averaging over the two different orientations is an important step, as illustrated in Fig. 4.7. The bias comes from the fact that hopping terms along (1,1) and (1,-1) (x and y) for the dimer (block construction) are no longer equivalent when solving the impurity problem since they are not treated equally in the cluster. The most dramatic effect is observed for the dimer cluster as is shown in Fig. 4.7(c), but the bias is removed by considering an averaging between (1,1) and (1,-1) oriented dimers. In comparison, the block construction suffers less from the bias, especially since in the block construction scheme, only one of the two the 4-site plaquettes, which are perfectly symmetric (if the coupling to the other sublattice is discarded), is solved in the impurity model, such that the asymmetry is incorporated at the non-interacting level only (though having some consequences on the self-energy via the self-consistency).

The remaining asymmetry originates from the antiferromagnetic stripes, as can be seen in Fig. 4.8. By changing the stripe orientation, we can perfectly switch the role of x and y. The most visible consequence is on the energy maps, while in the momentum resolved spectral functions only slight differences can be noticed at the *A* and *A'* high-symmetry points (in terms of spectral weight and energy position), which would probably be washed out by the experimental resolution. Unfortunately, even for the equal energy maps the extra modulation from matrix elements and the high probability that the incoming light spot covers zones with different stripes orientations prevent from the observation of this asymmetry in the experimental data.
## 4.5 TEMPERATURE DEPENDENCE

In this section, we first present the finite temperature results obtained with the CT-HYB [116, 118] solver, and show that two different insulating regimes can be identified: paramagnetic and antiferromagnetically ordered. At the magnetic ordering temperature, we observe an intriguing enhancement of the real part of the self-energy at $\omega_n \to 0$, that is explained with a simple effective model in the second part.

### 4.5.1 RESULTS

All the results we have presented so far were computed at $T = 0$K. However, there have been multiple predictions about the Néel temperature $T_N$ for the antiferromagnetic ordering of t-CuO in the literature [124, 128, 129, 137]. From an heuristic extrapolation, Siemons et al. predicted that the tetragonal phase of CuO would exhibit a high $T_N \approx 800$K (see Fig. 1(a) of Ref. 124). In CuO-Cu$_2$O heterostructures, an antiferromagnetic alignment was observed at $T_N = 600$K by Rabinovich et al. and attributed to the formation of cubic crystals of CuO. From ab initio, Peralta et al. predicted $T_N = 300-410$K with the self-interaction-corrected local density-function method [286], while Chen et al. estimated the ordering temperature to be $T_N = 865-1000$K with hybrid DFT [165] based on the Heyd-Scuseria-Ernzerhof method [287]. These theoretical predictions were obtained by estimating the exchange couplings, and evaluating $T_N$ within molecular-field theory [288, 289], to which an Anderson rescaling [17] is applied. Although these estimates cover a wide range of temperatures, they are all well above the monoclinic transition temperature $T_N \approx 200$K. This motivates for further investigation of the finite-temperature behavior of t-CuO. Our simple 2D model prevent from inferring on the value of the Néel temperature itself, however using CTQMC solvers working at finite temperature within antiferromagnetic CDMFT calculations, we can study the properties of the system at temperatures above and below the ordering temperature of the clusters $T_c$.

In Fig. 4.9(a) we show the staggered magnetization as a function of the inverse temperature $\beta$ as obtained with the dimer and the block construction clusters (Fig. 4.1(d),(g)). First, we note the asymptotic convergence of the staggered magnetization towards the $T = 0$ K value obtained with the MPS-based solver for $\beta \to \infty$. We identify an inverse temperature at which the order melts, namely $\beta_c \approx 18.5 \text{ eV}^{-1}$ (14.5 eV$^{-1}$), corresponding to a critical temperature of $T_c \approx 627$ K (800 K) obtained with the block-construction (dimer) cluster. $T_c$ is determined using a fit function of the type:

$$M(T) = \theta(T_c - T)\gamma\left(1 - \frac{T}{T_c}\right)^{\alpha},$$

where $T_c$, $\gamma$, $\alpha$ are fitting parameters and $\theta(T)$ is the Heaviside step function. By inspection of the self-energies close to the transition region, lower and upper critical temperatures are determined: the lower boundary is set by the temperature at which the spin splitting of the real-part of the self-energy vanishes, while the upper boundary is set such that the imaginary part of the diagonal components of the self-energy tends to 0 as $\omega_n \to 0$. By varying the fit boundaries within these two limits, we obtain a collection of fit from which the average critical value is extracted, as well
<figure><img src="image_52.png" /><figcaption>Figure 4.9: *Figure taken from Ref. 267.* (a) Staggered magnetization calculated using the dimer cluster and the block-construction. The dashed black lines indicate the $\beta = \infty$ result computed with the MPS-based impurity solver on the imaginary axis. The vertical lines depict the inverse critical temperature $\beta_c = (18.5\pm0.7) \text{ eV}^{-1}$ ((14.5\pm0.8) eV$^{-1}$) for the block-construction (dimer) cluster. The shaded area depicts the error bar for $\beta_c$. (b) Real part of the diagonal components of the self-energy for different inverse temperatures $\beta$ indicated in (a). The curves shown on the left correspond to the spin up (solid) and down (dashed) components on a cluster site. On the right, we show the self-energy at the two cluster momenta $K_1 = (0,0)$ (solid) and $K_2 = (0, \frac{\pi}{a})$ (dashed) respectively.</figcaption></figure>
as the error bars in Fig. 4.9(a) from the standard deviation of $\beta_c$ in the fit collection.

While the dimer cluster overestimates the magnetic order, the block-construction, which includes slightly longer-ranged magnetic fluctuations, leads to a smaller value of $T_c$. As mentioned above, here we study a simplified 2D model which does not include the inter-layer magnetic exchange coupling, and therefore for which the long-range antiferromagnetic order should not be stable at finite temperature [290]. The non-zero staggered magnetization of our CDMFT calculations is rather a consequence of choosing one of the two possible stripe directions within the mean-field scheme, than an actual hallmark of the long-range magnetic order. Despite prohibiting a direct determination of $T_N$, the reduction of $T_c$ upon extending the cluster size and its relatively high value nevertheless show the importance of including in-plane spin fluctuations.

At all temperatures, although there can be different magnetic orders, the system remains insulating. In order to understand the different mechanisms at play as a function of temperature, it is instructive to inspect the real part of the Matsubara self-energies, as shown in Fig. 4.9(b). For the sake of clarity, we discuss here the self-energies of the dimer cluster, but none of the conclusions change in the block-construction cluster. We compare the self-energies at three characteristic temperatures: in the paramagnetic regime ($\beta$ = 10eV$^{-1}$), within the transition region ($\beta$ = 20eV$^{-1}$), and in the ordered regime ($\beta$ = 40eV$^{-1}$). First, as the system enters the insulating ordered phase, we observe the asymptotic convergence towards the MPS results (see Fig. 4.9(b)). The frequency dependence of the self-energy gets strongly suppressed. This is well described in the atomic limit as derived in Ref. 291, or by the asymptotic development of the self-energy which becomes static in the antiferromagnetic ordered limit [292].

In the paramagnetic phase the insulating behavior can not be attributed to a freezing of dynamics due to large spin polarization, but rather to a momentum-selective level splitting (right panel of Fig. 4.9(b)). Close to the real axis, the $K_1 = (0, 0)$ orbital is very strongly favored with respect to the $K_2 = (0, \frac{\pi}{a})$ orbital. This is consistent with previous quantum cluster calculations performed for the dimer and larger cluster, and can be interpreted as a freezing of electron movement that is not generated by spin polarization but rather by penalizing electrons with non-zero momentum [48, 103, 108].

These different mechanisms have a drastic effect on the spectral function depending on the temperature. In Fig. 4.10 we show the momentum resolved spectral function obtained with the dimer cluster obtained at the three characteristic inverse temperature $\beta$ = 10, 20, 40eV$^{-1}$. These spectra were obtained by performing an analytic continuation with the Poor man's method from TRIQS [120, 215, 216]. We find good agreement with the experimental data only when the system is antiferromagnetic, thus showing the importance to take into account spin fluctuations in the model. At high temperature ($\beta$ = 10eV$^{-1}$), the spectral function is divided into a LHB and UHB, of which we show only the lower in Fig. 4.9, while at low-temperature ($\beta$ = 40eV$^{-1}$) the LHB shows a secondary splitting. It is the low-energy part of the LHB that agrees well with the experimental data. Moreover, we can reproduce with very good agreement this part of the spectrum
<figure><img src="image_53.png" /><figcaption>Figure 4.10: Momentum resolved spectral function obtained with the block construction. The k-path is the same as the one used in Fig. 4.6. The experimental dispersion was shifted by 0.4eV consistently with Fig. 4.6(b). The SCBA dispersion has also been shifted by $U/2$ to match the CDMFT one, and was extracted by taking the largest value in the low energy part of the SCBA spectrum. The size of the green markers represent the relative spectral weight intensity of the SCBA dispersion.</figcaption></figure>
by considering a hole propagating in an antiferromagnetic background. The author performed a SCBA calculation as described in Sec. 2.6 on one sublattice only, which accounts well for the low-energy part of the spectrum below the magnetic ordering temperature $T_c$. The agreement with ARPES and CDMFT is excellent, even if the SCBA dispersion is symmetric at the $B$ and $B'$, $A$ and $A'$ points since we considered only one sublattice. The consequence is therefore twofold: first it means that the low-energy band observed in ARPES and CDMFT physically corresponds to a hole propagating in an antiferromagnetic background. This is consistent with the necessity of including spin fluctuations in our effective model and of lowering the temperature so to obtain a good agreement with the experiments. Then, it is a further proof that thinking of t-CuO as two weakly connected sublattices is justified, since a hole propagating inside only one of the two sublattices is enough to get a very good agreement with experiments.

Overall, this underlines that there is a correlation-driven static level splitting present in the ordered phase whereas the paramagnetic phase is driven by dynamic splitting of momentum orbitals. One can note that despite the freezing of the dynamics of the self-energy as decreasing the temperature, at $\beta = 20\text{eV}^{-1}$ (i.e. at the magnetic transition), the spin-resolved self-energy shows an extra dynamic splitting at $w_n \to 0$, see the left panel of Fig. 4.9. The static part given by the high-frequency tail of the self-energy however is constantly increasing as decreasing the temperature. It is surprising to observe a splitting larger than U (the value that should be reached at 0K) at such high temperature, since from a simple screening picture we would expect the splitting to gradually decrease when increasing the temperature. In fact, we show in the following that with an effective single-site model incorporating non-local spin fluctuations in a mean-field fashion, such behavior can be recovered at the magnetic transition temperature.

### 4.5.2 EFFECTIVE SINGLE SITE MODEL

To identify the leading mechanism behind this extra dynamic splitting at the critical temperature, we derive here a simple single-site model including thermal fluctuations as well as antiferromagnetic fluctuations at the mean-field level. We take inspiration from the work of Stepanov et al. in Ref. 291 and consider a single Hubbard site subject to an external magnetic field representing the spin-exchange with the neighboring sites in a mean field fashion:

$$H = -\mu \sum_{\sigma} n_{\sigma} - h(n_{\uparrow} - n_{\downarrow}) + Un_{\uparrow}n_{\downarrow}, \qquad (4.4)$$

where $\sigma = \{\uparrow, \downarrow\}$ is the spin index, $\mu = U/2$ is the chemical potential set for half-filling, $h$ is the effective field, and $U$ the on-site Coulomb interaction. The model is illustrated in Fig. 4.11: being close to the $U \gg t$ limit, and knowing that the system remains insulating at all temperatures, we can consider our system as a collection localized spin-1/2 which interact via an antiferromagnetic exchange coupling (Fig. 4.11(a)). For simplicity, we restrain our model to a single CuO$_2$ sublattice, i.e. an antiferromagnetic square lattice. We then approximate these two-body interaction by a simple staggered magnetic field on each site (Fig. 4.11(b)), which originates from the spin environment, and which in principle could be determined self-consistently. Here, we do not aim at solving exactly this problem, but rather to extract asymptotic behaviors that could explain the
<figure><img src="image_54.png" /><figcaption>Figure 4.11: Sketch of the effective single site model restricted to one sublattice for simplicity. Since $U \gg t$ we can consider that the electrons are rather localized and antiferromagnetically aligned and interacting via the exchange term $J$. This many-body interacting system is approximated by the mean-field model shown in (b) in which the antiferromagnetic fluctuations with nearest-neighbors are integrated out as an effective magnetic field acting on each individual site.</figcaption></figure>

<table><caption>Table 4.1: Eigenvectors and eigenvalues of the Hamiltonian 4.4 for the interacting and non-interacting case ($U = 0$). Note that the chemical potential is $\mu = -U/2$.</caption><thead><tr><th>Eigenstate</th><th>Energy (interacting)</th><th>Energy (non-interacting)</th></tr></thead><tbody><tr><td>|0⟩</td><td>0</td><td>0</td></tr><tr><td>|↑⟩</td><td>-U/2 - h</td><td>-U/2 - h</td></tr><tr><td>|↓⟩</td><td>-U/2 + h</td><td>-U/2 + h</td></tr><tr><td>|↑↓⟩</td><td>0</td><td>-U</td></tr></tbody></table>

extra dynamic splitting in the self-energy.

In order to have access to the self-energy as a function of field $h$ and inverse temperature $\beta$, we determine the non-interacting and interacting Green's functions using the finite-temperature Lehmann's representation for a single site:

$$G_{\sigma}(i\omega_n) = \frac{1}{Z} \sum_{i,j} |\langle i|c_{\sigma}|j\rangle|^2 \frac{e^{-\beta E_i} + e^{-\beta E_j}}{i\omega_n + E_i - E_j},$$

and then use the Dyson equation:

$$\Sigma_{\sigma}(i\omega_n) = [G_{\sigma}^{0}(i\omega_n)]^{-1} - [G_{\sigma}(i\omega_n)]^{-1}. \qquad (4.5)$$
Since the on-site Coulomb interaction is the largest characteristic energy of the system, we assume $\beta U \gg 1$ and $\beta U \gg \beta h$. Using the eigenvectors and eigenvalues of Table 4.1, one can obtain the partition functions:

$$Z_0 = 2e^{\beta(\frac{U}{2})}(\cosh(\beta\frac{U}{2}) + \cosh(\beta h))$$
$$Z = 2(1 + e^{\beta\frac{U}{2}} \cosh(\beta h))$$

where the subscript 0 refers to the non-interacting system. The interacting and non-interacting Green's functions can then be written (using $\beta U \gg 1$ and $\beta U \gg \beta h$ for the interacting case):

$$G_{\uparrow}^{0}(i\omega_{n}) = \frac{1}{i\omega_{n} + h + \frac{U}{2}}$$
$$G_{\downarrow}^{0}(i\omega_{n}) = \frac{1}{i\omega_{n} - h + \frac{U}{2}}$$
$$G_{\uparrow}(i\omega_{n}) \approx \frac{1}{(i\omega_{n} + h)^2 - \frac{U^2}{4}} \left(i\omega_{n} + h - \frac{U}{2}\tanh(\beta h)\right)$$
$$G_{\downarrow}(i\omega_{n}) \approx \frac{1}{(i\omega_{n} - h)^2 - \frac{U^2}{4}} \left(i\omega_{n} - h + \frac{U}{2}\tanh(\beta h)\right).$$

Now, using the Dyson equation 4.5, we obtain the self-energy:

$$\Sigma_{\uparrow}(i\omega_n) = i\omega_n + h + \frac{U}{2} - \frac{(i\omega_n + h)^2 - \frac{U^2}{4}}{i\omega_n + h - \frac{U}{2}\tanh(\beta h)}$$
$$\Sigma_{\downarrow}(i\omega_n) = i\omega_n - h + \frac{U}{2} - \frac{(i\omega_n - h)^2 - \frac{U^2}{4}}{i\omega_n - h + \frac{U}{2}\tanh(\beta h)}.$$

We first need to check that this expression for the self-energy behaves correctly deep in the paramagnetic and in the antiferromagnetic regime. In the paramagnetic limit, i.e. $\beta h \to 0$, one gets:

$$\lim_{\beta h \to 0} \Sigma_{\uparrow}(i\omega_n) = \frac{U}{2} + \frac{U^2}{4(i\omega_n + h)}$$
$$\lim_{\beta h \to 0} \Sigma_{\downarrow}(i\omega_n) = \frac{U}{2} + \frac{U^2}{4(i\omega_n - h)}.$$

The Hubbard-I [35, 201] limit is recovered up to a constant $U/2$ shift, which is consistent with the $\beta = 10\text{eV}^{-1}$ data of Fig. 4.9(a): the two spin flavors become degenerate and the high frequency tail goes to zero ($U/2$ in our expression here due to the chemical potential shift). This is in good agreement with Fig. 4.12 where is shown the splitting in the real part of the self-energy between the two spin flavors as a function of inverse temperature. At high temperature, i.e. $\beta h \to 0$ in our effective model, the splitting vanishes. We also find a good agreement in the antiferromagnetic
<figure><img src="image_55.png" /><figcaption>Figure 4.12: Splitting in the real part of the local component of the Self-energy between the two spin flavors obtained with the dimer cluster as a function of the inverse temperature $\beta$. The shaded vertical region depicts the critical inverse temperature $\beta_c = 14.5 \pm 0.8\text{eV}^{-1}$. The horizontal black dotted line indicates the value of the Coulomb on-site interaction $U = 7\text{eV}$.</figcaption></figure>

limit $\beta h \to \infty$ ($h$ finite), since we recover the real static splitting $U$ between the two different spin flavors:

$$ \lim_{\beta h \to \infty} \Sigma_{\uparrow}(i\omega_n) = 0 $$

$$ \lim_{\beta h \to \infty} \Sigma_{\downarrow}(i\omega_n) = U. $$

This is in very good agreement with the CDMFT calculations at $\beta > 40\text{eV}^{-1}$ shown in Fig. 4.12.

Now that we have shown that the correct limits are recovered as a function of temperature, we place ourselves below the critical temperature $T_c$, so that $\beta h$ is finite, and compute the low- and high-frequency limits:

$$ \Sigma_{\uparrow}(i\omega_n \to \infty) = \frac{U}{2} - \frac{U}{2} \tanh(\beta h) $$

$$ \Sigma_{\downarrow}(i\omega_n \to \infty) = \frac{U}{2} + \frac{U}{2} \tanh(\beta h) $$

$$ \Sigma_{\uparrow}(i\omega_n \to 0) = \frac{U}{2} - \frac{U}{2} \frac{1}{\tanh(\beta h)} $$

$$ \Sigma_{\downarrow}(i\omega_n \to 0) = \frac{U}{2} + \frac{U}{2} \frac{1}{\tanh(\beta h)}. $$

One can immediately see that as the temperature increases, i.e. as $\beta h$ decreases, the splitting of the high-frequency tail for the two spin species decreases. The same trend is observed in the CDMFT calculations close to the transition, as can be seen in Fig. 4.12. Moreover, since in the low-frequency limit the $\tanh(\beta h)$ appears in the denominator, the splitting becomes larger than $U$ as the temperature increases until $h$ remains finite. This is again perfectly consistent with the trend observed in Fig. 4.12, where we find that at the magnetic transition the splitting at the $w_n \to 0$ limit diverges at the transition temperature $T_c$.

Therefore, this simple single-site effective mean-field model can capture the surprising enhancement of the spin splitting in the self-energy at low-frequency. This behavior in the ordered phase
<figure><img src="image_56.png" /><figcaption>Figure 4.13: *Figure taken from Ref. 267.* (a) Internal energy Ω as a function of filling *n* for different solutions within VCA: Antiferromagnetic stripe order (AFS), superconductivity (SC) of *d*<sub>xy</sub> or *d*<sub>x<sup>2</sup>-y<sup>2</sup></sub> symmetry as well as coexistence of all three. In (b) are shown the corresponding order parameters; the colors correspond to the solutions of (a). For all calculations we used the full 8-site diamond cluster of Fig. 4.1(g), i.e. *without* block-construction and optimized the functional in addition with respect to the (cluster) chemical potential *μ* (*μ*').</figcaption></figure>

can in large parts be traced back to an interplay between thermal and magnetic fluctuations at
the magnetic transition.

**4.6 SUPERCONDUCTIVITY**

So far in this chapter we have studied the correlations effects in undoped t-CuO. Although the parent compound already shows remarkable behavior such as the sublattice decoupling, a much richer physics may occur upon hole-doping. Experimentally achieving the chemical doping of t-CuO is extremely difficult, and has not been realized yet, since the standard chemical substitution used for cuprates is impossible due to the lack of charge reservoir layers. Theoretical doping is much simpler since we can tune at will the chemical potential in the effective 2D one-band model. In this section, we study the emergence of superconductivity upon hole-doping: its symmetry and possible interplay with the antiferromagnetic fluctuations. To this aim, B. Lenz performed VCA [68–70] calculations, a method particularly well suited to study the energetics of different symmetry breaking solutions of the model and their competition.

In usual cuprates, superconductivity occurs upon hole-doping and has an order parameter of *d*<sub>x<sup>2</sup>-y<sup>2</sup></sub> symmetry due to the geometry of the CuO<sub>2</sub> layers. In t-CuO however, we expect that the sublattice decoupling survives as doping the material, such that the order parameter would rather be of *d*<sub>xy</sub> symmetry (since the two CuO<sub>2</sub> sublattice are rotated by 45° with respect to the usual cuprates). To check for these different symmetry channels, a paramagnetic calculation is performed first, using the full 8-site cluster of Fig. 4.1(g) *without* block construction, and the internal energy Ω is calculated as a function of doping (ranging from *δ* = 0.0 to *δ* = 0.15). The result is shown as the gray curve in Fig. 4.13(a). Then, calculations were performed for different ordering fields separately: antiferromagnetic stripes, *d*<sub>x<sup>2</sup>-y<sup>2</sup></sub> and *d*<sub>xy</sub> superconductivity. All channels, separately, provide energetically favorable solutions for fillings *n* < 1 compared to the paramagnetic
solution, with the antiferromagnetic stripes being the most favorable at almost all doping until
$n \approx 0.9$. Interestingly, the magnetization is slightly enhanced at very low dopings ($0.95 < n < 1$)
with respect to the undoped solution, see Fig. 4.13(b). While both pairing channels also lower
the internal energy compared to the paramagnetic case, the $d_{xy}$ superconductivity shows a sig-
nificantly larger order parameter than $d_{x^2-y^2}$, consistently with our expectations based on the
sublattice decoupling scenario.

When allowing for competition, we find that the coexistence of superconductivity and antifer-
romagnetic stripes leads to the overall lowest energy solution at zero temperature, see red curve
in Fig. 4.13(a). This is consistent with previous studies finding coexistence of the two phases [50,
145–148]. In Fig. 4.13(b) we observe that the corresponding order parameters are reduced in the
coexistence solution as compared to the pure solutions mentioned above. This indicates a compe-
tition between magnetic and superconducting orders upon doping [145, 146]. Most interestingly,
the $d_{x^2-y^2}$ superconducting order parameter is strongly suppressed by the presence of the antifer-
romagnetic stripe order such that the Cooper pairing is mainly of $d_{xy}$-symmetry. Finally, we note
that superconductivity of $d_{xy}$ symmetry actually corresponds to $d_{x^2-y^2}$ symmetry within each of
the two sublattices (by placing ourselves in the 45° rotated reference frame). Therefore, in the
context of sublattice decoupling, the energetically most favorable solution could be interpreted
as the emergence of a $d_{x^2-y^2}$ superconducting state coexisting with an (Néel) antiferromagnetic
order on each sublattice.

4.7 CONCLUSION

In this chapter, we give a first formal justification of the description of t-CuO as two weakly inter-
acting CuO₂ sublattices [135] which explains the weak symmetry breaking found in ARPES [134].
Indeed, by inspecting selected moments of Matsubara self-energies, we show that the dynamic
correlations effectively weakens the sublattice coupling. Making use of this mechanism, we mo-
tivate the construction of a super-cluster which improves the momentum-resolved spectral func-
tions and equal energy maps computed with CDMFT on the real axis [58, 141–144], and is in ex-
cellent agreement with experiments [134]. We perform calculations at finite temperatures from
which we identify the driving mechanism for the insulating behavior found in the ordered as well
as the PM phases. We find an enhancement of the spin splitting in the Matsubara self-energy at
the magnetic transition temperature, which we trace back to be of intertwined magnetic and ther-
mal origin by constructing an effective single-site mean-field model.

Given the good agreement of our results with experiment, we believe that a minimal one-band
Hubbard model is sufficient to capture most electronic and magnetic properties of t-CuO as long
as dynamical local and short-range fluctuations are treated properly. Further, using VCA we are
able to make predictions about the presence and symmetry of superconductivity upon hole dop-
ing. We find that the decoupling of the sublattices carries through to the superconducting state,
in coexistence with antiferromagnetic stripe order. The $d_{xy}$ symmetry of the superconducting or-
der parameter can be interpreted as a superconducting state of $d_{x^2-y^2}$-type within each sublattice.
Due to its tetragonal symmetry, the lack of interstitial atoms between the well separated 2D CuO layers, and the fact that the electronic properties are mainly governed by its interlaced CuO₂ sublattices, we believe that t-CuO may be the ideal material to gain a more complete understanding of the physics behind cuprate superconductivity. Although doping of t-CuO by chemical substitution is probably not feasible experimentally, the study of doped t-CuO by other experimental techniques like space charge doping [293], which has been successfully applied to other cuprates [294], could be an option. Another interesting route to pursue experimentally consists in growing CuO layers on top of a different substrate. Recently, copper-oxide films have been grown on Bi₂Sr₂CaCu₂O₈+δ, which resulted in nodeless pairing in the superconducting state [295]. Whereas the monolayer was most likely of CuO₂ structure, the possibility of CuO could not be ruled out and is supported by *ab initio* calculations [296].
<figure><img src="image_57.png" /></figure>
# 5 INTERFACE PROPERTIES OF TETRAGONAL
CuO

## 5.1 PUZZLING SPIN AND ORBITAL MOMENTS

As shown in the previous chapter, using CDMFT, we are able to explain the experimentally observed sublattice decoupling [126, 134, 135], reproduce the ARPES spectrum of Ref. 134, and predict the nature of the superconducting order parameter [267]. We moreover confirm the antiferromagnetic stripe order observed in RIXS [135] and µ-SR [136] as stemming from the superposition of two usual antiferromagnetic square sublattices, in agreement with hybrid DFT calculations [128, 129, 131]. Along with the experimental and theoretical knowledge on the structure of t-CuO [124, 125, 128–131], we would now dispose of a coherent understanding of the bulk properties of t-CuO, if it was not for the recent finding of an anisotropic paramagnetic spin moment and a pinned orbital moment in C/t-CuO/STO samples using XMCD at the Cu $L_{2,3}$ edge [136].

The existence of such moments is puzzling since t-CuO is antiferromagnetically ordered and does not contain any element with strong spin-orbit coupling (SOC). Hernandez et al. proposed a simple model [136] (see Fig. 5.1), in which t-CuO would be composed of ferromagnetically ordered CuO layers (in-plane), which are antiferromagnetically stacked along the *c* axis (out-of-plane). The last CuO layer would be paramagnetic hence would follow the magnetic field, while another layer would be uncompensated and therefore would yield the pinned moments. Although this minimal model qualitatively accounts for the experimental observations, it is not satisfactory since:

(i) the ferromagnetic ordering inside each layer is in contradiction with all previous experimental [135] and theoretical [128, 129, 131] findings which support antiferromagnetic stripe ordering of the Cu spins inside the layers,

(ii) the uncompensated layer is composed of spin moments and not orbital ones, whereas the pinned moment is mainly of orbital character,

(iii) there is no reason why one layer should have a paramagnetic behavior, especially if all the Cu spins inside such layer are ferromagnetically ordered.

None of the experimental and theoretical work already published can provide a clear explanation. Our minimal model for bulk t-CuO in the previous chapter does not account for the presence of a paramagnetic moment, neither for an orbital one. More generally, the bulk properties of t-CuO seem incompatible with the XMCD measurements, what calls for a careful analysis of its interface with STO. Indeed, oxide heterostructures are known to display a large range of emergent phenomena [297, 298], and more especially STO-based junctions which can yield a
<figure><img src="image_58.png" /><figcaption>Figure 5.1: Sketch of the model proposed in Ref. 136. The paramagnetic top layer aligns with the external magnetic field, while the pinned and antiferromagnetic layers below are insensitive. The incoming x-ray circularly polarized is shown with the pink arrows (only one of the two circular polarization is illustrated).</figcaption></figure>

2DEG at their interface. The latter is observed in BiMnO₃/STO, LaAlO₃/STO, γ-Al₂O₃/STO and LaVO₃/STO [149, 150, 299, 300] for instance. The existence of a 2DEG at the t-CuO/STO interface is a promising prospect, for it was observed by Mardegan *et al*. that γ-Al₂O₃/STO samples display a paramagnetic spin moment [150], knowing that high electron mobility has been observed at the interface [301–303]. This scenario however does not provide a clear insight for the presence of a pinned orbital moment.

The pinned moment being an order of magnitude smaller than the spin one [136], and its origin being even more unclear, we will focus our modeling on the search of a 2DEG at the t-CuO/STO interface. In the case of BiMnO₃/STO, LaAlO₃/STO, and LaVO₃/STO, the metallic behavior is due to the polar discontinuity at the interface [149, 150, 299, 300]. However, the CuO layers are not polar, and we have seen in the previous chapter that it is a Mott insulator with a gap Δ > 2.35eV. Since STO is a band insulator, the presence of a 2DEG at the bare interface does not seem straightforward. Though, it is known for other heterostructures or bulk materials that the formation of O vacancies, to which STO is highly susceptible, can spark a 2DEG [304–308]. Considering that the STO substrate is fixed at 600K during the epitaxial growth, even if the pressure in oxygen is high and if the sample is annihilated, vacancies are very likely to be formed.

Although the 2DEG scenario is promising, Hernandez *et al.* argued against it since they did not observe any sign of Ti³⁺ in the XMCD signal at the Ti *L*₂,₃ edge [136], contrary to LaAlO₃/STO and γ-Al₂O₃/STO in which a spin moment is carried on the Ti atoms [149, 150, 299]. In this chapter, using first principle DFT+U calculations, we show that a 2DEG at the t-CuO/STO interface is in fact *not incompatible* with the XMCD measurements. More precisely, the formation of vacancies at the interface can lead to the appearance of a polarized 2DEG hosted in the CuO layer. Most importantly, even if the defect is located in the TiO₂ interface layer, the overall valence of the Ti atoms is unchanged with respect to the bulk, such that the 2DEG is invisible at the Ti *L*₂,₃ edge. In
our calculations, the interface metallic states have a non-zero spin polarization, although there is
no external magnetic field, due to the finite size of the unit cells and the specific arrangement of
the stripes and vacancies. Since the latter are randomly distributed in the real material, without
external magnetic field the 2DEG would be non-polarized, as expected for a paramagnetic spin
moment. We find that both the local and global C₄ symmetries are broken in the CuO layer above
the vacancy, which leads to the pinning of the in-plane component of the orbital moment, while
its out-of-plane counterpart remains collinear and proportional to the out-of-plane spin moment.
Our scenario therefore provides a mechanism for the paramagnetic spin moments, being a con-
sequence of a 2DEG, as well as for the pinned in-plane orbital moment, which results from the
breaking of the local and global C₄ symmetries when forming vacancies in the interface TiO₂
layer.

**5.2 MODEL AND METHOD**

We describe in the following how we model the *t-CuO/STO* heterostructure, with and without
vacancies, and how the calculations are carried out.

5.2.1 MODEL

To model the t-CuO/STO junction we use unit cells of the type shown in Fig. 5.2(a): the CuO layers are stacked onto a 2 unit cell-thick TiO₂-terminated STO substrate [125, 130, 133]. Cu atoms in t-CuO have a half-filled d<sub>xy</sub> orbital which would lead within LDA to a metallic ground state. We therefore apply DFT+U[151-153], as presented in Sec 1.4, using the semilocal Perdew-Burke-Ernzerhof (PBE) functional [163, 164], with a local U = 6eV on the Cu d orbitals. Since XMCD measurements are highly sensitive to the magnetic properties of the sample, we make sure that our unit cells enable the stripe ordering. We moreover double the unit cell in the (x,y) plane (see Fig. 5.4) when inserting an oxygen vacancy and, to keep the computations tractable, we restrict our model to a 4 CuO layer coverage. We refer to the different layers following the nomenclature of Fig. 5.2(a): CuO-{1, 2, 3, 4}, SrO-{1, 2} and TiO₂-{1, 2}. Note that the large unit cells considered also prevent the use of hybrid functionals, which have shown to yield excellent results for t-CuO [128-131]. In particular, Franchini et al. studied the thickness dependence of the t-CuO film on the crystal structure and the DOS, and found that the experimental results were best reproduced with a 8 CuO layer coverage [130]. To be consistent with such results, we decided to use as a starting point the inter-layer distances from Ref. 130. The STO lateral lattice parameter is set to *a* = 3.9Å, as obtained with hybrid DFT [165, 309] in excellent agreement with experiments [310].

5.2.2 METHODOLOGY

First an ionic relaxation is performed using the VASP [252–255] with a 6 × 6 × 1 Monkhorst-Pack **k**-grid, until the maximum force on each atom is smaller than 10meV. We use projector augmented wave-PBE pseudopotentials [311] with a cut-off energy of 400eV for the bare interface, which is increased up to 600eV for the larger unit cells with vacancies. The CuO layers are
<figure><img src="image_59.png" /><figcaption>Figure 5.2: (a) Layer-resolved DOS calculated for the t-CuO/STO junction without defect. The colored shaded area depicts the occupied states, and the grey shaded patches highlight the gap in each layer. We also illustrate the crystal structure of the t-CuO/STO junction (before ionic relaxation). (b) Projected DOS on the p<sub>x,y,z</sub> orbitals of the O atoms in the two TiO<sub>2</sub> layers corresponding to the two framed panels in (a).</figcaption></figure>

relaxed while the STO atoms are kept fixed except when considering an oxygen vacancy in the TiO₂-2 interface layer, in which case the whole structure is optimized. We did check that keeping the STO structure fixed did not change any conclusion. In fact, in the absence of vacancies in the STO layers the relaxation is small, in contrast to t-CuO which always shows a large buckling (see Fig. 5.11). We illustrate it for instance in the case of the bare interface in Fig. 5.3. Observing larger distortions in t-CuO seems quite natural since its crystal structure is sparser than that of its substrate, and we will show in the following that those distortions have a direct impact on the electronic structure and the magnetic moments.

In order to check the consistency of the VASP calculations, we also performed PBE+U calcu-
lations using Wien2k package [173, 174]. The results are qualitatively the same for all the cases
considered, however we notice that despite using the same PBE functional and U value, the gap
obtained in Wien2k is always larger (of about 0.5eV for the CuO layers) than in VASP (for U = 0
the two codes agree), as shown in Fig. 5.3 for the bare interface. Changing the smearing factor
and type in VASP does not change this discrepancy, as well as the wavefunction cut-off energy.
We then perform a Wien2k calculation on the relaxed structured resulting from VASP for each
of the cases considered, and from which all analysis and conclusions will be drawn. A 6 × 6 × 1
Monkhorst-Pack k-grid is used, and due to computational limitations we managed the size of the
planewave basis set by reducing the RKmax value when necessary (the smallest value being 5.14).
We did check for the bare interface that lowering the RKmax value did not change significantly
the DOS. To have access to the value of the orbital moment we also performed calculations with
SOC on the Sr, Ti and Cu atoms.

5.3 BARE INTERFACE

From previous studies on the *t-CuO/STO* junction we know that the bare (i.e. defect-free) inter-
face should not yield a *2DEG* [130, 133]. This is consistent with our results, see Fig. 5.2(a), in
<figure><img src="image_60.png" /><figcaption>Figure 5.3: Layered resolved DOS for the bare junction using Wien2k (**left**) and VASP (**center**) with relaxation of the CuO layers only, and VASP with relaxation of the whole structure (**right**). For the VASP calculations we used a Gaussian smearing $\gamma = 0.05$eV.</figcaption></figure>

which we plot the layer-projected DOS. The gap $\Delta_{DFT+U} \approx 1.5$eV in t-CuO is layer independent and in good agreement with previous DFT+U calculations [133], although smaller than the photoemission lower bound $\Delta \ge 2.4$eV [134]. The situation is different in the STO substrate where a clear difference can be seen between the interface and bulk TiO$_2$ layers: the gap is reduced from 2eV in the bulk (TiO$_2$-1) to 1.5eV at the interface (TiO$_2$-2). The states reducing the gap are traced back to be of O $p_x/p_y$ origin. Fig. 5.2(b) shows that the $p_x/p_y$ orbital states are pushed towards the Fermi level at the interface, whereas the $p_z$ orbital states are stable in energy and even slightly de-populated at the interface as evidenced by the loss of density around -3.5eV. This is consistent with previous first principles calculations performed for TiO$_2$-terminated STO and BaTiO$_3$, and is the consequence of the emergence of pure O $p_x/p_y$ states not hybridized with the Ti orbitals at the surface [312, 313]. Interestingly, adding a thin film of t-CuO above STO does not change this behavior. Due to the large original band gap of STO the system remains insulating since no 2DEG is observed at the bare interface, and the t-CuO antiferromagnetic stripe order is not disturbed.

As expected, from the bare interface we can not draw an explanation of the XMCD results since the system remains insulating and perfectly antiferromagnetic, ruling out any possibility of a paramagnetic spin moment. In the following we explore the effects of inserting O vacancies at the interface layers. Among the many possibilities to do so, we focus on five of them which should already give a clear trend of the effects on the electronic properties. We first study the case of a vacancy in the CuO-1 layer (Fig. 5.4(a)) in Sec. 5.4. Then, in Sec. 5.5 we focus on the most promising case in which a vacancy is considered in the TiO$_2$-2 layer (Fig. 5.4(b)). Since a 2DEG
<figure><img src="image_61.png" /><figcaption>Figure 5.4: (a)-(c) Zoom on the two interface CuO-1 and TiO₂-2 layers presenting the 5 vacancy configurations considered in this work. The black lines between the O and Cu/Ti atoms show that half of the O in the CuO-1 layer do not have a Ti neighbor.</figcaption></figure>

appears in the latter case, we increase the vacancy concentration further (Fig. 5.4(c)) trying two
different configurations.

**5.4 VACANCY IN THE CuO INTERFACE LAYER**

There are two geometrically different ways of creating an O vacancy in the interface CuO-1 layer (see Fig 5.4(a)): either by removing an O having a neighboring Ti in the layer below, or an O having no such neighbor. It is important to note that these two cases preserve the global C₄ rotational symmetry both in the CuO-1 and the TiO₂-2 layers. For the two configurations interestingly give the same results, we mostly focus here, in Fig. 5.5, on the ones obtained considering the vacancy without neighboring Ti atom. STO behaves similarly to the bare interface configuration, see Fig. 5.5(a). The gap is again reduced by 0.5eV at the interface, but with the addition of small extra contributions around ±0.4eV sparked by the large in-gap states appearing in the interface CuO-1 layer. The latter shows a drastic change compared to the bare interface case since a new set of in-gap states of mostly Cu *d* orbital character appears. Contrary to the usual effect of chemical doping in cuprates, where a ZRS [38] is formed between the Cu *d* and O *p*, here the O does not significantly contribute. Those states are localized around the vacancy, as can be seen in Fig. 5.5(b) where we show the projected DOS on the NN Cu atoms with respect to the vacancy, NNN and NNNN. The fact that only the four NN have a sizable contribution, that the oxygen contribution is weak, and that the surrounding layers DOS display a vanishing density at ±0.4eV show the localized nature of these states, both in- and out-of-plane. As expected, it is the *d*<sub>xy</sub> orbital of the NN Cu atoms surrounding the vacancy that mostly contributes to these in-gap states, see Fig. 5.5(c), although we notice a sizable contribution from the *d*<sub>z²</sub> orbital in the occupied region.
<figure><img src="image_62.png" /></figure>

Figure 5.5: (a) Layer-resolved DOS calculated for the case of an O vacancy in the CuO-1 layer. (b)-(c) Site and orbital projected DOS for the Cu atoms at the CuO-1 interface layer corresponding to the framed panel in (a). (d) Spin density map of the interface CuO-1 layer. (e) Band structure and projected DOS zoomed on the in-gap states along with the fitted Wannier bands (top) resulting in the wannier orbitals showed in the *bottom* around the vacancy, which are centered on the four NN Cu sites highlighted with the dotted frame in (d). The band color depicts the ratio $(\lambda_{xy} / \lambda_{z^2})^2$ between the contributions from $d_{xy}$ and $d_{z^2}$.
<table><caption>Table 5.1: Hopping parameters extracted from the Wannier fit showed in Fig. 5.5(e) for only two of the four NN Cu atoms (labeled 1, 2). The unit-cell translation indicates if the hopping is inside the first unit-cell (0, 0, 0), or outside (x, y, 0).</caption><thead><tr><th>Unit-Cell translation</th><th>t<sub>1,1</sub> (eV)</th><th>t<sub>1,2</sub> (eV)</th><th>t<sub>2,1</sub> (eV)</th><th>t<sub>2,2</sub> (eV)</th></tr></thead><tbody><tr><td>(0, 0, 0)</td><td>-0.002</td><td>0.364</td><td>0.364</td><td>-0.002</td></tr><tr><td>(1, 0, 0)</td><td>-0.006</td><td>-0.001</td><td>-0.001</td><td>-0.006</td></tr><tr><td>(-1, 0, 0)</td><td>-0.006</td><td>-0.001</td><td>-0.001</td><td>-0.006</td></tr><tr><td>(0, 1, 0)</td><td>-0.001</td><td>0.038</td><td>-0.003</td><td>-0.001</td></tr><tr><td>(0, -1, 0)</td><td>-0.001</td><td>-0.003</td><td>0.038</td><td>-0.001</td></tr></tbody></table>

The missing O leads to an increase of the electronic density on the surrounding NN Cu sites, which results, in combination with the $d_{z^2}/d_{xy}$ orbital mixing, in a peculiar effective orbital carrying the in-gap states on the NN sites. It is shown in the spin density plot of Fig. 5.5(d): the local $C_4$ symmetry on each NN Cu site is broken, but the global symmetry around the vacancy is conserved. The symmetry of the in-gap states around the Fermi level and the position of the Cu atoms around the vacancy advocate for the formation of bonding/anti-bonding localized states around the vacancy, thus preventing the formation of an electron gas. This can be seen by performing a fit of the in-gap bands with Wannier functions using wannier90 [183, 184], as shown in Fig. 5.5(e). The two bands are degenerate in spin, hence we present the fit performed only for the up spin part (the down spin giving the exact same results). We performed the fit starting from the two $d_{xy}$ orbitals of two of the NN Cu atoms, i.e. two Cu for the up spin flavor, and obtained a perfect matching with the two PBE+U bands. Consistently with the spin density, the Wannier orbitals break the local $C_4$ symmetry, and they expand towards the vacancy position in the middle of the four NN Cu. They interestingly have a $d_{z^2}$-like shape but rotated in-plane with the main lobe pointing towards the vacancy. This translates in the tight-binding Hamiltonian obtained from the Wannierization to a large hopping amplitude between the NN Cu sites inside the unit-cell $t_{in} \approx 0.364$eV, and small hopping outside $t_{out} \approx t_{in}/10$, with on-site energy $\epsilon \approx 0$. Hence, the splitting between the two bands $\Delta \approx 2t_{in}$ is consistent with a dimerization picture, and large comparatively to the small bandwidth $W$ since $t_{out}$ is one order of magnitude smaller than $t_{in}$, see Table 5.1.

Alike the previous case, the system remains insulating (although the gap is reduced at the interface CuO-1 layer) and antiferromagnetic, i.e. no paramagnetic spin moments would be observed. Still, it is interesting to study the orbital moments obtained from the calculations including SOC. Indeed, as advocated in Ref. 314, a crystalline defect could be at the origin of an orbital momentum pinning. If the magnetization axis is set out-of plane (0,0,1), the orbital moment is parallel and proportional to the local spin, such that on the NN Cu sites the absolute value of the orbital moment is reduced (see Fig. 5.6(a)). On every other site the magnetization is $m_{spin} \approx \pm 0.70\mu_B$ and $m_{orb} \approx \pm 0.14\mu_B$, while on the four NN sites it is weaker: $m_{spin} \approx \pm 0.45\mu_B$ and $m_{orb} \approx \pm 0.09\mu_B$. Moreover, as the four NN sites are equally dispatched on the up and down spin stripes the total orbital moment is zero. Most importantly, the out-of plane part of the orbital moment is collinear with the spin.
<figure><img src="image_63.png" /><figcaption>Figure 5.6: Schematic representation of the orbital momentum orientation for two different magnetization axis: (a) (0,0,1), (b) (1,0,0). The length of the vectors illustrates the variation of the absolute value of the moments.</figcaption></figure>

The situation is different when setting the magnetization axis in-plane. If it is along (1,0,0), or (0,1,0), the orbital moment is collinear with the spin on all the Cu sites but the four NN ones, on which it is pinned along (±1,±1,0), as illustrated in Fig. 5.6(b). Moreover, its value is not proportional anymore to the spin moment, it is even inversely proportional if the magnetization axis is set along (1,1,0) for the two NN atoms aligned along (1,1,0) around the vacancy. Indeed, in the case of (1,1,0), most of the sites display a spin moment $m_{spin} \approx \pm0.70\mu_B$ and $m_{orb} \approx \pm0.05\mu_B$ (along the magnetization direction), unless for the two aforementioned sites which carry $m_{spin} \approx \pm0.45\mu_B$ and $m_{orb} \approx \pm0.07\mu_B$. Therefore, although the total orbital moment averages out to zero because the global $C_4$ symmetry is not broken, we identify here a pinning mechanism which is intimately linked with the breaking of the on-site $C_4$ symmetry. On the four NN sites, the in-plane part of the orbital moment can be non-collinear and its absolute value even inversely proportional to that of the spin moment.

None of these conclusions change when considering an O vacancy lying above a Ti site instead of vacuum. We do not extensively describe this case here, as it is very similar, but for completeness we show the layered-DOS in Fig. 5.7. The same in-gap states appear, again confined in the CuO-1 layer only, and since this type of vacancy also preserves the global $C_4$ symmetry of the unit cell, the very same arguments and conclusions apply.

Hence, while this scenario is not suited for explaining the paramagnetic spin moment feature, it provides precious insights on the orbital moment behavior close to a vacancy. A key argument in favor of this pinning mechanism is that it would be consistent with the XMCD measurements in which only the in-plane component of the orbital moment is pinned [136]. What is missing in
<figure><img src="image_64.png" /><figcaption>Figure 5.7: (a) Layer-resolved DOS calculated for the case of an O vacancy in the CuO-1 layer above a Ti site. (b)-(c) Site and orbital projected DOS for the Cu atoms at the CuO-1 interface layer corresponding to the framed panel in (a).</figcaption></figure>

these two configurations is the breaking of the *global* $C_4$ symmetry to allow the emergence of a non-zero in-plane pinned orbital moment in the unit cell.

## 5.5 VACANCY IN THE TiO$_2$ INTERFACE LAYER

The physics is different if the O vacancy is created in the TiO$_2$-2 layer instead (Fig. 5.4(b)). The layer- and spin-resolved DOS are presented in Fig. 5.8(a): the CuO layer at the interface is metallic and spin polarized. Despite the global shift of the chemical potential, the system is only metallic at the interface since sizable DOS at the Fermi level only appears there. Similarly to the previous cases, those states are mostly of Cu *d* orbital character with almost no contribution from the *O* *p* orbitals. A striking result is that the DOS of the interface TiO$_2$-2 layer, where the O vacancy is located, is not affected by the defect: the valence of the Ti atoms is unchanged and they remain non-magnetic. This is an important result, since it shows that it is possible to generate a polarized 2DEG at the interface CuO-1 layer which is invisible at the Ti $L_{2,3}$ edge, since their *d* orbitals remain empty and non polarized. Moreover, note that the O vacancy does not intrinsically favor
<figure><img src="image_65.png" /><figcaption>Figure 5.8: (a) Layer-resolved DOS for up (left) and down (right) spins in the case of a vacancy in the TiO$_2$-2 layer. (b) Projected band structure and DOS on the Cu $d_{xy}$ and $d_{z^2}$ orbitals of the interface CuO-1 layer. The band color depicts the ratio between the contributions from $d_{xy}$ and $d_{z^2}$, while the opacity indicates their overall contribution to the bands. The black arrows indicate the bandwidth $W_{↑/↓}$.</figcaption></figure>

positive spin moments. The sign of the resulting moment depends on the defect and the stripes position so that if we shift the vacancy, or if we swap the up and down spin stripes we obtain a down-spin polarized electron gas as shown in Fig. 5.9, obtained by keeping the same structure but inverting up and down spins. Therefore, on average in the real system without external magnetic field the 2DEG should not be polarized, as expected for paramagnetic spin moments. These considerations contradict the claim in Ref. 136 that a 2DEG at the interface should be seen in the XMCD spectrum at the Ti $L_{2,3}$-edge, and restore it as a plausible explanation for the puzzling paramagnetic moments observed at the Cu $L_{2,3}$-edge. For the sake of clarity we will only discuss the configuration of Fig. 5.8 in the remainder of this section, keeping in mind that spins are ex-changeable.

In Fig. 5.8(b) we show the spin and orbital projected band structure. There is a striking difference between the two spin flavors: we find a single down spin $d_{xy}$ band crossing the Fermi level, whereas we observe for the up spin two split bands of a mixed Cu $d_{xy}/d_{z^2}$ character. The orbital mixing and the bands shape are analogous to the SOC-and-distortion-induced splitting of the Ir d bands in Sr$_2$IrO$_4$ [285], although it is remarkable that in our case only the up spin electrons are concerned. If the splitting of the up spin band is due to the structural distortions similarly to Sr$_2$IrO$_4$, the band mixing however is due to the $C_4$ symmetry breaking, not SOC. We note that the bands at the Fermi energy have a small bandwidth $W_{↑/↓} \approx 300$meV, hence being sign of a low-mobility 2DEG.

In order to better understand the spin up/down difference, we show in Fig. 5.10 the spin-density at the interface CuO-1 layer around its initial $Z_0$ position. We observe a large distortion: the spin
<figure><img src="image_66.png" /><figcaption>Figure 5.9: Layer-resolved DOS for up (left) and down (right) spins in the case of a vacancy in the TiO<sub>2</sub>-2 layer. We used the same structure as in Fig. 5.8(a) but up and down spin stripes have been swapped.</figcaption></figure>
<figure><img src="image_67.png" /><figcaption>Figure 5.10: Spin density map of the interface CuO-1 layer above (top), at (center), and below (bottom) the initial position Z₀. As a guide to the eyes, the position of the O vacancy and the Ti atoms in the TiO₂-2 layer below is shown. The black frame highlight the three Cu atom on which is focused Fig. 5.11(b).</figcaption></figure>
density for the same layer is non zero over more than 1Å along the c axis. Two Cu¹,² atoms be-
longing to the down spin stripe, see the top panel, show a breaking of the local C₄ symmetry due
to the extra electrons coming from the TiO₂-2 layer below, similarly to the four NN Cu sites in
the previous section. Interestingly, these charges are not located on the Cu³ site just above the
vacancy, but rather on the Cu¹,² sites. This is caused by the large buckling induced by the intro-
duction of the vacancy, as can be seen in Fig. 5.11(a). The average displacement along the c axis
per layer is evaluated as follow:

$$
d_{\text{layer}} = \frac{1}{N_{\text{layer}}} \sum_{i \in \text{layer}} |d_i| = \frac{1}{N_{\text{layer}}} \sum_{i \in \text{layer}} |z_i - z_{\text{avg}}| \quad (5.1)
$$

where Nlayer is the number of atom per layer, zᵢ the z position of the atom, and zavg the average
z position of the layer. It is larger at the CuO-1 interface layer, and is of the order of 0.5Å. It is
interesting to note that the distortions are larger than for the bare interface and the case of a va-
cancy in the CuO-1 layer, although here the defect is located in the TiO₂-2 layer below. Moreover,
the overall distortion follows the sublattice decoupling observed experimentally [126, 134, 135],
that we have explained in the previous chapter by means of CDMFT [267]. Note that the upturn
in the atomic displacement in STO as going towards the first layer is due to the vacuum below the
unit cell since the whole cell is optimized. It remains significantly smaller than the buckling in
the CuO-1 layer. Such displacement induces large changes in the Cu-O bonds, especially around
the Cu³ located above the vacancy, see Fig. 5.11(b). Since the latter has no O neighbor below and
its in-plane Cu-O distances are smaller than for the two surrounding Cu¹,² atoms, the splitting
between the dₓᵧ and d_z² orbitals is expected to be enhanced. The dₓᵧ energy in the two neighbor-
ing Cu¹,² sites is lower because their Cu-O distance is larger, i.e. the in-plane crystal field acting
on dₓᵧ is reduced, while the d_z² is higher in energy because of the O below. Therefore, the extra
charges from the TiO₂-2 layer will occupy the Cu¹,² sites belonging to the down-spin stripe since
they are more favorable in energy. Note that the Cu¹,² atoms are two of the four NNN of Cu³,
and it is not trivial why those two sites would prevail the two other NNN and even the four NN
Cu sites. The reason is that for all these sites, their Cu-O bonds are shorter than the ones of Cu¹,²
because they remain closer to their O neighbors, as can be seen in the *middle* and *bottom* panels
of Fig. 5.10, i.e. the same argument applies.

Since all the *d* orbitals but *d*<sub>xy</sub> are filled, and the *d*<sub>xy</sub> is already occupied with a down spin elec-
tron, the extra electrons have an up spin such that the electron gas is polarized. The breaking
of the local *C*<sub>4</sub> symmetry induces the mixing of the *d*<sub>xy</sub> and *d*<sub>z<sup>2</sup></sub> orbitals, which are furthermore
closer in energy on the *Cu*<sup>1,2</sup> sites. Finally, the splitting Δ of the bands in the up spin DOS is due
to a dimerization: the neighboring *Cu*<sup>1,2</sup> sites couple and contribute equally to the band above
and below the Fermi energy, resulting in the doubling of the unit cell for the up spin electrons, in
contrast to the *Cu*<sup>3</sup> site above the vacancy which solely contributes to the single down spin band
crossing the Fermi energy. Δ is however significantly smaller (Fig. 5.8(b)) that in the previous
case (Fig. 5.5(e)): a Wannierization of the up spin bands using only the two neighboring sites
gives an intra-unit-cell hopping *t*<sub>in</sub> ≈ 0.035eV, and a inter-unit-cell hopping of the same order of
magnitude *t*<sub>out</sub> ≈ *t*<sub>in</sub>.
<figure><img src="image_68.png" /></figure>

Figure 5.11: (a) Buckling on a portion of the crystal structure (**left**), illustration of how are evaluated the displacement $d_i$ (**top right**) and evaluation of the average buckling for each layer for the three first scenarios (**bottom right**) following Eq. 5.1. (b) View of the three Cu$^{1,2,3}$ atoms highlighted in Fig. 5.10 with a schematic representation of the crystal field splitting of the $d_{xy}$ and $d_{z^2}$ orbitals.

An important difference with the case of a vacancy in the CuO-1 layer is that in the latter case the $C_4$ symmetry around the vacancy is conserved (although locally it is broken on the **NN** Cu site). Here, however, the global $C_4$ symmetry is broken as can be seen in the **bottom panel** of Fig. 5.10: among the four O$^{1,2,3,4}$ atoms located in the CuO-1 layer around the vacancy position, only two of them (O$^{1,2}$) are located above Ti sites. Those two atoms are pulled down during the relaxation process, weakening the Cu-O-Cu bond between the three Cu$^{1,2,3}$ sites highlighted by the dotted frame in Fig. 5.10. Hence, in contrast with the previous case, the four **NNN** Cu sites around the vacancy position are not equivalent, and the extra electrons are hosted by only two of them (Cu$^{1,2}$), leading to a breaking of the global unit cell $C_4$ symmetry.

If the magnetization points out-of-plane, so is the orbital moment and it is *blind* to the in-plane symmetry breaking, such that the system is equivalent to the case where a vacancy is inserted in the CuO-1 layer. Indeed, the orbital moment remains proportional and collinear to the spin one, and is therefore reduced on the two Cu$^{1,2}$ sites: on most of the Cu sites $m_{spin} \approx \pm 0.70\mu_B$ and $m_{orb} \approx \pm 0.14\mu_B$, while on the two **NNN** ones $m_{spin} \approx -0.36\mu_B$ and $m_{orb} \approx -0.08\mu_B$. This leads to a non-zero out-of-plane orbital moment in the unit-cell, as illustrated in Fig. 5.12 where are also plotted the Wannier orbitals obtained by fitting the two split bands in the up spin band structure as well as the lonely band crossing the Fermi energy in the down spin sector. Hence, in the real material where we would have a random superposition of vacancies giving rise to up and down spin polarized **2DEG**, the out-of-plane component of the orbital moment should also average to zero.

If the magnetization is in-plane instead, we observe again a pinning of the in-plane orbital moment component on the three aforementioned sites, as shown in Fig. 5.12. In our configuration, the two Cu$^{1,2}$ sites are aligned along (1,1,0), such that the magnetization axis along (1,1,0) and (1,-
<figure><img src="image_69.png" /><figcaption>Figure 5.12: Schematic representation of the orbital momentum orientation for two different magnetization axis: (a) (0,0,1), (b) (1,0,0). The length of the vectors illustrates the variation of the absolute value of the moments. For the sake of clarity the isovalue of the Wannier orbitals was set to relatively high value (0.43), preventing from the appearance of surfaces elsewhere in the unit cell.</figcaption></figure>

1,0) are no longer equivalent (contrary to the case of a vacancy in the CuO layer, see Fig. 5.6). We observe that, if we choose the magnetization axis along (1,1,0), the orbital moment is increased in absolute value on these two Cu<sup>1,2</sup> sites ($m_{orb} \approx -0.09\mu_B$ instead of $m_{orb} \approx \pm 0.05\mu_B$ on the other sites) although they have a reduced spin-moment with respect to the other sites ($m_{spin} \approx -0.36\mu_B$ instead of $m_{spin} \approx \pm 0.70\mu_B$). Again, if the magnetization axis is along (1,0,0) or (0,1,0), the in-plane orbital moment on all the Cu atoms is collinear with the spin moment, unless for the two Cu<sup>1,2</sup> atoms for which it is pinned along (-1,-1,0), as well as the Cu<sup>3</sup> site above the vacancy for which it is pinned along (1,1,0). The sum of the three pinned in-plane components is non-zero and positive in this configuration.

Therefore, we have shown that a 2DEG appears in the CuO-1 interface layer when creating a vacancy in the TiO<sub>2</sub>-2 layer, which is invisible to XMCD at the Ti L<sub>2,3</sub> edge. Around the vacancy, the orbital moment is locally pinned along (1,1,0) direction and can be non-collinear with the spin moments, with an overall non-zero value. If we swap the stripes, we obtain a non-zero orbital moment of opposite sign, i.e. a random superposition should average to zero if the magnetization axis is the same. Also, the same effect could be seen in a (1,-1,0) configuration, in which case the orbital moment would be pinned along (1,-1,0) direction. If we consider a random distribution of these four configurations ((1,1,0) or (1,-1,0) and up and down spin polarized), and also consider that they all share the same magnetization axis, then the global orbital momentum should average to zero.

However, if the magnetization axis is different at the different vacancy sites, then it is possible to obtain a non-zero orbital moment because the of pinning, while the average spin moment remains zero. The striking difference here is that one vacancy gives rise to a non-zero pinned or-
bital momentum, and it is possible to consider a distribution of vacancies resulting in an average zero magnetization but non-zero orbital moment, whereas in the previous section it is impossible since the orbital momentum is locally averaged to zero because of the global $C_4$ symmetry conservation around the vacancy site. This is a key feature of this scenario that should apply on similar types of monoxide/dioxide junctions: inserting a vacancy in the dioxide interface layer induces a $C_4$ symmetry breaking in the monoxide layer above which could be at the origin of a non-zero pinned orbital moment. Most importantly, only the in-plane component of the orbital moment is pinned, whereas the out-of-plane part is collinear and proportional to the spin, in excellent agreement with the XMCD measurements finding no out-of-plane pinning [136].

Finally, one may have noticed that with our choice of on-site interaction $U = 6$eV, and in the absence of vacancies, we obtain a gap for t-CuO of roughly 1.5eV, much weaker than the ARPES lower bound $\Delta > 2.35$eV [134] and the STM upper bound $\Delta < 3.68$eV [126]. Since we get a gap of about 2eV for STO, then with our choice of $U$ value the first empty states available are in the t-CuO interface layer and we obtain the physics we have described above. By choosing a higher $U$ value, say to better reproduce the experiments, the t-CuO first available states would be pushed to higher energies, such that they could be replaced by the STO states at some critical value of the on-site interaction. In such a case, upon the formation of vacancies we could expect that the extra charges would remain in the TiO₂-2 interface layer, in contradiction with the results we have drawn in this chapter. However, PBE also underestimates the gap of STO, which is found experimentally to be $\Delta_{STO} \approx 3.3$ [315], i.e. larger than the ARPES lower bound and of the same order than the STM upper bound for t-CuO. Given that STM overestimates the gap values due to the band bending phenomenon [316], especially in the case of insulators with a large gap, it is reasonable to consider that the gap configuration we obtain in our calculations is similar to that of the real system, albeit an overall underestimation of both gap values. In other words, at the interface, the first available states to be filled upon inserting vacancies in the TiO₂-2 interface layer should be located in the CuO-1 layer, and our observations apply.

## 5.6 TWO VACANCIES IN THE TiO₂ INTERFACE LAYER

An important question is whether this scenario still holds at a higher concentration of vacancies in the TiO₂-2 layer. For instance in SrVO₃ two neighboring oxygen vacancies, instead of a single one, are needed to reproduce the additional spectral weight seen in photoemission [317]. Fig. 5.13 shows the layer- and spin-resolved DOS obtained for a larger concentration of oxygen vacancies with two different configurations (see Fig. 5.4(c)): either two vacancies still far from each other, or forming a pair around a single Ti atom. Both configurations induce an electron gas in the interface CuO-1 layer while the TiO₂-2 layer is still insulating. Since we doubled the number of vacancies, there is an additional charge transfer to the interface CuO-1 layer which translates into a larger DOS below the Fermi level. Similarly to the previous section, the additional states at the Fermi level are almost purely of Cu d character.

When the vacancies form a pair around the same Ti atom, the latter is not able to totally transfer
<figure><img src="image_70.png" /><figcaption>Double neighbour O vacancy in TiO<sub>2</sub></figcaption></figure>

<figure><img src="image_71.png" /><figcaption>Figure 5.13: Layer-resolved DOS for the two interface layers in case of higher concentration of vacancies.</figcaption></figure>
the excess charges and becomes slightly spin-polarized ($m_{spin} \approx -0.42\mu_{B}$), contrary to the case in which the vacancies are still apart from each other for which no spin moment is observed on the Ti atom. Such a spin moment on the Ti atoms is not consistent with the XMCD measurements since it would be detected at the Ti $L_{2,3}$ edge, and allows us to discard the scenario of oxygen vacancies pairing around the same Ti atom. In contrast, when the two vacancies do not combine around the same Ti atom, the 2DEG remains invisible at the Ti $L_{2,3}$ edge, and we observe an in-plane pinning of the orbital moment in the CuO-1 layer.

## 5.7 CONCLUSION

To summarize, performing DFT+U calculations on large supercells including oxygen vacancies at the t-CuO/STO interface, we are able to identify a scenario compatible with the puzzling XMCD measurements of Ref. 136. Namely, we show that oxygen vacancies in the TiO$_{2}$-2 layer induce a 2DEG confined in the CuO-1 layer at the interface, which would explain the experimental findings: we observe no valence change nor non-zero magnetic moment at the interface Ti sites, i.e. the 2DEG is invisible at the Ti $L_{2,3}$ edge, and show that a vacancy induces the appearance of paramagnetic metallic states with a mixed $d_{xy}/d_{z^2}$ character. We identify the breaking of both the local and global $C_{4}$ symmetries in the CuO-1 layer above the vacancy as being a mechanism for the emergence of an in-plane pinned orbital moment, in excellent agreement with the nature of the pinned moment observed in XMCD [136]. Such phenomenon is not restricted to the case of t-CuO/STO junctions since it simply stems from the stacking of a monoxide layer above a dioxide and should then generalize to other heterostructures. We moreover show that if the oxygen vacancies form pairs around a single Ti atom, the latter becomes spin polarized and is therefore not invisible anymore in XMCD.

This scenario also provides an understanding of the apparent inconsistency between the previous experimental observations and the XMCD measurements. Indeed, if the moments and the 2DEG live at the interface, as we argue in this chapter, then ARPES would not be able to resolve them since the probing depth is of the order of 2nm at the photon energy used in Ref. 134 and the t-CuO films were ~ 3.6nm-thick. In contrast, XMCD in total electron yield mode has a probing depth of 5 – 10nm [318] and would therefore be sensitive to the interface. The RIXS measurements were performed on t-CuO samples synthesized on Nb-doped STO substrates [135]. Such Nb-doped substrates are metallic and would most probably change drastically the interface properties as compared to the t-CuO/STO junctions measured in XMCD. Moreover, even in the case of an undoped STO substrate, RIXS is not suited for observing paramagnetic spin moments, neither pinned orbital moments, especially when the measurements are performed without external magnetic field.

It would therefore be crucial to further investigate the interface properties of t-CuO/STO both experimentally and theoretically. Local probes like STM would be good candidates in order to study specifically the interface layers. Another possibility would be to generate vacancies with irradiation, as was done for SrVO$_{3}$ [317], although one would need to create them specifically at
the *t*-CuO/STO interface. Since *t*-CuO was successfully synthesized up to 10 unit cells [126], it might be possible to use the small probing depth of *XMCD* in total electron yield to suppress, or enhance, the interface contribution to the spectrum, as was done in Ref. 149 for LaAlO$_{3}$/STO. Changing the termination of the *STO* substrate could also be a way to suppress the 2DEG, so the paramagnetic spin moment, as it was shown for instance that the AlO$_{2}$/STO interface becomes insulating when *STO* is SrO-terminated [306]. Finally, we show that the orientation of the magnetization axis should play an important role in the emergence of the pinned orbital moment so that its experimental determination is an important step in the understanding of the observed magnetic moments.
# PART V

## CONCLUSION
<figure><img src="image_72.png" /></figure>
Throughout this thesis, we studied the rich and complex physics of strongly correlated copper oxide systems by focusing on two different compounds: $Ca_{2-x}Na_xCuO_2Cl_2$ (NaCCOC) and tetragonal CuO (t-CuO). These two materials are particularly interesting since they have a simple crystal structure, even without charge reservoir layers for t-CuO, which makes the comparison between the cuprates' low-energy effective models [30–33, 38, 263] and experiments easier. We explored their spectral and magnetic properties by means of state-of-the art numerical techniques, density function theory (DFT)(+U) [158, 159] and cluster dynamical mean-field theory (CDMFT) [42, 43, 49–51] with 2 × 2 and 8-site clusters, that were systematically compared to experimental results from the literature and from our two angle-resolved photoemission spectroscopy (ARPES) measurement campaigns, the first performed at the Photon Factory (Tsukuba, Japan), and the second one at SOLEIL (Saint-Aubin, France).

For the low energy properties of both NaCCOC and t-CuO, we used a single-band Hubbard model [35–37] based on the Zhang-Rice singlets (ZRS) picture [38]. The framework of CDMFT naturally gives access to the spectral function which can be directly compared to the ARPES measurements. We found an excellent agreement for NaCCOC at $x = 0, x = 0.05$ and $x = 0.1$, despite a missing quasiparticle peak at the nodal point [92] in the ARPES measurements at $x = 0.05$ and $x = 0.1$. The theoretical spectrum captures the waterfall feature [24, 86] at the nodal point $(\pi/2, \pi/2)$ with both cluster sizes, with the kink at the right position in momentum and energy, both for the undoped and doped systems. We interpreted this spectral feature as a cross-over between a high-energy part where the correlations are purely local and which is well described by the Hubbard-I approximation [35], and a low-energy part where electrons are dressed by the antiferromagnetic spin correlations and form spin polarons [95–98]. This result already underlines the uttermost necessity of including non-local correlations when treating the cuprates' low energy models, since the purely local Hubbard-I picture fails to describe the physics close to the Fermi level. Moreover, especially with the 8-site cluster we were able to reproduce the two-peak structure observed at the nodal point [91, 92], which consists of a sharp nodal quasiparticle peak at the Fermi level emerging from a less coherent dispersive feature at around −0.2eV.

In the spectral function of NaCCOC, we observed the opening of a pseudogap (PG) upon hole doping, both in our calculations and ARPES measurements, in agreement with previous experimental [12, 14, 24, 26, 90–94] and theoretical works [104, 107, 111–114, 246]. In the latter, the PG was traced back to short-range spin fluctuations. Our analysis of the simple 2 × 2 cluster results is in line with this interpretation: the state histogram obtained from the impurity solver revealed that it is dominated by only three groups of states. The leading contribution is the singlet state $|N = 4, S = 0\rangle$, which is the most probable state from $x = 0$ until $x \sim 0.12$. Within our new choice of basis, this state involves a spin singlet between the D13, D42 orbitals of the plaquette cluster (as defined in Sec. 3.4), which induces the suppression of the zero-energy excitations when the system is connected to a bath of non-interacting electrons [121]. The flat dispersion, i.e. the van Hove singularity, at the nodal point, which would jump at the Fermi level upon hole-doping, is therefore split into two subbands below and above the Fermi level since the first remaining photoemission channel is a transition from the singlet state $|N = 4, S = 0\rangle$ to the triplet $|N = 4, S = 1\rangle$.
Interestingly, in the state $|N = 4, S = 1\rangle$ the D13, D42 orbitals are in a spin triplet state: we therefore interpreted the formation of the PG in the spectral function as stemming from the breaking of the singlets into triplets due to the cancelling of the zero-energy photoemission channel [121]. This is intimately related to the short-range spin fluctuations, which provides the non-degenerate singlet nature of $|N = 4, S = 0\rangle$. Furthermore, we showed that this is not a special feature of the 2 × 2 plaquette cluster since the 8-site cluster also opens a PG. While the $(\pm\pi/2, \pm\pi/2)$ orbitals included in this cluster could be expected to behave similarly to the D13, D42 orbitals, we argued that their hybridization to the bath is larger both at the non-interacting and interacting level. An isolated ground state on the cluster is no longer a good approximation for the $(\pm\pi/2, \pm\pi/2)$ orbitals, such that there is no suppression of the quasi-particle peak at the nodal point.

We therefore showed that the waterfall feature and the opening of the PG in NaCCOC, and cuprates in general, can be understood under the same framework of short-range spin fluctuations. Short-range correlations are also essential to find the origin of the sublattice decoupling in t-CuO, which can not be captured with purely local fluctuations. Indeed, we demonstrated with two different impurity solvers that it is induced by the suppression of the inter-sublattice self-energy component, which leaves the hopping $t_d$ bare while the NNN hopping term $t$ is renormalized by a factor ~ 2 by the correlations. Because it is energetically favorable to form bonds between sites linked with the leading order hopping term, and that all hopping processes involving $t$ only connect sites from the same sublattice, the CuO planes can be thought of as two weakly coupled CuO$_2$ sublattices in agreement with the experimental observations [134, 135]. We made use of this finding to motivate the construction of a super-cluster, allowing to improve the momentum resolution of our calculations. We obtained a spectral function reproducing astonishingly well the ARPES measurements [134], both using the continuous-time quantum Monte Carlo (CTQMC) [60, 62, 116] solver working at finite temperature used by the author, and the matrix-product state (MPS)-based [56–58, 141–144] solver working directly on the real frequency and at zero temperature used by M. Bramberger and M. Grundner. Combining these different solvers is extremely useful for it allows to connect and compare the finite to zero temperature results and properties.

From calculations at finite temperature, we showed that there exist two different insulating regimes in t-CuO: a paramagnetic regime at high temperature in which the self-energy shows a momentum dependent divergence at $\omega_n \to 0$, and a magnetically ordered regime at low temperature characterized by a static spin dependent shift in the self-energy. The transition between these different phases yields drastic transformations of the spectral functions: when decreasing the temperature, the lower Hubbard band (LHB) originally composed of a single broad dispersive feature turns into two split sub-bands. Only the spectral function calculated at low-temperature with the CTQMC solver corresponds to the one computed with the MPS-based solver at effectively zero temperature, and so to the ARPES measurements [134]. This underlines the effects of the spin fluctuations on the spectral properties of cuprates, which are further confirmed by the fact that the low-energy subband at low temperature is well accounted for by a self-consistent Born approximation (SCBA) calculation, similarly to NaCCOC. We moreover identified an intriguing
enhancement of the spin dependent shift in the real part of the self-energy at the transition temperature, that can be explained at the local mean-field level as a combination of spin and thermal fluctuations at the transition, based on a simple effective model derived by the author.

To further understand the behavior of one-particle quantities like the spectral function, it is useful to study the two-particle observables. From the analysis of the dynamical spin structure factor (DSSF) in NaCCOC, we confirmed the predominance of the spin-fluctuations at $Q = (\pi, \pi)$ at all dopings, in agreement with early quantum Monte Carlo (QMC) calculations [111], with the fluctuation diagnostic method applied to dynamical cluster approximation (DCA) [112] and diagrammatic Monte Carlo [113, 114] calculations, as well as the dual boson expansion around the DMFT solution [75]. Most interestingly, we identified a freezing of the spin-spin correlation function upon hole doping at long times $\tau = \beta/2$ and low temperatures $\beta \ge 40$, in line with recent findings [122, 123]. It can be explained by relating $S(\tau)$ to its real frequency counterpart, which directly shows that the finite value at $\tau = \beta/2$ is related to the presence of low-energy spin excitations. At $x \ne 0$, the DSSF’s gap closes: a low-energy and almost doping independent spectral weight appears, leading to a doping independent increase of $S(\tau = \beta/2)$ with respect to the undoped system. At too large doping, however, the freezing disappears due to the weakening of the spin fluctuations [123]. Moreover, we showed that the value of the correlation function at $\tau = \beta/2$ does not depend on the inter-site distance. While this was observed with the small plaquette cluster, we argued that it should hold for larger systems given that the DSSF peak at $(\pi, \pi)$ is broad in momentum and finite at low-frequency. This is an intriguing result that motivates for further investigations with larger system sizes.

Although most of the results in this thesis could not be obtained with DFT, it remains a necessary technique to be kept in the theoreticians’ toolkit for the study of strongly correlated systems’ ground state properties. The t-CuO/STO interface is a good example, since DFT+U is the only possible framework to investigate the emergence of paramagnetic spin and pinned orbital moments observed with x-ray magnetic circular dichroïsm (XMCD) at the Cu $L_{2,3}$ edge [136]. Using PBE+U [151–153], we showed that the formation of O vacancies in the TiO$_2$ interface layer can induce the formation of a 2D electron gas (2DEG), which can explain the measurement of a paramagnetic spin moment. Above all, this 2DEG would be invisible at the Ti $L_{2,3}$ edge since the valence of the Ti atoms in the TiO$_2$ interface layer remains unchanged, which is consistent with the experimental observations [136]. The large structural distortions, as well as the mismatch of the monoxide CuO and the dioxide TiO$_2$ layer induce the breaking of both local and global $C_4$ symmetries at the interface, which in turn leads to the in-plane pinning of the orbital moment. We emphasized that only the out-of-plane component is not pinned, in agreement with the XMCD measurements [136]. Therefore, using PBE+U we were able to construct a scenario that reconciles the XMCD observations with the physics of t-CuO/STO.

These findings remain to be confirmed experimentally, and would probably require the use of experimental techniques that can include or exclude at will the influence of the interface. Since XMCD is always performed in total electron yield for these kind of systems, a surface dependent
mode, one could grow thin films of t-CuO of varying thickness and study the dependence of the
magnetic moments, similarly to Ref. 149. If the sample is thick enough, the XMCD becomes
insensitive to the interface, and, according to our model, the spin and orbital moments should
disappear. Another way would be to change the TiO₂ termination of STO to SrO, which could for
instance prevent the formation of a 2DEG at the interface as is observed at the AlO₂/STO inter-
face [306, 319].

Further works are to be carried out on the doped phase of t-CuO, which has not yet been re-
alized experimentally, and has not yet been investigated theoretically apart from the variational
cluster approximation (VCA) calculations presented in this thesis [267]. It would be interesting
to further study the possible consequences of the sublattice decoupling, other than the unusual
dxy superconducting order parameter, on the expected phases of cuprates like the PG, or charge
density wave (CDW).

For the PG phase, based on the conclusions we drew from the study of NaCCOC, questions arise
about the behavior of the spin-freezing phenomenon: how does it evolve as a function of the
cluster size, how does it behave at the transition with superconductivity ? For the effect of clus-
ter size, it would be extremely interesting to confirm the observation of the spin-freezing in a
8-site DCA cluster [123] with CDMFT, especially since we would directly have access to the
inter-site correlation functions. This question opens a more general perspective on other tech-
niques than quantum cluster methods [52], such as diagrammatic expansions around the DMFT
solution [71-76]. A great advantage of these techniques is that they include long-range corre-
lations, and allow to obtain a much better momentum resolution for the one- and two-particle
quantities. On different grounds, computationally cheaper techniques like the two-particle self-
consistent method [320, 321], or the fluctuating local field [322, 323], would be useful tools to
study large systems while still including collective effects non-perturbatively.

In this thesis, we confronted the rich and non-trivial parent cuprate compound as well as the
PG phase. This is only a limited part of the zoo of phenomena observed in this class of materi-
als [14, 24]. While we have shown the importance of spin fluctuations, one should not forget that
they live and interact with other types of fluctuations, as the charge channel which is for instance
highlighted by the observation of CDW [25-29]. This complicates even more the understanding
of copper oxides, and most notably the origin of superconductivity which remains elusive. For
all these reasons, copper oxide systems are among the most interesting materials under study
nowadays. One should therefore think twice before throwing an old oxidized saucepan to the
trash, for it might yield an unexpected behavior.
# PART VI

## APPENDIX
<figure><img src="image_73.png" /></figure>
A PROBING THE ELECTRONIC STRUCTURE:
ANGLE-RESOLVED PHOTOEMISSION
SPECTROSCOPY

Angle-resolved photoemission spectroscopy (ARPES) is one of the most powerful experimental techniques to probe the momentum-resolved electronic structure of materials. It relies on the photoelectric effect, discovered by Hertz [324], Hallwachs [325] and Lenard [326], and explained by Einstein [327] in 1905, which was first used to experimentally measure atomic binding energies in 1957 by Siegbahn [328]. The principle is simple, and sketched in Fig. A.1: when a beam of light is shone on a sample, the electrons of the solid can be emitted into the vacuum by absorbing a photon. The measure of the emission angles ($\phi, \theta$) and the kinetic energy $E_{kin}$ of the photoelectrons in the analyzer allows to reconstruct the momentum-resolved electronic structure. In this appendix, we briefly present the ARPES cross-section and how it is linked to the Green's function that we calculate using the many-body methods presented in Chap. 2. We also show how to transform the raw experimental data, which is obtained as a function of angles ($\phi, \theta$), back to the momentum space. We do not aim at going into the details of the technique, for which we refer the interested reader to the reviews of Damascelli [24] and Moser [258].

<figure><img src="image_74.png" /><figcaption>Figure A.1: Sketch (**left**) of the ARPES experimental setup, and photography of the setup at SOLEIL (Saint-Aubin, France) synchrotron (**right**).</figcaption></figure>
## A.1 THE PHOTOEMISSION CROSS-SECTION AND THE GREEN'S FUNCTION

In a simplified notation, we can write the Hamiltonian of the material as:

$$H = \frac{1}{2m_e} \hat{\mathbf{p}}^2 + \hat{V}, \qquad (A.1)$$

where $\hat{\mathbf{p}}$ is the momentum operator, $\hat{V}$ the potential, we set $\hbar = 1$, and the bold characters denote the vector quantities. When light is shone on the sample, the momentum operator $\hat{\mathbf{p}}$ has to be replaced as follows:

$$\hat{\mathbf{p}} \rightarrow \hat{\mathbf{p}} - e\mathbf{A}, \qquad (A.2)$$

where $\mathbf{A}$ is the vector potential. The Hamiltonian is therefore transformed:

$$\begin{aligned} \tilde{H} &= \frac{1}{2m_e} (\hat{\mathbf{p}} - e\mathbf{A})^2 + \hat{V} \\ &= \underbrace{\frac{1}{2m_e} \hat{\mathbf{p}}^2 + \hat{V}}_{H} - \underbrace{\frac{e}{m_e} \mathbf{A} \cdot \hat{\mathbf{p}}}_{H_{int}} - \frac{e}{2m_e} (-i\nabla \mathbf{A}) + \frac{e^2}{2m_e} \mathbf{A}^2. \end{aligned} \qquad (A.3)$$

The third term is equal to zero except possibly at the surface, it is then neglected in the following. Similarly for the fourth term which describes two-photon processes, and can be ignored for the range of photon energies used in ARPES [24, 258]. We get to the approximate Hamiltonian

$$\tilde{H} \approx H + H_{int}, \qquad (A.4)$$

with $H_{int}$ encoding the light-matter interaction. It induces optical excitations from a $N$-electron ground state $|\psi_{GS}^N\rangle$ to a final state $|\psi_f^N\rangle$ containing $N-1$ electrons and one photoelectron. $H_{int}$ can be considered as a weak perturbation to $H$, such that Fermi's golden rule can be used to account for the transition probability:

$$w_f = 2\pi |\langle\psi_f^N| H_{int} |\psi_{GS}^N\rangle|^2 \delta(E_f^N - E_{GS}^N - \nu), \qquad (A.5)$$

where $E_f^N$ and $E_{GS}^N$ are the final and ground state energies, and $\nu$ is the photon energy. We assume to be at zero-temperature, such that the initial many-body state is conveniently restricted to the ground state $|\psi_{GS}^N\rangle$.

At this stage, there are two possibilities to calculate the transition probability. The most rigorous one is the one-step model, in which the photon absorption, electron removal and detection are treated as a single coherent process [329]. This treatment is extremely complex, and out of the scope of this appendix. We rather consider the phenomenological three-step model, in which the photoemission is divided into three separated steps:

(i) Creation of a photoelectron in the solid by the absorption of a photon,

(ii) Travel to the surface,

(iii) Emission into vacuum.