One- and Two-Particle Correlation Functions in the Cluster Perturbation Theory for Cuprates

Abstract

The physics of high-$T_c$ superconducting cuprates is obscured by the effect of strong electronic correlations. One way to overcome this problem is to seek an exact solution at least within a small cluster and expand it to the whole crystal. Such an approach is at the heart of cluster perturbation theory (CPT). Here, we developed CPT for the dynamic spin and charge susceptibilities (spin-CPT and charge-CPT), with the correlation effects explicitly taken into account by the exact diagonalization. We applied spin-CPT and charge-CPT to the effective two-band Hubbard model for the cuprates obtained from the three-band Emery model and calculated one- and two-particle correlation functions, namely, a spectral function and spin and charge susceptibilities. The doping dependence of the spin susceptibility was studied within spin-CPT and CPT-RPA, that is, the CPT generalization of the random phase approximation (RPA). In the underdoped region, both our methods resulted in the signatures of the upper branch of the spin excitation dispersion with the lowest excitation energy at the $(\pi ,\pi )$ wave vector and no presence of low-energy incommensurate excitations. In the high doping region, both methods produced a low energy response at four incommensurate wave vectors in qualitative agreement with the results of the inelastic neutron scattering experiments on overdoped cuprates.

1. Introduction

The electronic and magnetic subsystems of a correlated material are strongly coupled. Classic examples are the cuprate high-$T_c$ superconductors, whose phase diagram contains antiferromagnetic (AFM), metallic, charge-density wave, and pseudogap states. Undoped cuprates corresponding to the half-filled system are AFM insulators. In this case, the spin dynamics correspond to the excitation of magnons and could be described by linear spin-wave theory. The physics, however, become more complicated once the holes are doped in the system. The holes are delocalized, and the AFM long-range order quickly disappears, allowing the new cooperative phases to arise. Recent resonant inelastic X-ray scattering (RIXS) studies have revealed magnon-like excitations in cuprates over a wide doping range from underdoped to heavily overdoped systems [1,2,3,4]. Moreover, momentum-dependent charge excitations have been observed in RIXS in the pseudogap phase of cuprates [5]. The spin and charge dynamics are coupled to the doping-dependent changes in the electronic structure observed in angle-resolved photoemission spectroscopy (ARPES) [6,7]. Furthermore, electronic excitations can be described by one-particle Green functions, and the calculation of the magnetic and charge excitations relies on two-particle response functions. Two-particle correlation functions provide important information about the ordered phases of a strongly correlated system. On the experimental side, inelastic neutron scattering allows one to directly observe the dynamic spin susceptibility, i.e., the two-particle spin–spin correlation function [8,9,10,11,12,13,14].

Cuprates have a quasi-two-dimensional structure, and the conductivity is provided by the electrons in the copper–oxygen plane. This was the reason for intensive studies involving two-dimensional models and, in particular, the effective low-energy Hubbard model, which provided a simple approach that retained the essential physics. However, a detailed study of the electronic and magnetic properties of cuprates requires a more realistic model, such as the three-band Emery model that includes both Cu-$d_{x^2-y^2}$ and O-$p_{x,y}$ orbitals [15].

There are a number of methods for studying the physics of strongly correlated systems and cuprates in particular. One subclass of these methods was developed to progress beyond a mean-field theory by relying on the numerical solution of a somehow simplified problem. These approaches include quantum cluster theories [16], extensions of the dynamic mean-field theories [17,18], quantum Monte Carlo approaches [19,20,21], and variational Gutzwiller-type methods [22]. Here, we focus on cluster perturbation theory (CPT). CPT is a straightforward way of calculating one-particle correlation functions, i.e., spectral functions [23,24]. The latter can be directly compared to ARPES data. CPT is one of a number of quantum cluster theories [16], and it is the most economical cluster method in terms of the necessary computation power. In CPT, the first step is the exact diagonalization (ED) of a small cluster. Therefore, the short-range correlations are treated precisely. In the second step, the intercluster interactions are included according to some kind of perturbation theory. There have been few attempts to expand CPT for the calculation of two-particle correlation functions. The authors of [25] used variational cluster approximation (CPT modified by a self-consistent procedure) to solve the Bethe–Salpeter equation for a two-dimensional Hubbard model. In [26], the spin susceptibility was calculated using the determinant quantum Monte Carlo method and CPT for the one-band Hubbard model. The authors of [27] extended CPT to compute two-particle correlation functions by approximately solving the Bethe–Salpeter equation for the one-band one-dimensional Hubbard model.

Here, we developed CPT for dynamical spin susceptibility, an approach we call spin-CPT. This was based on an explicit calculation of correlation functions via exact diagonalization with a subsequent extraction of the two-particle spectral lattice function in a CPT-like manner, similar to [26,28,29]. Then, we applied spin-CPT to a two-dimensional model of cuprates (the effective Hubbard model for the CuO plane based on the Emery model). We compared the results of the spin-CPT and CPT-RPA approaches [27,30]. The latter is a straightforward generalization of the random phase approximation (RPA), where the bare-electron Green functions forming the electron–hole bubble are replaced by those obtained by CPT. We showed that spin-CPT produced a low-energy response at four incommensurate wave vectors, which qualitatively agreed with the results of the inelastic neutron scattering on overdoped cuprates. Additionally, it allowed us to obtain a spectral intensity distribution resembling the upper branch of an experimentally observed hourglass dispersion in a wide doping range.

This paper is organized as follows: In Section 2, we discuss the model and approximations used for the study. In Section 3, the main results are presented. In Section 4, the results are discussed. In Section 5, the concluding remarks are provided.

2. Model and Methods

2.1. Model

The Hamiltonian of the Emery model [15] describes the copper $d_{x^2-y^2}$ and oxygen $p_{x,y}$ orbitals:

(1)$\begin{array}{ccc}\hfill H_{pd}& =& \sum _{i,\sigma }{\epsilon }_d{n}_{i\sigma }^d+\sum _{j,\sigma }{\epsilon }_p{n}_{j\sigma }^p+\lambda U_d\sum _i{n}_{i\uparrow }^d{n}_{i\downarrow }^d+\lambda U_p\sum _j{n}_{j\uparrow }^p{n}_{j\downarrow }^p+\lambda V_{pd}\sum _{⟨i,j⟩}{n}_i^d{n}_j^p\hfill \\ & +& t_{pd}\sum _{⟨i,j⟩,\sigma }{\left(-1\right)}^{P_{ij}}\left(d_{i\sigma }^{†}{p}_{j\sigma }^{}+\mathrm{H}.\mathrm{c}.\right)+t_{pp}\sum _{⟨⟨j,j^{\prime }⟩⟩,\sigma }{\left(-1\right)}^{P_{j{j}^{\prime }}^{\prime }}\left(p_{j\sigma }^{†}{p}_{j^{\prime }\sigma }^{}+\mathrm{H}.\mathrm{c}.\right),\hfill \end{array}$

where i denotes the copper sites, j denotes the oxygen sites, $n_{i\left(j\right)\sigma }$ is the number operator of holes on site $i\left(j\right)$ with spin $\sigma$, $n_{i\left(j\right)}=n_{i\left(j\right)\uparrow }+n_{i\left(j\right)\downarrow }$, $P_{ij}$ and $P_{j{j}^{\prime }}^{\prime }$ are phase factors, and $d_{i\sigma }^{}$ ($p_{j\sigma }^{}$) destroys a hole with spin $\sigma$ on a $d_{x^2-y^2}$ ($p_{x\left(y\right)}$) orbital and site i (j). The Coulomb interaction normalization constant $\lambda$ is introduced for convenience to allow the variation of the interaction strength.

The following hopping parameters are used (here and below, all energy values are in eV): $t_{pd}=1.36$, $t_{pp}=0.86$. These are the hopping integrals of the multiband $p-d$ model calculated for La${}_2$CuO${}_4$ [31] by projecting [32] the ab initio local density approximation electronic structure onto the Wannier function basis. The Coulomb parameters used here are $U_d=9$, $U_p=4$, and $V_{pd}=1.5$. The charge-transfer energy parameter is taken to have a typical value, ${\epsilon }_p=3.6$ [33]. The one-electron energy of the copper d-orbital is set to zero, ${\epsilon }_d=0$.

Since CPT has a cluster at its heart, the choice of this cluster is important. In a CuO plane, one can naturally distinguish a single element consisting of a copper ion surrounded by four oxygen atoms. The oxygen, however, belongs to the two neighboring cells. To overcome this difficulty, we used the canonical fermions of Shastry [34], following [35,36,37,38,39,40,41,42], thereby orthogonalizing the $p_{x,y}$ orbitals and ending up with a representation of a square lattice of orthogonal CuO${}_4$ cells:

(2)$\begin{array}{ccc}\hfill H& =& \sum _{f,\sigma }\left[\sum _{\alpha }{\epsilon }_{\alpha }{n}_{f\sigma }^{\alpha }-2t_{pd}{\mu }_{00}\left(d_{f\sigma }^{†}{b}_{f\sigma }^{}+b_{f\sigma }^{†}{d}_{f\sigma }^{}\right)\right]+\sum _{f,\alpha }\lambda U_{\alpha }{n}_{f\uparrow }^{\alpha }{n}_{f\downarrow }^{\alpha }+\sum _{f,\alpha <\beta }\lambda V_{\alpha \beta }{n}_f^{\alpha }{n}_f^{\beta }\hfill \\ & +& \sum _{f\ne g,\sigma }\left[-2t_{pd}{\mu }_{fg}\left(d_{f\sigma }^{†}{b}_{g\sigma }^{}+b_{f\sigma }^{†}{d}_{g\sigma }^{}\right)+2t_{pp}{\nu }_{fg}\left(a_{f\sigma }^{†}{a}_{g\sigma }^{}-b_{f\sigma }^{†}{b}_{g\sigma }^{}\right)\right.\hfill \\ & -& \left.2t_{pp}{\chi }_{fg}\left(a_{f\sigma }^{†}{b}_{g\sigma }^{}+b_{f\sigma }^{†}{a}_{g\sigma }^{}\right)\right]+H_{int}^{cc},\hfill \end{array}$

where f indexes a cell position (site); ${\epsilon }_{\alpha }$ is the one-electron energy for an orbital index $\alpha$, which takes values $\left\{a,b,d\right\}$, where d is the copper orbital and a and b are the cell oxygen orbitals with the energies ${\epsilon }_a={\epsilon }_p+2t_{pp}{\nu }_{00}$ and ${\epsilon }_b={\epsilon }_p-2t_{pp}{\nu }_{00}$, respectively; the operator $d_{f\sigma }^{}$ annihilates a hole with spin $\sigma$ on a copper orbital at site f; operators $a_{f\sigma }^{}$ ($b_{f\sigma }^{}$) annihilate a hole with spin $\sigma$ on an oxygen orbital a (b) at site f; $n_{f\sigma }^{\alpha }$ are the number operators; and $n_f^{\alpha }=n_{f\uparrow }^{\alpha }+n_{f\downarrow }^{\alpha }$. The Wannier coefficients $\mu$, $\nu$, and $\chi$ and other constants related to the parameters of Equations (1) and (2) are the same as in [38,40,41,42]. In particular, the effective Coulomb interactions and the original interactions are related as follows: $U_a=U_b\approx 0.21U_p$, $V_{ad}=V_{bd}\approx 0.91V_{pd},U_{ab}\approx 0.17U_p$. Finally, $H_{int}^{cc}$ refers to all non-local interaction terms, including three-site and four-site ones, which did not have a significant influence on the results obtained using Equation (2) and are omitted hereafter.

Such a cell representation is especially useful when one is interested in the low-energy properties. Since only b orbitals are directly connected by the hopping processes with d orbitals, they should have a negligible effect on the low-energy part of the electronic structure. On the other hand, at a high doping level, a orbitals might become important [36]. In Figure 1, a comparison of the spectral functions obtained with and without the a orbital is shown for the overdoped case with a doping level of $p=0.25$ using a $2\times 2$ cluster consisting of four orthogonal cells. A distinguishable difference was observed only 2 eV below the Fermi level within the bands dominated by the oxygen spectral weight [43]. Since the needed computational power rose with the number of possible states, which increased with the number of orbitals considered, we could increase the cluster size once we removed one of the orbitals. Thus, in the main part of this paper, we omit the a orbitals. This allowed us to extend the cluster size within CPT to $3\times 3$, comprising nine orthogonal cells, with moderate computational effort.

2.2. Methods

The central object of our calculations is the transverse part of the dynamic spin susceptibility ${\chi }_{\alpha \beta ,\gamma \delta }$($\mathbf{q}$, $\Omega )$. It is defined as an analytical continuation of the Matsubara two-particle spin–spin correlation function

(3)${\chi }_{\alpha \beta ,{\alpha }^{\prime }{\beta }^{\prime }}(\mathbf{q},\mathrm{i}{\omega }_n)=\underset{0}{\overset{1/T}{\int }}d\tau e^{\mathrm{i}{\omega }_n\tau }⟨T_{\tau }{S}_{\alpha {\alpha }^{\prime }}^{+}(\mathbf{q},\tau )S_{{\beta }^{\prime }\beta }^{-}(-\mathbf{q},0)⟩.$

where $S_{\alpha {\alpha }^{\prime }}^{+}$ and $S_{{\beta }^{\prime }\beta }^{-}$ are the spin raising and lowering operators, respectively; $\mathbf{q}$ is the wave vector; ${\omega }_n$ is the Matsubara frequency; $\alpha$, $\beta$, ${\alpha }^{\prime }$, and ${\beta }^{\prime }$ are the orbital indices; $T_{\tau }$ is the ordering operator over Matsubara (imaginary) time $\tau$; and T is the temperature in the energy units. Averaging is carried out with the fully interacting Hamiltonian. By using the Wick’s theorem, we obtain the zeroth order approximation for the spin susceptibility in terms of the one-particle Green function $G_{\alpha \beta \sigma }(\mathbf{p},\mathrm{i}{\omega }_m)$ for an electron with spin $\sigma$:

(4)${\chi }_{\alpha \beta ,{\alpha }^{\prime }{\beta }^{\prime }}^{\left(0\right)}(\mathbf{q},\mathrm{i}{\omega }_n)=-T\sum _{{\omega }_m,\mathbf{p}}{G}_{\alpha \beta \uparrow }(\mathbf{p},\mathrm{i}{\omega }_m)G_{{\alpha }^{\prime }{\beta }^{\prime }\downarrow }(\mathbf{p}+\mathbf{q},\mathrm{i}{\omega }_n+\mathrm{i}{\omega }_m).$

The analytical continuation $\mathrm{i}{\omega }_n\to \Omega +\mathrm{i}\eta$ includes the positive infinitesimal parameter $\eta$ that we use to introduce Lorentzian broadening in the course of the numerical calculations. Since we considered the paramagnetic phase, the spin index was omitted.

The physical spin susceptibility is calculated on real frequencies $\Omega$ as the sum over all orbital indices coinciding at vertices:

(5)$\chi (\mathbf{q},\Omega )={\displaystyle \frac{1}{2}}\sum _{\alpha ,\beta }{\chi }_{\alpha \alpha ,\beta \beta }(\mathbf{q},\Omega ).$

Analogously to the spin correlation function (3), the charge correlation function is defined as

(6)${\Pi }_{\alpha \beta ,{\alpha }^{\prime }{\beta }^{\prime }}(\mathbf{q},\mathrm{i}{\omega }_n)=\underset{0}{\overset{1/T}{\int }}d\tau e^{\mathrm{i}{\omega }_n\tau }⟨T_{\tau }{\rho }_{\alpha {\alpha }^{\prime }}(\mathbf{q},\tau ){\rho }_{{\beta }^{\prime }\beta }(-\mathbf{q},0)⟩.$

where ${\rho }_{\alpha {\alpha }^{\prime }}$ is the particle density operator. The physical charge susceptibility on real frequencies $\Omega$ is calculated as

(7)$\Pi (\mathbf{q},\Omega )={\displaystyle \frac{1}{2}}\sum _{\alpha ,\beta }{\Pi }_{\alpha \alpha ,\beta \beta }(\mathbf{q},\Omega ).$

2.2.1. CPT

The electronic structure of a two-band model is studied using CPT, as discussed in detail in [23,24]. There are two steps in this approach. The first is the exact solution for the small cluster usually achieved by the exact diagonalization of the model Hamiltonian. Thus, the lattice is partitioned into a superlattice of clusters with a new translational order of an artificial origin. The second step is the reconstruction of the whole lattice by adding intercluster interactions, which is carried out in a perturbative manner.

For the Hamiltonian (2), we calculate the electronic spectral function

(8)$A_{\alpha \beta }(\mathbf{k},\omega )=-{\displaystyle \frac{1}{\pi }}\mathrm{Im}{G}_{\alpha \beta }(\mathbf{k},\omega +\mathrm{i}\eta ),$

where $G_{\alpha \beta }(\mathbf{k},\omega )$ is the electronic Green function on the real frequencies $\omega$, $\mathbf{k}$ is a wave vector, $\alpha$ and $\beta$ are the orbital indices, and $\eta \to 0^{+}$. The trace of (8) provides the total (physical) spectral function $A(\mathbf{k},\omega )={\displaystyle \sum _{\alpha }}{A}_{\alpha \alpha }(\mathbf{k},\omega )$. CPT allows one to calculate $A(\mathbf{k},\omega )$ within a small cluster and an arbitrary momentum resolution; it is also characterized by a fast convergence of results with an increasing cluster size [44]. Here, the intracluster interactions are taken into account by means of the ED of a cluster with open boundary conditions using the Lanczos method, while the intercluster hoppings are treated perturbatively. This approach allowed us to include the short-range (intracluster) correlations explicitly. The calculations are carried out at a zero temperature using a square cluster consisting of nine cells (sites).

The main CPT equation is

(9)$\widehat{G}{(\tilde{\mathbf{k}},\omega )}^{-1}={\widehat{G}}^c{\left(\omega \right)}^{-1}-\widehat{T}\left(\tilde{\mathbf{k}}\right),$

where all square matrices have dimensions $L\times N_c$; L is the number of orbitals, which in our case was equal to two; $N_c$ is the number of sites within the cluster; $\tilde{\mathbf{k}}$ is a wave vector defined within the cluster Brillouin zone; $\omega$ is a frequency; ${\widehat{G}}^c\left(\omega \right)$ is the exact local (intracluster) propagator; and $\widehat{T}\left(\tilde{\mathbf{k}}\right)$ is a Fourier transform of the hopping matrix. The translational invariance is restored in CPT as

(10)$G_{\alpha \beta }(\mathbf{k},\omega )=\frac{1}{N_c}\sum _{i,j}{e}^{-i\left({\mathbf{r}}_i-{\mathbf{r}}_j\right)\mathbf{k}}{G}_{\alpha \beta ,ij}(\mathbf{k},\omega ),$

where $i\left(j\right)$ is an intracluster site index, and ${\mathbf{r}}_{i\left(j\right)}$ is a corresponding radius vector.

2.2.2. CPT-RPA

CPT can be used to improve the calculation of the dynamic spin susceptibility $\chi (\mathbf{k},\Omega )$ within RPA by replacing the bare electron–hole bubble, ${\chi }^{\left(0\right)}(\mathbf{k},\Omega )$, with that calculated using CPT Green functions. This allows one to include self-energy corrections over the vertex renormalization of RPA. This methodology is discussed in detail in [27,30]. Thus, given an RPA vertex, the obtained magnetic structure is determined by the electronic spectrum. Instead of the “bare” susceptibility, the electron–hole bubble is calculated via the cluster electronic spectral functions [30]:

(11)${\chi }_{\alpha \beta ,\gamma \delta }^{c\left(0\right)}(\mathbf{k},\Omega )=-\sum _{\mathbf{q}}\int \int d{\omega }^{\prime }d{\omega }^{″}{A}_{\alpha \beta }(\mathbf{q},{\omega }^{\prime })A_{\gamma \delta }(\mathbf{k}+\mathbf{q},{\omega }^{″}){\displaystyle \frac{f\left({\omega }^{\prime }\right)-f\left({\omega }^{″}\right)}{{\omega }^{\prime }-{\omega }^{″}+\Omega +\mathrm{i}\eta }},$

where $f\left(\omega \right)$ is the Fermi function, and $\eta$ is the positive infinitesimal.

The transverse dynamic susceptibility within CPT-RPA is calculated as [45,46]

(12)${\chi }_{\alpha \beta ,\gamma \delta }(\mathbf{k},\Omega )={\chi }_{\alpha \beta ,\gamma \delta }^{c\left(0\right)}(\mathbf{k},\Omega )+\sum _{{\alpha }^{\prime },{\beta }^{\prime },{\gamma }^{\prime },{\delta }^{\prime }}{\chi }_{\alpha \beta ,{\alpha }^{\prime }{\beta }^{\prime }}^{c\left(0\right)}(\mathbf{k},\Omega ){\overline{U}}_{{\gamma }^{\prime }{\delta }^{\prime }}^{{\alpha }^{\prime }{\beta }^{\prime }}{\chi }_{{\gamma }^{\prime }{\delta }^{\prime },\gamma \delta }(\mathbf{k},\Omega ).$

A vertex $\overline{\mathbf{U}}$ including another renormalization constant ${\lambda }^{\prime }$ that will be discussed below is defined as $\overline{\mathbf{U}}={\lambda }^{\prime }\mathbf{U}$, where the nonzero components of the RPA vertex $\mathbf{U}$ are $U_{\alpha \alpha }^{\alpha \alpha }=U_{\alpha }$ and $U_{\alpha \beta }^{\alpha \beta }=V_{\alpha \beta }$, with $\alpha \ne \beta$ [45,46]. The particle–hole bubble ${\chi }^{c\left(0\right)}(\mathbf{k},\Omega )$ was calculated in this study using the CPT spectral function as in [30].

Dynamic charge susceptibility $\Pi (\mathbf{q},\Omega )$ is obtained through CPT-RPA by replacing the vertex $\mathbf{U}$ in the RPA Equation (12) with the vertex $-\mathbf{V}$, where $V_{\alpha \alpha }^{\alpha \alpha }=U_{\alpha }$, $V_{\beta \beta }^{\alpha \alpha }=2V_{\alpha \beta }$, $V_{\alpha \beta }^{\alpha \beta }=-V_{\alpha \beta }$ [45].

2.2.3. Explicit Calculation of the Two-Particle Correlation Function (Spin-CPT and Charge-CPT)

An alternative to the RPA-like approach is to calculate the two-particle correlation function directly without the intermediate steps of finding the one-particle Green function and building a spin or charge correlation function out of it. Here, we calculated the right side of Equation (3) in a CPT-like manner and refer to this approach as spin-CPT. One-particle CPT takes the exact intracluster correlations and treats the intercluster ones approximately, for example, within the Hubbard-I approximation for the Green functions built on the Fermi-type Hubbard operators [47,48,49,50]. A straightforward application of this ideology with respect to Bose-type Hubbard quasiparticles forming the spin excitations leads to an approximation in which intercluster interaction does not enter the Green function calculation procedure, since no spin–spin interaction terms to be taken into account by the generalized Hubbard-I approximation are present in the Hamiltonian (2). Moreover, the intercluster vertex corrections have been shown to have no crucial effect on such cluster calculations [26,28,29]. Thus, the whole procedure is completed in two steps and is identical to the approach implemented in [26,29]. At first, the transverse spin susceptibility is calculated by the ED within the cluster via the following equation:

(13)${\chi }_{\alpha \beta ,ij}^c(\Omega )=\sum _{\mu }\left[\frac{\left.⟨0\right|S_{\alpha ,i}^{+}\left|\mu \right.⟩\left.⟨\mu \right|S_{\beta ,j}^{-}\left|0\right.⟩}{\Omega -\left(E_{\mu }-E_0\right)+\mathrm{i}\eta }-\frac{\left.⟨0\right|S_{\beta ,j}^{-}\left|\mu \right.⟩\left.⟨\mu \right|S_{\alpha ,i}^{+}\left|0\right.⟩}{\Omega +\left(E_{\mu }-E_0\right)+\mathrm{i}\eta }\right],$

where $S^{-}$ and $S^{+}$ are the spin ladder operators, and $\left|0\right.⟩$ and $\left|\mu \right.⟩$ are the ground and excited states, which are obtained in the Lanczos approximations. Then, the translational invariance for the whole lattice is restored analogously to Equation (10):

(14)${\chi }_{\alpha \beta }(\mathbf{k},\Omega )=\frac{1}{N_c}\sum _{i,j}{e}^{-i\left({\mathbf{r}}_i-{\mathbf{r}}_j\right)\mathbf{k}}{\chi }_{\alpha \beta ,ij}^c(\Omega ).$

The dynamic charge susceptibility $\Pi (\mathbf{q},\Omega )$ can also be calculated in a similar manner on a cluster as follows:

(15)${\Pi }_{\alpha ,\beta ,i,j}^c(\Omega )=\sum _{\mu }\left[\frac{\left.⟨0\right|n_{\alpha ,i}\left|\mu \right.⟩\left.⟨\mu \right|n_{\beta ,j}\left|0\right.⟩}{\Omega -\left(E_{\mu }-E_0\right)+\mathrm{i}\eta }-\frac{\left.⟨0\right|n_{\beta ,j}\left|\mu \right.⟩\left.⟨\mu \right|n_{\alpha ,i}\left|0\right.⟩}{\Omega +\left(E_{\mu }-E_0\right)+\mathrm{i}\eta }\right],$

where $n_{\alpha ,i}$ is the number of the particle operator. Since such an approach allows one to calculate the charge correlation function directly, we refer to it as charge-CPT.

3. Results

To determine the changes in the low-energy electronic structure with doping and interactions, we started with the calculation of the Fermi surface for the CuO model (2). Two doping values corresponding to an integer number of electrons per cluster are shown in Figure 2, where we present the spectral function at the Fermi level. We later refer to such plots as the ‘Fermi surface plots’, although they shown neither the Fermi surface nor the Fermi contour itself. However, they are useful, since one can straightforwardly compare them to ARPES plots, which commonly show the distributions of the spectral weight in the narrow energy window near the Fermi level. The momentum–energy distribution of the spectral weight is shown in Figure 3. This distribution represents the band structure corresponding to the Fermi surface plots of Figure 2 for a doping level of $p=1/9$. Note the formation of the Hubbard subbands and the pseudogap at the Fermi level with the increasing interaction. Without Coulomb interactions ($\lambda =0$), the Fermi surface was electron-like. This was a direct consequence of the choice of band structure parameters. This shape did not agree with the ARPES results that showed a hole-like Fermi surface and even Fermi arcs at a low doping level p [6,7]. However, once we included the electronic correlations induced by the Coulomb interaction by increasing the parameter $\lambda$, the Fermi surface changed dramatically at a small p. For $\lambda =0.3$, it was already hole-like. The spectral weight decreased from the nodal to the antinodal direction, which was a signature of the pseudogap. Qualitatively, the obtained Fermi arc distribution of the spectral wight agreed with the results of the dynamic variational Monte Carlo approach [51]. On the contrary, the correlation effects were small at a high doping level. The Lifshitz transition occurred in between the presented values of doping at $p\approx 0.13$. The overall evolution according to the doping amount qualitatively agreed with the ARPES data [6,7] and the results of the mean-field strong coupling approach [52,53], where the appearance of two quantum phase transitions associated with the Fermi surface topology change from small hole pockets to a large hole-like and then electron-like Fermi surface was demonstrated. For a full-strength interaction, $\lambda =1$, the pseudogap was clearly visible at a low doping amount, while the electron-like Fermi surface at $p=2/9$ was affected only slightly. This time, the change in topology appeared at $p\approx 0.17$.

Now, we turn to the CPT-RPA results, which demonstrate how the spin spectrum was affected by the one-particle processes entering the particle–hole bubble ${\chi }^{c\left(0\right)}(\mathbf{k},\Omega )$ through the spectral functions and enhanced by the RPA vertex. The dynamic spin susceptibility under CPT-RPA with a Lorentzian broadening of $\eta =0.1$ is presented in Figure 4 and Figure 5. At a low doping level and relatively weak interaction strength (Figure 4), changes in the electronic structure such as the emergence of the pseudogap resulted in the formation of a response resembling a spin wave spectrum. The lowest energy region was dominated by a contribution of the spectral weight at the antiferromagnetic wave vector due to the presence of the short-range antiferromagnetic correlations dominating the electronic spectrum.

Figure 5 shows the doping dependence of the magnetic spectrum under CPT-RPA for a strong interaction. Here, the Coulomb interaction included in CPT and ${\chi }^{c\left(0\right)}(\mathbf{k},\Omega )$ was at full strength with $\lambda =1$. Parameter ${\lambda }^{\prime }$, which controlled the strength of the interactions entering the RPA vertex, however, was taken to be less than unity to avoid magnetic instability at finite doping, i.e., the divergence of the real part of the zero-frequency RPA spin susceptibility under CPT-RPA. We adjusted ${\lambda }^{\prime }$ so that the system was in the vicinity of instability, which appeared first at the Lifshitz transition, where the particle density at the Fermi level was maximal ($p\approx 0.17$ in the present case). In the pseudogap regime, the response was again somehow similar to the spin wave spectrum; however, it was highly damped. The energy of the $(\pi ,\pi )$ spectral weight maximum tended to zero at the Lifshitz transition point; see Figure 5b. A qualitative change with doping occurred in the overdoped case, where the low-energy spectral weight peaks shifted to the incommensurate wave vectors; see Figure 5c. This can be clearly seen in the corresponding constant-energy slices shown in Figure 6. At a low doping level, the maximal spectral weight was confined near the antiferromagnetic wave vector over a wide energy range. At a high doping level, on the contrary, low-energy excitations emerged at incommensurate wave vectors, forming a cross-like shape with the maximal spectral weight in the antinodal direction. An increase in the excitation energy led to a decrease in the incommensurability.

The doping dependence of the spin susceptibility according to another approach, spin-CPT, is shown in Figure 7. Here, the correlation effects within a cluster were treated explicitly by the ED. The results were in qualitative agreement with the quantum Monte Carlo results for the one-band Hubbard model [26]. At half-filling, the result was substantially different from the picture provided by CPT-RPA. The response function in Figure 7a clearly resembles a spin wave spectrum. However, since only short-range correlations were present within the cluster, the maximal intensity at the antiferromagnetic wave vector shifted to higher frequencies compared to the case of a long-range order. The corresponding constant-energy cuts are shown in Figure 8 with energy values chosen to demonstrate the most prominent features for each doping. At higher energy values, the cuts revealed a dip near the commensurate wave vector $(\pi ,\pi )$, as expected for damped spin-waves and not observed for CPT-RPA. At a low doping level, the results were quite similar, but the spin response peaked less sharply at $(\pi ,\pi )$. The energy evolution of constant-energy cuts in Figure 8b,f,i shows the formation of a cross-like feature at higher energy values with the peaks of the spectral weight shifted from the center to the antinodal direction. In the overdoped regime with $p=2/9$, the spectral peak at $(\pi ,\pi )$ can be observed only at a high frequency in Figure 7c. All low-energy responses were at incommensurate wave vectors and were more pronounced in the antinodal direction; see Figure 8i.

Note that some typical artifacts of the cluster calculations such as size effects were present in the results. These effects should not be taken as physical results. For example, the repeating pattern barely seen in Figure 8i is typical for cluster spectral function calculations. Additionally, the momentum resolution was coarse and provided only qualitative results due to the small cluster size and the absence of intercluster interactions in the computation scheme.

Apart from the spin susceptibility, we also calculated the dynamic charge susceptibility. This was obtained via CPT-RPA and charge-CPT; see Figure 9. For both approaches, when the doping level was increased, the spectral weight appeared at low energies. For CPT-RPA, a feature reminiscent of noninteracting susceptibility was observed, while for charge-CPT, the spectral weight below 2 eV resembled a dispersion law with the bandwidth significantly larger than that of the spin excitations. The small cluster size used herein did not allow us to elucidate complicated charge ordering effects such as those obtained via the quantum Monte Carlo method on large clusters [54]. The doping-dependent evolution of the dynamic charge susceptibility according to charge-CPT was in agreement with the most general tendencies observed using the quantum Monte Carlo approach [55,56], but the strong size effects in our calculations did not allow us to obtain more subtle features than those discussed above.

4. Discussion

Now, we turn to the qualitative comparison of the results obtained in the form of experimental spectra from inelastic neutron scattering (INS) on cuprates. At half-filling, cuprates are antiferromagnetic insulators due to a strong Coulomb interaction, and their spin spectrum has a spin-wave character [57]. Unfortunately, CPT-RPA failed to reproduce a spin-wave spectrum as observed in cuprates at half-filling (the spectrum obtained in this case was uninformative and is not shown). On the other hand, spin-CPT showed signatures of such a spectrum, which could be explained as stemming from the explicit inclusion of the short-range antiferromagnetic correlations within a cluster.

The addition of charge carriers into CuO${}_2$ planes leads to the formation of a more complicated spectrum, whose nature is still under debate [8,9]. Often, the lower downward dispersing and higher upwards dispersing branches of an hourglass spectrum are noted. These are reported to have different temperature dependencies [10] and presumably are of different natures. In the underdoped pseudogap region, several types of behavior have been found in different materials. For example, the spin spectrum of La${}_{2-x}$Sr${}_x$CuO${}_4$ has an hourglass shape in both superconducting and normal states [11], while there are no signatures of the lower downward dispersing branch in the normal state of YBa${}_2$Cu${}_3$O${}_{6+x}$ [12]. For underdoped HgBa${}_2$CuO${}_{4+\delta }$, there is no lower branch in either the superconducting or the normal states [13]; however, the lower branch has been observed near optimal doping in the superconducting state [14]. Both our methods in the underdoped region resulted in the signatures of the upper branch with the lowest excitation energy at $(\pi ,\pi )$ and no presence of low-energy incommensurate excitations. In particular, the absence of the low-energy incommensurate excitations for CPT-RPA, where the susceptibility was calculated from the electronic structure exhibiting pseudogap behaviour, might point to the suppression of such excitations by the pseudogap.

In the overdoped regime, both our methods produced a response at the lowest energy values around the four incommensurate wave vectors located at qualitatively similar positions to those in the INS data for overdoped La${}_{2-x}$Sr${}_x$CuO${}_4$ [58,59]. At high energy values, CPT-RPA presented a picture somewhat similar to that of INS, i.e., a weaker broad feature around $(\pi ,\pi )$. For spin-CPT, the response at and close to the $(\pi ,\pi )$ point appeared only at very high energies, although the general fact of the weakening of the high-energy response close to the $(\pi ,\pi )$ point was similar to experimental observations. A possible explanation for the observed discrepancies was that the RPA part of the CPT-RPA method underestimated the spin–spin correlations in the cuprates and thus could not completely reproduce the spin-wave-like high-energy excitations such as those seen in INS. Spin-CPT, in turn, had coarse momentum and energy resolutions and thus failed to reproduce the fine dispersive features present in the CPT-RPA results.

RIXS experiments have been successful in studying high-energy magnetic excitations, though with some restrictions on the probed area of wave vectors $\mathbf{q}$. Via RIXS, magnon-like excitations in cuprates were detected, dispersing upward to energy values as high as ∼400 meV in a wide doping range from underdoped to heavily overdoped samples [1,2,3,4]. In all the spin spectra presented herein for the wave vectors in the $[0,0.5\pi ]$ range, mainly available in RIXS, an overall magnon-like shape of dispersion could be traced, where the highest response intensity was observed at energy values increasing as the wave vector increased from $(0,0)$ to $(0,\pi )$ and from $(0,0)$ to $(\pi /2,\pi /2)$. The energy width of the response was larger in the first direction than in the second. Momentum-dependent charge excitations were also found in RIXS [5]. These were dispersing upward to an energy value twice the maximal value of spin excitations, which is in general consistent with the charge susceptibility spectra obtained herein.

5. Conclusions

We formulated a cluster perturbation theory for two-particle correlation functions: namely, spin-CPT and charge-CPT for spin and charge dynamic susceptibilities, respectively. Both quantities were calculated by spin/charge-CPT and CPT-RPA in a wide doping range for the effective two-orbital Hubbard model obtained from the Emery $p-d$ model. In the underdoped case, the magnetic response showed signatures of the upper branch of the hourglass dispersion and no presence of low-energy incommensurate spin excitations. Similar spectral weight distributions have been obtained using perturbation theory for strongly-correlated fermions in the normal phase [60,61,62,63]. In the overdoped case, both spin-CPT and CPT-RPA produced low-energy spin excitations located near the four incommensurate wave vectors. These results agreed with the spectra for doped cuprates obtained in INS and RIXS. The calculated momentum-dependent charge excitations were in general agreement with RIXS.

We demonstrated the crucial role of Coulomb interactions for the formation of the Fermi surface shape. A ‘bare’ band structure for the Emery model with the parameters obtained from the ab initio approach resulted in an electron-like Fermi surface, in disagreement with the ARPES data. On the contrary, the electronic correlations induced by Coulomb interactions changed the Fermi surface shape to be hole-like at a low doping level and led to the Fermi arc distribution of the spectral weight (corresponding to the pseudogap) in the underdoped region. At a high doping level, the correlation effects were shown to be small, resulting in an electron-like Fermi surface with the Lifshitz transition occurring at $p\approx 0.13$. The overall changes in the Fermi surface with doping qualitatively agreed with the ARPES data and the results of other mean-field strong coupling approaches.

One may expect that CPT-RPA would overestimate the contribution of the itinerant electrons to a two-particle quantity such as spin susceptibility. However, it is important that it considers the effect of the electronic structure, which is usually obtained quite accurately in CPT. Additionally, the electronic structure affects the spin susceptibility through the particle–hole bubble enhanced by the RPA vertex. Spin-CPT, on the other hand, converged to the exact result with an increase in the system size but underestimated the long-range correlations. In general, CPT-based approaches are simple and economic methods to calculate momentum- and energy-resolved one- and two-particle correlation functions and allow the exact consideration of short-range correlations. The latter play a key role in doped high-$T_c$ cuprates.

Author Contributions

Conceptualization, M.M.K., S.G.O. and V.I.K.; calculations, V.I.K.; formal analysis, S.G.O.; writing, V.I.K., S.V.N. and M.M.K. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data is contained within this article.

Acknowledgments

We acknowledge useful discussions with E.I. Shneyder.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1. The electronic spectral function calculated via CPT using a [Forumla omitted. See PDF.] cluster for the Hamiltonian (2) including (a) and excluding (b) the a orbital. Lorentzian broadening at [Forumla omitted. See PDF.] was used here and below.

Figure 2. The electronic spectral function at the Fermi level for doping with [Forumla omitted. See PDF.] (a,c,e) and [Forumla omitted. See PDF.] (b,d,f) obtained using the Coulomb interaction normalization constants [Forumla omitted. See PDF.] (a,b), [Forumla omitted. See PDF.] (c,d), and [Forumla omitted. See PDF.] (e,f).

Figure 3. Energy–momentum distribution of the electronic spectral function calculated for doping at [Forumla omitted. See PDF.] using the Coulomb interaction normalization constants [Forumla omitted. See PDF.] (a), [Forumla omitted. See PDF.] (b), and [Forumla omitted. See PDF.] (c).

Figure 4. Dynamic spin susceptibility [Forumla omitted. See PDF.] calculated via CPT-RPA with a [Forumla omitted. See PDF.] cluster for doping at [Forumla omitted. See PDF.] and Coulomb interaction renormalization constants [Forumla omitted. See PDF.] (a) and [Forumla omitted. See PDF.] (b).

Figure 5. Dynamic spin susceptibility [Forumla omitted. See PDF.] calculated via CPT-RPA with [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] for doping at [Forumla omitted. See PDF.] (a), [Forumla omitted. See PDF.] (b), and [Forumla omitted. See PDF.] (c).

Figure 6. Constant energy cuts of the dynamical spin susceptibility calculated under CPT-RPA with [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] for doping at [Forumla omitted. See PDF.] (a,c,e) and [Forumla omitted. See PDF.] (b,d,f) and an energy of [Forumla omitted. See PDF.] (corresponding value is shown in each panel).

Figure 7. Dynamic spin susceptibility [Forumla omitted. See PDF.] calculated via spin-CPT for doping at [Forumla omitted. See PDF.] (a), [Forumla omitted. See PDF.] (b), and [Forumla omitted. See PDF.] (c).

Figure 8. Constant energy cuts of the dynamic spin susceptibility under spin-CPT for doping at [Forumla omitted. See PDF.] (a–c), [Forumla omitted. See PDF.] (d–f), and [Forumla omitted. See PDF.] (g–i) and an energy of [Forumla omitted. See PDF.] (corresponding value is shown in each panel).

Figure 9. Dynamic charge susceptibility [Forumla omitted. See PDF.] according to charge-CPT (a,c,e) and CPT-RPA (b,d,f) for doping at [Forumla omitted. See PDF.] (a,b), [Forumla omitted. See PDF.] (c,d), and [Forumla omitted. See PDF.] (e,f).

References

1. Le Tacon, M.; Ghiringhelli, G.; Chaloupka, J.; Sala, M.M.; Hinkov, V.; Haverkort, M.; Minola, M.; Bakr, M.; Zhou, K.; Blanco-Canosa, S. et al. Intense paramagnon excitations in a large family of high-temperature superconductors. Nat. Phys.; 2011; 7, pp. 725-730. [DOI: https://dx.doi.org/10.1038/nphys2041]

2. Le Tacon, M.; Minola, M.; Peets, D.C.; Moretti Sala, M.; Blanco-Canosa, S.; Hinkov, V.; Liang, R.; Bonn, D.A.; Hardy, W.N.; Lin, C.T. et al. Dispersive spin excitations in highly overdoped cuprates revealed by resonant inelastic X-ray scattering. Phys. Rev. B; 2013; 88, 020501. [DOI: https://dx.doi.org/10.1103/PhysRevB.88.020501]

3. Wang, L.; He, G.; Yang, Z.; Garcia-Fernandez, M.; Nag, A.; Zhou, K.; Minola, M.; Tacon, M.L.; Keimer, B.; Peng, Y. et al. Paramagnons and high-temperature superconductivity in a model family of cuprates. Nat. Commun.; 2022; 13, 3163. [DOI: https://dx.doi.org/10.1038/s41467-022-30918-z] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35672416]

4. Robarts, H.C.; Barthélemy, M.; Kummer, K.; García-Fernández, M.; Li, J.; Nag, A.; Walters, A.C.; Zhou, K.J.; Hayden, S.M. Anisotropic damping and wave vector dependent susceptibility of the spin fluctuations in La_2−xSr_xCuO₄ studied by resonant inelastic X-ray scattering. Phys. Rev. B; 2019; 100, 214510. [DOI: https://dx.doi.org/10.1103/PhysRevB.100.214510]

5. Ishii, K.; Tohyama, T.; Asano, S.; Sato, K.; Fujita, M.; Wakimoto, S.; Tustsui, K.; Sota, S.; Miyawaki, J.; Niwa, H. et al. Observation of momentum-dependent charge excitations in hole-doped cuprates using resonant inelastic X-ray scattering at the oxygen K edge. Phys. Rev. B; 2017; 96, 115148. [DOI: https://dx.doi.org/10.1103/PhysRevB.96.115148]

6. Damascelli, A.; Hussain, Z.; Shen, Z.X. Angle-resolved photoemission studies of the cuprate superconductors. Rev. Mod. Phys.; 2003; 75, pp. 473-541. [DOI: https://dx.doi.org/10.1103/RevModPhys.75.473]

7. Sobota, J.A.; He, Y.; Shen, Z.X. Angle-resolved photoemission studies of quantum materials. Rev. Mod. Phys.; 2021; 93, 025006. [DOI: https://dx.doi.org/10.1103/RevModPhys.93.025006]

8. Tranquada, J.M.; Xu, G.; Zaliznyak, I.A. Superconductivity, antiferromagnetism, and neutron scattering. J. Magn. Magn. Mater.; 2014; 350, pp. 148-160. [DOI: https://dx.doi.org/10.1016/j.jmmm.2013.09.029]

9. Fujita, M.; Hiraka, H.; Matsuda, M.; Matsuura, M.; M. Tranquada, J.; Wakimoto, S.; Xu, G.; Yamada, K. Progress in Neutron Scattering Studies of Spin Excitations in High-Tc Cuprates. J. Phys. Soc. Jpn.; 2012; 81, 011007. [DOI: https://dx.doi.org/10.1143/JPSJ.81.011007]

10. Sato, K.; Ikeuchi, K.; Kajimoto, R.; Wakimoto, S.; Arai, M.; Fujita, M. Coexistence of Two Components in Magnetic Excitations of La_2−xSr_xCuO₄ (x = 0.10 and 0.16). J. Phys. Soc. Jpn.; 2020; 89, 114703. [DOI: https://dx.doi.org/10.7566/JPSJ.89.114703]

11. Lipscombe, O.J.; Vignolle, B.; Perring, T.G.; Frost, C.D.; Hayden, S.M. Emergence of Coherent Magnetic Excitations in the High Temperature Underdoped La_2−xSr_xCuO₄ Superconductor at Low Temperatures. Phys. Rev. Lett.; 2009; 102, 167002. [DOI: https://dx.doi.org/10.1103/PhysRevLett.102.167002]

12. Hinkov, V.; Bourges, P.; Pailhes, S.; Sidis, Y.; Ivanov, A.; Frost, C.; Perring, T.; Lin, C.; Chen, D.; Keimer, B. Spin dynamics in the pseudogap state of a high-temperature superconductor. Nat. Phys.; 2007; 3, pp. 780-785. [DOI: https://dx.doi.org/10.1038/nphys720]

13. Chan, M.; Dorow, C.; Mangin-Thro, L.; Tang, Y.; Ge, Y.; Veit, M.; Yu, G.; Zhao, X.; Christianson, A.; Park, J. et al. Commensurate antiferromagnetic excitations as a signature of the pseudogap in the tetragonal high-Tc cuprate HgBa₂CuO_4+δ. Nat. Commun.; 2016; 7, 10819. [DOI: https://dx.doi.org/10.1038/ncomms10819]

14. Chan, M.K.; Tang, Y.; Dorow, C.J.; Jeong, J.; Mangin-Thro, L.; Veit, M.J.; Ge, Y.; Abernathy, D.L.; Sidis, Y.; Bourges, P. et al. Hourglass Dispersion and Resonance of Magnetic Excitations in the Superconducting State of the Single-Layer Cuprate HgBa₂CuO_4+δ Near Optimal Doping. Phys. Rev. Lett.; 2016; 117, 277002. [DOI: https://dx.doi.org/10.1103/PhysRevLett.117.277002]

15. Emery, V.J. Theory of high-T_c superconductivity in oxides. Phys. Rev. Lett.; 1987; 58, pp. 2794-2797. [DOI: https://dx.doi.org/10.1103/PhysRevLett.58.2794]

16. Maier, T.; Jarrell, M.; Pruschke, T.; Hettler, M.H. Quantum cluster theories. Rev. Mod. Phys.; 2005; 77, pp. 1027-1080. [DOI: https://dx.doi.org/10.1103/RevModPhys.77.1027]

17. Kotliar, G.; Savrasov, S.Y.; Haule, K.; Oudovenko, V.S.; Parcollet, O.; Marianetti, C.A. Electronic structure calculations with dynamical mean-field theory. Rev. Mod. Phys.; 2006; 78, pp. 865-951. [DOI: https://dx.doi.org/10.1103/RevModPhys.78.865]

18. Rohringer, G.; Hafermann, H.; Toschi, A.; Katanin, A.A.; Antipov, A.E.; Katsnelson, M.I.; Lichtenstein, A.I.; Rubtsov, A.N.; Held, K. Diagrammatic routes to nonlocal correlations beyond dynamical mean field theory. Rev. Mod. Phys.; 2018; 90, 025003. [DOI: https://dx.doi.org/10.1103/RevModPhys.90.025003]

19. Prokof’ev, N.V.; Svistunov, B.V.; Tupitsyn, I.S. Exact quantum Monte Carlo process for the statistics of discrete systems. J. Exp. Theor. Phys. Lett.; 1996; 64, pp. 911-916. [DOI: https://dx.doi.org/10.1134/1.567243]

20. Kashurnikov, V.A.; Krasavin, A.V. Effective quantum Monte Carlo algorithm for modeling strongly correlated systems. J. Exp. Theor. Phys.; 2007; 105, pp. 69-78. [DOI: https://dx.doi.org/10.1134/S1063776107070163]

21. Gull, E.; Millis, A.J.; Lichtenstein, A.I.; Rubtsov, A.N.; Troyer, M.; Werner, P. Continuous-time Monte Carlo methods for quantum impurity models. Rev. Mod. Phys.; 2011; 83, pp. 349-404. [DOI: https://dx.doi.org/10.1103/RevModPhys.83.349]

22. Spałek, J.; Fidrysiak, M.; Zegrodnik, M.; Biborski, A. Superconductivity in high-Tc and related strongly correlated systems from variational perspective: Beyond mean field theory. Phys. Rep.; 2022; 959, pp. 1-117. [DOI: https://dx.doi.org/10.1016/j.physrep.2022.02.003]

23. Sénéchal, D.; Perez, D.; Pioro-Ladrière, M. Spectral Weight of the Hubbard Model through Cluster Perturbation Theory. Phys. Rev. Lett.; 2000; 84, pp. 522-525. [DOI: https://dx.doi.org/10.1103/PhysRevLett.84.522] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/11015954]

24. Sénéchal, D.; Perez, D.; Plouffe, D. Cluster perturbation theory for Hubbard models. Phys. Rev. B; 2002; 66, 075129. [DOI: https://dx.doi.org/10.1103/PhysRevB.66.075129]

25. Brehm, S.; Arrigoni, E.; Aichhorn, M.; Hanke, W. Theory of two-particle excitations and the magnetic susceptibility in high-Tc cuprate superconductors. Europhys. Lett.; 2010; 89, 27005. [DOI: https://dx.doi.org/10.1209/0295-5075/89/27005]

26. Kung, Y.F.; Bazin, C.; Wohlfeld, K.; Wang, Y.; Chen, C.C.; Jia, C.J.; Johnston, S.; Moritz, B.; Mila, F.; Devereaux, T.P. Numerically exploring the 1D-2D dimensional crossover on spin dynamics in the doped Hubbard model. Phys. Rev. B; 2017; 96, 195106. [DOI: https://dx.doi.org/10.1103/PhysRevB.96.195106]

27. Raum, P.T.; Alvarez, G.; Maier, T.; Scarola, V.W. Two-particle correlation functions in cluster perturbation theory: Hubbard spin susceptibilities. Phys. Rev. B; 2020; 101, 075122. [DOI: https://dx.doi.org/10.1103/PhysRevB.101.075122]

28. Chen, C.C.; van Veenendaal, M.; Devereaux, T.P.; Wohlfeld, K. Fractionalization, entanglement, and separation: Understanding the collective excitations in a spin-orbital chain. Phys. Rev. B; 2015; 91, 165102. [DOI: https://dx.doi.org/10.1103/PhysRevB.91.165102]

29. Pärschke, E.M.; Wang, Y.; Moritz, B.; Devereaux, T.P.; Chen, C.C.; Wohlfeld, K. Numerical investigation of spin excitations in a doped spin chain. Phys. Rev. B; 2019; 99, 205102. [DOI: https://dx.doi.org/10.1103/PhysRevB.99.205102]

30. Nikolaev, S.V.; Korshunov, M.M. Spin and Charge Susceptibilities of the Two-Orbital Model within the Cluster Perturbation Theory for Fe-Based Materials. J. Supercond. Nov. Magn.; 2016; 29, 3093. [DOI: https://dx.doi.org/10.1007/s10948-016-3784-8]

31. Korshunov, M.M.; Gavrichkov, V.A.; Ovchinnikov, S.G.; Nekrasov, I.A.; Pchelkina, Z.V.; Anisimov, V.I. Hybrid LDA and generalized tight-binding method for electronic structure calculations of strongly correlated electron systems. Phys. Rev. B; 2005; 72, 165104. [DOI: https://dx.doi.org/10.1103/PhysRevB.72.165104]

32. Anisimov, V.I.; Kondakov, D.E.; Kozhevnikov, A.V.; Nekrasov, I.A.; Pchelkina, Z.V.; Allen, J.W.; Mo, S.K.; Kim, H.D.; Metcalf, P.; Suga, S. et al. Full orbital calculation scheme for materials with strongly correlated electrons. Phys. Rev. B; 2005; 71, 125119. [DOI: https://dx.doi.org/10.1103/PhysRevB.71.125119]

33. Hybertsen, M.S.; Schlüter, M.; Christensen, N.E. Calculation of Coulomb-interaction parameters for La₂CuO₄ using a constrained-density-functional approach. Phys. Rev. B; 1989; 39, pp. 9028-9041. [DOI: https://dx.doi.org/10.1103/PhysRevB.39.9028]

34. Shastry, B.S. t-J model and nuclear magnetic relaxation in high-T_c materials. Phys. Rev. Lett.; 1989; 63, pp. 1288-1291. [DOI: https://dx.doi.org/10.1103/PhysRevLett.63.1288]

35. Lovtsov, S.; Yushankhai, V. Effective singlet-triplet model for CuO₂ plane in oxide superconductors: The change fluctuation regime. Phys. C Supercond.; 1991; 179, pp. 159-166. [DOI: https://dx.doi.org/10.1016/0921-4534(91)90024-S]

36. Jefferson, J.H.; Eskes, H.; Feiner, L.F. Derivation of a single-band model for CuO₂ planes by a cell-perturbation method. Phys. Rev. B; 1992; 45, pp. 7959-7972. [DOI: https://dx.doi.org/10.1103/PhysRevB.45.7959]

37. Schüttler, H.B.; Fedro, A.J. Copper-oxygen charge excitations and the effective-single-band theory of cuprate superconductors. Phys. Rev. B; 1992; 45, pp. 7588-7591. [DOI: https://dx.doi.org/10.1103/PhysRevB.45.7588]

38. Feiner, L.F.; Jefferson, J.H.; Raimondi, R. Effective single-band models for the high-T_c cuprates. I. Coulomb interactions. Phys. Rev. B; 1996; 53, pp. 8751-8773. [DOI: https://dx.doi.org/10.1103/PhysRevB.53.8751]

39. Raimondi, R.; Jefferson, J.H.; Feiner, L.F. Effective single-band models for the high-T_c cuprates. II. Role of apical oxygen. Phys. Rev. B; 1996; 53, pp. 8774-8788. [DOI: https://dx.doi.org/10.1103/PhysRevB.53.8774]

40. Gavrichkov, V.A.; Ovchinnikov, S.G.; Borisov, A.A.; Goryachev, E.G. Evolution of the band structure of quasiparticles with doping in copper oxides on the basis of a generalized tight-binding method. J. Exp. Theor. Phys.; 2000; 91, pp. 369-383. [DOI: https://dx.doi.org/10.1134/1.1311997]

41. Makarov, I.A.; Shneyder, E.I.; Kozlov, P.A.; Ovchinnikov, S.G. Polaronic approach to strongly correlated electron systems with strong electron-phonon interaction. Phys. Rev. B; 2015; 92, 155143. [DOI: https://dx.doi.org/10.1103/PhysRevB.92.155143]

42. Shneyder, E.I.; Nikolaev, S.V.; Zotova, M.V.; Kaldin, R.A.; Ovchinnikov, S.G. Polaron transformations in the realistic model of the strongly correlated electron system. Phys. Rev. B; 2020; 101, 235114. [DOI: https://dx.doi.org/10.1103/PhysRevB.101.235114]

43. Kung, Y.F.; Chen, C.C.; Wang, Y.; Huang, E.W.; Nowadnick, E.A.; Moritz, B.; Scalettar, R.T.; Johnston, S.; Devereaux, T.P. Characterizing the three-orbital Hubbard model with determinant quantum Monte Carlo. Phys. Rev. B; 2016; 93, 155166. [DOI: https://dx.doi.org/10.1103/PhysRevB.93.155166]

44. Huang, E.W.; Ding, S.; Liu, J.; Wang, Y. Determinantal quantum Monte Carlo solver for cluster perturbation theory. Phys. Rev. Res.; 2022; 4, L042015. [DOI: https://dx.doi.org/10.1103/PhysRevResearch.4.L042015]

45. Graser, S.; Maier, T.A.; Hirschfeld, P.J.; Scalapino, D.J. Near-degeneracy of several pairing channels in multiorbital models for the Fe pnictides. New J. Phys.; 2009; 11, 025016. [DOI: https://dx.doi.org/10.1088/1367-2630/11/2/025016]

46. Korshunov, M.M. Itinerant spin fluctuations in iron-based superconductors. Perturbation Theory: Advances in Research and Applications; Pirogov, Z. Nova Science Publishers Inc.: New York, NY, USA, 2018; Chapter 2 pp. 61-138.

47. Nikolaev, S.V.; Ovchinnikov, S.G. Cluster perturbation theory in Hubbard model exactly taking into account the short-range magnetic order in 2 × 2 cluster. JETP; 2010; 111, pp. 635-644. [DOI: https://dx.doi.org/10.1134/S1063776110100146]

48. Nikolaev, S.V.; Ovchinnikov, S.G. Effect of hole doping on the electronic structure and the Fermi surface in the Hubbard model within norm-conserving cluster pertubation theory. J. Exp. Theor. Phys.; 2012; 114, pp. 118-131. [DOI: https://dx.doi.org/10.1134/S1063776111150143]

49. Kuz’min, V.I.; Nikolaev, S.V.; Ovchinnikov, S.G. Comparison of the electronic structure of the Hubbard and t − J models within the cluster perturbation theory. Phys. Rev. B; 2014; 90, 245104. [DOI: https://dx.doi.org/10.1103/PhysRevB.90.245104]

50. Kuz’min, V.I.; Visotin, M.A.; Nikolaev, S.V.; Ovchinnikov, S.G. Doping and temperature evolution of pseudogap and spin-spin correlations in the two-dimensional Hubbard model. Phys. Rev. B; 2020; 101, 115141. [DOI: https://dx.doi.org/10.1103/PhysRevB.101.115141]

51. Rosenberg, P.; Sénéchal, D.; Tremblay, A.M.S.; Charlebois, M. Fermi arcs from dynamical variational Monte Carlo. Phys. Rev. B; 2022; 106, 245132. [DOI: https://dx.doi.org/10.1103/PhysRevB.106.245132]

52. Korshunov, M.M.; Ovchinnikov, S.G. Doping-dependent evolution of low-energy excitations and quantum phase transitions within an effective model for high-Tc copper oxides. Eur. Phys. J. B; 2007; 57, pp. 271-278. [DOI: https://dx.doi.org/10.1140/epjb/e2007-00179-2]

53. Ovchinnikov, S.G.; Shneyder, E.I.; Korshunov, M.M. From underdoped to overdoped cuprates: Two quantum phase transitions. J. Physics Condens. Matter.; 2011; 23, 045701. [DOI: https://dx.doi.org/10.1088/0953-8984/23/4/045701]

54. Mai, P.; Nichols, N.S.; Karakuzu, S.; Bao, F.; Del Maestro, A.; Maier, T.A.; Johnston, S. Robust charge-density wave correlations in the electron-doped single-band Hubbard model. arXiv; 2022; arXiv: 2210.14930[DOI: https://dx.doi.org/10.1038/s41467-023-38566-7]

55. Gröber, C.; Eder, R.; Hanke, W. Anomalous low-doping phase of the Hubbard model. Phys. Rev. B; 2000; 62, pp. 4336-4352. [DOI: https://dx.doi.org/10.1103/PhysRevB.62.4336]

56. Dong, X.; Chen, X.; Gull, E. Dynamical charge susceptibility in the Hubbard model. Phys. Rev. B; 2019; 100, 235107. [DOI: https://dx.doi.org/10.1103/PhysRevB.100.235107]

57. Coldea, R.; Hayden, S.M.; Aeppli, G.; Perring, T.G.; Frost, C.D.; Mason, T.E.; Cheong, S.W.; Fisk, Z. Spin Waves and Electronic Interactions in La₂CuO₄. Phys. Rev. Lett.; 2001; 86, pp. 5377-5380. [DOI: https://dx.doi.org/10.1103/PhysRevLett.86.5377]

58. Wakimoto, S.; Yamada, K.; Tranquada, J.M.; Frost, C.D.; Birgeneau, R.J.; Zhang, H. Disappearance of Antiferromagnetic Spin Excitations in Overdoped La_2−xSr_xCuO₄. Phys. Rev. Lett.; 2007; 98, 247003. [DOI: https://dx.doi.org/10.1103/PhysRevLett.98.247003]

59. Lipscombe, O.J.; Hayden, S.M.; Vignolle, B.; McMorrow, D.F.; Perring, T.G. Persistence of High-Frequency Spin Fluctuations in Overdoped Superconducting La_2−xSr_xCuO₄ (x = 0.22). Phys. Rev. Lett.; 2007; 99, 067002. [DOI: https://dx.doi.org/10.1103/PhysRevLett.99.067002]

60. Eremin, M.V.; Aleev, A.A.; Eremin, I.M. On interlay between the magnetic susceptibilities of localized and itinerant electrons in hole-doped HTSCs. JETP Lett.; 2006; 84, pp. 167-170. [DOI: https://dx.doi.org/10.1134/S0021364006150148]

61. Eremin, M.V.; Aleev, A.A.; Eremin, I.M. Dynamic spin susceptibility of hole-doped high-temperature superconductors in a singlet-correlated conduction band model. J. Exp. Theor. Phys.; 2008; 106, pp. 752-764. [DOI: https://dx.doi.org/10.1134/S1063776108040146]

62. Eremin, M.V.; Shigapov, I.M.; Eremin, I.M. Dual features of magnetic susceptibility in superconducting cuprates: A comparison to inelastic neutron scattering. Eur. Phys. J. B; 2012; 85, 131. [DOI: https://dx.doi.org/10.1140/epjb/e2012-20539-y]

63. Eremin, M.V.; Shigapov, I.M.; Thuy, H.T.D. Collective spin excitations in the singlet-correlated band model: A comparison with resonant inelastic X-ray scattering. J. Phys. Condens. Matter.; 2013; 25, 345701. [DOI: https://dx.doi.org/10.1088/0953-8984/25/34/345701] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23896734]