Exact analysis of the interplay of charge order and unconventional pairings in the 2D Hatsugai-Kohmoto model

Abstract

We provide here a study of some competing ordering tendencies exhibited by the exactly solvable 2D Hatsugai–Kohmoto (HK) model on a square lattice. To this end, we investigate the interplay between superconductivity, charge-density wave (CDW) and pair-density wave (PDW) orders as a function of interaction, doping parameter, magnetic field, and uniaxial strain. As a result, we confirm the intertwined nature of CDW and PDW fluctuating orders for intermediate-to-strong couplings. We also verify that, while an applied magnetic field favors the formation of a CDW and allows the subsequent emergence of a PDW as a secondary order, strain effects favor unidirectional PDW as a primary order over the subdominant appearance of a stripe-like CDW. These results underscore the value of the HK model as an interesting platform in order to investigate (via an exactly solvable framework) the emergence of charge order and unconventional superconductivity in fermionic systems with strong interactions. Finally, we briefly discuss an orbital generalization of the HK model, which has been recently argued to be relevant to describe the properties of realistic strongly correlated systems.

I Introduction

The mechanism of high-temperature superconductivity in the cuprates is probably one of the most profound enigmas in condensed matter physics [1] and remains largely uncomprehended to this date. In this respect, it is fair to say that there is a reasonably broad consensus in the community that a minimal model that potentially describes the physics of these materials microscopically is the 2D Hubbard model [2, 3] on a square lattice with the on-site interaction $U$ being of the order of the bandwidth $W$. On the other hand, this intermediate interaction regime of the model is widely acknowledged to be one of the hardest problems to be addressed theoretically (both analytically and numerically) [4, 5] and, for this reason, simpler models that are easier to solve and can potentially capture some important aspects of the physics of these compounds are extremely interesting in the field. One such model is the so-called 2D Hatsugai-Kohmoto (HK) model, proposed in the early 1990s [6] (a similar model was proposed in Ref. [7]), which is an exactly solvable interacting fermionic model that describes both a Mott insulating phase at half-filling and a non-Fermi liquid (NFL) metallic phase that emerges upon doping [8, 9, 10]. It was later demonstrated that this NFL metallic phase also turns out to be unstable towards the formation of superconductivity at low temperatures [11, 12], which is qualitatively consistent with the behavior observed, e.g., in the cuprate compounds [2, 13, 14], among other systems. For this reason, our main aim in this paper will be to continue investigating some additional properties of this model in an attempt to potentially gain some insight into the remarkable physics of such strongly correlated materials.

An important aspect that we will be especially concerned with here will be the interplay of some competing electronic phases that appear in the 2D HK model at low temperatures from an exact point of view. As will become clear shortly, we will obtain that in addition to the Mott insulating and NFL phases [6, 11, 12], this model can describe a variety of other ordering tendencies such as, e.g., unidirectional and bidirectional charge-density waves (CDW), spin-singlet and spin-triplet superconductivity (SC), and also pair-density-wave (PDW) phases, which demonstrate the richness of the underlying physics already captured by this simple effective model. However, we point out that not all these emergent phases will necessarily develop long-range order at low temperatures in the model, but the main message that we want to convey here is that we will be able to identify the intertwined nature of some ordering tendencies and also obtain a hierarchy between those phases that may give some insight of what might happen in more realistic microscopic descriptions of some strongly correlated superconductors.

Furthermore, we will study, e.g., the effects of uniaxial strain and applied magnetic field in the 2D HK model as a means to enhance (or suppress) specific phases that appear in the model in order to disentangle the different ordering tendencies that turn out to be very close in energy in the corresponding phase diagram. Altering the non-interacting band structure will also be explored here in order to understand how robust are some emergent phases described here with respect to these small modifications in the model. Finally, we will discuss a possible generalization of the HK model, which might eventually be of relevance in order to compare its low-temperature physics with more realistic strongly correlated systems, such as, e.g., the 2D Hubbard model.

Therefore, this paper will be organized as follows. In Sec. II, we will define the 2D HK Hamiltonian. In this part, we will be very concise and will refrain from discussing all previous results obtained for this model available in the literature (for this reason, we also refer here to some previous papers that discuss several aspects of its properties - see, e.g., [6, 11, 12, 15, 16, 17, 18]). We point out that our present paper represents a significant extension of our earlier work [17] in that we now investigate charge order with a modulation $\mathbf{q}\neq 0$ associated with $s$- and $d$-wave symmetries and also pairing orders (i.e., both SC and PDW) associated with $p$- and $d$-wave symmetries and, in addition, we provide a comprehensive analysis of several variations of the HK model such as including an applied external magnetic field, uniaxial strain effects and modifications of the band structure. In Sec. III, we will discuss the methodology that will later be used to determine the corresponding phase diagrams of the model. In Sec. IV, we will present the results obtained within our present approach. In Sec. V, we end with a summary of our work and also provide an outlook on possible future directions. Finally, in Appendix A, we provide a demonstration of the susceptibility of the CDW phases that emerge in the 2D HK model, while, in Appendix B, we briefly discuss an orbital generalization of the HK model.

II HK Model

Our starting point is the 2D HK Hamiltonian [6, 11] (for a recent review about this model, see also Ref. [18]) defined on a square lattice with an additional Zeeman term, which is given by:

$$H_{HK} =\sum_{\mathbf{k},\sigma}\left(\xi_{\mathbf{k}}-\sigma B\right)c_{\mathbf{k}\sigma}^{\dagger}c_{\mathbf{k}\sigma}+U\sum_{\mathbf{k}}n_{\mathbf{k}\uparrow}n_{\mathbf{k}\downarrow},$$ (1)

where $c^{\dagger}_{\mathbf{k}\sigma}$ and $c_{\mathbf{k}\sigma}$ are, respectively, the fermionic creation and annihilation operators for electrons with momentum $\mathbf{k}$ and spin projection $\sigma=\uparrow,\downarrow$. The term $\xi_{\mathbf{k}}=-2(t_{x}\cos k_{x}+t_{y}\cos k_{y})-\mu$ describes the single-particle band dispersion, with $t_{x}=t-\delta$ and $t_{y}=t+\delta$ denoting the hopping amplitudes along the $x$- and $y$-directions, respectively. We point out that these hopping amplitudes are allowed to be modified by the most direct effect of the strain $\delta$ (which is the increase in the orbital overlap integral by compression, generating anisotropy in the hopping in $x$ and $y$ directions [19]) and $\mu$ is the chemical potential. The term proportional to $B$ corresponds to the Zeeman term due to an external magnetic field. The number operator is defined as $n_{\mathbf{k}\sigma}=c^{\dagger}_{\mathbf{k}\sigma}c_{\mathbf{k}\sigma}$, and $U>0$ denotes the local repulsive interaction in momentum space. All external perturbations (magnetic field, strain, etc.) included in the 2D HK model preserve the exact solvability of the model. This enables us to investigate these modifications of such a strongly correlated model while maintaining analytical tractability. In order to analyze the emergence of some correlated phases in this model, we will add to Eq. (1) a coupling in the charge-density wave and pairing channels encapsulated in the general term denoted by $H_{V}=H_{SC}+H_{CDW}$, whose terms in the r.h.s. of the equation will be explained below. Here, we will focus only on the aspects of the physics of the model that will be our main concern in this work, since other properties and additional technical details can also be found in several papers available in the literature (see, e.g., Refs. [6, 11, 12, 15, 16, 18]) and also in a previous paper by the authors [17].

There are various possible microscopic mechanisms for Cooper pairing (phonon-mediated, spin-fluctuation-mediated, etc.). In what follows, we will be agnostic with regard to the origin of pairing and will simply assume an attractive BCS-type interaction [20, 21], which is written in momentum space as

$$H_{SC}=-V_{p}\sum_{{\mathbf{k},\mathbf{k}^{\prime},\sigma,\sigma^{\prime}}}\gamma_{\mathbf{k}}\gamma_{\mathbf{k^{\prime}}}\,c^{\dagger}_{\mathbf{k}+\mathbf{q}/2,\sigma}c^{\dagger}_{\mathbf{-k}+\mathbf{q}/2,\sigma^{\prime}}c_{\mathbf{-k^{\prime}}+\mathbf{q}/2,\sigma^{\prime}}c_{\mathbf{k^{\prime}}+\mathbf{q}/2,\sigma},$$ (2)

where $V_{p}>0$ (we point out that we also consider here another simplifying assumption that the pairing interaction for both spin-singlet and spin-triplet channels is equal). For spin-singlet pairing, the pair operator associated with momentum $\mathbf{q}$ – which can be zero for a uniform SC phase or finite for a PDW phase – is defined as $\Delta_{\mathbf{k}}(\mathbf{q})=\gamma_{\mathbf{k}}\,c_{-\mathbf{k}+\frac{\mathbf{q}}{2},\downarrow}\,c_{\mathbf{k}+\frac{\mathbf{q}}{2},\uparrow}$, with $\gamma_{\mathbf{k}}$ being a form factor that encodes even-parity symmetry ($s$- and $d$-wave) [22]. For spin-triplet pairing, the corresponding pair operator is defined as $\Delta_{\mathbf{k}}(\mathbf{q})=\gamma_{\mathbf{k}}\,c_{-\mathbf{k}+\frac{\mathbf{q}}{2},\uparrow}\,c_{\mathbf{k}+\frac{\mathbf{q}}{2},\uparrow}$, with $\gamma_{\mathbf{k}}$ being an odd-parity ($p$-wave) form factor [23]. Therefore, we will consider the following form factors associated with each of the above-mentioned symmetries:

$$\gamma_{\mathbf{k}}=\begin{cases}1,&\text{($s$-wave)};\\ \cos k_{x}-\cos k_{y},&\text{($d$-wave)};\\ \sin k_{x}+i\sin k_{y},&\text{($p$-wave)}.\end{cases}$$ (3)

By following a similar strategy explained before, we will also assume an attractive coupling in the CDW channel, which is given as follows:

$$H_{\text{CDW}}=-V_{c}\sum_{\mathbf{k},\mathbf{k^{\prime}},\sigma,\sigma^{\prime}}\gamma_{\mathbf{k}}\gamma_{\mathbf{k^{\prime}}}\,{c}^{\dagger}_{\mathbf{k}+{\mathbf{q}/2},\sigma}\,{c}^{\dagger}_{\mathbf{k^{\prime}}-\mathbf{q}/2,\sigma^{\prime}}\,{c}_{\mathbf{k^{\prime}}+{\mathbf{q}/2},\sigma^{\prime}}\,{c}_{\mathbf{k}-\mathbf{q}/2,\sigma},$$ (4)

where $V_{c}>0$ (for a fixed $\mathbf{q}\neq 0$), and the CDW operator associated with a modulation $\mathbf{q}$ is defined as $\rho_{\mathbf{k}}(\mathbf{q})=\gamma_{\mathbf{k}}\sum_{\sigma}{c}^{\dagger}_{\mathbf{k}+{\mathbf{q}}/{2},\sigma}\,{c}_{\mathbf{k}-{\mathbf{q}}/{2},\sigma}$ [we will consider here both $s$- and $d$-wave symmetries for the CDW, whose form factors are defined in Eq. (3)].

The exact Green’s function in the Hatsugai–Kohmoto model can be written in closed form as [6, 11]

$$G_{\mathbf{k},\sigma}(i\omega) =\frac{1-\langle n_{\mathbf{k},\bar{\sigma}}\rangle}{i\omega-\xi_{\mathbf{k},\sigma}}+\frac{\langle n_{\mathbf{k},\bar{\sigma}}\rangle}{i\omega-\xi_{\mathbf{k},\sigma}-U},$$ (5)

where $\xi_{\mathbf{k},\sigma}=\xi_{\mathbf{k}}-\sigma B$ and $\langle n_{\mathbf{k},\bar{\sigma}}\rangle$ denotes the average occupation of momentum $\mathbf{k}$ with opposite spin $\bar{\sigma}=-\sigma$. This occupation is given by the expression:

$$\langle n_{\mathbf{k},\sigma}\rangle=\frac{e^{-\beta\xi_{\mathbf{k},\sigma}}+e^{-\beta(2\xi_{\mathbf{k}}+U)}}{1+e^{-\beta\xi_{\mathbf{k},\uparrow}}+e^{-\beta\xi_{\mathbf{k},\downarrow}}+e^{-\beta(2\xi_{\mathbf{k}}+U)}}.$$ (6)

Let us analyze the above result in the absence of external perturbations in the zero-temperature limit. In this case, the average occupation is given by the following piecewise expression [6, 11, 12]:

$$\langle n_{\mathbf{k},\sigma}\rangle=\begin{cases}0,&\mu<\epsilon_{\mathbf{k}}<\frac{W}{2}\hskip 42.67912pt(\mathbf{k}\in\Omega_{0}),\\ \frac{1}{2},&\mu-U<\epsilon_{\mathbf{k}}<\mu\hskip 27.5992pt(\mathbf{k}\in\Omega_{1}),\\ 1,&-\frac{W}{2}<\epsilon_{\mathbf{k}}<\mu-U\hskip 15.36429pt(\mathbf{k}\in\Omega_{2}),\end{cases}$$ (7)

where $W$ is the bandwidth of the model and the momentum regions $\Omega_{0}$, $\Omega_{1}$, and $\Omega_{2}$ correspond, respectively, to empty, single-particle-occupied, and double-particle-occupied momentum states. We note that due to the existence of the $\Omega_{1}$ region, the ground state displays a massive degeneracy, which is lifted by the external perturbation denoted by $H_{V}$. Finally, in order to connect the above division of the momentum space with the hole doping parameter $x$ in the 2D HK model, we obtain by demanding that the $\Omega_{1}$ region lies within the bandwidth the following conditions at $T=0$ [17]:

$$x=\begin{cases}\frac{1}{2}-\tilde{\mu},&\text{if }\tilde{\mu}-u<-\frac{1}{2}\text{ and }-\frac{1}{2}<\tilde{\mu}\leq\frac{1}{2},\\ u-2\tilde{\mu},&\text{if }\tilde{\mu}-u>-\frac{1}{2}\text{ and }-\frac{1}{2}<\tilde{\mu}\leq\frac{1}{2},\end{cases}$$ (8)

where we defined the dimensionless quantities $\tilde{\mu}=\mu/W$ and $u=U/W$.

III Methodology

To determine the ordering tendencies in the 2D HK model, we begin by computing several order-parameter susceptibilities. The pair susceptibility for the spin-singlet case is defined as

$$\chi_{p,(s)}^{(0)}(i\nu_{n},\mathbf{q}) =\int_{0}^{\beta}d\tau\,e^{i\nu_{n}\tau}\big\langle T\,\Delta_{\mathbf{q}}(\tau)\,\Delta_{\mathbf{q}}^{\dagger}(0)\big\rangle_{0}$$

$$=-\frac{1}{\beta}\sum_{\omega_{m}}\sum_{\mathbf{k}}\gamma_{\mathbf{k}}\gamma_{\mathbf{k}}\,G_{\mathbf{k}+\frac{\mathbf{q}}{2},\,\uparrow}(i\omega_{m})\,G_{-\mathbf{k}+\frac{\mathbf{q}}{2},\,\downarrow}(i\nu_{n}-i\omega_{m}),$$ (9)

where the form factors can be either $s$- or $d$-wave. The product of the Green’s functions yields the expression for the “bare” pair susceptibility:

$$\chi_{p,(s)}^{(0)}(i\nu_{n},\mathbf{q})=\sum_{\mathbf{k},a,b}\gamma_{\mathbf{k}}\gamma_{\mathbf{k}}\,n_{\mathbf{k}+\frac{\mathbf{q}}{2},\uparrow}^{a}\;n_{-\mathbf{k}+\frac{\mathbf{q}}{2},\downarrow}^{b}\;\frac{f\!\bigl(\xi_{\mathbf{k}+\frac{\mathbf{q}}{2},\uparrow}^{a}\bigr)+f\!\bigl(\xi_{-\mathbf{k}+\frac{\mathbf{q}}{2},\downarrow}^{b}\bigr)-1}{i\nu_{n}-\xi_{\mathbf{k}+\frac{\mathbf{q}}{2}}^{a}-\xi_{-\mathbf{k}+\frac{\mathbf{q}}{2}}^{b}},$$ (10)

where $a$ and $b$ label the upper (denoted by $u$) and lower (denoted by $l$) bands, defined by $\xi_{\mathbf{k}}^{u}=\xi_{\mathbf{k}}+U$ and $\xi_{\mathbf{k}}^{l}=\xi_{\mathbf{k}}$. Note that the energies in the denominator do not carry spin indices. This is because the particles being paired have opposite spins. As a result, the magnetic contributions enter with opposite signs and cancel each other out.

For the spin-triplet case, the “bare” pair susceptibility becomes

$$\chi_{p,(t)}^{(0)}(i\nu_{n},\mathbf{q})=\sum_{\mathbf{k},a,b}\gamma_{\mathbf{k}}\gamma_{\mathbf{k}}\,n_{\mathbf{k}+\frac{\mathbf{q}}{2},\downarrow}^{a}\;n_{-\mathbf{k}+\frac{\mathbf{q}}{2},\downarrow}^{b}\;\frac{f\!\bigl(\xi_{\mathbf{k}+\frac{\mathbf{q}}{2},\uparrow}^{a}\bigr)+f\!\bigl(\xi_{-\mathbf{k}+\frac{\mathbf{q}}{2},\uparrow}^{b}\bigr)-1}{i\nu_{n}-\xi_{\mathbf{k}+\frac{\mathbf{q}}{2},\uparrow}^{a}-\xi_{-\mathbf{k}+\frac{\mathbf{q}}{2},\uparrow}^{b}},$$ (11)

where now the form factor is $p$-wave.

The “bare” charge susceptibility is obtained from the density–density correlation function:

$$\chi_{c}^{(0)}(i\nu_{n},\mathbf{q}) =\int_{0}^{\beta}d\tau\,e^{i\nu_{n}\tau}\big\langle T\,\rho_{c}(\mathbf{q},\tau)\,\rho_{c}(-\mathbf{q},0)\big\rangle_{0}$$

$$=-\frac{1}{\beta}\sum_{i\omega_{m}}\sum_{\mathbf{k},\sigma}\gamma_{\mathbf{k}}\gamma_{\mathbf{k}}\,G_{\mathbf{k}-\frac{\mathbf{q}}{2},\,\sigma}(i\omega_{m})\;G_{\mathbf{k}+\frac{\mathbf{q}}{2},\,\sigma}(i\omega_{m}+i\nu_{n}),$$ (12)

which leads to the following expression for the charge susceptibility:

$$\chi_{c}^{(0)}(i\nu_{n},\mathbf{q}) =-\sum_{\mathbf{k},\sigma}\sum_{a,b}\gamma_{\mathbf{k}}\gamma_{\mathbf{k}}\,n_{\mathbf{k}-\frac{\mathbf{q}}{2},\bar{\sigma}}^{a}\,n_{\mathbf{k}+\frac{\mathbf{q}}{2},\bar{\sigma}}^{b}\frac{f\bigl(\xi_{\mathbf{k}-\frac{\mathbf{q}}{2},\sigma}^{a}\bigr)-f\bigl(\xi_{\mathbf{k}+\frac{\mathbf{q}}{2},\sigma}^{b}\bigr)}{i\nu_{n}+\xi_{\mathbf{k}-\frac{\mathbf{q}}{2}}^{a}-\xi_{\mathbf{k}+\frac{\mathbf{q}}{2}}^{b}}.$$ (13)

Note that once again the energies in the denominator do not include spin indices. In this case, although the Green’s functions are spin dependent, the energy dispersions appear with opposite signs for the same spin projection, and the spin-dependent contributions cancel out.

It is important to note that we have computed only the “bare” susceptibility in each case. As shown in Ref. [11], the full susceptibility $\chi_{p}(i\nu_{n},\mathbf{q})$ is related to the bare one via a Dyson-like equation:

$$\chi_{p}(i\nu_{n},\mathbf{q})=\frac{\chi_{p}^{(0)}(i\nu_{n},\mathbf{q})}{1-V_{p}\,\chi_{p}^{(0)}(i\nu_{n},\mathbf{q})}.$$ (14)

A similar result also holds for the CDW susceptibilities, as shown in our derivation in Appendix A. Therefore:

$$\chi_{c}(i\nu_{n},\mathbf{q})=\frac{\chi_{c}^{(0)}(i\nu_{n},\mathbf{q})}{1-V_{c}\,\chi_{c}^{(0)}(i\nu_{n},\mathbf{q})}.$$ (15)

An immediate consequence of these relations is that the critical temperature associated with the instability of each phase is determined by the following condition:

$$\chi_{i}^{(0)}(i\nu_{n}=0,\mathbf{q})\Big|_{T=T_{c}^{(i)}}=\frac{1}{V_{i}},$$ (16)

where the index $i$ refers to the type of the order-parameter susceptibility [e.g., pair (denoted by “$p$”) and CDW (denoted by “$c$”)] and the corresponding critical temperature $T_{c}^{(i)}$. To keep the parameter space of the present model manageable, we will also assume from now on that $V_{p}=V_{c}=V$, for simplicity.

Finally, another important quantity to consider is the ferromagnetic susceptibility. As shown in Ref. [24], in the calculation of both the uniform charge and spin susceptibilities (i.e., for $\mathbf{q}=0$), the linear response theory applied to the HK model breaks down. However, this problem can be circumvented by employing the following thermodynamic expression for the uniform spin susceptibility, which is given by:

$$\chi_{s}=\lim_{B\rightarrow 0}\frac{\partial M}{\partial B}=\lim_{B\rightarrow 0}\sum_{\mathbf{k}}\left(\frac{\partial n_{\mathbf{k}\uparrow}}{\partial B}-\frac{\partial n_{\mathbf{k}\downarrow}}{\partial B}\right)\mu_{B},$$ (17)

where $M$ is the magnetization, $\mu_{B}$ is the Bohr magneton, and $n_{\mathbf{k}\sigma}$ is defined in Eq. (6).

Figure 1: Phase diagram of dimensionless interaction $u=U/W$ versus doping parameter $x$ of the competing ordering tendencies that emerge in the 2D HK model with $V=0.65$, with no magnetic field ($B=0$) and no strain ($\delta=0$). The label ‘CDW + PDW’ means that both CDW and PDW at $\mathbf{q}=(\pi,\pi)$ with $s$-wave symmetry are found within this region. The beige color denotes a region where a CDW at $\mathbf{q}=(\pi,\pi)$ with $s$-wave symmetry (with a subdominant PDW) appears in the model. The red part denotes $s$-wave singlet SC. At intermediate-to-strong coupling, CDW orders allow the emergence of PDW as a secondary order with the same symmetry and the same wavevector. The dashed lines are only a guide to the eye.

IV Results and Discussion

To begin with, we discuss how the competition among the different ordering tendencies in the 2D HK model takes place in the absence of strain ($\delta=0$) and external magnetic field ($B=0$). As will become clear soon, from the analysis of the order-parameter susceptibilities (see, e.g., Refs. [25, 26, 27, 28] in the context of other 2D correlated models), we obtain strong evidence of multiple short-order tendencies of comparable strengths within the present model. For this reason, we will focus henceforth on the dominant orders only, for each regime of the dimensionless interaction $u$ and doping parameter $x$.

For infinitesimal $V\rightarrow 0$, we confirm that the only instability of NFL phase in the model away from half-filling is towards the formation of superconductivity, which agrees with previous results regarding this model that are available in the literature [11, 12, 16, 17]. However, for a finite $V$, we obtain a variety of additional short-range orders that emerge in the corresponding phase diagram, which are summarized in our Fig. 1. We observe that at weak coupling near half-filling and beyond, a bidirectional fluctuating $d$-wave CDW order appears near $\mathbf{q}\approx(\pi,\pi)$. At weak-coupling and near half-filling, the essential driver for this order is the presence of approximate nesting of the Fermi surface.

As the local interaction $u$ becomes intermediate-to-strong, an intertwined CDW with a PDW near $\mathbf{q}\approx(\pi,\pi)$ dominates and this behavior takes place from a small doping up to a critical doping parameter $x_{c}$, in agreement with our previous results in Ref. [17]. This is labeled as ‘CDW + PDW’ in Fig. 1 and both phases possess $s$-wave symmetry. This feature of the PDW being intertwined with the CDW region displayed by the 2D HK model is very interesting, since it indicates a close connection between these two phases [29, 30, 31, 32, 33, 34]. Therefore, the 2D HK model provides good evidence supporting this interdependent nature of CDW and PDW phases within the underdoped regime and intermediate $u$. A similar mechanism was proposed to describe some aspects of the physics of the pseudogap regime in the cuprates in the underdoped regime (see, e.g., Refs. [29, 30, 31, 32, 33, 34]).

Moreover, for an even higher interaction $u$, another charge ordered phase emerges in the model: a unidirectional (stripe-like) CDW near $\mathbf{q}=(\pi,0)$ with $d$-wave symmetry (with an accompanying subdominant PDW order with the same symmetry and associated with the same wavevector). For such strong couplings, the CDW order is not necessarily responsive to details of the Fermi surface structure, and depends to a great extent only on the strength of $V$. Lastly, far from half-filling, first a $p$-wave spin-triplet SC phase emerges and then, with further doping, an $s$-wave spin-singlet SC appears. The former SC phase is related to the presence of ferromagnetic fluctuations in the model, which mediate the pairing interaction (cf. the discussion about the ferromagnetic fluctuations for the 2D HK model in the case of $n=1$ orbital in Appendix B). As soon as these fluctuations become extremely weak as a result of strong overdoping, a conventional SC phase (i.e., $s$-wave) emerges at low temperatures.

Figure 2: Phase diagram of dimensionless interaction $u=U/W$ versus doping parameter $x$ of the competing orders in the 2D HK model for $B=0.05$, $V=0.4$ and $\delta=0$. The light blue regions correspond to $d$-wave CDW with $\mathbf{q}=(\pi,0)$, and the beige color denotes the $s$-wave CDW with $\mathbf{q}=(\pi,\pi)$. The dashed lines are only a guide to the eye.

IV.1 Magnetic field effect

The next step is to analyze the effect of a finite external magnetic field applied to the model ($B\neq 0$). Overall, we observe that both the spin-singlet SC and PDW phases are suppressed by an external magnetic field (see Fig. 2). For the spin-triplet SC and CDW cases, the situation becomes more interesting, since these phases strongly compete with each other, depending on the interaction $u$ and the doping parameter $x$. For very small interaction potentials ($V\to 0$), the spin-triplet SC phase becomes the dominant instability of the model throughout the phase diagram. This is related to the presence of ferromagnetic fluctuations in the model, as was already discussed in the previous section, which become even stronger as a result of the magnetic field. However, for slightly larger $V$, the CDW phase becomes increasingly relevant and starts to dominate the $p$-wave SC phase for some regimes, as shown in Fig. 2.

It is important to note that this result holds for small magnetic fields. As the strength of $B$ increases, a stronger interaction $V$ is required for the CDW phase to occupy a significant portion of the phase diagram, otherwise the spin-triplet SC phase prevails. Moreover, the $p$-wave SC phase is now clearly enhanced near the Mott-insulating regime, where this pairing state becomes the dominant order. The CDW fluctuations also exhibit a subtler response to the application of a finite $B$: although it is modestly enhanced at intermediate dopings, the main effect of the magnetic field is, in fact, to shift the regions where this phase is already favored due to an approximate nesting condition. The only exception to this rule is the $d$-wave unidirectional (i.e., stripe-like) CDW at the wave vector $\mathbf{q}=(\pi,0)$, which now appears in a region of the phase diagram of weak-to-intermediate $u$ and at moderate-to-large doping $x$.

Figure 3: Phase diagram of dimensionless interaction $u=U/W$ versus doping parameter $x$ of the 2D HK model for the strain anistropy given by $\delta_{c}=W/8$ and interaction $V=0.15$. Both PDW orders represented above are unidirectional and have ordering vector near $\mathbf{q}=(\pi,0)$. The yellow region denotes the $d$-wave CDW phase near $\textbf{q}=(\pi,\pi)$, as the beige color corresponds the $s$-wave CDW near $\mathbf{q}=(\pi,\pi)$. The dashed lines are only a guide to the eye.

IV.2 Strain effect

In the previous section, we have obtained that the SC correlations appear stronger whenever the CDW correlations are weaker and vice-versa (i.e., they both compete with each other). Since the application of strain $\delta$ is well known to reduce charge ordering tendencies, we now turn to a discussion of these effects in the present model. We show the corresponding results in Fig. 3. From this plot, we observe that the spin-triplet SC phase at intermediate-to-large doping remains to some extent robust throughout the phase diagram. Despite this, we note that the $d$-wave spin-singlet SC phase now appears in some regions, and, in addition, it becomes displaced toward the edges of the phase diagram as the strain anisotropy increases. For even larger $\delta$, the $d$-wave SC region continues this trend, and eventually, for higher values ($\delta>0.15$), it disappears completely.

The CDW phase follows the aforementioned expected behavior. Although some regions (where CDW dominates) remain qualitatively similar as before, it now becomes clearly more fragile, disappearing almost entirely for smaller $V$. Its persistence along some characteristic parts of the phase diagram arises because the CDW with $\mathbf{q}=(\pi,\pi)$ is intrinsically a bit more robust in these regions; thus, even when the application of strain tends to suppress it, some remnants of the phase are maintained. Consequently, if the interaction $V$ is fixed and $\delta$ is increased further, the CDW phase gradually vanishes altogether from the phase diagram.

It is also interesting to point out that a unidirectional PDW phase associated with the ordering vector $\mathbf{q}=(\pi,0)$ is strongly favored under strain application for large $u$ and within the underdoped regime. In fact, when the strain increases to a critical $\delta_{c}=W/8$, a particularly interesting phenomenon emerges: the PDW phase becomes exactly degenerate with the SC phase and spans almost the whole phase diagram. In Fig. 3, we denote this fact by the label ‘SC + PDW’.

Now, focusing specifically on the interconnection between the spin-singlet PDW phase and the CDW, the hierarchy between these two phases becomes notably different from the previously discussed cases. The application of strain now changes the balance between those intertwined phases, with PDW now becoming the primary order for stronger interactions with the subsequent appearance of a stripe-like CDW order as the secondary phase. Nevertheless, both phases remain associated with the same symmetry and possess the same modulation described by the wavevector $\mathbf{q}$. The exception to this rule is the regime of weak-to-intermediate interaction and intermediate doping, where the previous behavior remains true for the CDW with $\mathbf{q}=(\pi,\pi)$ wavevector, since such a bidirectional CDW is still stabilized in regions of the phase diagram with approximate nesting of the underlying Fermi surface of the 2D HK model.

From an experimental point of view, it is still very challenging to determine whether unidirectional PDW can be favored by strain effects in the context of the cuprate superconductors. The experimental measurements of a PDW state in these compounds have been limited so far to surface-sensitive scanning tunneling microscopy (STM) techniques (see, e.g., Ref. [35]). Nevertheless, we point out that nuclear magnetic resonance (NMR) experiments are a promising pathway in order to demonstrate the PDW in the bulk of the cuprate compounds under strain and these measurements are currently underway. Therefore, we believe that the results of the HK model could potentially provide a guide for future research into the PDW phase in these systems.

Figure 4: Phase diagram of dimensionless interaction $u=U/W$ versus doping parameter $x$ of the 2D HK model for $t^{\prime}=0.45t$ and $V=0.15$ (with $\delta=0$ and $B=0$). The grey region corresponds to the NFL metallic phase. The dashed lines are only a guide to the eye.

IV.3 Inclusion of next-nearest neighbor hopping $t^{\prime}$

Another interesting case that will be considered here refers to the inclusion of next-nearest-neighbor hopping $t^{\prime}$ in the band structure in the model (assuming, for simplicity, no uniaxial strain and no applied magnetic field). In this case, the energy dispersion of the 2D HK takes the form

$$\xi_{\mathbf{k}}=-2t(\cos k_{x}+\cos k_{y})+4t^{\prime}\cos k_{x}\cos k_{y}-\mu.$$ (18)

Our results obtained for this case are shown in Fig. 4. From this plot, we observe that increasing $t^{\prime}$ significantly weakens CDW correlations throughout the phase diagram. Indeed, as emphasized previously, since CDW ordering tendencies compete with pairing correlations with momentum $\mathbf{q}=0$ in the model, SC becomes stronger in this case. We note that for intermediate-to-strong interaction $u$, spin-triplet SC emerges upon doping the Mott insulator away from half-filling up to critical doping. This is related to the strong ferromagnetic fluctuations that exist in the model within this regime already mentioned before, which mediate the Cooper pairing in this case. This is expected to be an artifact of the 2D HK model, which can be corrected by generalizing the model and making it more realistic. Recently, several works in the literature have appeared [36, 37, 38, 39] that put forward a generalization of the 2D HK model via an introduction of a large number $n$ of orbitals, where it has been suggested that such an orbital HK (OHK) model, with some assumptions and in a specific limit, would approach the behavior, e.g., of the paradigmatic 2D Hubbard model defined with $N/n$ clusters ($N$ being the total number of sites) [40]. In agreement with the latter work, we show in Appendix B by reproducing successfully the expected behavior for such systems that, for strong interaction (i.e., $u>1$), ferromagnetic fluctuations turn out to be indeed suppressed in the model for $n\geq 4$. This further supports the hypothesis that, for a more realistic OHK model, spin-triplet SC order is not expected to emerge within this regime for a finite $U$. Lastly, we point out that a $d$-wave SC phase emerges in the 2D HK model for large $t^{\prime}$ with further doping or for small $u$.

V Summary and Outlook

In this work, we have examined the interplay between several competing ordering tendencies in the 2D HK model (that displays a Mott insulating phase at half-filling and an NFL metallic phase upon doping) under the influence of some external perturbations, including magnetic fields, uniaxial strain, and next-nearest-neighbor hopping. By means of an exact analytical calculation of the corresponding order-parameter susceptibilities, we have constructed detailed phase diagrams that capture the emergence and competition of SC, CDW, and PDW orders in the model. However, we point out that the 2D HK model has some limitations compared to more realistic strongly correlated models (such as, e.g., the 2D Hubbard model) in view of the presence of the local-in-momentum coupling $U$, whereas in fundamental models to describe, e.g., the cuprate compounds, a strong momentum dependence in the interaction is necessary.

A central result of our analysis is that such a relatively simple model can reproduce some key physical features, potentially associated with more complex strongly correlated systems. These include the competition between SC and CDW orders, the emergence of an intertwined connection between PDW and CDW orders at intermediate-to-strong coupling in the underdoped regime, and even a correspondence with the physics of the 2D Hubbard model by means of a generalization of the 2D HK model by including a large number of orbitals. This is both remarkable and significant: the exact solvability of the orbital HK model may provide an analytical scheme to investigate the physics of the 2D Hubbard model at intermediate-to-large couplings without the need for heavy numerical simulations [40].

Therefore, the HK model offers a uniquely transparent and tractable platform for exploring fluctuation-induced strongly correlated superconductivity, and also its interplay with both CDW and PDW orders. Future directions include investigating the correlated phases analyzed here by incorporating additional orbital degrees of freedom to the HK model and also extending it to include spin–orbit-coupling (see, e.g., Refs. [41, 42, 43, 44]) potentially bridging the gap between the HK model and real materials that exhibit strong correlations, topological properties, and unconventional pairing phases. Another interesting direction of research refers to the calculation of the transport properties associated with the NFL phase displayed by the 2D HK model and also its possible generalizations (see, e.g., Refs. [45, 46, 47] for an analysis of such properties in the context of another NFL phase).

Appendix A Derivation for the CDW susceptibilities at low temperatures

Let us now derive exactly the equation for the charge susceptibility $\chi_{c}(i\nu_{n})$ at $T\ll U$, $W$ in the presence of a CDW interaction (see Eq. (4)). Our goal is to evaluate the Matsubara charge susceptibility at finite coupling, i.e.

$$\chi_{c}=\int_{0}^{\beta}d\tau\,e^{i\nu_{n}\tau}\langle\hat{T}_{\tau}\rho_{c}(\tau,\mathbf{q})\rho_{c}(0,-\mathbf{q})\rangle_{V_{c}},$$ (19)

and relate it to the bare susceptibility $\chi_{c}^{(0)}(i\nu_{n})$, defined by the same expression evaluated at zero interaction (i.e., $V_{c}=0$).

We proceed by expanding the thermal expectation value in the interaction picture. The time-ordered correlator is given by

$$\langle\hat{T}_{\tau}\rho_{c}(\mathbf{q},\tau)\rho_{c}(-\mathbf{q},0)\rangle=\frac{\langle\hat{T}_{\tau}S(\beta)\rho_{c}(\mathbf{q},\tau)\rho_{c}(-\mathbf{q},0)\rangle_{0}}{\langle\hat{T}_{\tau}S(\beta)\rangle_{0}},$$ (20)

with the subscript $0$ denoting expectation values taken with respect to the non-interacting Hamiltonian and the $S$-matrix is defined as

$$S(\beta)=\hat{T}_{\tau}\exp\left(-\int_{0}^{\beta}d\tau^{\prime}\,H_{\text{CDW}}(\tau^{\prime})\right).$$ (21)

Expanding the numerator, we obtain the following

$$\langle\hat{T}_{\tau}S(\beta,0)\,\rho_{c}(\mathbf{q},\tau)\rho_{c}(-\mathbf{q},0)\rangle_{0}=\\ \sum_{m=0}^{\infty}\frac{(V_{c})^{m}}{m!}\left\langle\hat{T}_{\tau}\left(\int_{0}^{\beta}d\tau_{1}\,H_{\text{CDW}}(\tau_{1})\right)^{m}\rho_{c}(\mathbf{q},\tau)\rho_{c}(-\mathbf{q},0)\right\rangle_{0},$$

which can be usefully interpreted in terms of connected and disconnected diagrams, as typically obtained via Wick’s theorem [48, 49]. All disconnected-diagram-like terms are canceled by the denominator. To make this explicit, let us consider the first-order contribution, where the variable $\mathbf{q}$ is treated as an index to avoid cluttering the notation:

$$\left\langle\hat{T}_{\tau}\,\rho_{-q}(\tau_{1})\rho_{q}(\tau_{1})\rho_{q}(\tau)\rho_{-q}(0)\right\rangle=\quad\>\left\langle\hat{T}_{\tau}\,\rho_{-q}(\tau_{1})\rho_{q}(\tau_{1})\right\rangle\left\langle\hat{T}_{\tau}\,\rho_{q}(\tau)\rho_{-q}(0)\right\rangle$$

$$+\left\langle\hat{T}_{\tau}\,\rho_{-q}(\tau_{1})\rho_{q}(\tau)\right\rangle\left\langle\hat{T}_{\tau}\,\rho_{q}(\tau_{1})\rho_{-q}(0)\right\rangle$$

$$+\left\langle\hat{T}_{\tau}\,\rho_{-q}(\tau_{1})\rho_{-q}(0)\right\rangle\left\langle\hat{T}_{\tau}\,\rho_{q}(\tau_{1})\rho_{q}(\tau)\right\rangle.$$

Here, we note two important points. First, the term

$$\left\langle\hat{T}_{\tau}\,\rho_{-q}(\tau_{1})\rho_{-q}(0)\right\rangle\left\langle\hat{T}_{\tau}\,\rho_{q}(\tau_{1})\rho_{q}(\tau)\right\rangle$$

can be interpreted in terms of a disconnected diagram, and contributes only when $\mathbf{q}=0$. By contrast, if $\mathbf{q}\neq 0$ (as considered in the main text), this term vanishes. For the two remaining terms, the product

$$\left\langle\hat{T}_{\tau}\,\rho_{-q}(\tau_{1})\rho_{q}(\tau)\right\rangle\left\langle\hat{T}_{\tau}\,\rho_{q}(\tau_{1})\rho_{-q}(0)\right\rangle$$

can be interpreted in terms of a connected-diagram term and yields a contribution of the form $\chi(\tau_{1}-\tau)\chi(\tau_{1})$. The other term is again disconnected. A similar structure persists at higher orders. For instance, at the next order, six non-zero terms arise (four disconnected and two connected diagram terms), resulting in the series

$$\chi_{c}=\chi_{0}+V_{c}\chi_{0}^{2}+(V_{c})^{2}\chi_{0}^{3}+\cdots,$$ (22)

where we now change the notation slightly and set $\chi_{0}\equiv\chi^{(0)}_{c}$. Summing all terms, we obtain the following result for the geometric series:

$$\chi_{c}=\sum_{m=0}^{\infty}V_{c}^{m}\chi_{0}^{m+1}(i\nu_{n})=\frac{\chi_{0}(i\nu_{n})}{1-V_{c}\chi_{0}(i\nu_{n})}.$$ (23)

This is the Dyson-like equation for the charge susceptibility of the present model, which is discussed in the main text.

Figure 5: Ferromagnetic susceptibility for $x=0.1$ and $u=1.5$ for different numbers of orbitals ($n$) included in the 2D OHK model. It is clear that beyond four orbitals ($n\geq 4$), the ferromagnetic susceptibility of the model becomes suppressed at low temperatures. Here, we assume units such that $\mu_{B}=1$.

Appendix B Ferromagnetic susceptibility in the 2D orbital HK model

The orbital extension of the 2D HK model (henceforth, denoted by OHK model) is particularly interesting. As shown in some recent works [36, 37, 38, 39, 40], adding such degrees of freedom makes the 2D HK model more realistic since it can potentially capture the behavior displayed by the 2D Hubbard model defined with $N/n$ clusters ($N$ being the total number of sites). This offers a promising route for understanding some key features of the Hubbard model via an exactly solvable framework [40], while also brings the corresponding effective HK model potentially closer to the realistic behavior of some strongly correlated materials.

The form of the extended Hamiltonian, now including $n$ orbitals, is given by [36]

$$H_{\text{OHK}}=\sum_{\mathbf{k},\alpha,\beta,\sigma}t_{\alpha\beta}(\mathbf{k})\,c^{\dagger}_{\mathbf{k}\alpha\sigma}c_{\mathbf{k}\beta\sigma}-\mu\sum_{\mathbf{k},\alpha,\sigma}n_{\mathbf{k}\alpha\sigma}+\sum_{\mathbf{k},\alpha,\beta}U_{\alpha\beta}\,n_{\mathbf{k}\alpha\uparrow}\,n_{\mathbf{k}\beta\downarrow}.$$ (24)

As can be seen, the structure of this Hamiltonian is similar to the single-orbital case [Eq. (1)], but with the addition of the indices $\alpha$ and $\beta$, which label the different orbitals. The term $t_{\alpha\beta}(\mathbf{k})$ denotes the momentum-dependent hopping matrix between orbitals, while $U_{\alpha\beta}$ represents the interaction matrix.

Although this extension of the model may appear uneventful at first sight, it has a significant impact. A clear example is the Green’s function in the OHK model, which can no longer be written in the same form as in the single-orbital case, thereby requiring the use of the Lehmann representation. This is because the occupation number operator $n_{\mathbf{k}\alpha,\sigma}$ no longer commutes with the Hamiltonian, and the structure of the creation and annihilation operators becomes non-trivial. Consequently, these operators must be evaluated in the eigenbasis of the interacting system, which can be understood more clearly by examining the following structure of the Green’s function:

$$\mathcal{G}_{\alpha\beta}^{\text{OHK}}(\mathbf{k},\omega+i0^{+})=\frac{1}{Z}\sum_{m,n}\frac{e^{-\beta E_{m}}+e^{-\beta E_{n}}}{\omega+i0^{+}+E_{m}-E_{n}}\,\langle m|\hat{c}_{\mathbf{k}\alpha\sigma}|n\rangle\langle n|\hat{c}^{\dagger}_{\mathbf{k}\beta\sigma}|m\rangle,$$

where the matrix elements represent the creation and annihilation operators in the eigenbasis of the full Hamiltonian. Similarly, $E_{m}$ and $E_{n}$ are the corresponding many-body eigenenergies of the system. The partition function is given by ${Z}=\mathrm{Tr}\,e^{-\beta H_{\text{OHK}}}$, from which quantities such as the average occupation can be computed.

Observables such as the uniform spin susceptibility or pair susceptibility (both addressed in the present work) can be calculated using the same formalism, but now with a modified structure for the Green’s function and partition function that incorporates the orbital degrees of freedom. As mentioned in the main text, an important result that we want to emphasize here is that by increasing the number of orbitals in the 2D OHK model, this leads to a suppression of the ferromagnetic susceptibility near half-filling. This behavior is clearly observed in Fig. 5. Since the 2D Hubbard model is expected to be described by the OHK model in the large-$n$ limit, this result is consistent with the expected behavior of the Hubbard model in that regime.