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

HAL
open science

# Charge, bond, and pair density wave orders in a strongly correlated system

Anurag Banerjee, Catherine Pépin, Amit Ghosal

► To cite this version:

Anurag Banerjee, Catherine Pépin, Amit Ghosal. Charge, bond, and pair density wave orders in a strongly correlated system. Physical Review B, 2022, 105 (13), pp.134505. 10.1103/Phys-RevB.105.134505. cea-04212156

HAL Id: cea-04212156

https://cea.hal.science/cea-04212156

Submitted on 20 Sep 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Charge, bond, and pair density wave orders in a strongly correlated system

Anurag Banerjee,¹ Catherine Pépin,¹ and Amit Ghosal²

¹Institut de Physique Théorique, Université Paris-Saclay, CEA, CNRS, F-91191 Gif-sur-Yvette, France

²Indian Institute of Science Education and Research Kolkata, Mohanpur, India-741246

The coexistence of multiple quasi-degenerate orders is the hallmark of the strongly correlated materials. Experiments often reveal several spatially modulated orders in the underdoped cuprates. This has come to the forefront with the possible detection of the pair density wave states. However, microscopic calculations often struggle to stabilize such spatially modulating orders as the ground state in the strong correlation limit. This work uses the $t-t'$-J-model with an additional nearest-neighbor repulsion to stabilize spatially oscillating charge, bond, and pairing orders in the underdoped regime. We employ the standard Gutzwiller approach while treating the inhomogeneity for the spatial orders using the self-consistent Hartree-Fock-Bogoliubov methodology. Our calculations reveal that unidirectional bond density states coexisting with charge and pairing modulations can have lower energy than the uniform superconducting state over an extensive doping range. These modulating states vanish monotonically as the modulation wavevector becomes shorter with increased dopings. The finite momentum orders give way to a vestigial nematic phase upon increasing doping which only breaks the rotational symmetry of the system. The nematic order vanishes upon further increasing doping, and only uniform superconductivity survives. The spatial features of the ground state at each doping reveal multiple wavevectors, which potentially drives the incommensuration of charge orders. Interestingly, the spatially modulating states are absent when the strong correlations criteria are relaxed, suggesting that the removal of double occupancy aids the stabilization of density wave orders.

## I. INTRODUCTION

The presence of multiple broken symmetry phases and competing orders has become the new paradigm for describing correlated quantum matter featuring complex phase diagram [1-4]. Given the abundance of rich phase diagrams in different materials suggests the intertwining of several orders [5, 6]. In this scenario, a precursor state of all potentially degenerate states forms at a high energy scale which cannot be linked with one particular order. At low energies, slight energy differences stabilize one order over the other.

A unique realization for intertwined orders is the vestigial or composite phase, where a higher-order component of the likely symmetry broken phase condenses although the individual components vanish [5, 6]. Therefore, a vestigial state has a smaller subset of broken symmetries than the parent state. For instance, a nematic phase that breaks the rotational symmetry is a vestige of the primary checkerboard and stripe phase where both translation and rotational symmetries are broken [7]. Recently, the sublattice phase-resolved study of electronic structure in Bi-based cuprates compound finds a vestigial nematic phase around 14-17% hole doping [8]. Moreover, the nematic order coexists with the charge density wave (CDW) orders at lower dopings [8].

The signatures of charge modulation are routinely observed at moderate hole-doped materials in all cuprate families [9-11]. Usually, the modulation wavevector, $Q$, reduces as the doping increases [12, 13] except for the $La_2CuO_4$ where the wavevector increases with doping [14-16]. Moreover, detailed investigations have also revealed a bond-centered charge order with a primarily d-wave form factor [17-19]. Furthermore, disorder severely affects the charge modulations such that the unidirectional stripe and the bidirectional checkerboard patterns become indistinguishable [20]. The local spectroscopic probes indicate short-ranged domains with unidirectional charge modulation [21, 22], albeit quantum os-

cillation measurements favor the bidirectional picture [23].

Recently the study of the modulated orders in cuprates has gained impetus with the observations of pair density wave (PDW) state coexisting with usual d-wave superconductivity (SC) using Josephson scanning tunneling microscopy (STM) [24, 25]. In the same regime, d-wave CDW oscillations are also detected [25]. The importance of detecting the PDW state is underlined by the several exotic proposals of the celebrated pseudogap state that hinges on the quantum disordering [26] or fractionalization of the PDW state [27, 28]. However, it is extremely difficult to stabilize such PDW states as the ground state of a microscopic Hamiltonian [29-33].

The $tJ$-model and its parent Hubbard model at large repulsion $U$ are generally studied to extract the physics of spatially modulated states [34-36]. The variational approaches on the two-dimensional (2D) $tJ$ model find a striking near-degeneracy among the unidirectional charge order with stripe SC, a PDW state with spin density wave component, and the uniform d-wave SC state [37]. A few recent calculations on the same model observe a pair density wave state with vanishing mean SC pairing energetically very close to the uniform d-wave SC [38, 39]. Such quasi-degeneracy among different states suggests that a small perturbation to the $tJ$ model will favor one state over the other, thus giving rise to a wide variety of phenomenology.

Using commonly accepted minimal models to describe strongly correlated quantum matters, i.e., the $t-t'-J$ model with an additional nearest-neighbor repulsion, we study the interplay of unconventional superconductivity with spatially modulated charge-ordered states in the ground state. Early investigations on such model manifest promising signatures of density wave orders [40-44]. A careful analysis of the variation of the ordering wavevector with doping, however, is still missing. In addition, the real-space pattern of the charge, bond order, and the possibility of finding spatially modulated pairing amplitude (similar to the Fulde-Ferrell-
Larkin-Ovchinnikov phase but in the absence of magnetic field) need to be studied. Furthermore, the mechanism of the demise of modulated orders with increasing doping remains unexplored. In this work, we focus on the unidirectional charge orders that are routinely observed in the underdoped cuprates [21, 22, 45]. Furthermore, we neglect the magnetically ordered states to keep the numerical calculation tractable. In the next section, we describe the model and the method in detail.

## II. MODEL AND METHODS

We work with the standard tJ-model on a two-dimensional square lattice with all the double occupancy removed, with an additional nearest neighbor repulsion. The model is given by

$$ \mathcal{H} = \sum_{i,j,\sigma} (t_{ij} \tilde{c}_{i\sigma}^{\dagger} \tilde{c}_{j\sigma} + \text{h.c.}) + J \sum_{\langle i,j \rangle} \left( \tilde{S}_i \cdot \tilde{S}_j - \frac{1}{4} \tilde{n}_i \tilde{n}_j \right) \\ + W \sum_{\langle i,j \rangle, \sigma, \sigma'} \hat{n}_{i\sigma} \hat{n}_{j\sigma'} - \mu \sum_{i,\sigma} \hat{n}_{i,\sigma}. \quad (1) $$

Here, $t_{ij} = -t$ if $i, j$ are nearest neighbors and $t_{ij} = t'$ if $i, j$ are next nearest neighbor sites. We fix the $t = 1$ and $t' = 0.25t$, and all the energy scales are in the units of $t$. To compare with the experimental energy scales in cuprates, one can use $t \approx 300$ meV.

The tJ Hamiltonian is an effective low-energy model obtained from a Schrieffer-Wolf transformation of the Hubbard model for large on-site repulsion $U$ [46]. The $S_i$ are the electron spin operator and the exchange interaction $J \approx (4t^2)/U$, is set to $J = 0.3t$ in this work. Furthermore, $W$ is the nearest neighbor repulsion term, which we fix to $W = 0.6t$ for the rest of the analysis.

The local number operators for a particular spin is given by $\hat{n}_{i,\sigma} = c_{i,\sigma}^{\dagger}c_{i,\sigma}$ and the local density is the expectation value of the number operator. The 'tilde' appearing on the creation (annihilation) operators at site $i$ with spin $\sigma = \uparrow, \downarrow$, $c_{i,\sigma}^{\dagger}$ ($c_{i,\sigma}$) in Eq. (1) is explained below. The local density is thus given by $\rho_i = \sum_{\sigma} \langle c_{i\sigma}^{\dagger}c_{i\sigma} \rangle$. $\mu$ is the chemical potential that fixes the average density of carriers $\rho$ in the system. The spatial average of the local density gives the average density of electrons $\rho$. The (under) doping level from the half-filling is defined via $\delta = (1 - \rho)$.

Strong on-site repulsions reduce the Hilbert space of the tJ-model by eliminating all double occupancies. The 'tilde' on operators appearing in the Eq. (1) signifies that they operate on a truncated Hilbert space which projects out double occupation on any site due to strong repulsive interactions. The operators without 'tilde' live in the usual unrestricted Hilbert space.

$$ \tilde{c}_{i\sigma} = c_{i\sigma}(1 - \hat{n}_{i\sigma}), \quad (2) $$

where $c$ is the normal annihilation operator. A straightforward implementation of such a restriction is given by the

Gutzwiller approximation, which transforms the strongly interacting problem to an unrestricted Hilbert space, at the expense of introducing “Gutzwiller factors”, which renormalize all operators [47–49].

$$ \langle \tilde{c}_{i\sigma}^{\dagger} \tilde{c}_{j\sigma} \rangle \approx g_{ij}^{t} \langle c_{i\sigma}^{\dagger} c_{j\sigma} \rangle_{0}, \quad (3) $$

$$ \langle \tilde{S}_i \cdot \tilde{S}_j \rangle \approx g_{ij}^{xy} \langle \mathbf{S}_i \cdot \mathbf{S}_j \rangle_0, \quad (4) $$

$$ \langle n_{i\sigma} n_{j\sigma} \rangle \approx \langle n_{i\sigma} n_{j\sigma} \rangle_0, \quad (5) $$

where $\langle \dots \rangle$ is the expectation in the restricted Hilbert space whereas $\langle \dots \rangle_0$ is the expectation of operators in the unrestricted Hilbert space. The Gutzwiller factors in the absence of any magnetic orders are given by

$$ g_{ij}^{t} = \sqrt{\frac{4x_i x_j}{(1+x_i)(1+x_j)}}, \quad (6) $$

$$ g_{ij}^{xy} = \frac{4}{(1+x_i)(1+x_j)}, \quad (7) $$

where $x_i = 1 - \rho_i$. After removing the double-occupancy, we perform the mean field decomposition of the interaction terms in Hartree, Fock and Cooper channels as given below,

$$ \sum_{\sigma} \langle c_{i\sigma}^{\dagger} c_{i\sigma} \rangle_0 = \rho_i, \quad (8) $$

$$ \langle c_{i\sigma}^{\dagger} c_{j\sigma} \rangle_0 = \chi_{ij}, \quad (9) $$

$$ \langle c_{j\downarrow} c_{i\uparrow} + c_{i\downarrow} c_{j\uparrow} \rangle_0 = \Delta_{ij}. \quad (10) $$

Note that we are not allowing for any magnetic orderings in this system by focusing on superconducting with charge and bond density wave (BDW) orders only. Therefore, we eliminate the spin-flip Fock terms from the mean-field Hamiltonian, to preserve the spin rotational symmetry. To perform the mean field analysis we closely follow Ref. [38, 50, 51]. The explicit expression for the mean-field Hamiltonian is given by

$$ \mathcal{H}_{MF} = \sum_{i,a,\sigma} (t g_{i,i+a}^t - W_{i,i+a}^{FS}) c_{i,\sigma}^{\dagger} c_{i+a,\sigma} \\ + \sum_{i,b,\sigma} t' g_{i,i+b}^t c_{i,\sigma}^{\dagger} c_{i+b,\sigma} + \sum_{i,\sigma} (-\mu + \mu_i^{HS}) \hat{n}_{i\sigma} \\ + \sum_{i,a} \left[ G_1^{i,a} \Delta_{i,i+a} \left( c_{i,\uparrow}^{\dagger} c_{i+a,\downarrow}^{\dagger} + c_{i+a,\uparrow}^{\dagger} c_{i,\downarrow}^{\dagger} \right) + \text{h.c.} \right]. \quad (11) $$

Here, the sum over $\mathbf{a} = \pm \hat{x}, \pm \hat{y}$ are the nearest neighbor vectors and $\mathbf{b} = \pm (\hat{x} \pm \hat{y})$ are the next neighbor vector. The other
variables are given by,

$$W_{i,a}^{FS} = \frac{J}{2} \left( \frac{3g_{i,i+a}^{xy}}{2} - \frac{1}{2} \right) \chi_{i,i+a} + W \chi_{i,i+a}, \quad (12)$$

$$G_1^{i,a} = -\frac{J}{4} \left( \frac{3g_{i,i+a}^{xy} + 1}{4} \right), \quad (13)$$

$$\mu_i^{HS} = -4t \sum_{a,\sigma} \frac{\partial g_{i,i+a}^t}{\partial n_i} \chi_{i,i+a} + 4t' \sum_{b,\sigma} \frac{\partial g_{i,i+b}^t}{\partial n_i} \chi_{i,i+b} \\ - \frac{3J}{2} \sum_{a,\sigma} \frac{\partial g_{i,i+a}^{xy}}{\partial n_i} \left( \frac{\Delta_{i,i+a}^2}{4} + \chi_{i,i+a}^2 \right) \\ - \frac{J-2W}{4} \sum_{i,a,\sigma} \rho_{i+a} \quad (14)$$

Next, we diagonalize the Hamiltonian using Bogoliubov-deGennes(BdG) transformation [52–54],

$$c_{i\sigma} = \sum_{n} \left( \gamma_{n\sigma} u_{i,n} - \sigma \gamma_{n,-\sigma}^{\dagger} v_{i,n}^{*} \right), \quad (15)$$

where $\gamma^\dagger(\gamma)$ are Bogoliubov creation (annihilation) operator. The resulting eigensystem is then solved self-consistently for all the local Hartree, Fock and Bogoliubov order parameters and $\mu$. If there are $N$ lattice sites, then there are $(5N+1)$ variables to be solved self-consistently. We initialize the SC pairing $\Delta_{ij}$ and the bond density $\chi_{ij}$, and local density to be the modulating around a mean value with a desired wavevector $Q = (1/\lambda)2\pi/a_0\hat{x}_i$. Here $\lambda$ is the period of oscillations and $a_0$ is the lattice-spacing. We focus on the unidirectional density wave orders in this paper. The initial guess is given by

$$\Delta_{i,i+a} = \Delta_0 + \Delta_Q \cos(Q \cdot r_i). \quad (16)$$

$\Delta_0$ is the guess mean SC pairing amplitude, and $\Delta_Q$ is the initial guess for SC pairing modulation amplitude. We initialize the bond and local density similarly and use such a trail solution to achieve self-consistency. We perform further analysis on those solutions for the cases we achieve self-consistency with the desired wavevector.

One can evaluate the self-consistent d-wave and extended s-wave components from the SC pairing $\Delta_{ij}$ by using the following definition

$$\Delta_d(i) = \frac{1}{4} (\Delta_{i,i+\hat{x}} + \Delta_{i,i-\hat{x}} - \Delta_{i,i+\hat{y}} - \Delta_{i,i-\hat{y}}), \quad (17)$$

$$\Delta_s(i) = \frac{1}{4} (\Delta_{i,i+\hat{x}} + \Delta_{i,i-\hat{x}} + \Delta_{i,i+\hat{y}} + \Delta_{i,i-\hat{y}}). \quad (18)$$

Similar definition can be made for the d-wave and extended s-wave component for the bond density order $\chi_{ij}$.

$$\chi_d(i) = \frac{1}{4} (\chi_{i,i+\hat{x}} + \chi_{i,i-\hat{x}} - \chi_{i,i+\hat{y}} - \chi_{i,i-\hat{y}}), \quad (19)$$

$$\chi_s(i) = \frac{1}{4} (\chi_{i,i+\hat{x}} + \chi_{i,i-\hat{x}} + \chi_{i,i+\hat{y}} + \chi_{i,i-\hat{y}}). \quad (20)$$

To extract the finite-momentum component of the charge, bond and pairing, we first calculate the Fourier components

of the variation of an ordering. For example, the Fourier component for the local density, is given by

$$\tilde{\rho}(\mathbf{q}) = \frac{1}{N} \sum_i \exp(i\mathbf{q} \cdot \mathbf{r}_i)(\rho_i - \rho). \quad (21)$$

Note that we remove the average density from the local density to consider only the oscillating part instead of the uniform part. The principal peak from the Fourier component gives us the modulation wavevector $Q$ for the desired quantity. The modulation amplitude from the local density is evaluated by the strength of $\tilde{\rho}(q)$ at this wavevector such that $\chi_{\rho}^Q \equiv \tilde{\rho}(Q)$. Similarly, we can define the other modulation amplitude by using the Fourier transform of the respective order parameters, as

$$\tilde{\chi}_d(\mathbf{q}) = \frac{1}{N} \sum_i \exp(i\mathbf{q} \cdot \mathbf{r}_i)(\chi_d(i) - \chi_d), \quad (22)$$

$$\tilde{\chi}_s(\mathbf{q}) = \frac{1}{N} \sum_i \exp(i\mathbf{q} \cdot \mathbf{r}_i)(\chi_s(i) - \chi_s), \quad (23)$$

$$\tilde{\Delta}_d(\mathbf{q}) = \frac{1}{N} \sum_i \exp(i\mathbf{q} \cdot \mathbf{r}_i)(\Delta_d(i) - \Delta_d), \quad (24)$$

$$\tilde{\Delta}_s(\mathbf{q}) = \frac{1}{N} \sum_i \exp(i\mathbf{q} \cdot \mathbf{r}_i)(\Delta_s(i) - \Delta_s). \quad (25)$$

Here $\chi_d$ are the spatially averaged component d-wave bond density. The finite-momentum component at the dominant wavevector is again defined as $\chi_d^Q \equiv \tilde{\chi}_d(Q)$. Similar definitions are also used for the other variables.

We can also track the broken rotational symmetry of the bond density and cooper pair density operators by using [7]

$$\mathcal{N}_{\chi} = \frac{1}{2} \sqrt{\langle \chi_x^2 \rangle - \langle \chi_y^2 \rangle} \quad (26)$$

$$\mathcal{N}_{\Delta} = \frac{1}{2} \sqrt{\langle \Delta_x^2 \rangle - \langle \Delta_y^2 \rangle} \quad (27)$$

where the average $\langle...\rangle$ is over the lattice sites. Notice that this quantity only contributes when there is a nematicity and will vanish for purely extended s-wave and d-wave component of bond or cooper pair density.

In order to describe our results in the subsequent sections, our strategy is the following. First, we extract the optimal wavevector by comparing the variational ground state energy of the self-consistent BdG solutions at each doping. Subsequently, we utilize the optimal self-consistent solution to study the variation of different order parameters with doping. Finally, we also show the local features of the optimal solutions for a few dopings.

### III. RESULTS

Our calculations are done on a two-dimensional square lattice unit-cell of linear dimension $40a_0$ to $60a_0$. The lattice size was tuned withing this range to keep the lattice commensurate with the period of modulations. We found that the self-consistency becomes easier if we accommodate even number
<figure><img src="image_3.png" /></figure>

FIG. 1. Compares the ground state energy of the self-consistent solutions with different wavevectors for a few dopings. The thick lines are the guide to the eyes. The dotted lines display the $E_{GS}$ for the uniform and the nematic states. (a) For $\delta = 0.075$ there is clearly a minima for $\lambda = 4$. As the doping rises to $\delta = 0.1$, the optimal modulation wavelength increases to $\lambda = 5$. The wavelength keeps increasing further with reduced doping, becoming $\lambda = 8$ for $\delta = 0.125$. At $\delta = 0.14$, the nematic state becomes the optimal ground state. The inset of (f) demonstrates the decrease of optimal wavevector with doping.

of full periodic density wave oscillations within the unit-cell. For instance, for a system with period of $8a_0$, we use a system size of $48a_0 \times 48a_0$. We have also checked that our qualitative results are independent of the system size.

A. Optimal modulation wavevector with doping

We begin our analysis by studying the ground state energy for different self-consistent solutions of the BdG equations. Initially, we start with a trail state with a single modulation wavevector at $Q = (1/\lambda)2\pi/a_0$ for the pairing amplitude, bond and charge density and try to achieve self-consistency. For each case, we achieve self-consistent solutions for the desired wavevector; we calculate the ground state energy $E_{GS} = \langle H \rangle_{MF}$, which is presented in Fig. (1) for different wavevectors and for several dopings. The minimum of $E_{GS}(Q)$ gives us the optimal solution at each doping level.

We also present the estimate of the $E_{GS}$ for the uniform and

the nematic orders ($Q = 0$) with the dashed lines. The nematic order breaks the $C_4$-rotational symmetry of the square lattice, although it is spatially uniform. Therefore, this is the charge (and pairing) nematic state where the x and the y-component of the bond (and pairing) density are different(See Sec (III C 3) for details).

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

FIG. 2. Presents the evolution of different order parameters with doping. (a) Shows the evolution of the amplitude of oscillations for the different quantities at the optimal ground state. The blue lines exhibit the oscillation amplitude for the local density at the optimal wavevector. Similarly, the orange and green curves show the evolution of finite-momentum amplitude for extended s-wave and d-wave bond density. Moreover, the same for the superconducting pairing in extended s-wave and d-wave channels is presented in red and violet traces. The lines are the guide to the eyes. All these spatially modulated orders vanish near the $\delta \approx 0.14$. (b) Displays the mean d-wave superconducting pairing in the blue trace, compared with the case where no charge or bond density orders are allowed $\Delta_d^0$ in cyan. The average $\Delta_d$ significantly reduces in the regions where the charge and bond orders are stable. Moreover, we show the renormalized d-wave pairing $\tilde{\Delta}_d$, which also reduces in the regime where charge and pairing oscillations are dominant i.e., $\delta < 0.16$. (c) Presents the nematicity in the $\chi$ and $\Delta$ with red and black traces respectively. The nematicity vanishes around the doping $\delta = 0.16$.

At low-doping, the evolution of $E_{GS}$ has a minima at $\lambda = 4$, as seen in Fig. (1a) and Fig. (1b). The purely nematic and the uniform states also have much higher energy than the optimal solution. As the doping is increased to $\delta = 0.1$, the energy for a system with density waves of primarily $Q = (1/5)2\pi/a_0$ minimizes the system. The optimal modulation wavevector shortens further to $Q_0 = 1/6 - 1/8$ in r.l.u. with increasing doping as seen in Fig. (1d) and Fig. (1e). The decreasing optimal wavevector with doping is also presented in the inset of Fig. (1f). Increasing the doping further makes the nematic state energetically favorable than the spatially modulating ones, as shown in Fig. (1f). For $\delta > 0.16$, only the uniform
SC states are sustained.

Interestingly, $E_{GS}$ for the doping range $\delta = 0.1 - 0.14$ reveal that states with different wavevectors are quasi-degenerate as marked by the shallowing minima. Such energetically close-lying orders open up the possibility of a state with an admixture of different wavevectors. This can lead to an incommensuration of the density wave [55, 56]. We will show in the subsequent sections, that although we start with a single wavevector, higher harmonics of the same wavevector appears in the self-consistent solutions.

Finally, we note that the demise of the modulated orders with doping takes an intriguing pathway. The sharp minima at a particular wavevector at low-dopings indicate that the modulated state with the optimal $Q$ is stable. However, around $\delta = 0.14$ as shown in Fig. (1f), the modulated states with different wavevectors as well as the nematic and the uniform SC states strikingly become quasi-degenerate. In this regime, minute energy differences can tilt the balance of these orders locally. Beyond $\delta > 0.16$, only the uniform d-wave superconducting self-consistent solutions are obtained, indicating the demise of modulated and nematic phases.

## B. Variation of order parameters with doping

Above variational analysis identifies the modulated charge patterns in the ground state, but now we proceed to explore the strength of the corresponding spatially modulated orders. In Fig. (2a), we show the strength of different modulation amplitude at the dominant wavevector defined using Eq. (21-25).

Interestingly, over the range of doping, the finite momentum components of BDW, CDW, and PDW orders are significant for $\delta < 0.14$ as shown in Fig. (2a). However, the d-wave BDW order is the dominant almost over the whole doping range studied. Focusing, on the modulation amplitude for dSC pairing, $\Delta_d^Q$, we observe a non-monotonic behavior with doping. For low doping, the d-wave SC pairing modulations are weaker than other energy scales. However, $\Delta_d^Q$ increase rapidly $\delta = 0.1$. Similarly, there remains a significant extended s-wave pairing modulation at all dopings below $\delta = 0.14$. We emphasize that the wavevector of the modulated orders reduces around this doping as established in the Sec. (III A). Recently, STM experiments observed signatures of unidirectional PDW state with $Q = (1/8)2\pi/a_0$ coexisting with uniform superconducting order in underdoped BSSCO at zero magnetic fields [24, 25].

We have also presented the mean d-wave SC pairing amplitude in Fig. (2b) and compared it with the situation when all the modulated orders are suppressed in $\Delta_d^0$. In the regime where the spatially varying states exist, the average $\Delta_d$ reduces significantly from the uniform case. The renormalized d-wave pairing amplitude $\tilde{\Delta}_d$ also diminishes due to the interplay with the density wave orders.

We also find the signatures of nematicity throughout the doping range $\delta < 0.16$, which is indicated by the non-zero $N_\Delta$ and $N_\chi$ in Fig. (2c). The mean bond density and the Cooper pair density is different along the x and y-directions, thus breaking the $C_4$ rotational symmetry of the square lattice.

For the present calculation, the nematic order coexists only with uniform dSC order in the doping range $\delta = 0.14 - 0.16$. However, it also coexists with all spatially modulating orders in the doping range $\delta = 0.06 - 0.14$ as shown in Fig. (2).

## C. Spatial profile for the optimal solution

In this subsection, we explore the spatial structure of the optimal solution at a few doping levels. STM experiments routinely map the spatial variation of diverse orders. Therefore, we start by concentrating on the low-doping regime where we find a unidirectional BDW, PDW, and CDW state with a periodic modulation of four lattice spacing.

$$1. \text{ CDW, BDW and PDW state with } Q = (1/4, 0)$$

The Fig. (3a) presents the bond density $\chi_{ij}$ in the color plot. As mentioned before, we consider unidirectional modulated states, and therefore, the bond order only modulates along the x-direction while it is uniform along the y-direction. Interestingly, the $\chi_{ij}$ at each site displays different x and y components as shown in Fig. (3a), leading to a charge nematicity coexisting with the density wave orders.

Next, we extract the d-wave and the extended s-wave components of the self-consistent Fock shift using Eq. (19) and Eq. (20). Both the $\chi_d$ and $\chi_s$ modulates in space with a period of four lattice spacings. We also present the variation of the local density $\rho_i$ in the inset of Fig. (3c). We observe spatial modulation of the local density, although it has a much rapid variation instead of regular behavior over a single period.

Fig. (3c) exhibits the Fourier transform of bond and charge density as calculated from Eq. (21-23). Our goal is to capture the dominant oscillations for each component of bond and charge modulations. Hence, we removed the mean to eliminate the peak coming from the uniform values of these quantities. The principal peak for the d-wave BDW and s-wave BDW is observed at the expected value of $Q = (1/4)2\pi/a_0$, and both $\tilde{\chi}_s(Q)$ and $\tilde{\chi}_d(Q)$ have the same strength at this doping. Moreover, there is no other additional peak for both $\tilde{\chi}_d(q)$ and $\tilde{\chi}_s(q)$. The amplitude of the CDW oscillation is much weaker than the BDW ones. Also, there is a sub-dominant peak at $Q = (1/2)2\pi/a_0$, arising from higher-order harmonics of the local density oscillation.

We present the SC pairing $\Delta_{ij}$ on bonds in Fig. (3d) and the extracted d-wave and extended s-wave in Fig. (3e). Notice that for a uniform $tJ$ model, the expectation value for the SC pairing with a d-wave form-factor is finite, but the same for extended s-wave vanishes [50, 51]. Therefore, the $\Delta_s$ pairing modulation transpires around zero, whereas pairing modulation of the $\Delta_d$ occurs around a non-zero mean value [57]. A PDW state is cooperative pairing modulations such that the average pairing vanishes. For the usual pairing modulation due to the interplay of charge or bond density wave order, the mean pairing is significantly larger than the modulation amplitude. Therefore, we observe a PDW oscillation for the extended s-wave pairing, whereas a normal pairing modulation in the d-
<figure><img src="image_5.png" /></figure>

FIG. 3. Spatial profile for the optimal self-consistent solution at $\delta = 0.075$. Here the periodicity of unidirectional oscillation is $4a_0$. (a) Shows the bond density $\chi_{ij} = \sum_{\sigma} \langle c_{i\sigma}^{\dagger} c_{j\sigma} \rangle$, where i, j are nearest neighbouring sites. (b) The d-wave and extended s-wave components of the $\chi$ are extracted in (b). Both clearly show pure oscillations with the periodicity of four lattice spacings. The inset of (c) shows the spatial profile for the local density, which also oscillates. (c) Exhibits the Fourier transform of these quantities. The dominant wavevector for all the oscillations is at $Q = (1/4)2\pi/a_0$ as expected. (d) Shows the similar bond plots for local SC pairing $\Delta_{ij}$. (e) Shows the corresponding d-wave and the extended s-wave component of SC pairing. (f) Illustrates the Fourier transform for the oscillating part of $\Delta_d$ and $\Delta_s$. The pairing modulations also manifest a dominant peak at $Q = (1/4)2\pi/a_0$ at this dopings.

wave channel. While the component of extended s-wave pairing is essentially absent in the uniform superconductor, the presence of the modulated charge orders brings extended s-wave pairing modulation into life. The Fourier transform for the pairing gives an equally strong pairing modulation amplitude for both the form factor with no subdominant peak, as presented in Fig. (3d).

## 2. CDW, BDW and PDW state with $Q = (1/8, 0)$

For $\delta = 0.125$, the primary wavevector of the modulated orders in our microscopic model reduces to $Q = (1/8)2\pi/a_0$. Fig. (4a) shows the bond density and Fig. (4b) the d-wave and s'-wave component of the same. The modulation of the bond density is not smooth over the period and therefore accommodates higher harmonics. This is also true for the density modulations presented as an inset of Fig. (4c). The presence of higher harmonics of oscillations is verified by the strong subdominant peaks of the Fourier transform of the BDW and CDW orders as exhibited in Fig. (4c). Interestingly, the d-wave BDW becomes the dominant form factor at $q = 0.125$ in

reciprocal lattice unit (rlu). However, the $\tilde{\chi}_s$ has an intriguing evolution. For the extended s-wave BDW modulations, the dominant peak shifts at $q = (1/4)2\pi/a_0$, which is at $q = 2Q$.

Similarly, the SC pairing amplitude is presented in Fig. (4d), and the extended s-wave and d-wave components are presented in Fig. (4e). The pairing modulation is smoother contrasted to the BDW or CDW oscillations. In Fig. (4f) the Fourier transform of $\Delta_d$ and $\Delta_s$ confirms the dominant ordering at $q = (1/8)2\pi/a_0$. The secondary peaks are also much weaker than the principal ones.

## 3. Nematic state

Finally, for $\delta = 0.14$, the optimal self-consistent solution supports a $Q = 0$ nematic state. Fig. (5a) shows the bond density amplitude. The color difference indicates the breaking of the rotational symmetry of the square lattice. However, there are no modulations of the bond density either along the x or y-directions. Moreover, the lack of spatially modulated order is confirmed by the SC pairing in Fig. (5b) and the local density in Fig. (5c).
<figure><img src="image_6.png" /><figcaption>FIG. 4. Spatial profile for the optimal self-consistent solution at $\delta = 0.125$. Here the periodicity of unidirectional oscillation is $8a_0$. (a) Shows the bond density $\chi_{ij}$, where i, j are nearest neighboring sites. (b) The d-wave and extended s-wave components of the $\chi$ are extracted in (b). Both show some jerky oscillations with the periodicity of eight lattice spacings. The inset of (c) shows the spatial profile for the local density, which oscillates. (c) Exhibits the Fourier transform of the quantities mentioned above. The dominant wavevector for $\chi_d$ and $\chi_\rho$ are at $Q = (1/8)2\pi/a_0$, whereas the same for $\chi_s$ is at $Q = (1/4)2\pi/a_0$. However, several peaks at the multiples of the ordering wavevector, implies that density modulation is not regular over a period at this doping. (d) Shows the similar bond plots for local SC pairing $\Delta_{ij}$. (e) Shows the corresponding d-wave and the extended s-wave component of SC pairing. (f) Displays the Fourier transform for the oscillating part of $\Delta_d$ and $\Delta_s$. The pairing modulations also exhibit a dominant peak at $Q = (1/8)2\pi/a_0$ at this doping.</figcaption></figure>

# IV. DISCUSSIONS

We study the $t-t'-J$ model with additional nearest-neighbor repulsion and find modulating bond density wave states coexisting with charge and pairing modulated states in the underdoped regime. The modulated states have lower energy than the uniform d-wave superconducting state in the doping range $\delta = 0.06-0.14$. The optimal wavevector for all the modulated quantities is commensurate with the lattice with higher harmonics of oscillations. Over the doping range $\delta = 0.06-0.16$, there remains a significant nematicity such that the $C_4$-rotational symmetry is broken. Moreover, for $\delta = 0.14-0.16$, the system supports a purely nematic state without any translation symmetry breaking. Finally, beyond this doping, only a uniform d-wave superconducting order is stabilized. In the following subsections, we discuss how this work compares with other results in the literature in the context of the cuprates.

## A. Role of strong correlation in stabilizing modulated orders

We attempted to study the modulated states by relaxing the strong correlation constraint by setting the Gutzwiller factors to unity. In the past, studies have been performed for d-wave SC in the $tJ$ model, and this procedure is known as Inhomogeneous Mean Field Theory (IMT) [50, 51]. In sharp contrast to the strongly correlated results presented in this paper, we failed to stabilize any modulated solutions across the parameter regime studied. Of course, an absence of modulated self-consistent solutions might be due to a lack of exhaustive search or a bad initial guess of the order parameters. However, any initial trial state with modulated charge, pairing, and bond densities flowed strictly towards a uniform d-wave SC solution. Furthermore, the ease of achieving modulated self-consistent solutions for multiple wavevectors over a wide doping range in the highly correlated regime indicates a lack of modulated states in the absence of a strong correlation. The propensity to form modulated orders near the half-filling where the interactions are stronger also hints at the same possibility. Recent explorations of the pair density of states in the
<figure><img src="image_7.png" /><figcaption>FIG. 5. Shows the spatial dependence for the optimal self-consistent solution at $\delta = 0.14$. Although the optimal self-consistent solution has no modulation in this parameter regime, it breaks the $C_4$-rotational symmetry of the square lattice leading to a charge nematic state. (a) Shows the bond density $\chi_{ij}$, which is different for the x-bonds and y-bonds. Similarly, the local SC pairing is not purely d-wave as shown in (b) as expected in the uniform tJ-model. Interestingly, there is no spatial variation of these quantities or the local density as presented in (c).</figcaption></figure>

tJ-model also observed a lack of modulated orders when the Gutzwiller factors are set to unity [58].

## B. Evolution of optimal wavevector with doping

In the $t-t'-J$ model with additional nearest-neighbor repulsion, the optimal modulation wavevector reduces with increasing doping. For doping around $\delta = 0.06-0.09$ the optimal wavevector is around $Q = (1/4)2\pi/a_0$. As the doping is increased further, the periodicity increases to $8a_0$ around $\delta = 0.125$ as shown in the inset of Fig. (1f). Generally, in the cuprates, the wavevectors also diminish with increased hole doping [12, 13] except in La-based material [14–16]. However, in $La_2CuO_4$ the additional spin density wave play a vital role in dictating the modulation wavevector [59, 60]. It would be interesting to discern the role of magnetic orders on the modulation wavevector in the future.

Additionally, in cuprates the wavevector varies from $Q = 0.33-0.2 (2\pi/a_0)$ leading to a periodicity of 3–5 lattice spacing up to $\delta = 0.2$. However, this contrasts with theoretical calculations that often predict longer modulation wavelength in the tJ-model and the associated Hubbard model at

large $U$ [35]. The evolution of the ground state energy in Fig. (1) hints towards an explanation of this apparent anomaly. Notice that below $\delta < 0.1$, the difference between the $E_{GS}$ at the optimal wavevector with other wavevector is significant. Note that we only studied the commensurate wavelengths with the lattice spacings. However, the energy difference between the states with different wavelengths becomes minute, $\delta E \sim 0.001t$ as doping increases beyond $\delta > 0.1$. Since these states are near degenerate, various fluctuations absent in our study can modify the delicate balance of energy. Consequently, the state with a smaller wavelength can become the ground state.

Additionally, there can be another important consequence for near-degeneracy of modulated states with different wavevectors as shown in Fig. (1c)-Fig. (1f). This can possibly support a state with an admixture of different wavevectors. Furthermore, the irregular evolution of the $\chi_d$ and $\chi_s$ over a period in Fig. (4b) also hints towards the incipient incommensuration. Such patterns generates strong higher harmonics in the Fourier transform of the BDW order in Fig. (4c). A splitting or broadening of the primary peak of the $\tilde{\chi}(q)$ indicates an incommensuration of any density wave pattern. Here, strong higher harmonics with the presence of nearby quasi-degenerate states can drive a multi-peak incommensuration around $\delta \approx 0.125$.

## C. Nearest neighbor repulsion term in the tJ-model

The $t-J$ model and its “parent” – Hubbard model with large on-site repulsion $U$ are generally considered the minimal model for cuprate superconductors. However, the interactions are not entirely screened in realistic materials to be completely local. A residual nearest-neighbor repulsion is estimated to be one order of magnitude smaller than the on-site repulsion in cuprates [61–63]. Such interactions are often neglected due to their reduced magnitude. However, they become crucial when the ground state without them shows near-degeneracy among several broken symmetry states, as in the underdoped regime of the $t-J$ model.

The modulated orders in the $t-t'-J$ model have been studied extensively in the past [38, 44, 58, 64]. These studies found different modulating bond, charge, and pair density wave states energetically very close yet higher than the uniform superconducting ground state. Therefore, additional ingredients absent in the $t-t'-J$ must be responsible for stabilizing the spatially modulating density wave orders as the ground state in cuprates. Ref. [44] suggests studying a tJ-model with additional on-site repulsion stabilizing nematic phases with superconductivity. There are other studies suggesting that adding a nearest-neighbor repulsion term to prefer the density wave orders as the ground state [40, 42]. In particular, the nearest neighbor repulsion is expected to be detrimental to the d-wave superconductivity and can thus stabilize the nearby density waves orders as the ground state. Similar to these studies, we also observed that a critical value of nearest neighbor repulsion is essential to stabilize the density wave orders [40, 42].
Furthermore, increasing the nearest neighbor repulsion term keeping the other terms constant allows for the modulated orders to survive up to higher dopings. We have confirmed that for $W = 0.75t$, the spatially modulated orders remain the ground state up to $\delta = 0.16$ whereas the ground state becomes purely nematic in $\delta = 0.16 - 0.19$. The energy difference between nematic and the uniform state remains significantly more prominent for larger $W$ leading to unambiguous identification of the nematic regime.

### D. Pairing modulation

We observe modulations of the Cooper pair density both in the extended $s$-wave and the $d$-wave channel over the doping range $\delta = 0.06 - 0.14$. A modulation of Cooper pairing can be a direct consequence of the coexistence with the BDW order. On one hand, the modulation of the bond density leads to an oscillation of the local density of states which dictates the propensity for developing local SC pair and drives the modulation of the pairing. On the other hand, a more exotic pair density wave is a state of spatially oscillating Cooper pair such that the average pairing vanishes or is very small. In our analysis, the $\Delta_s$ shows a vanishingly small mean with a robust modulation amplitude over the whole doping range as presented in Fig. (3e) and Fig. (4e). On the other hand, the $d$-wave SC pairing modulates around a finite non-zero value. Consequently, this indicates a PDW with $s'$ form factor coexisting with other modulated states over a wide doping range.

Recent experiments on BSSCO [24, 25] claims to observe PDW modulations at $Q_p = (1/8)2\pi/a_0$ of $\sim 6$ meV coexisting with a uniform SC at $\sim 36$ meV. The observed PDW state has either extended $s$-wave or $s$-wave form factor. These studies also observe CDW oscillations at $Q_c = Q_p$ and at $Q_c = 2Q_p$. The model calculations around $\delta = 0.125$ reveal a similar behavior – an extended $s$-wave PDW undulation with a non-zero mean $d$-wave superconductor. Furthermore, the optimal self-consistent solution displays a PDW order $q = Q$ with the CDW orders at both $q = Q$ and $q = 2Q$, very similar to what we observe in Fig. (4c). It is encouraging that our simplified model for modulated orders in strongly correlated systems can capture some features of the experimental signatures.

### E. Vestigial nematic order

Our model calculation reveals that a translationally invariant nematic state stabilizes in the doping range $\delta = 0.14 - 0.16$. Note that the nematic state is a higher order component of the superconducting and bond order parameter. If the bond order $\chi$ is written as a vector then the nematic order is given by $N_\chi = 1/2(\langle\chi^\dagger\sigma_z\chi\rangle$, where $\sigma_z$ is the z-component of the Pauli matrix. In the nematic state the combined order parameter is finite although the $\langle\chi_x^Q\rangle = 0$ and $\langle\chi_y^Q\rangle = 0$. Consequently, nematic state breaks only a subset of the all possible symmetry that can be broken in this regime leading to a vestigial nematic state.

<figure><img src="image_8.png" /><figcaption>FIG. 6. Presents a schematic symmetry broken phases with doping from our model calculations at zero temperature. The doping regime presented below the broken line is beyond the scope of this study. Below this doping we expect magnetic and other symmetries to be broken. From the low to intermediate dopings, both the translation symmetry, $\mathcal{T}$ and the $C_4$ rotational symmetry is broken, leading to the observation of the bond, charge, and pair density wave orders. Beyond a doping level, only $C_4$ rotational symmetry is broken, leading to a nematic state. Note that fewer symmetries are broken as the doping is increased. Beyond critical doping, only translational and rotational invariant $d$-wave SC remains robust. Uniform superconductors only break the particle number conservation symmetry. The uniform $d$-wave SC remains robust over the whole doping regime studied, only decaying at much higher doping.</figcaption></figure>

Recently sublattice resolved electronic structure experiments in cuprates [8] reveal a $Q = 0$ nematic state over the doping range $\delta = 0.06 - 0.17$. The vestigial nematic state coexists with finite-$Q$ $d$-wave charge order and a uniform $d$-wave in the intermediate dopings. A pure nematic state without any translational symmetry breaking is also found in the doping range $\delta = 0.15 - 0.17$. We note that our optimal solution also sustains a pure $C_4$ rotational symmetry breaking order around doping. Beyond this doping, only uniform $d$-wave superconductivity stabilizes in our model calculation.

### F. Phase diagram

In Fig. (6), we provide a schematic map of different symmetry breaking orders based on our model calculations at zero temperature. The doping regime presented below the broken line is beyond the scope of this study. Uniform superconductors only break the particle number conservation symmetry. At zero temperature, our model calculation suggests that uniform $d$-wave SC remains robust over the whole doping regime studied, only decaying at much higher doping. In contrast, the finite-$Q$ density wave states become relevant only in the lightly doped regime. The density wave order breaks both the translation, $\mathcal{T}$ and the $C_4$ rotational symmetry of the square lattice. Note that Fig. (2b) shows that the mean $d$-wave pairing amplitude reduces drastically in the regime where other symmetries are broken. Consequently, the suppression of SC pairing is expected from our calculations at $T = 0$. A similar drop in superconducting correlations is observed in cuprates
around the hole-doping of 12%, where the amplitude of the modulated orders is strong [65, 66].

Between the finite-Q orders and the uniform SC, there is a regime where the rotational symmetry remains broken although the translational invariance is restored. Such a state is a vestigial nematic state which breaks a smaller subset of symmetries than the primary state of density wave orders. The series of cascading broken symmetries as the doping reduces hints at the possible entanglement between multiple orders near the highly correlated regime.

The charge density waves [67] and the associated vestigial nematicity [6–8] are considered an essential ingredient of the pseudogap state. Our study suggests that ground state of cuprates in the underdoped regime indeed breaks certain symmetries generating various density wave states, which manifests in experiments as the pseudogap state in cuprates at finite temperature. Therefore, to capture the pseudogap state, it is

essential to analyze the effects of different fluctuations on all different broken symmetry states. Furthermore, thermal and quantum fluctuations can lead to spatial reorganization into domains of finite-Q and Q = 0 orders [68–70]. Additionally, the Mott insulator and magnetically ordered states can significantly alter the states at low doping. A rigorous treatment of the magnetic orders with these broken symmetry states, along with the effects of thermal fluctuations, remains to be explored in the future.

V. ACKNOWLEDGEMENT

The authors thank Yvan Sidis, Maxence Grandadam, and Debmalya Chakraborty for valuable discussions. This work has received financial support from the ERC, under grant agreement AdG694651-CHAMPAGNE.

[1] B. Keimer and J. Moore, Nature Physics 13, 1045 (2017).

[2] J. P. Sun, K. Matsuura, G. Z. Ye, Y. Mizukami, M. Shimozawa, K. Matsubayashi, M. Yamashita, T. Watashige, S. Kasahara, Y. Matsuda, J.-Q. Yan, B. C. Sales, Y. Uwatoko, J.-G. Cheng, and T. Shibauchi, Nature Communications 7, 12146 (2016).

[3] E. Dagotto, Science 309, 257 (2005).

[4] D. Bounoua, L. Mangin-Thro, J. Jeong, R. Saint-Martin, L. Pinsard-Gaudart, Y. Sidis, and P. Bourges, Communications Physics 3, 1 (2020).

[5] E. Fradkin, S. A. Kivelson, and J. M. Tranquada, Rev. Mod. Phys. 87, 457 (2015).

[6] R. M. Fernandes, P. P. Orth, and J. Schmalian, Annual Review of Condensed Matter Physics 10, 133 (2019).

[7] L. Nie, G. Tarjus, and S. A. Kivelson, Proc. Natl. Acad. Sci. 111, 7980 (2014).

[8] S. Mukhopadhyay, R. Sharma, C. K. Kim, S. D. Edkins, M. H. Hamidian, H. Eisaki, S.-i. Uchida, E.-A. Kim, M. J. Lawler, A. P. Mackenzie, et al., Proceedings of the National Academy of Sciences 116, 13249 (2019).

[9] J. E. Hoffman, E. W. Hudson, K. M. Lang, V. Madhavan, H. Eisaki, S. Uchida, and J. C. Davis, Science 295, 466 (2002).

[10] W. D. Wise, M. C. Boyer, K. Chatterjee, T. Kondo, T. Takeuchi, H. Ikuta, Y. Wang, and E. W. Hudson, Nat. Phys. 4, 696 (2008).

[11] A. Frano, S. Blanco-Canosa, B. Keimer, and R. J. Birgeneau, Journal of Physics: Condensed Matter 32, 374005 (2020).

[12] R. Comin, A. Frano, M. M. Yee, Y. Yoshida, H. Eisaki, E. Schierle, E. Weschke, R. Sutarto, F. He, A. Soumyanarayanan, Y. He, M. Le Tacon, I. S. Elfimov, J. E. Hoffman, G. A. Sawatzky, B. Keimer, and A. Damascelli, Science 343, 390 (2014).

[13] E. H. da Silva Neto, P. Aynajian, A. Frano, R. Comin, E. Schierle, E. Weschke, A. Gyenis, J. Wen, J. Schneeloch, Z. Xu, S. Ono, G. Gu, M. Le Tacon, and A. Yazdani, Science 343, 393 (2014).

[14] J. Fink, E. Schierle, E. Weschke, J. Geck, D. Hawthorn, V. Soltwisch, H. Wadati, H.-H. Wu, H. A. Dürr, N. Wizent, B. Büchner, and G. A. Sawatzky, Phys. Rev. B 79, 100502(R) (2009).

[15] J. Fink, V. Soltwisch, J. Geck, E. Schierle, E. Weschke, and B. Büchner, Phys. Rev. B 83, 092503 (2011).

[16] M. Hücker, M. v. Zimmermann, G. D. Gu, Z. J. Xu, J. S. Wen, G. Xu, H. J. Kang, A. Zheludev, and J. M. Tranquada, Phys. Rev. B 83, 104506 (2011).

[17] E. M. Forgan, E. Blackburn, A. T. Holmes, A. K. R. Briffa, J. Chang, L. Bouchenoire, S. D. Brown, R. Liang, D. Bonn, W. N. Hardy, N. B. Christensen, M. V. Zimmermann, M. Hücker, and S. M. Hayden, Nature Communications 6, 10064 (2015).

[18] R. Comin, R. Sutarto, F. He, E. H. da Silva Neto, L. Chauviere, A. Frano, R. Liang, W. N. Hardy, D. A. Bonn, Y. Yoshida, H. Eisaki, A. J. Achkar, D. G. Hawthorn, B. Keimer, G. A. Sawatzky, and A. Damascelli, Nature Materials 14, 796 (2015).

[19] R. Comin, R. Sutarto, E. H. da Silva Neto, L. Chauviere, R. Liang, W. N. Hardy, D. A. Bonn, F. He, G. A. Sawatzky, and A. Damascelli, Science 347, 1335 (2015).

[20] G. Campi, A. Bianconi, N. Poccia, G. Bianconi, L. Barba, G. Arrighetti, D. Innocenti, J. Karpinski, N. D. Zhigadlo, S. M. Kazakov, M. Burghammer, M. V. Zimmermann, M. Sprung, and A. Ricci, Nature (London) 525, 359 (2015).

[21] Y. Kohsaka, C. Taylor, K. Fujita, A. Schmidt, C. Lupien, T. Hanaguri, M. Azuma, M. Takano, H. Eisaki, H. Takagi, S. Uchida, and J. C. Davis, Science 315, 1380 (2007).

[22] C. V. C. Parker, P. P. Aynajian, E. H. E. da Silva Neto, A. A. Pushp, S. S. Ono, J. J. Wen, Z. Z. Xu, G. G. Gu, and A. A. Yazdani, Nature 468, 677 (2010).

[23] M. K. Chan, N. Harrison, R. D. McDonald, B. J. Ramshaw, K. A. Modic, N. Barišić, and M. Greven, Nature Communications 7, 12244 (2016).

[24] M. H. Hamidian, S. D. Edkins, S. H. Joo, A. Kostin, H. Eisaki, S. Uchida, M. J. Lawler, E.-A. Kim, A. P. Mackenzie, K. Fujita, J. Lee, and J. C. S. Davis, Nature 532, 343 (2016).

[25] Z. Du, H. Li, S. H. Joo, E. P. Donoway, J. Lee, J. S. Davis, G. Gu, P. D. Johnson, and K. Fujita, Nature 580, 65 (2020).

[26] Z. Dai, T. Senthil, and P. A. Lee, Phys. Rev. B 101, 064502 (2020).

[27] D. Chakraborty, M. Grandadam, M. H. Hamidian, J. C. S. Davis, Y. Sidis, and C. Pépin, Phys. Rev. B 100, 224511 (2019).

[28] M. Grandadam, D. Chakraborty, and C. Pépin, Journal of Superconductivity and Novel Magnetism 33, 2361 (2020).

[29] P. A. Lee, Phys. Rev. X 4, 031017 (2014).
[30] X. Y. Xu, K. T. Law, and P. A. Lee, Physical review letters **122**, 167001 (2019).

[31] E. Berg, E. Fradkin, and S. A. Kivelson, Physical review letters **105**, 146403 (2010).

[32] C. Peng, Y.-F. Jiang, T. P. Devereaux, and H.-C. Jiang, *npj Quantum Materials* **6**, 64 (2021).

[33] S. S. Dash and D. Sénéchal, Physical Review B **103**, 045142 (2021).

[34] S. Sachdev and R. La Placa, Phys. Rev. Lett. **111**, 027202 (2013).

[35] B.-X. Zheng, C.-M. Chung, P. Corboz, G. Ehlers, M.-P. Qin, R. M. Noack, H. Shi, S. R. White, S. Zhang, and G. K.-L. Chan, Science **358**, 1155 (2017).

[36] A. Himeda, T. Kato, and M. Ogata, Phys. Rev. Lett. **88**, 117001 (2002).

[37] P. Corboz, T. M. Rice, and M. Troyer, Phys. Rev. Lett. **113**, 046402 (2014).

[38] P. Choubey, W.-L. Tu, T.-K. Lee, and P. Hirschfeld, New Journal of Physics **19**, 013028 (2017).

[39] P. Choubey, S. H. Joo, K. Fujita, Z. Du, S. Edkins, M. Hamidian, H. Eisaki, S. Uchida, A. Mackenzie, J. Lee, et al., Proceedings of the National Academy of Sciences **117**, 14805 (2020).

[40] A. Allais, J. Bauer, and S. Sachdev, Phys. Rev. B **90**, 155114 (2014).

[41] A. Allais, J. Bauer, and S. Sachdev, Indian J. Phys. **88**, 905 (2014).

[42] J. Spałek, M. Zegrodnik, and J. Kaczmarczyk, Physical Review B **95**, 024506 (2017).

[43] J. D. Sau and S. Sachdev, Physical Review B **89**, 075129 (2014).

[44] M. Zegrodnik and J. Spałek, Phys. Rev. B **98**, 155144 (2018).

[45] H. Zhao, Z. Ren, B. Rachmilowitz, J. Schneeloch, R. Zhong, G. Gu, Z. Wang, and I. Zeljkovic, Nature materials **18**, 103 (2019).

[46] K. Chao, J. Spałek, and A. Oles, Journal of Physics C: Solid State Physics **10**, L271 (1977).

[47] P. Fazekas, *Lecture notes on electron correlation and magnetism*, Vol. 5 (World Scientific, 1999).

[48] N. Fukushima, Physical Review B **78**, 115105 (2008).

[49] W.-H. Ko, C. P. Nave, and P. A. Lee, Physical Review B **76**, 245113 (2007).

[50] D. Chakraborty and A. Ghosal, New Journal of Physics **16**, 103018 (2014).

[51] A. Garg, M. Randeria, and N. Trivedi, Nature Physics **4**, 762 (2008).

[52] P.-G. De Gennes and P. A. Pincus, *Superconductivity of metals and alloys* (CRC Press, 2018).

[53] A. Ghosal, M. Randeria, and N. Trivedi, Phys. Rev. Lett. **81**, 3940 (1998).

[54] A. Ghosal, M. Randeria, and N. Trivedi, Phys. Rev. B **65**, 014501 (2001).

[55] M. Vershinin, Science **303**, 1995 (2004).

[56] I. El Baggari, B. H. Savitzky, A. S. Admasu, J. Kim, S.-W. Cheong, R. Hovden, and L. F. Kourkoutis, Proceedings of the National Academy of Sciences **115**, 1445 (2018).

[57] Note that the presence of the weak nematicity generates a small value of mean $\Delta_s$.

[58] W.-L. Tu and T.-K. Lee, Scientific Reports **9**, 1719 (2019).

[59] L. Nie, A. V. Maharaj, E. Fradkin, and S. A. Kivelson, Phys. Rev. B **96**, 085142 (2017).

[60] H. Miao, J. Lorenzana, G. Seibold, Y. Peng, A. Amorese, F. Yakhou-Harris, K. Kummer, N. Brookes, R. Konik, V. Thampy, et al., Proceedings of the National Academy of Sciences **114**, 12430 (2017).

[61] A. Reymbaut, M. Charlebois, M. F. Asiani, L. Fratino, P. Sémon, G. Sordi, and A.-M. Tremblay, Physical Review B **94**, 155146 (2016).

[62] D. Sénéchal, A. G. R. Day, V. Bouliane, and A.-M. S. Tremblay, Phys. Rev. B **87**, 075123 (2013).

[63] M. Jiang, U. R. Hähner, T. C. Schulthess, and T. A. Maier, Physical Review B **97**, 184507 (2018).

[64] K.-Y. Yang, W. Q. Chen, T. M. Rice, M. Sigrist, and F.-C. Zhang, New Journal of Physics **11**, 055053 (2009).

[65] M. Hücker, N. B. Christensen, A. T. Holmes, E. Blackburn, E. M. Forgan, R. Liang, D. A. Bonn, W. N. Hardy, O. Gutowski, M. v. Zimmermann, S. M. Hayden, and J. Chang, Phys. Rev. B **90**, 054514 (2014).

[66] S. Blanco-Canosa, A. Frano, E. Schierle, J. Porras, T. Loew, M. Minola, M. Bluschke, E. Weschke, B. Keimer, and M. Le Tacon, Phys. Rev. B **90**, 054513 (2014).

[67] B. Loret, N. Auvray, Y. Gallais, M. Cazayous, A. Forget, D. Colson, M.-H. Julien, I. Paul, M. Civelli, and A. Sacuto, Nature Physics **15**, 771 (2019).

[68] S. Lee, J. Collini, S. X.-L. Sun, M. Mitrano, X. Guo, C. Eckberg, J. Paglione, E. Fradkin, and P. Abbamonte, Phys. Rev. Lett. **127**, 027602 (2021).

[69] A. Banerjee, A. Garg, and A. Ghosal, Phys. Rev. B **98**, 104206 (2018).

[70] D. Agterberg and H. Tsunetsugu, Nature Physics **4**, 639 (2008).