Multistage Electronic Nematic Transitions in Cuprate Superconductors: Functional-Renormalization-Group Analysis

TL;DR

This paper uses a functional-renormalization-group approach to analyze complex multistage nematic phase transitions in cuprate superconductors, revealing uniform and density-wave nematic orders driven by spin fluctuations.

Contribution

It uncovers the dominant uniform nematic instability and subsequent density-wave order in cuprates without assuming specific wavevectors, highlighting multistage transitions driven by fluctuation processes.

Findings

Uniform ($q=0$) d-symmetry charge modulation is dominant.

A nematic p-orbital density wave appears at $Q_a = (rac{ ext{pi}}{2}, 0)$.

Multistage nematic transitions are driven by fluctuation-exchange processes.

Abstract

Recently, complex phase transitions accompanied by the rotational symmetry breaking have been discovered experimentally in cuprate superconductors. To find the realized order parameters, we study various charge susceptibilities in an unbiased way, by applying the functional-renormalization-group method to the realistic $d$-$p$ Hubbard model. Without assuming the wavevector of the order parameter, we reveal that the most dominant instability is the uniform ($q = 0$) charge modulation on the $p_x$ and $p_y$ orbitals, which possesses the d-symmetry. This uniform nematic order triggers another nematic p-orbital density wave along the axial (Cu-Cu) direction at $Q_a = (\pi/2,0)$. It is predicted that uniform nematic order is driven by the spin fluctuations in the pseudogap region, and another nematic density-wave order at $q = Q_a$ is triggered by the uniform order. The predicted multistage nematic transitions are caused by the Aslamazov-Larkin-type fluctuation-exchange processes.

Tables

Table 1. Table 1: Hopping integrals for the N = 0 𝑁 0 N=0 and N = 1 𝑁 1 N=1 models given in Ref. Hansmann  et al. ( 2014 ) . The unit is eV.

NMTO	${ϵ}_d-{ϵ}_p$	$t_{dd}$	$t_{pd}$	$t_{pd}^{\prime }$	$t_{pp}$	$t_{pp}^{\prime }$	$t_{pp}^{\prime \prime }$	$t_{pp}^{\prime \prime \prime }$
$N=0$	$0.43$	$-0.10$	$0.96$	$-0.10$	$0.15$	$-0.24$	$0.02$	$0.11$
$N=1$	$0.95$	$0.15$	$1.48$	$0.08$	$0.91$	$0.03$	$0.15$	$0.03$

Equations

$${\chi }^{\mathrm{spin}}(\mathbf{q})=\frac{1}{2}{\int }_0^{1/T}d\tau ⟨S_d(\mathbf{q},\tau )S_d(-\mathbf{q},0)⟩,$$

$${\chi }_d^{p\text{mbox}-\mathrm{orb}}(\mathbf{q})=\frac{1}{2}{\int }_0^{1/T}d\tau ⟨n_d^{p\text{mbox}-\mathrm{orb}}(\mathbf{q},\tau )n_d^{p\text{mbox}-\mathrm{orb}}(-\mathbf{q},0)⟩,$$

$${\chi }^{d\text{mbox}-\mathrm{orb}}(\mathbf{q})=\frac{1}{2}{\int }_0^{1/T}d\tau ⟨B(\mathbf{q},\tau )B(-\mathbf{q},0)⟩,$$

$$H_U=\frac{1}{4}\sum _i\sum _{\sigma {\sigma }^{\prime }\rho {\rho }^{\prime }}{U}^{0;\sigma {\sigma }^{\prime }\rho {\rho }^{\prime }}{d}_{i\sigma }^{†}{d}_{i{\sigma }^{\prime }}{d}_{i{\rho }^{\prime }}^{†}{d}_{i\rho },$$

$$U^{0;\sigma {\sigma }^{\prime }\rho {\rho }^{\prime }}=\frac{1}{2}{U}^{0;s}{\vec{\sigma }}_{\sigma {\sigma }^{\prime }}\cdot {\vec{\sigma }}_{{\rho }^{\prime }\rho }+\frac{1}{2}{U}^{0;c}{\delta }_{\sigma ,{\sigma }^{\prime }}{\delta }_{{\rho }^{\prime },\rho },$$

$${\Gamma }^{\sigma {\sigma }^{\prime }\rho {\rho }^{\prime }}(\mathbf{k}+\mathbf{q},\mathbf{k};{\mathbf{k}}^{\prime }+\mathbf{q},{\mathbf{k}}^{\prime })$$

$$=\frac{1}{2}{\Gamma }^s(\mathbf{k}+\mathbf{q},\mathbf{k};{\mathbf{k}}^{\prime }+\mathbf{q},{\mathbf{k}}^{\prime }){\vec{\sigma }}_{\sigma {\sigma }^{\prime }}\cdot {\vec{\sigma }}_{{\rho }^{\prime }\rho }$$

$$+\frac{1}{2}{\Gamma }^c(\mathbf{k}+\mathbf{q},\mathbf{k};{\mathbf{k}}^{\prime }+\mathbf{q},{\mathbf{k}}^{\prime }){\delta }_{\sigma ,{\sigma }^{\prime }}{\delta }_{{\rho }^{\prime },\rho },$$

$$\frac{d}{d\Lambda }{\Gamma }_{\mathrm{RG}}(k_1,k_2;k_3,k_4)=-\frac{T}{N}\sum _{k,k^{\prime }}\left[\frac{d}{d\Lambda }G(k)G(k^{\prime })\right]$$

$$\displaystyle\ \ \ \ \ \ \times\Bigl{[}\Gamma_{\rm RG}(k_{1},k_{2};k,k^{\prime})\,\Gamma_{\rm RG}(k,k^{\prime};k_{3},k_{4})$$

$$-{\Gamma }_{\mathrm{RG}}(k_1,k_3;k,k^{\prime }){\Gamma }_{\mathrm{RG}}(k,k^{\prime };k_2,k_4)$$

$$\displaystyle\ \ \ \ \ \ -\frac{1}{2}\Gamma_{\rm RG}(k_{1},k;k^{\prime},k_{4})\,\Gamma_{\rm RG}(k,k_{2};k_{3},k^{\prime})\Bigr{]},$$

$$\Delta \Sigma (k)=(1-P_{{\mathrm{A}}_{1g}})T\sum _{q}V(q)G(k+q),$$

$${\lambda }_{\mathbf{q}}\Delta {\Sigma }_{\mathbf{q}}(k)=T\sum _{k^{\prime }}K(\mathbf{q};k,k^{\prime })\Delta {\Sigma }_{\mathbf{q}}(k^{\prime }),$$

$$K(\mathbf{q};k,k^{\prime })=\left(\frac{3}{2}{V}_0^s(k-k^{\prime })+\frac{1}{2}{V}_0^c(k-k^{\prime })\right)$$

$$\times G_0(k^{\prime }+\mathbf{q}/2)G_0(k^{\prime }-\mathbf{q}/2)$$

$$-T\sum _p(\frac{3}{2}{V}_0^s(p+\mathbf{q}/2)V_0^s(p-\mathbf{q}/2)$$

$$+\frac{1}{2}{V}_0^c(p+\mathbf{q}/2)V_0^c(p-\mathbf{q}/2))$$

$$\times G_0(k-p)\left({\Lambda }_{\mathbf{q}}(k^{\prime };p)+{\Lambda }_{\mathbf{q}}(k^{\prime };-p)\right),$$

Full text

I Introduction

In the normal state of high-$T_{c}$ cuprate superconductors, interesting unconventional order parameters emerge due to the strong interference among the spin, charge, and orbital degrees of freedom. These phenomena should be directly related to the fundamental electronic states in the pseudogap region. The emergence of the charge-density-wave (CDW) states inside the pseudogap region has been confirmed by the x-ray and STM measurements Ghiringhelli et al. (2012); Chang et al. (2012); Fujita et al. (2014); Wu et al. (2015); Hamidian et al. (2015); Comin and Damascelli (2016), as schematically shown in Fig. 1(a). The observed CDW pattern is shown in Fig. 1(b), in which the density modulations mainly occur on the oxygen $p_{x}$ and $p_{y}$ orbitals with antiphase ($d$-symmetry) form factor. The discovery of the CDW has promoted significant progress in the theoretical studies, such as the spin-fluctuation-driven density-wave scenarios Davis and Lee (2013); Metlitski and Sachdev (2010); Husemann and Metzner (2012); Efetov et al. (2013); Sachdev and La Placa (2013); Mishra and Norman (2015); Yamakawa and Kontani (2015); Orth:2017 and the superconducting-fluctuation scenarios Berg et al. (2009); Fradkin et al. (2015); Wang et al. (2015); Lee (2014).

The origin and nature of the pseudogap phase below $T^{*}$ remain unsolved. For example, it is unclear whether the pseudogap is a distinct phase or a continuous crossover. The short-range spin-fluctuations at $T\sim T^{*}$ induce the large quasiparticle damping TPSC ; Kotliar ; Maier , which causes the pseudogap in the density-of-states. On the other hand, the phase transition around $T^{*}$ have been reported by the resonant ultrasound spectroscopy Shekhter et al. (2013), ARPES analysis ARPES-Science2011 , and magnetic torque measurement Matsuda:2016 . Especially, Ref. Matsuda:2016 discovered the $C_{4}$ symmetry breaking (nematic) transition, and its natural candidate is the uniform CDW with $d$-symmetry schematically shown in Fig. 1(c). Then, a fundamental question is what mechanism can account for such unconventional multistage CDW transitions. No CDW instabilities are given by the mean-field-level approximations, like the random-phase-approximation (RPA), unless large inter-site interactions are introduced Bulut ; Yamakawa and Kontani (2015). Therefore, higher-order many body effects, called the vertex corrections (VCs), should be essential for the CDW formation Onari and Kontani (2012); Davis and Lee (2013); Metlitski and Sachdev (2010); Husemann and Metzner (2012); Efetov et al. (2013); Sachdev and La Placa (2013); Mishra and Norman (2015); Yamakawa and Kontani (2015).

In many spin-fluctuation-driven CDW scenarios, the CDW wavevector is given by the minor nesting vector $\bm{Q}_{\mathrm{a}}$ or $\bm{Q}_{\mathrm{d}}$ in Fig. 1(d); $\bm{Q}_{\mathrm{a}}$ is the “axial-wavevector” parallel to the nearest Cu-Cu direction, and $\bm{Q}_{\mathrm{d}}$ is the “diagonal-wavevector”. The experimental axial CDW is obtained if the Aslamazov-Larkin VCs (AL-VCs) are taken into account Yamakawa and Kontani (2015). In addition, the uniform (${{\bm{q}}}={\bm{0}}$) CDW instability has been studied intensively based on the Hubbard models Halboth and Metzner (2000); Honerkamp (2005); Husemann and Metzner (2012); Su:2011co ; Kitatani:2016ur . In these studies, however, it was difficult to exclude the possibility that the CDW susceptibility has the maximum at finite ${{\bm{q}}}$.

Theoretically, it is difficult to analyze the spin and charge susceptibilities with general wavevector $\bm{q}$ on equal footing, by including the VCs in an unbiased way. For this purpose, in principle, the functional renormalization-group (fRG) method would be the best theoretical method. The pioneering fRG studies Honerkamp (2005); Husemann and Metzner (2012) were performed only in the weak-coupling region, so the obtained CDW instability is small and its ${{\bm{q}}}$-dependence is not clear. In order to overcome this problem, we have to improve the numerical accuracy of the fRG method, and apply it to the two-dimensional Hubbard model in the strong-coupling region.

In this paper, we study the orbital-dependent spin and charge susceptibilities for various symmetries on equal footing, by analyzing the higher-order VCs in an unbiased way using the improved fRG method. We find that the uniform CDW accompanied by the $p$-orbital polarization ($n_{x}\neq n_{y}$), shown in Fig. 1(c), is driven by the antiferro spin fluctuations. In this uniform nematic CDW phase, another nematic CDW instability emerges at the wavevector ${{\bm{q}}}={{\bm{Q}}}_{\rm a}$ as shown in Fig. 1(b). The present study indicates that the uniform $p$-orbital polarization appears in the pseudogap region, and the axial ${{\bm{q}}}={{\bm{Q}}}_{\rm a}$ CDW is induced at $T_{\rm CDW}<T^{*}$. These multistage CDW transitions in under-doped cuprates originate from the higher-order AL-type VCs.

In the present study, we use the functional RG + constrained RPA (RG+cRPA) method. The advantage of this method had been explained in Refs. Tsuchiizu et al. (2013, 2015, 2016); Tazai et al. (2016) and Appendix A in detail.

II Model and Theoretical Method

Here, we study a standard three-orbital $d$-$p$ Hubbard model Yamakawa and Kontani (2015); Tsuchiizu et al. (2016); Hansmann et al. (2014) expressed as $H=\sum_{\bm{k},\sigma}\bm{c}_{\bm{k},\sigma}^{\dagger}\,\hat{h}_{0}(\bm{k})\,\bm{c}_{\bm{k},\sigma}+U\sum_{\bm{j}}n_{d,\bm{j},\uparrow}n_{d,\bm{j},\downarrow}$, where $\bm{c}_{\bm{k},\sigma}^{\dagger}=(d_{\bm{k},\sigma}^{\dagger},p_{x,\bm{k},\sigma}^{\dagger},p_{y,\bm{k},\sigma}^{\dagger})$ is the creation operator for the electron on $d$, $p_{x}$, and $p_{y}$ orbitals, and $\hat{h}_{0}(\bm{k})$ is the kinetic term given as the 0MTO model in Refs. Hansmann et al. (2014); note . (The numerical results are unchanged if another realistic 1MTO model is used; see Appendix B.) $U$ is the Hubbard-type on-site Coulomb interaction for the $d$ orbital, and $n_{d,\bm{j},\sigma}=d^{\dagger}_{\bm{j},\sigma}d_{\bm{j},\sigma}$ at site $\bm{j}$. Hereafter, we study the $10\%$ hole doping case. The Fermi surface (FS) and the band structure are shown in Figs. 1(d) and (e), respectively.

By using the RG+cRPA theory in Ref. Tsuchiizu et al. (2016), we find that the spin susceptibility for $d$-electrons,

[TABLE]

and the $B_{1g}$-symmetry ($d$-symmetry) charge susceptibility for $p$-electrons,

[TABLE]

are the most enhanced susceptibilities Tsuchiizu et al. (2016). Here, $S_{d}(\bm{q},\tau)$ is the $d$-electron spin operator, and $n^{p\mbox{-}\mathrm{orb}}_{d}(\bm{q})\equiv n_{x}(\bm{q})-n_{y}(\bm{q})$ ($n_{x(y)}(\bm{q})=\sum_{\bm{k},\sigma}p_{x(y),\bm{k},\sigma}^{\dagger}p_{x(y),\bm{k+q},\sigma}$) is the $p$-orbital charge-density-wave ($p$O-CDW) operator with $B_{1g}$ symmetry. If $\chi^{p\mbox{-}\mathrm{orb}}_{d}(\bm{q})$ diverges at $\bm{q}=\bm{Q}_{\rm a}$ [$\bm{q}=\bm{0}$], the $p$O-CDW order shown in Fig. 1(b) [Fig. 1(c)], which is the CDW on $p$-orbitals, is realized. We verified that the charge susceptibilities with non-$B_{1g}$-symmetries, such as the $A_{1g}$-symmetry total charge susceptibility for $n\equiv n_{d}+n_{x}+n_{y}$, remains small even in the strong-coupling region Tsuchiizu et al. (2016).

In the RG+cRPA method, we calculate the scattering processes involving higher-energy states $|E_{\bm{k},\nu}|>\Lambda_{0}$ [$\nu$ being the band index; see Fig. 1(e)] using the RPA with the energy-constraint, and incorporate their contributions into the initial vertex functions of the RG equations Tsuchiizu et al. (2013, 2015, 2016); Tazai et al. (2016). Using the RPA, the higher-energy processes are calculated accurately by dropping the VCs, which are less important for higher-energy processes. The lower-energy scattering processes for $|E_{{{\bm{k}}},\nu}|<\Lambda_{0}$ are calculated by solving the RG equations, based on the $N_{p}$-patch RG scheme Halboth and Metzner (2000); Metzner et al. (2012). Hereafter, we put $N_{p}=128$ and $\Lambda_{0}=0.5$ eV. In the RG+cRPA method, the numerical accuracy of the susceptibilities is greatly improved even in the weak-coupling region since the cRPA is used for the higher-energy processes, for which the $N_{p}$-patch RG scheme is less accurate. We verified that the numerical results are essentially independent of the choice of $\Lambda_{0}$ when $E_{\rm F}\gtrsim\Lambda_{0}\gg T$.

By solving the RG equations, many-body vertices are gradually renormalized as reducing the energy scale $\Lambda_{l}=\Lambda_{0}e^{-l}$ with increasing $l\ (\geq 0)$. In principle, the renormalization of the vertex saturates when $\Lambda_{l}$ reaches $\sim T$ Metzner et al. (2012); Bourbonnais et al. (2004). Here, we introduce the lower-energy cutoff $\Lambda_{\rm low}\ (\sim T)$ in the RG equations for the four-point vertex $\Gamma_{l}^{s(c)}$, and stop the renormalization at $\Lambda_{l}=\Lambda_{\rm low}$; see Appendix A and Ref. Husemann and Metzner (2012). (We do not introduce the lower-energy cutoff in the RG equations for $\chi^{s,c}({{\bm{q}}})$.) In the previous study Tsuchiizu et al. (2016), we set a large cutoff $\Lambda_{\rm low}=\pi T$ to achieve stable numerical results. When $\Lambda_{\rm low}\gg T$, the uniform (${{\bm{q}}}={\bm{0}}$) nematic susceptibility is especially underestimated compared to ${{\bm{q}}}\neq{\bm{0}}$ instabilities, as we will discuss later. Since we have improved the numerical accuracy in solving the RG equations, we can use a smaller natural cutoff $\Lambda_{\rm low}=T$. For this reason, we can obtain the ${{\bm{q}}}$-dependence of the susceptibility accurately, including ${{\bm{q}}}\approx{\bm{0}}$.

We find that the numerical accuracy and stability are improved by employing the Wick-ordered scheme of the fRG formalism, in which the cutoff function $\Theta_{<}^{\Lambda}({\epsilon})=\Theta(\Lambda-|{\epsilon}|)$ is used for the Green function Metzner et al. (2012). In this scheme, in principle, the VCs due to the higher-energy processes are included more accurately, compared to using another cutoff function $\Theta_{>}^{\Lambda}({\epsilon})=\Theta(|{\epsilon}|-\Lambda)$ based on the Kadanoff-Wilson scheme used in Ref. Tsuchiizu et al. (2016).

III Multistage Electronic Nematic Transitions

In Figs. 2(a) and (b), we show the $p$O-CDW susceptibility $\chi^{p\mbox{-}\mathrm{orb}}_{d}(\bm{q})$ given by the RG+cRPA method for $U=4.32$ eV at $T=0.1$ eV. The obtained large peaks at $\bm{q}=\bm{0}$, $\bm{Q}_{\mathrm{a}}$, and $\bm{Q}_{\mathrm{d}}$ originate from the VCs, since the RPA result is less singular as seen in Fig. 2(b). As shown in Figs. 2(a)-2(c), the most dominant peak locates at $\bm{q}=\bm{0}$. This is consistent with the experimental uniform nematic transition at $T^{*}\ (>T_{\rm CDW})$ Matsuda:2016 . We also obtain the peak structures at $\bm{q}=\bm{Q}_{\mathrm{a}}$ and $\bm{Q}_{\mathrm{d}}$, consistently with our previous fRG study Tsuchiizu et al. (2016). Figure 2(c) shows that $\chi^{p\mbox{-}\mathrm{orb}}_{d}(\bm{0})$ monotonically increases with decreasing $T$, consistently with the recent electronic nematic susceptibility measurement Shibauchi-nem . At low temperatures, $\chi^{p\mbox{-}\mathrm{orb}}_{d}(\bm{Q}_{\mathrm{a}})$ increases steeply and becomes larger than $\chi^{p\mbox{-}\mathrm{orb}}_{d}(\bm{Q}_{\mathrm{d}})$, shown in the inset of Fig. 2(c). Note that the temperature $T=0.1$ eV is comparable to $T^{*}\sim 300$ K if the mass-enhancement factor $m^{*}/m_{\rm band}\sim 3$ is taken into account.

The enhancement of $\chi^{p\mbox{-}\mathrm{orb}}_{d}(\bm{q})$ is caused by the spin fluctuations, due to the strong charge-spin interplay given by the VCs. The moderate peak at $\bm{Q}_{\mathrm{d}}$ is caused by the Maki-Thompson (MT)-VCs, given by the series of the single-fluctuation-exchange processes shown in Fig. 2(d) Sachdev and La Placa (2013); Mishra and Norman (2015). However, the MT-VCs cannot account for the dominant peaks at $\bm{q}=\bm{0}$ and $\bm{Q}_{\mathrm{a}}$. Recently, it was found that the uniform nematic order in the Fe-based superconductors Onari and Kontani (2012); Kontani and Yamakawa (2014) and Sr3Ru2O7 Ohno et al. (2013); Tsuchiizu et al. (2013) is driven by the AL-VC, given by the series of the double-fluctuation-exchange processes shown in Fig. 2(e). In fact, the first term in Fig. 2(e) is proportional to $\sum_{{\bm{k}}}\chi^{\rm spin}({{\bm{k}}}+{{\bm{q}}})\chi^{\rm spin}({{\bm{k}}})$, which takes large value for ${{\bm{q}}}={\bm{0}}$ when $\chi^{\rm spin}_{\rm max}\gg 1$ Onari and Kontani (2012); Grilli . Later, we will demonstrate that the AL-VC causes the uniform and axial CDW instabilities in the present $d$-$p$ model.

Next, we investigate the $U$-dependences of the susceptibilities. In the inset of Fig. 3(a), we show the $U$ dependence of $\chi^{\mathrm{spin}}_{\rm max}\equiv{\rm max}_{{\bm{q}}}\{\chi^{\mathrm{spin}}({{\bm{q}}})\}$. Thanks to the numerical accuracy of the RG+cRPA method, $\chi^{\mathrm{spin}}_{\mathrm{max}}$ perfectly follows the RPA result for wide weak-coupling region ($U<4$ eV). To clarify the close interplay between spin and orbital fluctuations, we plot the peak values of $\chi^{p\mbox{-}\mathrm{orb}}_{d}(\bm{q})$ as a function of $\chi^{\mathrm{spin}}_{\rm max}$ in Fig. 3(a). In contrast to $\chi^{\rm spin}_{\rm max}$, $\chi^{p\mbox{-}\mathrm{orb}}_{d}(\bm{q})$ strongly deviates from the RPA result, indicating the significance of the VCs. With increasing $U$, the peak position of $\chi^{p\mbox{-}\mathrm{orb}}_{d}(\bm{q})$ shifts to ${{\bm{q}}}={\bm{0}}$ at $\chi^{\mathrm{spin}}_{\rm max}\sim 2.5$, and $\chi^{p\mbox{-}\mathrm{orb}}_{d}(\bm{0})$ exceeds the spin susceptibility for $\chi^{\mathrm{spin}}_{\rm max}\gtrsim 10$ eV*-1*.

To understand the origin of the enhancement of $\chi^{p\mbox{-}\mathrm{orb}}_{d}(\bm{q})$, we analyze the scaling flow of the effective interaction for the $p$O-CDW introduced as $\Gamma_{d}^{p\mbox{-}\mathrm{orb}}(\bm{q})\equiv\Gamma^{c}_{x;x}(\bm{q})+\Gamma^{c}_{y;y}(\bm{q})-\Gamma^{c}_{x;y}(\bm{q})-\Gamma^{c}_{y;x}(\bm{q})$ with $\Gamma^{c}_{{\alpha};{\beta}}(\bm{q})\equiv\sum_{\bm{k},\bm{k}^{\prime}}\Gamma^{c}_{l}(\bm{k}+\bm{q},\bm{k};\bm{k}^{\prime}+\bm{q},\bm{k}^{\prime})\,u^{*}_{\alpha}(\bm{k}+{{\bm{q}}})\,u_{\alpha}(\bm{k})\,\cdot u_{\beta}(\bm{k}^{\prime}+{{\bm{q}}})\,u_{\beta}^{*}(\bm{k}^{\prime})$. Here $\Gamma^{c}_{l}$ is the charge-channel four-point vertex, which is a moderate function of the Fermi momenta in the parameter range of the present numerical study. $u_{\alpha}(\bm{k})$ is the matrix element connecting the $p$-orbitals (${\alpha}=x,y$) and the conduction band Tsuchiizu et al. (2016). The scaling flow of $\Gamma_{d}^{p\mbox{-}\mathrm{orb}}(\bm{q})$ is shown in Fig. 3(b), with the scaling parameter $l=\ln(\Lambda_{0}/\Lambda_{l})$. The negative effective interaction drives the enhancement of the corresponding instability. We also plot the effective interaction for the spin channel, $\Gamma^{\mathrm{spin}}(\bm{Q}_{\mathrm{s}})$. For the spin-channel, $\Gamma^{\mathrm{spin}}(\bm{Q}_{\mathrm{s}})\sim-U$ at $l=0$, and it is renormalized like the RPA as $\Gamma^{\mathrm{spin}}_{l}=\Gamma^{\mathrm{spin}}_{0}/(1-c|\Gamma^{\mathrm{spin}}_{0}|l)$ for $l\lesssim\ln(\Lambda_{0}/T)=1.6$, where $c$ is the density of states. For the charge-channel, although $\Gamma_{d}^{p\mbox{-}\mathrm{orb}}(\bm{q})$ at $l=0$ is quite small, it is strongly renormalized to be a large negative value. This result means that the CDW instability originates from the VC going beyond the RPA.

We also calculate the $d$-electron charge susceptibility with form factor $f_{{\bm{q}}}(\bm{k})$, which is given as

[TABLE]

where $B(\bm{q})=\sum_{\bm{k},\sigma}f_{{\bm{q}}}(\bm{k})\,d_{\bm{k}-{{\bm{q}}}/2,\sigma}^{\dagger}d_{\bm{k}+{{\bm{q}}}/2,\sigma}$. The numerically optimized $f_{{\bm{q}}}(\bm{k})$ at ${{\bm{q}}}=\bm{0}$ is shown in Fig. 3(c), which has the $B_{1g}$-symmetry. Its Fourier transformation gives the modulation of the effective hopping integrals, called the $d_{x^{2}-y^{2}}$-wave bond-order. Since ${{\bm{k}}}$-dependence of $f_{\bm{0}}(\bm{k})$ in Fig. 3(c) is similar to that of $|u_{x}({{\bm{k}}})|^{2}-|u_{y}({{\bm{k}}})|^{2}$, the obtained $\chi^{\mathrm{spin}}_{\mathrm{max}}$ dependence of $\chi^{d\mbox{-}\mathrm{orb}}({\bm{0}})$ with the optimized form factor is similar to that of $\chi^{p\mbox{-}\mathrm{orb}}_{d}(\bm{0})$ shown in Fig. 3(a). In the Appendix C, we analyze the single $d$-orbital Hubbard model, and find the strong enhancement of $\chi^{d\mbox{-}\mathrm{orb}}(\bm{q})$ with $B_{1g}$ form factor at ${{\bm{q}}}={\rm 0}$, ${{\bm{Q}}}_{\rm a}$ and ${{\bm{Q}}}_{\rm d}$, very similarly to the $p$O-CDW susceptibilities shown in Figs. 2 and 3.

As shown in Fig. 2 (c), $\chi^{p\mbox{-}\mathrm{orb}}_{d}({\bm{0}})$ increases divergently at $T\sim 0.1$ eV, and the uniform $p$-orbital polarization with $n_{x}\neq n_{y}$ depicted in Fig. 1(c) appears below the transition temperature. In order to discuss the CDW instabilities inside the nematic phase, we perform the RG+cRPA analysis in the presence of the uniform $p$O-CDW order $H^{\prime}=-\frac{1}{2}\Delta E_{p}[n_{x}(\bm{0})-n_{y}(\bm{0})]$. In Fig. 4(a), we plot the peak values of $\chi^{p\mbox{-}\mathrm{orb}}_{d}(\bm{q})$ in the uniform $p$O-CDW state with $\Delta E_{p}=0.3$ eV. Due to small $\Delta E_{p}>0$, $\chi^{p\mbox{-}\mathrm{orb}}_{d}(\bm{q})$ at $\bm{q}=\bm{Q}_{\mathrm{a}}^{x}$ (along $x$-axis) strongly increases whereas that at $\bm{q}=\bm{Q}_{\mathrm{a}}^{y}$ (along $y$-axis) decreases. Thus, the $p$O-CDW at $\bm{q}=\bm{Q}_{\mathrm{a}}^{x}$ is expected to emerge below $T_{\mathrm{\rm CDW}}$, consistently with the phase diagram in Fig. 1(a).

IV Origin of Nematic Orders

To confirm the mechanism of the nematic transition, we perform the diagrammatic calculation for the MT- and AL-VCs. These VCs can be obtained by solving the CDW equation introduced in Ref. Onari:2016gs in the study of Fe-based superconductors. We analyze the linearized CDW equation introduced in Appendix D and in Ref. Kawaguchi-JPSJ . By solving the CDW equation, both MT- and AL-VCs shown in Figs. 2(d) and 2(e) with the optimized form factors are systematically generated. Figure 4(b) shows the eigenvalue of the linearized equation, $\lambda_{{\bm{q}}}$, for $\Delta E_{p}=0.1$ eV. Here, ${\alpha}_{S}$ is the spin Stoner factor, and the horizontal axis is proportional to $\chi^{\rm spin}_{\rm max}$. The CDW susceptibility increases with the increase of $\lambda_{{\bm{q}}}$. In Fig. 4(b), we set the quasiparticle damping $\gamma=0.3$ eV. Note that $\lambda_{{\bm{q}}}$ decreases with $\gamma$, whereas its overall ${{\bm{q}}}$-dependence is independent of $\gamma$, as shown in Appendix D. The obtained results are qualitatively consistent with the RG+cRPA results in Fig. 4(a). In Appendix D, we reveal that the CDW instabilities at ${{\bm{q}}}={\bm{0}}$ and ${{\bm{q}}}={{\bm{Q}}}_{\rm a}$ are given by the higher-order AL-type VCs.

Finally, we explain why $\Lambda_{\rm low}$ should be set small. The RG equation for the ${{\bm{q}}}={\bm{0}}$ vertex, $\bar{\Gamma}^{c}_{l}({{\bm{k}}};{{\bm{k}}}^{\prime})\equiv\Gamma^{c}_{l}({{\bm{k}}},{{\bm{k}}};{{\bm{k}}}^{\prime},{{\bm{k}}}^{\prime})$, is given as $d\bar{\Gamma}^{c}({{\bm{k}}};{{\bm{k}}}^{\prime})/dl\propto\Lambda_{l}\dot{f}(\Lambda_{l})\sum_{{{\bm{k}}}^{\prime\prime}}\delta(|E_{{{\bm{k}}}^{\prime\prime}}|-\Lambda_{l})\bar{\Gamma}^{c}_{l}({{\bm{k}}};{{\bm{k}}}^{\prime\prime})\bar{\Gamma}^{c}_{l}({{\bm{k}}}^{\prime\prime};{{\bm{k}}}^{\prime})$ ​+(other two terms), where $\dot{f}({\epsilon})$ is the derivative of the Fermi distribution function. Since the factor $|\dot{f}(\Lambda_{l})|$ is small for $\Lambda_{l}\gtrsim 4T$, the obtained $\bar{\Gamma}^{c}_{l}$ is strongly reduced if $\Lambda_{\mathrm{low}}\gg T$. In contrast, $\Gamma^{s,c}$ for ${{\bm{q}}}\neq\bm{0}$ is not so sensitive to $\Lambda_{\mathrm{low}}$. For this reason, the renormalization effect of $\Gamma_{l}^{c}$ is underestimated for ${{\bm{q}}}\approx\bm{0}$ if $\Lambda_{\rm low}\gg T$. Then, the obtained $\chi^{p\mbox{-}\mathrm{orb}}_{d}(\bm{0})$ is suppressed to be smaller than the peak values at ${{\bm{q}}}\neq\bm{0}$ if a large cutoff $\Lambda_{\mathrm{low}}\gg T$ is used, similarly to the previous results for $\Lambda_{\mathrm{low}}=\pi T$ Tsuchiizu et al. (2016). In the present study, large $\chi^{p\mbox{-}\mathrm{orb}}_{d}(\bm{0})$ is correctly obtained thanks to the use of the small cutoff $\Lambda_{\mathrm{low}}=T$.

V Doping-dependence of the CDW susceptibilities

In above sections, we studied the spin- and charge-susceptibilities in the $d$-$p$ Hubbard model only for 10% hole-doped case ($p=0.10$). To understand the experimental phase diagram in Fig. 1 (a), however, we have to study the doping dependence of susceptibilities. This issue is a very important but difficult goal for theorists. We find that the CDW is driven by the strong spin fluctuations, which strongly develop when the hole-density $p$ approaches zero experimentally.

In the RG+cRPA theory, the CDW is driven by the strong spin fluctuations, and spin fluctuations develop as the hole-carrier $p$ approaches zero experimentally. For this reason, as shown in Fig. 5, the uniform CDW susceptibility $\chi_{d}^{p\mbox{-}\mathrm{orb}}(0)$ linearly increases as $p$ decreases accompanied by the increment of $\chi_{\rm max}^{\rm spin}$ for $p\sim 0$. This result is consistent with the experimental $p$-dependence of $T^{*}$ in Fig. 1 (a) in the main text. In Fig. 5, we modify $U$ slightly so that the experimental approximate relation $\chi_{\rm max}^{\rm spin}\propto 1/p$ is satisfied, as performed in our previous study Yamakawa and Kontani (2015). We put $U=4.31$, $4.25$ and $4.06$ eV for $p=0.05$, $0.10$ and $0.20$, respectively.

In contrast to $T^{*}$, $T_{\rm CDW}$ decreases near the half-filling for $p<0.1$ as depicted in Fig. 1 (a). This behavior is also understood qualitatively based on the spin-fluctuation-driven mechanism. In fact, the axial CDW wavelength ${{\bm{Q}}}_{\rm a}$ is given by the nesting vector between the neighboring hot spots, and $|{{\bm{Q}}}_{\rm a}|$ increases as $p$ approaches zero. The CDW instability driven by the Aslamazov-Larkin vertex correction, which is qualitatively proportional to $\sum_{{\bm{q}}}\chi^{s}({{\bm{q}}})\chi^{s}({{\bm{q}}}+{{\bm{Q}}}_{\rm a})$, is suppressed if $|{{\bm{Q}}}_{\rm a}|$ is very large as we explained in Ref. Yamakawa and Kontani (2015). Therefore, the difference in the doping-dependences of $T^{*}$ and $T_{\rm CDW}$ is qualitatively understood. It is an important future issue to reproduce the experimental phase diagram in Fig. 1 (a) more completely, which is one of the greatest goals in this field.

VI Summary

In summary, we studied various unconventional CDW instabilities in the $d$-$p$ Hubbard model by using the RG+cRPA method, and predicted the multistage CDW transitions in cuprate superconductors. Based on the proposed spin-fluctuation-driven CDW mechanism, the following understanding has been reached: The short-range spin-fluctuations drive the uniform nematic CDW around $T^{*}$, and it triggers the axial ${{\bm{q}}}={{\bm{Q}}}_{\rm a}$ CDW at $T_{\rm CDW}$ successively. We also explained the doping-dependence of $T^{*}$ based on the RG+cRPA theory. These results naturally explain the phase diagram in Fig. 1 (a), except for heavily under-doped region. Although the uniform CDW order cannot simply explain the pseudogap formation, the large quasiparticle damping TPSC ; Kotliar ; Maier due to the short-range spin-fluctuations may induce the pseudogap for $T\lesssim T^{*}$.

Acknowledgements.

We are grateful to Y. Matsuda, T. Hanaguri, T Shibauchi, Y. Kasahara, Y. Gallais, W. Metzner, T. Enss, L. Classen, and S. Onari for fruitful discussions. This work was supported by Grant-in-Aid for Scientific Research from the Ministry of Education, Culture, Sports, Science, and Technology, Japan, and in part by Nara Women’s University Intramural Grant for Project Research.

Appendix A RG equations for the four-point vertex

In the main text, we analyzed the $d$-$p$ Hubbard model by using the RG+cRPA method Tsuchiizu et al. (2015). This method is the combination of the fRG theory and the cRPA. The RG+cRPA method enables us to perform reliable numerical study in the unbiased way. In this method, we introduce the original cutoff energy $\Lambda_{0}$ in order to divide each band into the higher and the lower energy regions: The higher-energy scattering processes are calculated by using the cRPA: The lower-energy scattering processes are analyzed by solving the RG equations, in which the initial vertices in the differential equation are given by the cRPA.

In the present model, the bare Coulomb interaction term on $d$-electrons is given as

[TABLE]

where $U^{0;c}=U$ and $U^{0;s}=-U$.

The antisymmetrized full four-point vertex ${\Gamma}({{\bm{k}}}+{{\bm{q}}},{{\bm{k}}};{{\bm{k}}}^{\prime}+{{\bm{q}}},{{\bm{k}}}^{\prime})$, which is the dressed vertex of the bare vertex ${\hat{U}}$ in Eq. (5) in the microscopic Fermi liquid theory, is depicted in Fig. 6(a). Reflecting the SU(2) symmetry of the present model, ${\Gamma}$ is uniquely decomposed into the spin-channel and charge-channel four-point vertices by using the following relation:

[TABLE]

where $\sigma,\sigma^{\prime},\rho,\rho^{\prime}$ are spin indices, and $\vec{\bf{\sigma}}$ is the Pauli matrix vector. We stress that ${\Gamma}^{c,s}$ are fully antisymmetrized, so the requirement by the Pauli principle is satisfied. We note that ${\Gamma}^{\uparrow\uparrow\uparrow\uparrow}=\frac{1}{2}{\Gamma}^{c}+\frac{1}{2}{\Gamma}^{s}$, ${\Gamma}^{\uparrow\uparrow\downarrow\downarrow}=\frac{1}{2}{\Gamma}^{c}-\frac{1}{2}{\Gamma}^{s}$, and ${\Gamma}^{\uparrow\downarrow\uparrow\downarrow}={\Gamma}^{s}$.

In the RG formalism, the four-point vertex function is determined by solving the differential equations, called the RG equations. In the band representation basis, the explicit form of the RG equations is given by Tazai et al. (2016):

[TABLE]

where $G(k)$ is the Green function multiplied by the Heaviside step function $\Theta(\Lambda-|E_{{{\bm{k}}},\nu}|)$, and $k$ is the compact notation of the momentum, band, and spin index: $k=({{\bm{k}}},{\epsilon}_{n},\nu,\sigma)$. The diagrammatic representation of the RG equations is shown in Fig. 6(b). The first two contributions in the right-hand-side represent the particle-hole channels and the last contribution is the particle-particle channel.

In a conventional fRG method, $\Lambda_{0}$ is set larger than the bandwidth $W_{\rm band}$, and the initial value is given by the bare Coulomb interaction in Eq. (5). In the RG+cRPA method, we set $\Lambda_{0}<W_{\rm band}$, and the initial value is given by the constrained RPA to include the higher-energy processes without over-counting of diagrams Tsuchiizu et al. (2015).

In the main text, we introduced the lower-energy cutoff $\Lambda_{\rm low}\ (\sim T)$ in the RG equation for the four-point vertex: Eq. (7). For this purpose, we multiply the cutoff function $((\Lambda_{\rm low}/\Lambda)^{\zeta}+1)^{-1}$ to the right-hand-side of Eq. (7). Here, $\zeta$ is a parameter determining the width of this smooth cutoff, and we set $\zeta=10$ in the main text. We do not introduce the lower-energy cutoff in the RG equation for the susceptibilities.

Appendix B RG+cRPA analysis of the $d$-$p$ model with

different hopping parameters

In Ref. Hansmann et al. (2014), the authors derived the realistic $d$-$p$ models for La-based cuprate, by evaluating the hopping parameters on the basis of the $N$th order muffin tin orbitals ($N$MTO). The model parameters for $N=0$ and $N=1$ are given in Table 1. The bandstructure of the $N=0$ basis model (0MTO model) is very close to the LDA bandstructure near Fermi energy. For this reason, we have analyzed the 0MTO model in the main text. On the other hand, the $N=1$ basis model (1MTO model) appropriately reproduces the overall oxygen bonding band structure with deep bottom energy $E\simeq-8$ eV. In order to check the reliability of our RG+cRPA results, we analyze the $d$-$p$ model with the 1MTO model parameters.

Figure 7(a) shows the band structure of the 1MTO model. Here, we introduced the third-nearest $d$-$d$ hopping $t_{dd}^{\mathrm{3rd}}=-0.1$ eV to make the FS closer to Y-based cuprates. The FS of this model is shown in Fig. 7(b). Now, we analyze this model by using the RG+cRPA method. The parameters are the same as in the main text except for $U$. The number of patches is $N_{p}=128$ and the initial cutoff is $\Lambda_{0}=0.5$ eV. The temperature is fixed at $T=0.1$ eV.

In Fig. 7(c), we show the obtained $\chi^{p\mbox{-}\mathrm{orb}}_{d}(\bm{q})$ for $U=5.72$ eV. The RPA results are also shown for comparison. It has the largest peak at $\bm{q}=\bm{0}$ and the second largest peak at $\bm{q}=\bm{Q}_{\mathrm{a}}$, respectively. The obtained $\bm{q}$-dependence of $\chi^{p\mbox{-}\mathrm{orb}}_{d}(\bm{q})$ is similar to Fig. 2(b). We also investigate the $U$-dependences of the spin and charge susceptibilities. As shown in the inset of Fig. 7(d), relatively large $U$ is required for the enhancement of $\chi^{\mathrm{spin}}_{\rm max}$ in the 1MTO model, since the density of states of the $d$ orbital at the Fermi energy in the 1MTO model is smaller than that in the 0MTO model Hansmann et al. (2014). In Fig. 7(d), we plot the peak values of $\chi^{p\mbox{-}\mathrm{orb}}_{d}(\bm{q})$ as functions of $\chi^{\mathrm{spin}}_{\rm max}$. The obtained results are quite similar to Fig. 3(a) in the main text. Thus the spin-fluctuation-driven CDW instabilities are universal phenomena in both the 0MTO model and 1MTO model.

In summary, we investigate the $d$-$p$ model with 1MTO model parameters. We found that the results are very similar to those for the 0MTO model given in the main text. Therefore the mechanism of the spin-fluctuation-driven CDW instabilities revealed in the main text is universal, independently of the details of the model parameters.

Appendix C RG+cRPA analysis for the single $d$-orbital Hubbard model

In the main text, we studied the 0MTO $d$-$p$ Hubbard model based on the RG+cRPA theory, and found that the $p$O-CDW susceptibilities develop strongly in the strong spin-fluctuation region. Similar results are obtained in the 1MTO model in which ${\epsilon}_{d}-{\epsilon}_{p}$ is 0.53eV larger than that in the 0MTO model, as we show in Appendix B. In these $d$-$p$ models, any Coulomb interactions on $p$-orbitals are not taken into account. Therefore, spin-fluctuation-driven CDW formation is also expected to be realized in the single $d$-orbital Hubbard model with on-site Coulomb interaction.

Here, we study the single $d$-orbital Hubbard model with the first-, the second-, and the third-nearest hopping integrals as $t=-0.50$ eV, $t^{\prime}=0.083$ eV, and $t^{\prime\prime\prime}=-0.10$ eV, respectively. The bandstructure and Fermi surface for $n=0.90$ are shown in Figs. 8 (a) and (b), respectively. We calculate the $d$-electron charge susceptibility $\chi^{d\mbox{-}\mathrm{orb}}(\bm{q})$ with the $B_{1g}$ form factor $f_{\bm{q}}(\bm{k})=\cos(k_{x})-\cos(k_{y})$ introduced in Eq. (3) in the main text, The obtained results are summarized in Fig. 8 (c): Both $\chi^{d\mbox{-}\mathrm{orb}}({\bm{0}})$ and $\chi^{d\mbox{-}\mathrm{orb}}({{\bm{Q}}}_{\rm a})$ are strongly enlarged in the strong spin-fluctuation region, very similarly to the $p$O-CDW susceptibility shown in Fig. 3 (a) in the main text.

Therefore, it was verified that our main numerical results in the main text are unchanged even in the single $d$-orbital model, once the $B_{1g}$ form factor is taken into account. We also analyzed the CDW equation for the single $d$-orbital model, and obtained the strong CDW instability. The obtained form factor $\Delta\Sigma_{{\bm{0}}}({{\bm{k}}})$ has the $B_{1g}$-symmetry. In real space, this is the bond-order (=modulation of the hopping integrals) given by the Fourier transformation of the symmetry-breaking self-energy $\Delta\Sigma_{{\bm{0}}}({{\bm{k}}})$. Thus, the robustness of the spin-fluctuation-driven CDW mechanism has been clearly confirmed.

Appendix D Analysis of the linearized CDW equation

In the main text, we analyzed the $d$-$p$ Hubbard model for cuprate superconductors in an unbiased manner using the RG+cRPA method. We find that the nematic CDW with $d$-form factor is the leading instability. The axial nematic CDW instability at ${{\bm{q}}}={{\bm{Q}}}_{\rm a}$ is the second strongest, and its strength increases under the static uniform CDW order. This result leads to the prediction that uniform nematic CDW occurs at the pseudogap temperature $T^{*}$, and the axial CDW at wavevector ${{\bm{q}}}={{\bm{Q}}}_{\rm a}$ is induced under $T^{*}$.

In this section, we study the CDW formation mechanism in cuprate superconductors based on the diagrammatic method, in order to find what many-body processes cause the CDW order. Theoretically, the CDW order is given as the symmetry-breaking in the self-energy $\Delta\Sigma(k)$. According to Refs. Onari and Kontani (2012); Kawaguchi-JPSJ , the self-consistent CDW equation is given as

[TABLE]

where $P_{\rm A_{1g}}$ is the $A_{\rm 1g}$-symmetry projection operator, and $G(k)=(G_{0}^{-1}(k)-\Delta\Sigma(k))^{-1}$ is the $d$-electron Green function with the symmetry-breaking term $\Delta\Sigma$. $V(q)=U^{2}(\frac{3}{2}\chi^{s}(q)+\frac{1}{2}\chi^{c}(q)-\chi^{0}(q))+U$, where ${\chi}^{s(c)}(q)={\chi}^{0}(q)/({1}-(+){U}{\chi}^{0}(q))$ and $\chi^{0}(q)=-T\sum_{k}G(k+q)G(k)$.

In order to analyze the CDW state with arbitrary wavevector ${{\bm{q}}}$, we linearize Eq. (8) with respect to $\Delta\Sigma$:

[TABLE]

where $\lambda_{{\bm{q}}}$ is the eigenvalue for the CDW for the wavevector ${{\bm{q}}}$. The CDW with wavevector ${{\bm{q}}}$ appears when $\lambda_{{\bm{q}}}=1$, and the eigenvector $\Delta\Sigma_{{\bm{q}}}(k)$ gives the CDW form factor. The kernel $K({{\bm{q}}},k,k^{\prime})$ is given in Fig. 9(a), and its analytic expression is given as Kawaguchi-JPSJ :

[TABLE]

where $V_{0}^{s}(q)=U+U^{2}\chi_{0}^{s}(q)$, $V_{0}^{c}(q)=-U+U^{2}\chi_{0}^{c}(q)$, and $\Lambda_{{\bm{q}}}(k;p)\equiv G_{0}(k+\frac{{{\bm{q}}}}{2})G_{0}(k-\frac{{{\bm{q}}}}{2})G_{0}(k-p)$. The subscript [math] in Eq. (10) represents the functions with $\Delta\Sigma=0$.

By solving the linearized CDW equation (9), many higher-order vertex corrections (VCs) are systematically generated. Some examples of the generated VCs are shown in Fig. 9(b). If we drop the Hartree term and MT term in $K({{\bm{q}}};k,k^{\prime})$, we obtain the series of higher-order AL-VCs shown in Fig. 9(c). The AL terms drive the ${{\bm{q}}}={\bm{0}}$ CDW instability since its functional form $\propto\sum_{{\bm{k}}}\chi^{s}({{\bm{k}}}+{{\bm{q}}})\chi^{s}({{\bm{k}}})$ is large for ${{\bm{q}}}\approx{\bm{0}}$ Onari:2016gs .

Figure 10(a) shows the obtained ${{\bm{q}}}$-dependence of $\lambda_{{\bm{q}}}$ for ${\alpha}_{S}=0.995$ at $T=50$ meV. Here, we introduced the quasiparticle damping $\gamma=0.3$ eV into $G_{0}(k)$. Here, $\lambda_{{\bm{q}}}$ is the largest at ${{\bm{q}}}={\bm{0}}$, and the second largest maximum is at ${{\bm{q}}}={{\bm{Q}}}_{\rm a}$. ${\alpha}_{S}\equiv U{\max}_{{\bm{q}}}\{\chi^{0}_{0}({{\bm{q}}})\}$ is the spin Stoner factor. We also show the eigenvalue $\lambda^{\rm AL}_{{\bm{q}}}$ (and the second-largest eigenvalue $\lambda^{\rm AL,2nd}_{{\bm{q}}}$), which is obtained by dropping the Hartree and MT terms in the kernel. That is, $\lambda^{\rm AL}_{{\bm{q}}}$ is given by the higher-order AL processes shown in Fig. 9(c). At ${{\bm{q}}}={\bm{0}}$ and ${{\bm{Q}}}_{\rm a}$, $\lambda^{\rm AL}_{{\bm{q}}}$ is almost equal to the true eigenvalue $\lambda_{{\bm{q}}}$.

In the present analysis, we dropped the ${\epsilon}_{n}$-dependence of $\Delta\Sigma_{{\bm{q}}}(k)$ by performing the analytic continuation $(i{\epsilon}_{n}\rightarrow{\epsilon})$ and putting ${\epsilon}=0$. We also dropped the ${\epsilon}_{n}$-dependence of the quasiparticle damping $\gamma$. Due to these simplifications, the obtained $\lambda_{{\bm{q}}}$ is overestimated. Therefore, we do not put the constraint $\lambda_{{\bm{q}}}<1$ here.

In Fig. 10(b), we show the eigenvalue $\lambda^{\rm MT}_{{\bm{q}}}$, which is obtained by dropping the Hartree and AL terms in the kernel. It is much smaller than $\lambda_{{\bm{q}}}$ at ${{\bm{q}}}={\bm{0}}$ and ${{\bm{Q}}}_{\rm a}$, whereas $\lambda_{{\bm{q}}}$ at ${{\bm{q}}}={{\bm{Q}}}_{\rm d}$ is comparable to the true eigenvalue. Therefore, the origin of the CDW instability at ${{\bm{q}}}={\bm{0}}$ and ${{\bm{Q}}}_{\rm a}$ is the AL process, whereas that at ${{\bm{q}}}={{\bm{Q}}}_{\rm d}$ is mainly the MT process.

Figure 10 (c) shows the eigenvalues at ${{\bm{q}}}={\bm{0}}$, ${{\bm{Q}}}_{\rm a}$, and ${{\bm{Q}}}_{\rm d}$ as function of ${\alpha}_{S}$. As the spin susceptibility increases (${\alpha}_{S}\gtrsim 0.98$), $\lambda_{{\bm{q}}}$ is drastically enlarged by the VCs, and $\lambda_{{{\bm{q}}}={\bm{0}}}$ becomes the largest due to the AL processes. The form factor at ${{\bm{q}}}={\bm{0}}$, $\Delta\Sigma_{{\bm{0}}}({{\bm{k}}})$, has the $d$-wave symmetry as shown in Fig. 10 (d).

We stress that the eigenvalue $\lambda_{{\bm{q}}}$ is quickly suppressed by increasing $\gamma$, which is actually large in cuprates. Figures 10 (e) and (f) show the CDW susceptibilities for larger damping rate $\gamma=0.6$ eV, in the cases of $\Delta E_{p}=0$ and $\Delta E_{p}=0.1$ eV, respectively. (Note that the damping rate is renormalized to be $\sim\gamma/5$ in cuprates.) In Fig. 10 (e), $\lambda_{{\bm{q}}}$ reaches unity first at ${{\bm{q}}}={\bm{0}}$ with increasing $\alpha_{S}$. In the nematic state with $\Delta E_{p}=0.1$ eV shown in Fig. 10 (f), $\lambda_{{\bm{q}}}$ at ${{\bm{q}}}={{\bm{Q}}}_{\rm a}^{x}$ exceeds $\lambda_{{{\bm{Q}}}_{\rm d}}$ for ${\alpha}_{S}>0.996$. The corresponding eigenvalue is $\sim 1.4$, which decreases with increasing $\gamma$. This result supports the main result of the present RG+cRPA study shown in Fig. 4 (a) in the main text.

In summary, we analyzed the linearized CDW equation based on the $d$-$p$ Hubbard model, by including both the MT and AL VCs into the kernel. When the spin fluctuations are strong (${\alpha}_{S}\gtrsim 0.98$), the uniform nematic CDW has the strongest instability. The axial CDW instability is strongly magnified under the uniform CDW order, as we explain the main text. The obtained results are consistent with the results by the RG+cRPA in the main text. Thus, it is concluded that the higher-order AL processes give the CDW orders at ${{\bm{q}}}={\bm{0}}$ and ${{\bm{Q}}}_{\rm a}$.