Theory of d-wave high temperature superconductivity in the cuprates involving non-linear lattice modes

TL;DR

This paper proposes a d-wave superconductivity mechanism in cuprates involving nonlinear lattice modes, deriving a gap equation within an extended Hubbard model that explains key experimental features like the superconducting dome.

Contribution

It introduces a detailed derivation of the superconducting gap equation considering nonlinear lattice interactions and demonstrates how the extended Hubbard model captures experimental phenomena.

Findings

Explicit d-wave symmetry in electron-lattice interaction

Transition temperature with k-dependence in the Brillouin zone

Model reproduces the superconducting dome and phase diagram

Abstract

The transition mechanism in high temperature cuprate superconductors is an outstanding puzzle. A previous suggestion on the role of non-linear local lattice instability modes on the microscopic pairing mechanism in high temperature cuprate superconductors \cite{Lee:JSNM09} is re-examined to provide a viable mechanism for superconductivity in these cuprates via an unusual lattice vibration in which an electron is predominantly interacting with a nonlinear $Q_2$ mode of the oxygen clusters in the CuO$_2$ planes. It is shown that the interaction has explicit d-wave symmetry and leads to an indirect coupling of d-wave symmetry between electrons. As a follow-up of \cite{Lee:JSNM09}, in this paper, we report detailed derivation of the superconducting gap equation and numerical solutions for the transition temperature as inherently integrated into the so-called Extended Hubbard Model (EHM). A unique feature in the EHM is that the transition temperature has an inherent k-dependence. In addition, superconducting gap solutions are restrained to specific regions in the first Brillouin zone (1BZ). It is very feasible to expect that the EHM naturally inherits a huge parameter space in which experimentally measured results, such as the well-known superconducting dome and the phase diagram from electronic Raman scattering \cite{Sacuto:RPP13} can be accommodated. The EHM model hence offers a viable venue to search for or confirm any signature in k-point-sensitive experimental measurements.

Figures

Equations

$$H$$

$$H_{\mathrm{ep}}=-\sum _{i,\sigma }g{Q}_i{c}_{i,\sigma }^{†}{c}_{i,\sigma }.$$

$$H_{\mathrm{ep}}=-\frac{1}{\sqrt{N}}\sum _{\mathrm{k}}\sum _{i,\sigma }g\left(\mathrm{k}\right)c_{i,\sigma }^{†}{c}_{i,\sigma }\left(a_{\mathrm{k},\sigma }+a_{-\mathrm{k},-\sigma }^{†}\right)e^{i\mathrm{k}\cdot {\mathrm{R}}_i},$$

$$g\left(\mathrm{k}\right)=2g{\left(\frac{\hslash }{8m{\omega }_{\mathrm{k}}}\right)}^{\frac{1}{2}}{e}_{\alpha }\left(\mathrm{k}\right)\left[\cos \left(a{k}_x\right)-\cos \left(a{k}_y\right)\right].$$

$$g\left(\mathrm{k}\right)$$

$$g\left(0\right)$$

$$H_4=\frac{1}{4}A\sum _i{Q}_i^{†},$$

$$H^{\prime }$$

$$S$$

$$H^{\prime }$$

$$H_4^{\prime }\backsimeq -\sum _i{G}_i{c}_{i,\sigma }^{†}{c}_{i,\sigma },$$

$$G_i$$

$$G_{\mathrm{k}}$$

$$⟨Q_{\mathrm{k}}⟩=\frac{1}{\sqrt{N}}\sum _i⟨Q_i^3⟩e^{i\mathrm{k}\cdot {\mathrm{R}}_i}.$$

$$H$$

$$H_t$$

$$C$$

$${\lambda }_{\mathrm{k}}=\sqrt{E_{\mathrm{k}}^2+B_{\mathrm{k}}^2},$$

$$E_{\mathrm{k}}={ϵ}_{\mathrm{k}}-\mu -A_{\mathrm{k}}.$$

$${ϵ}_{\mathrm{k}}=-2t\left[\cos \left(a{k}_x\right)+\cos \left(a{k}_y\right)\right],$$

$$F\le F_v$$

$$Z_t$$

$$F_t$$

$${⟨O⟩}_t=\frac{\mathrm{Tr}\left(O{e}^{-\beta H_t}\right)}{\mathrm{Tr}\left(e^{-\beta H_t}\right)}.$$

$$F_v$$

$$\frac{\partial F_v}{\partial B_{\mathrm{k}}}=0,$$

$$B_{\mathrm{k}}=\frac{1}{4N}\sum _{\mathrm{q}}\left(V_{\mathrm{k}-\mathrm{q}}-2U\right)B_{\mathrm{q}}\frac{\tanh \left(\frac{\beta {\lambda }_{\mathrm{q}}}{2}\right)}{{\lambda }_{\mathrm{q}}}.$$

$${\mathrm{\Delta }}_{\mathrm{k}}=\sum _{\mathrm{q}}\left(V_{\mathrm{k}-\mathrm{q}}-2U\right){\mathrm{\Delta }}_{\mathrm{q}}{F}_{\mathrm{q}}\left(\beta ,\mu \right),$$

$$F_{\mathrm{q}}\left(\beta ,\mu \right)=\frac{\tanh \left(\frac{\beta {\lambda }_{\mathrm{q}}\left(\mu \right)}{2}\right)}{4N{\lambda }_{\mathrm{q}}},$$

$${\lambda }_{\mathrm{q}}\left(\mu \right)=\sqrt{E_{\mathrm{q}}^2\left(\mu \right)+{\mathrm{\Delta }}_{\mathrm{q}}^2},$$

$$E_{\mathrm{q}}\left(\mu \right)={ϵ}_{\mathrm{q}}-\mu .$$

$$V_{\mathrm{k}}=\frac{{g\left(\mathrm{k}\right)}^2}{\hslash {\omega }_{\mathrm{k}}}.$$

Full text

1 Introduction

Since the discovery of high temperature superconductivity in the cuprates in 1986 Bednorz:ZP86 , there is still no consensus on the microscopic mechanism despite the numerous theories and ideas put forward Norman:RPP03 , Lee:RPP08 (and references therein). The purpose of this paper is to report numerical computations of a theory of $d$-wave superconductivity involving non-linear lattice modes. The non-linear modes are 2D vibrations of the oxygens in the ${\mathrm{CuO}}_{2}$ plane known as $Q_{2}$ modes. In these modes, the vibrations of a pair of oxygens along the $x$-axis of a ${\mathrm{CuO}}_{2}$ cluster are anti-phase to the vibrations of a pair of oxygens along the $y$-axis. Detailed descriptions have been discussed in the earlier paper Lee:JSNM09 .

We have also shown Lee_JSNM11 that such planar non-linear vibrations can produce anomalous small magnetic moments Lee:JSNM14 . The broken rotational symmetry in each ${\mathrm{CuO}}_{2}$ unit cell due to the generation of small magnetic fields makes our model a candidate to explain intra-unit-cell electronic nematicity measured recently Lawler:Nature10 . In summary, in addition to an explanation of the pseudogap and Fermi Arcs Hufner_RPP08 , our electron lattice model is able to explain, compared to the RVB ANDERSON:Science87 and t-J model theories Ogata:RPP08 , quantum critical point fluctuation theories Varma:RPP16 , the following experimental observations:

The electron-lattice model involving a predominant 2-D lattice mode is consistent with the unusual or surprising observation of a ubiquitous single phonon energy in the photo-emission measurements of a large number of samples Lanzara:Nature01 . We note that the interpretation of that paper tends to favour electron-lattice interaction. 2. 2.

It has been shown Lee_JSNM11 that such planar non-linear vibrations can produce anomalous small magnetic moments Lee:JSNM14 . 3. 3.

The broken $p/2$ rotational symmetry in each ${\mathrm{CuO}}_{2}$ unit cell due to the generation of small magnetic fields makes our model a candidate to explain intra-unit-cell electronic nematicity in the pseudogap phase as has been measured recently Lawler:Nature10 . A model with fully symmetric electron-lattice interaction will not be able to have this feature.

Another puzzling feature of high-$T_{c}$ cuprates comes from electronic Raman scattering measurement as reported in Sacuto:RPP13 , which suggests that the density of Cooper pairs is strongly dependent on doping, and for low doping levels it is confined in k-space where they form “islands” around the nodes. We shall show that the model Lee:JSNM09 by construct provides a natural setting in which k-space electronic distribution of the Cooper pairs, such as that observed in Sacuto:RPP13 , can be understood.

The main purpose of this paper is to report on the superconducting behavior as inherently embedded in the effective model Hamiltonian of Lee:JSNM09 . We focus attention on this effective Hamiltonian and show that it leads to a $d$-wave gap equation. Numerical solutions of the $d$-wave gap equation have been studied and we report some of the results here.

2 Model Hamiltonian

The basis of the model has been discussed in detail in an earlier paper Lee:JSNM09 . We shall refer the model as the “Extended Hubbard Model” (EHM). We suggested an operating mechanism which we termed as “hole induced suppression of the Jahn-Teller effect” for a CuO cluster centered around a Cu ion. This means that when the cuprate is doped to produce holes located in the planar Cu-O bond, the presence of the hole breaks the symmetry of the oxygen coordination which leads to the well-known Jahn-Teller (JT) instability for Cu2+. This causes a hole induced suppression of the JT distortion. When an electron moves to the hole, the conditions for the JT instability are restored and the oxygen clusters will distort to a more stable configuration with lower potential energy. The relaxation energy is quite large and we suggested that this large energy can excite a non-linear mode involving the vibrations of the JT active oxygen modes. It should be noted that our suggestion was based on a qualitative argument involving broken symmetry. It has been shown by other workers studying the problem computationally that hole induced suppression of JT effect does indeed occur. They termed this suppression as “anti-Jahn-Teller effect” Kamimura_JSNM12 .

The total model Hamiltonian of the EHM is

[TABLE]

$n_{i,\sigma}=c^{\dagger}_{i,\sigma}c_{i,\sigma}$ is the electron number operator. The total Hamiltonian is made up of a Hubbard Hamiltonian Lee:RPP08 with an anharmonic lattice background $H_{4}$ and the electron-lattice interaction $H_{\mathrm{ep}}$. The last three terms in Eq. (1) represent the standard Hubbard Hamiltonian. The electron-lattice interaction term is written as

[TABLE]

The coupling parameter is $g$. It is convenient to express $Q_{i}$ as defined in Eq. (2) in terms of lattice phonon creation (destruction) operators $a^{\dagger}_{\boldsymbol{\mathrm{k}}}\ (a_{\boldsymbol{\mathrm{k}}})$,

[TABLE]

where

[TABLE]

${\omega}_{\boldsymbol{\mathrm{k}}}$ is the phonon frequency, $e_{\alpha}\left(\boldsymbol{\mathrm{k}}\right)$ is the phonon polarization vector along the $\alpha$-axis. ${\bf{\mathrm{R}}}_{i}$ is the vector position at site $i$. $g(\boldsymbol{\mathrm{k}})$ has the classic $d$-wave symmetry property,

[TABLE]

$t_{ij}$ is the hopping or tunneling parameter for an electron from site $i$ to site $j$, $U$ is the Hubbard on-site repulsion energy and $\mu$ is the chemical potential. The harmonic lattice background is given by $H_{\mathrm{harmonic}}$ and $H_{4}$ is the fourth order anharmonic term. We choose for simplicity

[TABLE]

where $A$ is the fourth order anharmonicity parameter. Detailed arguments for an effective Hamiltonian have been discussed in our earlier paper Lee:JSNM09 . The following effective extended Hubbard Hamiltonian $H^{\prime}$ is obtained from Eq. (1) via canonical transformation

[TABLE]

[TABLE]

Note that we have dropped the background harmonic term $H_{\mathrm{harmonic}}$ from $H^{\prime}$ as it is not directly involved in the dynamics of superconductivity in the EHM. We treat the strong electron correlation problem as in the Hubbard Model. In the Hubbard model, all repulsive interactions between electrons are projected out except the on-site repulsion for double occupation. The third term represents effective electron-electron interaction. We restrict the interaction $V_{ij}$ to nearest neighbours for simplicity. Explicit expression of $V_{ij}$ is referred to Eq. 65 in Section A.1 of the Appendix. The fourth term,

[TABLE]

represents phenomenologically the energy of polaron formation caused by the interaction of an electron with the non-linear mode. It is assumed that $G_{i}$ in $H^{\prime}_{4}$and its Fourier transform take the following form

[TABLE]

Here, $\left\langle Q_{\boldsymbol{\mathrm{k}}}\right\rangle\mathrm{\ }$is a thermal average and defined as

[TABLE]

We treat this thermal average as a parameter at this stage of approximation. We have approximated the contributions from anharmonicity by retaining only terms linear in $c^{\dagger}c$. Higher order terms in $c^{\dagger}c$ will be neglected in our first attempt to evaluate the anharmonic effects. Hence it must be noted that the parameter for effective electron-electron interaction includes both harmonic and anharmonic interactions. In the fifth term in Eq. (6), the chemical potential $\mu$ is determined by the doping level.

If we examine the terms in the canonically transformed Hamiltonian, Eq. (6), and compare them with the original terms added to a standard Hubbard Hamiltonian, Eq. (1), we can justify writing an extended Hubbard Hamiltonian to investigate the thermodynamic properties of the system. We note that there has been a considerable amount of effort Micnas:RMP90 in the use of extended Hubbard Hamiltonians to study superconductivity in narrow band systems with local attractive interactions. What is new in our study is the role of non-linear lattice modes in the polaron formation. We have shown that this enables the model to produce anomalous magnetic moments and show nematicity in the electron properties in the pseudogap phase Lee_JSNM11 .

$G_{\boldsymbol{\mathrm{k}}}$ in Eq. (2) possesses $d$-symmetry properties. It could be the JT stabilization or relaxation energy as described above which is quite large. Physically, we can expect the JT polaron to be trapped if thermal excitation $(\sim kT)$ is not sufficient to excite the non-linear mode to free the electron. The electron can move to another site by tunneling. Alternatively, the JT relaxation energy released by a nearby site could possibly excite a non-linear vibration to release the electron. This suggested physical scenario of charge carriers interacting with non-linear modes to form a polaron is different from the well-developed theory of polarons and bipolarons in a harmonic lattice Mott:WS96 .

The EHM by construct accommodates a mechanism (the JT $Q_{2}$ mode) that explains pseudogap phenomena in a natural setting without invoking exotic physics. Our first paper on EHM Lee:JSNM09 places emphasis on elucidating the pseudogap in the light of the EHM. The model can also generate superconducting gaps, as in the case of generic Hubbard models and variants, due to the presence of on-site hoping and repulsive mechanism intrinsic to this class of models. In addition, as will be explained in the following sections, transition temperature solutions are integrally tied to the dopant concentration (which creates holes in the cuprates) in the EHM. The solutions for transition temperature and supercondicting gap embedded in the model have non-trivial and rich structures. In the following we shall explore the characteristics of high-$T_{c}$ superconductivity phenomena embedded in the EHM.

3 Derivation of superconducting gap equation

We note that the structure of EHM is quite similar (but not entirely so) to the work by Micnas:PRB88 . To begin with, Eq. (6) is cast into the following form,

[TABLE]

We have dropped the prime superscript on $H^{\prime}$. Details of casting Eq. (6) into Eq. (9) are referred to subsection A.2 in the Appendix. $\epsilon_{\mathrm{\bf k}}$ is the Fourier transform of $t_{ij}$. To facilitate the derivation of superconducting gap equation from the model Hamiltonian, we begin with a trial Hamiltonian $H_{t}$,

[TABLE]

where two variational parameters are introduced, $A_{\boldsymbol{\mathrm{k}}}$ for pseudogap and $B_{\boldsymbol{\mathrm{k}}}$ for singlet superconductivity. To derive superconducting gap equation, we first diagonalize the trial Hamiltonian, then carry out variational calculation on the trial free energy with respect to the superconductivity variational parameter $B_{\boldsymbol{\mathrm{k}}}$ by neglecting $A_{\boldsymbol{\mathrm{k}}}$ (which is not related to superconductivity). The resultant diagonalized trial Hamiltonian is written in the form

[TABLE]

${\alpha}_{\boldsymbol{\mathrm{k}},\sigma},{\alpha}^{\dagger}_{\boldsymbol{\mathrm{k}},\sigma}$ are Bogolyubov-Valatin transformation operators of $c_{\boldsymbol{\mathrm{k}},\sigma},c^{\dagger}_{\boldsymbol{\mathrm{k}},\sigma}$, as defined in Eq. (A.3) in the Appendix. The eigenvalue ${\lambda}_{\boldsymbol{\mathrm{k}}}$ of $H_{t}$ is given by

[TABLE]

[TABLE]

Diagonalization of the trial Hamiltonian into the form Eq. (3), and derivation of the corresponding eigenvalue [Eq. (12)], are detailed in subsection A.3 in the Appendix. ${\epsilon}_{\boldsymbol{\mathrm{k}}}$ in Eq. (13), which originates from the kinetic energy (hopping term) $t_{ij}$ in the Hubbard Hamiltonian, can be approximated by tight-binding band energy,

[TABLE]

where $t$ is the nearest neighbor hopping parameter, while $a$ is the lattice parameter of the CuO lattice. This type of notation is used for Hubbard model, t-J model, etc., and reflects the physics of electron hopping from site to site.

To obtain the free energy $F$, we apply variational principle and write

[TABLE]

where $N_{e}=nN$, $n$ charge carrier density, and $N$ is the total number of lattice sites. The trial free energy $F_{t}$ is defined in terms of the partition function $Z_{t}$ for the trial Hamiltonian $H_{t}$,

[TABLE]

where $\beta$ is the inverse temperature $\beta=\frac{1}{k_{B}T}$, and $k_{B}$ Boltzmann constant. ${\left\langle H\right\rangle}_{t},{\left\langle H_{t}\right\rangle}_{t}$ which are defined as expectation values in the trial basis set, can be derived using the expectation value of an arbitrary operator in the same basis, as prescribed in statistical mechanics,

[TABLE]

The full expressions of $Z_{t},F_{t},\ {\left\langle H\right\rangle}_{t},{\left\langle H_{t}\right\rangle}_{t}$ and $F_{v}$ are derived in Appendix A.4. In particular, we have [from Eq. (A.4)]

[TABLE]

Minimizing $F_{v}$ with respect to $B_{\boldsymbol{\mathrm{k}}}$,

[TABLE]

results in the following equality for $B_{\boldsymbol{\mathrm{k}}}$,

[TABLE]

Details of the minimization procedure to arrive at Eq. (3) is reported in Appendix A.5. Note that in arriving at Eq. (3), the variable $A_{\boldsymbol{\mathrm{k}}}$ is ignored because we wish to concentrate only on the superconductivity sector of the current model. Replacing the symbol $B_{\boldsymbol{\mathrm{k}}}\to{\mathrm{\Delta}}_{\boldsymbol{\mathrm{k}}}$, where ${\mathrm{\Delta}}_{\boldsymbol{\mathrm{k}}}$ denotes the superconducting energy gap, we arrive at the gap equation for $d$-wave superconductivity,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The $V_{\boldsymbol{\mathrm{k}}}$ term in Eq. (20) is given as (see the derivation of Eq. (A.1) in Appendix A.1),

[TABLE]

It is the Fourier conjugate for $V_{ij}$ that appears in Eq. (6). In the EHM, $V_{\boldsymbol{\mathrm{k}}}$ carries a $d$-wave symmetry which is traced back to the $Q_{2}$-mode as embedded in $G_{\boldsymbol{\mathrm{k}}}$ term in Eq. (2). Explicitly, $V_{\boldsymbol{\mathrm{k}}-\boldsymbol{\mathrm{q}}}$, as appeared in Eq. (20), reads

[TABLE]

where we have defined

[TABLE]

In Eq. (25), we have assumed that $V$ is frequency-independent, i.e.,

[TABLE]

$V$, along with $U$, are treated as free parameters of the model.

The superconducting gap equation, Eq. (20), essentially agrees with Eq. (2.8) in Micnas:PRB88 and Eqs. (2), (3) in Micnas:JPCM02 , who employed Hartree-Fock factoring method to derive a mean field Hamiltonian. With nearest-neighbor interaction and no specific symmetry, Micnas:PRB88 used an ansatz for the gap function which has on-site extended $s$-wave and $d$-wave symmetries. The same results are arrived at in present work using variational method. However, in our case, we start with an interaction $V_{\boldsymbol{\mathrm{k}}}$ [Eq. (25)] which has $d$-wave symmetry.

The number of electrons (charge carriers) per ${\mathrm{Cu}}^{2+}$ lattice site, $n$, by definition, is given by

[TABLE]

Evaluation of the RHS of Eq. (28) leads to (for derivation, see Eq. (124) of subsection A.4 in the Appendix)

[TABLE]

Putting Eq. (28) and Eq. (29) together results in a constraint condition imposed by the charge carried concentration $n$,

[TABLE]

The concentration of charge carriers $n$ can be experimentally correlated to the hole doping concentration in the cuprate, $x$, so that

[TABLE]

Eq. (20), subjected to the constraint from Eq. (30), determines the superconductivity transition temperature in the present model. The (inverse) transition temperature, denoted ${\beta}_{c}(=\frac{1}{k_{B}T_{c}})$, at which superconductivity is switched on when temperature is reduced from above, can be predicted by the model Hamiltonian by solving Eq. (20).

The method of approximation adopted in present work to solve the superconducting gap equation, Eq. (20), is inspired by Micnas:PRB88 . To this end, we use an ansatz of the form

[TABLE]

where ${\mathrm{\Delta}}_{0},{\mathrm{\Delta}}_{\eta}\mathrm{\ }$are k-point independent. ${\mathrm{\Delta}}_{\boldsymbol{\mathrm{q}}}$ has a $d$-wave symmetry in the limit $U\to 0$. This could be shown as follows. Substitute Eq. (32) into Eq. (20), and employ the fact that in the limit $U\to 0,$ ${\mathrm{\Delta}}_{0}\to 0$, we obtain

[TABLE]

The ansatz Eq. (32) is then inserted into the RHS of Eq. (33) to compute a value of ${\mathrm{\Delta}}_{\boldsymbol{\mathrm{k}}}.$ The symmetry property of the computed ${\mathrm{\Delta}}_{\boldsymbol{\mathrm{k}}}$ can be revealed by referring to Fig. 1 where a pair of symmetric points, ${\boldsymbol{\mathrm{q}}}_{1}$ and ${\boldsymbol{\mathrm{q}}}_{2}$, symmetric about the diagonal line $k_{x}=k_{y}$ are considered. We consider a point $\boldsymbol{\mathrm{k}}$ which lies on the diagonal line such that $\boldsymbol{\mathrm{k}}=\left(k{\mathrm{sin}\frac{\pi}{4\ }\ },k{\mathrm{cos}\frac{\pi}{4\ }\ }\right)=\left(\frac{k}{\sqrt{}2\ },\frac{k}{\sqrt{}2\ }\right)$. Hence,

[TABLE]

Due to the reflection symmetry of ${\boldsymbol{\mathrm{q}}}_{1}$ and ${\boldsymbol{\mathrm{q}}}_{2}$ about the diagonal line, it can be deduced that

[TABLE]

Replacing ${\boldsymbol{\mathrm{q}}}_{i}$ by $\boldsymbol{\mathrm{k}}-{\boldsymbol{\mathrm{q}}}_{i}$,

[TABLE]

From Eqs. (21), (22), (33), (35) and  (36),

[TABLE]

The summation over $\boldsymbol{\mathrm{q}}$ in  (33), namely $\sum_{\boldsymbol{\mathrm{q}}}{{\eta}^{2}_{\boldsymbol{\mathrm{k}}-\boldsymbol{\mathrm{q}}}{\eta}_{\boldsymbol{\mathrm{q}}}F_{\boldsymbol{\mathrm{q}}}\left(\beta,\mu\right)}$, will vanish due to pair-wise cancellation in the form of

[TABLE]

due to Eqs. (35), (36), (37). The result is a gap ${\mathrm{\Delta}}_{\boldsymbol{\mathrm{k}}}$ that has $d$-wave symmetry in the limit $U\to 0$,

[TABLE]

Eq. (38) is true because in the microscopic model, $V_{\boldsymbol{\mathrm{k}}}$ is proportional to ${\eta}^{2}_{\boldsymbol{\mathrm{k}}}$. If $V_{\boldsymbol{\mathrm{k}}}$ were to depend linearly on ${\eta}_{\boldsymbol{\mathrm{k}}}$, Eq. (38) will not hold. In comparison, $V_{\boldsymbol{\mathrm{k}}}$ in the phenomenological model by Micnas:PRB88 is not known, and the above result does not follow directly, despite it is still possible to search for $s$-wave and $d$-wave gap solutions in such models. The present model agrees with that of Micnas:PRB88 for pure $d$-wave gap solutions, $U=0\Rightarrow{\mathrm{\Delta}}_{0}=0.$ As far as we are aware of, this is the first time an e-lattice model can show explicitly $d$-wave gap.

4 Solutions for transition temperature and superconducting gap

Superconductivity transition temperature, $T_{c}=\frac{1}{k_{B}{\beta}_{c}}$, can be deduced from the gap equation in the present model. The existence of the solution to ${\beta}_{c}$ in the limit $U\to 0$ (hence ${\mathrm{\Delta}}_{0}=0$) can be proven in a limiting case as discussed below. In the limit $U\to 0$, for $T$ close to $T_{c}$, ${\mathrm{\Delta}}_{\eta}\ll 1$, $\Rightarrow{\mathrm{\Delta}}_{\boldsymbol{\mathrm{k}}}\approx{\mathrm{\Delta}}_{\eta}{\eta}_{\mathrm{\bf k}}\ll 1.$ We expand $F_{\boldsymbol{\mathrm{q}}}$ in Eq. (33) in powers of ${\mathrm{\Delta}}_{\boldsymbol{\mathrm{q}}}$,

[TABLE]

From Eq. (33) and Eq. (4),

[TABLE]

Solving Eq. (4) for ${\mathrm{\Delta}}_{\eta}$, we get ${\mathrm{\Delta}}_{\ \eta}=0$, or

[TABLE]

As $T\to T_{c},{\mathrm{\Delta}}_{\eta}\to 0$, Eq. (41) reduces to

[TABLE]

Consider the special case of high transition temperature, i.e., ${\beta}_{c}=\frac{1}{k_{B}T_{c}}\ll 1$. This allows us to expand

[TABLE]

Eq. (42) becomes

[TABLE]

We must also expand in $\beta$ for the constraint Eq. (31) in conjunction with Eq. (44),

[TABLE]

Eq. (4) solves $\mu$ directly for a given value of $x$. Note from Eq. (44) that ${\beta}_{c}$ does not depend on $\mu$ as it gets dropped out during the process of approximation. Eq. (44) clearly shows that (inverse) transition temperature ${\beta}_{c}$ exists in high temperature limit. Note that the summation in the denominator of Eq. (44), $\sum_{\boldsymbol{\mathrm{q}}}{{\eta}^{2}_{\boldsymbol{\mathrm{k}}-\boldsymbol{\mathrm{q}}}{\eta}_{\boldsymbol{\mathrm{q}}}\leadstoext O\left(N\right)}$ so that the value of ${\beta}_{c}\sim O\left(1\right)$, which is an agreeable expectation based on order of approximation estimate. Note also that the k-dependence of the transition temperature in Eq. (44) in general does not get cancelled off. This is a distinctive feature of the EHM compared to other phenomenological models: that each k-point contributes differently to the transition temperature. $T_{c}$ measured in experiments, as interpreted in the EHM, is an effective value contributed by various k-point that are being probed.

In a more general scenario, switching on the $U$ term will break the pure $d$-wave symmetry of ${\mathrm{\Delta}}_{\boldsymbol{\mathrm{k}}}$, as is obvious from the ansatz for ${\mathrm{\Delta}}_{\boldsymbol{\mathrm{k}}}$. The solutions for ${\beta}_{c}$ in such a scenario are now discussed. Inclusion of the $U$ term would modify the simplified transition temperature Eq. (44), in addition to coupling the constraint from the dopant concentration $x$ to ${\beta}_{c}$. The gap equation Eq. (20) is highly non-linear and self-iterative. Full numerical solutions could only be obtained by deploying sophisticated computational approach (e.g., genetic algorithm) but we shall not pursue along this direction here. Instead, we shall only attempt to solve ${\beta}_{c}$ with some simplifying assumptions.

We begin with our assumed ansatz, Eq. (32). Slotting it into the LHS of the gap function Eq. (20), using Eq. (25), and distinguishing the gap equation into k-dependent and k-independent parts, the following two independent relations are obtained:

[TABLE]

Slotting ${\mathrm{\Delta}}_{\boldsymbol{\mathrm{q}}}={\mathrm{\Delta}}_{0}+{\mathrm{\Delta}}_{\eta}{\eta}_{\boldsymbol{\mathrm{q}}}$ back into the RHS of Eq. (4),

[TABLE]

[TABLE]

which are then combined into a matrix form,

[TABLE]

where

[TABLE]

To determine the critical temperature ${\beta}_{c}$, we require the determinant of the 2 by 2 matrix $M$ in Eq. (54) to vanish. This is a sufficient condition for the existence of non-trivial solutions to Eq. (4):

[TABLE]

To proceed with the abstraction of the transition temperature, the function $F_{\boldsymbol{\mathrm{q}}}\left(\beta\right)$ (as defined in Eq. (21)) near the transition temperature ${\beta}_{c}$ at a k-point q is assumed to take the approximated form

[TABLE]

where, as$\ \beta$ approaches ${\beta}_{c}$ from above ($T$ approaches $T_{c}$ from below), ${\mathrm{\Delta}}_{\boldsymbol{\mathrm{q}}}$ approximately vanishes, so that ${\lambda}_{\boldsymbol{\mathrm{q}}}\to E_{\boldsymbol{\mathrm{q}}}$. In our numerical scheme, ${\beta}_{c}$ at a given k-point k is obtained by seeking the root of Eq. (58),

[TABLE]

with the constraint from Eq. (30) considered concurrently. To this end, we make the replacement ${\lambda}_{\boldsymbol{\mathrm{q}}}\to E_{\boldsymbol{\mathrm{q}}}$ in the constraint Eq. (30), so that it now reads

[TABLE]

where

[TABLE]

The first Brillouin zone (1BZ) is discretized into $N$ equally spaced sites, ${\boldsymbol{\mathrm{k}}}_{j}\in\left\{{\boldsymbol{\mathrm{k}}}_{1},{\boldsymbol{\mathrm{k}}}_{2},\dots,{\boldsymbol{\mathrm{k}}}_{N}\right\}$. ${\beta}_{c}$, if existed, also depends on the level of dopant $x$, which is an independent variable of the model, and an experimentally controllable parameter. In our numerical scheme, a value of $N=N_{k}=3721$ was adopted. This corresponds to setting the k-points in the range $-\frac{\pi}{2a}\leq\ k_{x},k_{y}\leq\ \frac{\pi}{2a}$ with interval $\mathrm{\Delta}k_{x}=\mathrm{\Delta}k_{y}=\frac{\pi}{N_{1\mathrm{D}}a}$, $N_{1\mathrm{D}}=60$. A ${\beta}_{c}\ \mathrm{vs}.\ \ x$ curve is traced out for each ${\boldsymbol{\mathrm{k}}}_{j}$, where $x_{0}\leq x\leq x_{m}$, with $x_{0}=0.01,x_{m}=0.27$. Since the transition temperature is k-point specific and $x$-dependent, we denote ${\beta}_{c}={\beta}_{c}\left({\boldsymbol{\mathrm{k}}}_{j},x\right).$ The corresponding value for the superconducting gap ${\mathrm{\Delta}}_{{\boldsymbol{\mathrm{k}}}_{\boldsymbol{j}}}(x)$ at the k-point ${\boldsymbol{\mathrm{k}}}_{j}$ and dopant concentration $x$ exists only if the solution ${\beta}_{c}\left({\boldsymbol{\mathrm{k}}}_{j},x\right)$ exists. If ${\beta}_{c}\left({\boldsymbol{\mathrm{k}}}_{j},x\right)$ existed, the corresponding ${\mathrm{\Delta}}_{{\boldsymbol{\mathrm{k}}}_{j}}(x)$ is given by

[TABLE]

where

[TABLE]

Derivation of Eq. (62) and Eq. (63) is detailed in A.6 of the Appendix. ${\mathrm{\Delta}}_{0}$ is an arbitrary constant. To be quantitative in our numerical calculation, we assigned ${\mathrm{\Delta}}_{0}=t$, where $t$ is the hopping parameter in Eq. (14). Note that ${\mathrm{\Delta}}_{0},{\mathrm{\Delta}}_{\eta}$ are both k-independent per the ansatz Eq. (32).

A remark is in place for the numerical solution of ${\mathrm{\Delta}}_{{\boldsymbol{\mathrm{k}}}_{j}}(x)$, Eq. (62). The numerical values of $r_{\eta}\left({\boldsymbol{\mathrm{k}}}_{j},x\right)$ as determined via Eq. (63) will become infinite if $U$ is identically 0, despite the corresponding transition temperature (if existed) is finite. The apparent singularity in the gap value at $U=0$ is merely a numerical artifact due to the approximations made in current numerical scheme.

5 Numerical Results

The numerical computation for ${\beta}_{c}\left({\boldsymbol{\mathrm{k}}}_{j},x\right)$ and ${\mathrm{\Delta}}_{{\boldsymbol{\mathrm{k}}}_{\boldsymbol{j}}}\left(x\right)$ as predicted by the EHM is carried out in atomic units, where the following quantities assume numerical values as depicted:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

$a_{0}$ is the Bohr radius. Two parameters used for numerical calculation are the lattice constant $a$ for the CuO plane in the cuprates, and the tight-biding parameter $t$. In atomic unit, their values are respectively assigned as

[TABLE]

[TABLE]

Inverse temperature parameter $\beta$ is related to temperature $T$ via

[TABLE]

We shall express the value of $V$ and $U$ in unit of the hopping parameter $t$, $V=\ vt$, $U=\ ut$, where the reduced parameters are such that $u>0,v<0$. The numerical values of the gap function ${\mathrm{\Delta}}_{\boldsymbol{\mathrm{k}}}$ will also be expressed in unit of $t$.

Transition temperature ${\beta}_{c}\left({\boldsymbol{\mathrm{k}}}_{j},x\right)$ and superconducting gap ${\mathrm{\Delta}}_{{\boldsymbol{\mathrm{k}}}_{j}}\left(x\right)$, if existed, could be picked up by the numerical procedure (within the limits of the numerical resolution implemented) for any fixed values of $U,V$. However, only transition temperature of less than $T_{c\mathrm{last}}=500\ \mathrm{K}$ will be covered. Even if solutions at larger than 500 K exist, they will not be shown. Each k-point in the 1BZ can possibly generate a full $T_{c}$ vs. $x$ curve of their own, depending on the location of ${\boldsymbol{\mathrm{k}}}_{j}$, and the values of the free parameters, $U$ and $V$ assumed. These curves in general varies from k-point to k-point. In the special case of $U=0$, there shall be no solution for ${\beta}_{c}$ along the $k_{x}=\pm k_{y}$ line in the 1BZ due to the $d$-wave symmetry.

As an illustration, we show in Fig. 2 and Fig. 3 a few $T_{c}\ \mathrm{vs}.\ x$ curves at selected k-points for fixed parameters set $\left\{u=0.0,v=-20.0\right\}$ and $\left\{u=0.25,v=-1.0\right\}$ respectively. These figures provide a vivid illustration that transition temperature, as numerically solved from the superconducting gap equation, varies with dopant level $x$ that enters the solution via the constraint Eq. (30). Moreover, the shape of the $T_{c}\ \mathrm{vs}.\ x$ curve varies from k-point to k-point in the 1BZ. At a lower level, these $T_{c}\ \mathrm{vs}\ x$ solutions in turn are determined by two free parameters in the EHM, namely, $U$ and $V$.

For a fixed value of $x$ and parameter set {$u,v$}, some k-points in the 1BZ may develop ${\beta}_{c}$ while others don’t. The collection of all solutions for ${\beta}_{c}$ for all k-points (including those k-points where no solutions are developed) at a fixed value of $x$ and parameter set {$u,v$} can be displayed in the form of a contour plot. The series of figures illustrated in Fig. 4, Fig. 5 are the contour plots in the 1BZ for the parameters sets $\left\{u=0.0,v=\right.$ $\left.-20.0\right\}$ and $\left\{u=0.25,v=-1.0\right\}$ respectively. A collection of 12 plots, each corresponds to a specific value of $x$ ($0.0\leq x\leq 0.27$) are shown in each of Fig. 4 and Fig. 5. Each pixel in the contour plot is color coded to represent the value of $T_{c}({\boldsymbol{\mathrm{k}}}_{j},x)$ [which is equivalent to $\beta_{c}({\boldsymbol{\mathrm{k}}}_{j},x)$ via $T_{c}({\boldsymbol{\mathrm{k}}}_{j},x)=1/{k_{B}\beta_{c}({\boldsymbol{\mathrm{k}}}_{j},x)}$] at that k-point, ${\boldsymbol{\mathrm{k}}}_{{j}}$. The patterns of these plots, which are $x$- and $\left\{u,v\right\}$-dependent, show that solutions for ${\beta}_{c}$ only get developed in restricted regions in the 1BZ. In the illustrative example of the case {$u=0.0,v=-20.0$} in Fig. 4, regions which develop solution for ${\beta}_{c}$ concentrate in the approximately semi-circular rings, with small radii and variable widths, attached to the edges of $k_{x}=\pm\frac{2\pi}{a}$ in the 1BZ. Meanwhile, for $\left\{u=0.25,v=\right.$ $\left.-1.0\right\}$ (Fig. 5), the radii of these rings become larger in general, but their perimeters never cross the quarter segments confined by the diagonal and anti-diagonal lines, $k_{x}={\pm k}_{y}$. In addition, as $x$ tends towards the 0.27 limit, the semi-circles are deformed into triangular wedges with two of their common lengths pushed towards the boarders of the quarter segments (i.e., the $k_{x}={\pm k}_{y}$ lines). Each k-point (pixel) in the contour plot also has a gap value ${\mathrm{\Delta}}_{{\boldsymbol{\mathrm{k}}}_{\boldsymbol{j}}}\left({\beta}_{c},x\right)$ if a solution for ${\beta}_{c}({\boldsymbol{\mathrm{k}}}_{j},x)$ also exists there. Hence, the pattern of the contour plots for ${\mathrm{\Delta}}_{\boldsymbol{\mathrm{k}}}\left({\beta}_{c},x\right)$ are exactly similar as that for ${\beta}_{c}(\boldsymbol{\mathrm{k}},x)$. As such we have omitted the individual figures for the contour plots of ${\mathrm{\Delta}}_{\boldsymbol{\mathrm{k}}}\left({\beta}_{c},x\right)$. In general, the pattern of the solutions in the contour plot of the 1BZ not only varies with $x$, it is also sensitive to the choice of {$u,v$}. The contour plots selected in Fig. 4 and Fig. 5 illustrate the point that in the EHM, solutions to the superconducting gap function are restrained to specific regions in the 1BZ. This is a unique prediction of the EHM.

Another illustrative way to report the solutions obtained with our numerical scheme is to present them in the form of non color-coded point plot in $T_{c}-x$ space. This is illustrated for the case $\left\{u=0.0,v=-20.0\right\}$ in Fig. 6. It graphically presents the collection of all solutions ${\beta}_{c}$ found in the 1BZ at all values of $x$ scanned by our numerical scheme for a fixed parameter set {$u,v$}. The coordinates of each pixels, $\left\{x,T_{c}\right\}$, indicate the value of transition temperature ${\beta}_{c}$ obtained at the corresponding value $x$. Void pixels are those where no solution for ${\beta}_{c}$ (hence ${\mathrm{\Delta}}_{\boldsymbol{\mathrm{k}}}$) exists. Be noted that there is no values available for ${\mathrm{\Delta}}_{\boldsymbol{\mathrm{k}}}(\beta_{c},x)$ in the case with $u=0$ because it is numerically singular when $u=0$, as explained earlier. Color-coded point plot in $T_{c}-x$ space is used to illustrate the numerical results of the gap functions graphically. The collection of the gaps ${\mathrm{\Delta}}_{\boldsymbol{\mathrm{k}}}\left({\beta}_{c},x\right)\mathrm{/}t\mathrm{\ }$ contributed by all k-points in the 1BZ for the illustrative case {$u$=0.25, $v$=-1.0} is displayed in Fig. 7. Each pixel is color coded to reflect the value of ${\mathrm{\Delta}}_{\boldsymbol{\mathrm{k}}}\left({\beta}_{c},x\right)\mathrm{/}t$ at the coordinates $\left\{x,T_{c}\right\}$. Both color-coded and non color-coded point plots provide a sense of how solutions for ${\beta}_{c}$ is distributed in the $T_{c}-x$ space. In general, it is found that the solution patterns illustrated in Fig.(2) – Fig.(7) are numerically sensitive to the choice of the parameters, $u,v$. Rigorous and exhaustive exploration of the parameter space of the model and matching the solution patterns against existing experimental findings are not within the scope of this paper but reserved for future work.

6 Conclusion

Our theoretical and numerical works indicate that the EHM is phenomenologically rich. Superconducting gap solutions are restrained to specific regions in the 1BZ, and transitional temperature is $k$-dependent. It is very feasible to expect that the EHM naturally inherits a huge parameter space in which experimentally measured results, such as the well-known superconducting dome and the phase diagram from electronic Raman scattering Sacuto:RPP13 can be accommodated. The EHM model hence offers a viable venue to search for or confirm any signature in k-point-sensitive experimental measurements.

Acknowledgements

This work was begun at UKM and BSL thanks RAS for hospitality. BSL acknowledges financial support from Ministry of Science, Technology and Innovation and Malaysian Academy of Science under the Brain Gain Program Malaysia (MOSTI/BGM/B2L/

2/5063). YTL wishes to acknowledge the support of Universiti Sains Malaysia RU grant (grant number 1001/PFIZIK/811240). RAS thanks the Ministry of Higher Education of Malaysia for the support (grant number ERGS/1/2011/STG/UKM/01/25). All numerical works were exclusively conducted in USM by YTL. YTL gladfully acknowledges Dr. Chan Huah Yong from the School of Computer Science, USM, for providing us computing resources to carry out part of the calculations done in this paper.

Appendix A Appendix

A.1

Coefficients and their Fourier transforms of the model Hamiltonian

The total Hamiltonian in real space $H$ is as given in Eq. (1). The Hamiltonian $H$ is transformed into the canonical form $H^{\prime}$ via a standard displaced oscillator transformation as per Eq. (6). The primed terms $H^{\prime}_{1}$ and $H^{\prime}_{3}$ in Eq. (6) are:

[TABLE]

Referring to $H^{\prime}_{3}$ in (A.1), we identify

[TABLE]

The coefficients in Eq. (6), namely $t_{ij},V_{ij},c_{i,\sigma},c^{\dagger}_{i,\sigma},G_{i}$, are related to their Fourier counterparts via

[TABLE]

Comparing the expressions of $V_{ij}$ in Eq. (65) and Eq. (A.1),

[TABLE]

This explains how Eq. (24) is arrived at.

A.2

Derivation of Eq. (9) from Eq. (6)

Each term in Eq. (6) can be cast into Bloch representation with the aid of Fourier transformations from Eq. (A.1).

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Adding up Eq. (68) - Eq. (72), $H^{\prime}_{1}+H^{\prime}_{2}+H^{\prime}_{3}+H^{\prime}_{4}+H^{\prime}_{5}$ produces the Bloch representation model Hamiltonian, Eq. (9).

A.3

Diagonalization of trial Hamiltonian of Eq. (10) into the form Eq. (3)

The trial Hamiltonian of Eq. (10) is diagonalized by following a textbook approach, such as that discussed in pg. 164 in Kittel:QTOS87 . To this end, we also require the anti-commutative relations for the operators $c_{\boldsymbol{\mathrm{k}},\sigma},c^{\dagger}_{\boldsymbol{\mathrm{k}},\sigma}$. We write down the Heisenberg equation of motion for $c_{\boldsymbol{\mathrm{k}},\sigma}$ in terms of the trial Hamiltonian of Eq. (10),

[TABLE]

Expanding and simplifying the commutator in Eq. (73),

[TABLE]

We now look at the $I,II,III$ terms in Eq. (74) in turn. The term $I$ in Eq. (74) can be simplified to

[TABLE]

by using the following anti-commutative relations

[TABLE]

The term $II$ in Eq. (74) can be simplified to

[TABLE]

by using the following anti-commutative relation

[TABLE]

The term $III$ in Eq. (74) can be simplified to

[TABLE]

by using the following anti-commutative relation,

[TABLE]

Putting $I,II,III$ into Eq. (74),

[TABLE]

Similarly,

[TABLE]

In deriving Eq. (78) and Eq. (79) we have defined

[TABLE]

and assumed

[TABLE]

For the next step, we shall perform Bogolyubov-Valatin transformation (see, e.g., pg. 307 of Haken, Haken:QFTOS83 ) on $c^{\dagger}_{-\boldsymbol{\mathrm{k}}\boldsymbol{\mathrm{,-}}\sigma}$,$c_{\boldsymbol{\mathrm{-}}\boldsymbol{\mathrm{k}}\boldsymbol{\mathrm{,-}}\sigma},c_{\boldsymbol{\mathrm{k}},\sigma}{,c}^{\dagger}_{\boldsymbol{\mathrm{k}},\sigma}$ to obtain

[TABLE]

[TABLE]

[TABLE]

[TABLE]

with the condition $u^{2}_{\boldsymbol{\mathrm{k}}}+v^{2}_{\boldsymbol{\mathrm{k}}}=1$ and $u_{\boldsymbol{\mathrm{k}}},v_{\boldsymbol{\mathrm{k}}}$ are real numbers. In addition, we shall also assume

[TABLE]

Unless otherwise specify, as a short hand notation, we shall suppress the spin subscript in the$\ c_{\boldsymbol{\mathrm{k}}},c^{\boldsymbol{\dagger}}_{\boldsymbol{\mathrm{k}}}$ and ${\alpha}^{\dagger}_{\boldsymbol{\mathrm{k}}},{\alpha}_{\boldsymbol{\mathrm{k}}}$ operators with the understanding that they are implicitly $\sigma$-bearing, i.e.,

[TABLE]

[TABLE]

The operators in Eq. (84) fulfill the following anti-commutative relations:

[TABLE]

Inversing Eq. (85) we have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Insert Eq. (A.3) into the trial Hamiltonian $H_{t}$ of Eq. (10),

[TABLE]

If $H_{t}$ were diagonal in the $\left\{{{\alpha}^{\dagger}_{\boldsymbol{\mathrm{k}}},\alpha}_{\boldsymbol{\mathrm{k}}}\right\}$ representation, it can be expressed in the form

[TABLE]

where ${\lambda}_{\boldsymbol{\mathrm{k}}}$ the eigenvalue of $H_{t}$. The constant term, which does not contain any $\alpha$ operator, can be neglected when applying variation calculation on $H_{t}$. The equations of motion for ${\alpha}_{\boldsymbol{\mathrm{k}}}$ and ${,\ \alpha}^{\dagger}_{\boldsymbol{\mathrm{k}}}$ are respectively given by

[TABLE]

Putting $H_{t}$ in Eq. (87) together with ${\alpha}^{\dagger}_{\boldsymbol{\mathrm{k}}}$ in a commutator yields

[TABLE]

The terms $\left[{\alpha}^{\dagger}_{\boldsymbol{\mathrm{k}}},H_{t,III}\right]\ \mathrm{and\ }\left[{\alpha}^{\dagger}_{\boldsymbol{\mathrm{k}}},H_{t,IV}\right]$ that should otherwise appear in the RHS of Eq. (90) vanish due to the following anticommutative relations,

[TABLE]

Using the following results,

[TABLE]

the commutators in the RHS of Eq. (90) are reduced to

[TABLE]

and

[TABLE]

By Eq. (90), Eq. (92) and Eq. (93), the commutator $\left[{\alpha}^{\dagger}_{\boldsymbol{\mathrm{k}}},H_{t}\right]$ in Eq. (A.3) becomes

[TABLE]

Comparing Eq. (A.3) and Eq. (94), the eigenvalue ${\lambda}_{\boldsymbol{\mathrm{k}}}$ is identified, namely,

[TABLE]

We next look at $u^{2}_{\boldsymbol{\mathrm{k}}},v^{2}_{\boldsymbol{\mathrm{k}}}$ in Eq. (95). Upon diagonalization, the coefficients for the off-diagonal terms of $H_{t}$ in Eq. (87) (i.e. $H_{t,III}$, $H_{t,IV})$, which involve ${\alpha}^{\dagger}_{\boldsymbol{\mathrm{k}}}{\alpha}^{\dagger}_{\boldsymbol{\mathrm{-}}\boldsymbol{\mathrm{k}}}$,$\ {\alpha}_{-\boldsymbol{\mathrm{k}}}{\alpha}_{\boldsymbol{\mathrm{k}}}$, should vanish, i.e.,

[TABLE]

Since $u^{2}_{\boldsymbol{\mathrm{k}}}+v^{2}_{\boldsymbol{\mathrm{k}}}=1$, we parametrise

[TABLE]

It is also assumed that the parameter ${\theta}_{\boldsymbol{\mathrm{k}}}={\theta}_{\boldsymbol{\mathrm{-}}\boldsymbol{\mathrm{k}}}$ so that Eq. (97) is consistent with Eq. (83). The will cast Eq. (96) into

[TABLE]

Hence, from Eq. (96), Eq. (97) and Eq. (A.3), the eigenvalue ${\lambda}_{\boldsymbol{\mathrm{k}}}$ in Eq. (95) is now simplified to

[TABLE]

This completes the derivation of Eq. (12). We note that the result Eq. (99) does not have 4 branches as obtained in Micnas:PRB88 because no charge ordered states are considered in the present model. The term $H_{t,II}$ in Eq. (87) can be cast into the form

[TABLE]

where we have used the relation ${{\alpha}_{-\boldsymbol{\mathrm{k}}}\alpha}^{\dagger}_{-\boldsymbol{\mathrm{k}}}=1-{\alpha}^{\dagger}_{-\boldsymbol{\mathrm{k}}}{\alpha}_{-\boldsymbol{\mathrm{k}}}$ and replaced $-\boldsymbol{\mathrm{k}}\to\boldsymbol{\mathrm{k}}$. The constant term in Eq. (100) is evaluated by summing over the spin degree of freedom,

[TABLE]

In arriving at Eq. (101) the following identities, which can be derived from Eq. (97) and Eq. (A.3), have been used, namely,

[TABLE]

Be noted that the expression for the constant in Eq. (101), $C=\sum_{\boldsymbol{\mathrm{k}}}{E_{\boldsymbol{\mathrm{k}}}}-\sum_{\boldsymbol{\mathrm{k}}}{{\lambda}_{\boldsymbol{\mathrm{k}}}}$, is defined with summation over the variable $\boldsymbol{\mathrm{k}}$ only. Upon diagonalization, $H_{t,III}$, $H_{t,IV}$ in $H_{t}$ in Eq. (87) vanish. Combining $H_{t,II}$ [from Eq. (100)] with $H_{t,I}$ in Eq. (87) gives Eq. (3),

[TABLE]

This completes the diagonalization of the trial Hamiltonian*.*

A.4

Derivation of ${{F}}_{{t}},\ {\left\langle{H}\right\rangle}_{{t}},{\left\langle{{H}}_{{t}}\right\rangle}_{{t}}$ as appear in Eq. (15)

Derivation of ${\boldsymbol{F}}_{\boldsymbol{t}}$:

Partition function $Z_{t}$ for trial Hamiltonian $H_{t}$ is defined as

[TABLE]

where the constant $C$ is as defined in Eq. (101). The trial free energy is

[TABLE]

Derivation of ${\left\langle\boldsymbol{H}\right\rangle}_{\boldsymbol{t}}$:

The derivation ${\left\langle H\right\rangle}_{t}$ proceeds by beginning with the definition of the model Hamiltonian $H$, Eq. (9), which is separated into three parts,

[TABLE]

[TABLE]

In terms of the operators ${\alpha}_{\boldsymbol{\mathrm{k}},\sigma},\ {\alpha}^{\dagger}_{\boldsymbol{\mathrm{k}},\sigma}$, ${\left\langle h_{1}\right\rangle}_{t}$ reads

[TABLE]

The expectation values ${\left\langle{\alpha}^{\dagger}_{\boldsymbol{\mathrm{k}},\sigma}{\alpha}^{\dagger}_{\boldsymbol{\mathrm{-}}\boldsymbol{\mathrm{k}},-\sigma}\right\rangle}_{t},\ {\left\langle{\alpha}_{\boldsymbol{\mathrm{-}}\boldsymbol{\mathrm{k}},-\sigma}{\alpha}_{\boldsymbol{\mathrm{-}}\boldsymbol{\mathrm{k}},-\sigma}\right\rangle}_{t}$ in Eq. (107) vanish in representation in which ${\alpha}^{\dagger}_{\boldsymbol{\mathrm{k}},\sigma}{,\alpha}_{\boldsymbol{\mathrm{k}},\sigma}$ are diagonal. There are two expectation values left to be evaluated in Eq. (107), i.e.,

${\left\langle{\alpha}^{\dagger}_{\boldsymbol{\mathrm{k}},\sigma}{\alpha}_{\boldsymbol{\mathrm{k}},\sigma}\right\rangle}_{t}$ and ${\left\langle{\alpha}^{\dagger}_{\boldsymbol{\mathrm{-}}\boldsymbol{\mathrm{k}},-\sigma}{\alpha}_{\boldsymbol{\mathrm{-}}\boldsymbol{\mathrm{k}},-\sigma}\right\rangle}_{t}$. ${\left\langle{\alpha}_{\boldsymbol{\mathrm{k}},\sigma}{\alpha}^{\dagger}_{\boldsymbol{\mathrm{k}},\sigma}\right\rangle}_{t}$ is evaluated via

[TABLE]

In arriving at A 45 we have ignored the contribution from constant term $C$ in $H_{t}$ of A 41.

The subscrpts in ${\left\langle{\alpha}^{\dagger}_{\boldsymbol{\mathrm{-}}\boldsymbol{\mathrm{k}},-\sigma}{\alpha}_{\boldsymbol{\mathrm{-}}\boldsymbol{\mathrm{k}},-\sigma}\right\rangle}_{t}$ can be manipulated so that it can be expressed in terms of ${\left\langle{{\alpha}^{\dagger}_{\boldsymbol{\mathrm{k}},\sigma}\alpha}_{\boldsymbol{\mathrm{k}},\sigma}\right\rangle}_{t}$,

[TABLE]

where we have made use of ${\epsilon}_{-\boldsymbol{\mathrm{k}}}={\epsilon}_{\boldsymbol{\mathrm{k}}}$. Using Eq. (A.3), Eq. (108), Eq. (109), the square bracket term in Eq. (107) can be simplified to

[TABLE]

[TABLE]

To evaluate ${\left\langle h_{2}\right\rangle}_{t}$, we rearrange ${\left\langle c^{\dagger}_{\boldsymbol{\mathrm{k}}+\boldsymbol{\mathrm{q}},\sigma}c_{\boldsymbol{\mathrm{k}},\sigma}c^{\dagger}_{{\boldsymbol{\mathrm{k}}}^{\prime}-\boldsymbol{\mathrm{q}}{,\sigma}^{\prime}}c_{{\boldsymbol{\mathrm{k}}}^{\prime},{\sigma}^{\prime}}\right\rangle}_{t}$ by reverting the position of the $c^{\dagger}$operator twice and using anti-commutative relations for $c^{\dagger}_{\boldsymbol{\mathrm{k}}},c_{\boldsymbol{\mathrm{k}}}$ to arrive at

[TABLE]

In view of Eq. (112), $h_{2}$ reads (after a few manipulative steps in the subscript symbols)

[TABLE]

It can be shown that ${\left\langle c^{\dagger}_{\boldsymbol{\mathrm{q}},\sigma}{c^{\dagger}_{\mathrm{-}\boldsymbol{\mathrm{q}}\boldsymbol{\mathrm{,-}}\sigma}c_{\mathrm{-}\boldsymbol{\mathrm{k}}\mathrm{\ },-\sigma}c}_{\boldsymbol{\mathrm{k}},\sigma}\right\rangle}_{t}={\left\langle c^{\dagger}_{\boldsymbol{\mathrm{q}},\sigma}c^{\dagger}_{\mathrm{-}\boldsymbol{\mathrm{q}}\boldsymbol{\mathrm{,-}}\sigma}\right\rangle}_{t}{\left\langle{c_{\mathrm{-}\boldsymbol{\mathrm{k}}\mathrm{\ },-\sigma}c}_{\boldsymbol{\mathrm{k}},\sigma}\right\rangle}_{t}$ provided $\boldsymbol{\mathrm{q}}{\mathrm{\neq}}\boldsymbol{\mathrm{k}}$, or equivalently, $V_{0}$ is excluded in the summation in Eq. (113). ${\left\langle c^{\dagger}_{\boldsymbol{\mathrm{q}},\sigma}c^{\dagger}_{\mathrm{-}\boldsymbol{\mathrm{q}}\boldsymbol{\mathrm{,-}}\sigma}\right\rangle}_{t}$ is given by the following expression [after some algebra, and making use of Eq. (A.3), Eq. (108), Eq. (109), Eq. (114)],

[TABLE]

As for the derivation of Eq. (114), we can similarly obtain

[TABLE]

Using Eq. (114) and Eq. (115), the expectation value of ${\left\langle h_{2}\right\rangle}_{t}$ is now obtained,

[TABLE]

The $G_{\boldsymbol{\mathrm{q}}}$ term that enters ${\left\langle h_{3}\right\rangle}_{t}$ has $d$-wave symmetry and is physically significant pertaining to the pseudo-gap Lee:JSNM09 . However, as ${\left\langle h_{3}\right\rangle}_{t}$ does not contribute to the process of minimization of $F_{v}$ with respect to $B_{\boldsymbol{\mathrm{k}}}$ for deriving superconducting gap equation, we leave as it is without explicitly evaluating it. It is relabeled as ${\left\langle G\ \mathrm{term}\right\rangle}_{t}$,

[TABLE]

Putting the final expression of ${\left\langle h_{1}\right\rangle}_{t},{\left\langle h_{2}\right\rangle}_{t}$ and ${\left\langle h_{3}\right\rangle}_{t}$ together,

[TABLE]

Derivation of ${\left\langle{\boldsymbol{H}}_{\boldsymbol{t}}\right\rangle}_{\boldsymbol{t}}$:

Referring to Eq. (73),

[TABLE]

Using Eq. (A.3), Eq. (114), Eq. (115) we have

[TABLE]

This completes the derivation of ${\left\langle H_{t}\right\rangle}_{t}.$

Derivation of ${\boldsymbol{F}}_{\boldsymbol{v}}$:

Finally, putting everything together, i.e., Eq. (105) (for $F_{t}$), Eq. (118) (for ${\left\langle H\right\rangle}_{t}$) Eq. (120) (for ${\left\langle H_{t}\right\rangle}_{t}$), the variational free energy is written as

[TABLE]

Note that the last three terms in (A.4) are independent of $B_{\boldsymbol{\mathrm{k}}}$. Also, we have used the relation $N_{e}=Nn$, where $n$ denotes concentration of charge carries, defined as

[TABLE]

From Eq. (111), we can identify $n_{\boldsymbol{\mathrm{k}},\sigma}\equiv{\left\langle c^{\dagger}_{\boldsymbol{\mathrm{k}},\sigma}c_{\boldsymbol{\mathrm{k}},\sigma}\right\rangle}_{t}$ in terms of the variables $E_{\boldsymbol{\mathrm{k}}}$ and ${\lambda}_{\boldsymbol{\mathrm{k}}}$, namely,

[TABLE]

In view of Eq. (122) and Eq. (123), we have

[TABLE]

Note that in the RHS of Eq. (124), the spin degree of freedom has already been summed over.

A.5

Minimization of ${\boldsymbol{F}}_{\boldsymbol{v}}$ with respect to $B_{\boldsymbol{\mathrm{k}}}$ to derive the recurrent gap equation.

$F_{v}$ in Eq. (A.4) can be minimized term by term. Differentiating the first term in $F_{v}$ by making use of the relation

[TABLE]

we obtain

[TABLE]

where

[TABLE]

Differentiating the second term in $F_{v}$ results in

[TABLE]

Differentiating the third term in $F_{v}$ results in

[TABLE]

Putting (126), (127), (128), (129) together into (A.4),

[TABLE]

A.6

Derivation of and Eq. (62) and Eq. (63)

In linear algebra, a homogeneous system of Eq. (4) has infinitely many non-trivial solutions if the determinant vanishes, or equivalently, the matrix $M$ in Eq. (4) is a rank 1 matrix. The non-trivial solutions \left(\begin{array}[]{c}{\mathrm{\Delta}}_{0}\\ {\mathrm{\Delta}}_{\eta}\end{array}\right) has only one linearly independent arbitrary variable, which we will choose as ${\mathrm{\Delta}}_{0}$. ${\mathrm{\Delta}}_{\eta}$is parametrized in terms of ${\mathrm{\Delta}}_{0}$ via $r_{\eta}=\frac{{\mathrm{\Delta}}_{\eta}}{{\mathrm{\Delta}}_{0}}$. If the root for $D_{\boldsymbol{\mathrm{k}}}\left(\beta\right)=\left|M\right|=T_{1}T_{2\boldsymbol{\mathrm{k}}}+T_{3}T_{4\boldsymbol{\mathrm{k}}}\ =0$ [i.e., Eq. (58)] exists, the matrix M\is reduced into row-reduced echelon form (which is always possible because it is a rank 1 matrix with vanishing determinant),

[TABLE]

Thus,

[TABLE]

The numerical value of the gap ${\mathrm{\Delta}}_{\boldsymbol{\mathrm{k}}}$** **at k-point $\boldsymbol{\mathrm{k}}$ is given as per

[TABLE]

This completes the derivation of Eq. (62) and Eq. (63).