1. Introductions
The Hubbard model, though it is simple, is in the central position for understanding strongly correlated electron systems [1, 2]. The single-band Hubbard model [3] was considered as the starting point to explain the high-Tc superconductivity (SC) [4]. The strong Hubbard repulsion limit of the Hubbard model tends to the t–J model, which also was derived from a more realistic model for cuprates [5–7]. Numerous subsequent studies on these two models were done either analytically or numerically. Many numerical simulation results are very impressive but they are basically subject to the computational resources and so are far from conclusive ones. Useful analytical approaches include the Gutzwiller approximation [8], mean field (MF) theories [9–17], and the gauge theory [18–24] based on the slave particle formalism [25–31]. Analog to the slave particle models, a large class of models, e.g. spin-fermion models, were developed to study the strongly correlated systems such as cuprates based on the spin fluctuations [32–36].
Recently, the renormalized MF theory based on the Gutzwiller projection [14–17] has been generalized to that in the form of statistically-consistent Gutzwiller approximation [37], which was proved to be equivalent to the slave boson theory. The results in terms of the further subsequent generalization, i.e. a systematic diagrammatic expansion of the variational Gutzwiller-type wave function may be quantitatively compared with the experimental properties of cuprates [38].
We are not going to focus on the results obtained by these methods because they are too fruitful to be summarized. We will try to improve the slave boson method and fix some shortcomings of the theory. For example, the slave particle theory looks very powerful because it exactly maps a strongly correlated electron system to a weakly coupled slave particle one but things become difficult when dealing with local constraint conditions Ti = 0 (see equation (3 )), i.e. only one type of the single slave particle can occupy a lattice site i. The temporal component of the gauge field, λi(t) and the spatial components of the gauge field which are introduced to compensate for the gauge symmetry breaking by the MF approximation are not dynamic so the conventional perturbation theory is not applicable. In this paper, we try to solve these problems.
In terms of Dirac’s approach to solve the first-class constraint systems, a term −∑iλiTi with λi being the Lagrange multiplier is added to the Hamiltonian H. Since there are no temporal or spatial derivatives of λi in the Hamiltonian, [H, λi] = 0 and then λi will not evolve with time. In literature, λi was simply relaxed to a time-dependent field λi(t) and as the temporal component of the gauge field. This introduces the redundant degrees of freedom because λi should be kept static, i.e. an additional constraint ${\partial }_{t}{\lambda }_{i}(t)\equiv {\dot{\lambda }}_{i}(t)=0$ must be enforced. This point was missed before. Instead, in the MF approximation, a conventional approximation ${\lambda }_{i}(t)=\bar{\lambda }$, a constant with no spatial and temporal dependence, was taken. Although $\dot{\bar{\lambda }}$ is zero, obviously this brings many unphysical degrees of freedom so that the MF theory after this approximation is not reliable or controllable. Many further improvements are proposed to deal with this issue but they do not bring conclusive results [18–21]. Recent development in the statistical Gutzwiller approximation sheds light to systematically relieving the difficulties that original renormalization MF theory meets [37, 38]. In this paper, we make efforts to improve the slave particle theory by considering the additional constraint ${\dot{\lambda }}_{i}(t)=0$ instead of ${\lambda }_{i}=\bar{\lambda }$.
When an electron operator is decomposed into slave particles, a gauge symmetry is induced and λi(t) behaves as a gauge potential in the temporal direction. To remove the redundant gauge degrees of freedom, one has to introduce a GFC while keeping the physical observables is gauge invariant. Simply setting ${\lambda }_{i}(t)=\bar{\lambda }$ is a GFC but it is not a good GFC because it violates the constraint Ti = 0 and brings unphysical degrees of freedom. For a gauge theory with constraints, the GFC must be consistent with the constraints.
In this paper, we show ${\dot{\lambda }}_{i}(t)=0$ is a good GFC because it keeps Ti = 0 unchanging. Introducing this GFC may be thought of as an application of the general method of the gauge fixing in the gauge theory with Dirac’s first-class constraints [39, 40]. After gauge fixing, the gauge invariance of the system is equivalent to the Becchi-Rouet-Stora-Tyutin (BRST) global symmetry. For Dirac’s first-class constraints, the BRST symmetry is the criterion whether the GFC is correct because the BRST invariance requirement to the physical states is exactly equivalent to Dirac’s first-class constraints.
The constraint ${\dot{\lambda }}_{i}(t)=0$ endows λi(t) with the dynamics. The fluctuation from λi(t) through the interaction between λi(t) and the slave particles enables us to examine the stability of the variational wave function of a given phase. In this sense, the variational or MF states of the phases become controllable and we then solve the local constraint condition problem of the slave particle theory. For different purposes, we may deal with λi(t) in different ways. In the gauge theory of the slave boson, in order to study the normal state behavior of the systems, the spatial components of the gauge field are introduced [18–21]. The physical properties of the system were studied by integrating away the matter fields. The spatial components of the gauge field actually also play the role of the Lagrange multiplier of the counterflow constraint on the slave particles’ currents. The similar additional constraint to ${\dot{\lambda }}_{i}(t)=0$ is also needed but was not considered in the previous research. In this paper, we will not consider the spatial components of the gauge fields because there might be further complicated symmetry analysis more than the gauge symmetry in that case. We will leave it for a coming soon work.
In this paper, we focus on the pairing states of the spinons and study the BCS-type MF states. In this case, the spatial components of the gauge field are not necessary to be introduced. We examine the spinion pairing gap of the Hubbard and t–J models. For the Hubbard model, it was found that the MF state of the slave bosons at half-filling in a small U is an s-wave SC state [41, 42]. This is obviously wrong because the Hubbard model in the weak coupling limit is a conventional metal. We show that after considering the dynamics of λi(t) induced by the constraint ${\dot{\lambda }}_{i}(t)=0$, the SC state is not stable because integrating over λi(t) induces an unusual pairing instability of spinon’s Fermi surface and the spinon is gapped. This destroys the SC at half-filling.
Similarly, for the t–J model at the spinon pairing gap state or the SC state, integrating over λi(t) contributes an additional unusual term to the spinon pairing. This additional contribution does not destroy the MF SC gap but may substantially suppress it. Numerically, the gap will be smaller than at least a factor of one-fifth of the MF SC gap. Thus, one can expect the d-wave SC critical temperature Tc, whose MF value is of the order 1000 K, is substantially lowered and might be comparable with that of cuprates.
This paper is organized as follows: in section 2 , we explain why the additional constraint ${\dot{\lambda }}_{i}(t)=0$ is necessary for the slave particle theory. We relate this additional constraint to the BRST quantization of the gauge theory. In section 3 , we apply our theory to the Hubbard model at half-filling and show the SC state in the small U is unstable. In section 4 , we study the spinon gap suppression by λi(t)'s fluctuation in the t–J model and discuss the implication to Tc of curpates. In section 5 , we conclude this work and schematically look forward to the prospect of the applications of our theory to the strongly correlated systems.
2. Constraint to the Lagrange multiplier and BRST quantization
For a strongly correlated electron many-body system, a conventional perturbation theory based on the Fermi liquid theory does not work. In order to turn the strong interacting electron model to an equivalent weak coupling theory, a powerful method called the slave boson/fermion theory is applied [25–28]. For the electron operator ciσ at a lattice site i, the local quantum space is {∣0⟩, ∣↑⟩, ∣↓⟩, ∣↑ ↓⟩}. The completeness condition reads1 ) maps to a local constraint
$\begin{eqnarray}| 0\rangle \langle 0| +| \uparrow \rangle \langle \uparrow | +| \downarrow \rangle \langle \downarrow | +| \uparrow \downarrow \rangle \langle \downarrow \uparrow | =1.\end{eqnarray}$
The slave boson representation of the electron is mapping $| 0\rangle \to {h}^{\dagger },| \sigma \rangle \to {f}_{\sigma }^{\dagger }$ and ∣↑ ↓⟩ → d†. The operators h, d, fσ are called the holon, doublon, and spinon which destroy some vacuum ∣vac⟩. For the slave boson, h and d are the bosonic operators while fσ are the fermionic operators $\begin{eqnarray}\left[{d}_{i},{d}_{j}^{\dagger }\right]=[{h}_{i},{h}_{j}^{\dagger }]={\delta }_{{ij}},\{{f}_{i\sigma },{f}_{j\sigma ^{\prime} }^{\dagger }\}={\delta }_{{ij}}{\delta }_{\sigma ,\sigma ^{\prime} },\end{eqnarray}$
and so on. The completeness condition ( $\begin{eqnarray}{T}_{i}={h}_{i}^{\dagger }{h}_{i}+{f}_{i\uparrow }^{\dagger }{f}_{i\uparrow }+{f}_{i\downarrow }^{\dagger }{f}_{i\downarrow }+{d}_{i}^{\dagger }{d}_{i}-1=0,\end{eqnarray}$
i.e. a given lattice site can only be occupied by one given type of particles. The electron operator is decomposed into ${c}_{i\sigma }^{\dagger }={f}_{i\sigma }^{\dagger }{h}_{i}+\sigma {f}_{i,-\sigma }{d}_{i}^{\dagger }$ with σ = {↑, ↓} ≡ {+, −}. With the local constraint, the anti-commutation of the electron operators are equivalent to h, d s’ are bosonic while fσ s’ are fermionic. This equivalence also holds if fσ are bosonic and h, d are fermionic, which is called the slave fermion representation. In this work, we focus on the slave boson one although they are equivalent before MF approximations. Using the slave boson representation, a strongly correlated many-body electron Hamiltonian He can be mapped to a Hamiltonian H in the slave boson representation.2.1. Additional constraint
For a Hamilton system with constraints, we follow Dirac’s method to solve the constrained system and introduce a Lagrange multiplier λi. The Hamiltonian for the constrained problem is given by
$\begin{eqnarray}{H}_{\lambda }=H-\displaystyle \sum _{i}{\lambda }_{i}{T}_{i}.\end{eqnarray}$
In Schrödinger’s picture, H, λi and Ti are all time-independent. Notice that [Hλ, λi] = 0 and then λi does not evolve as time.Going to Heisenberg’s picture, all operators and fields Φi become time-dependent, ${{\rm{\Phi }}}_{i}(t)={{\rm{e}}}^{{\rm{i}}{H}_{\lambda }t}{{\rm{\Phi }}}_{i}{{\rm{e}}}^{-{\rm{i}}{H}_{\lambda }t}$ except λi since [Hλ, λi] = 0, where Φi stand for the ‘matter’ fields (holon, doublon, and spinon). On the other hand, because the constraint Ti = 0 has to be enforced in any space-time location, the Lagrange multiplier has to be relaxed to time-dependent. Thus, one has to add an additional constraint in order to be consistent with no time evolution of λi(t), namely
$\begin{eqnarray}{\dot{\lambda }}_{i}(t)=0.\end{eqnarray}$
By introducing a new Lagrange multiplier ${\pi }_{{\lambda }_{i}}(t)$ to force ${\dot{\lambda }}_{i}(t)=0$, the Lagrangian is then given by $\begin{eqnarray}\begin{array}{rcl}{L}_{\lambda } & = & \displaystyle \sum _{i}{\pi }_{{\lambda }_{i}}(t){\dot{\lambda }}_{i}(t)+\displaystyle \sum _{i\sigma }{f}_{i\sigma }^{\dagger }({\rm{i}}{\partial }_{t}+{\lambda }_{i}(t)){f}_{i\sigma }\\ & & +\displaystyle \sum _{i}({h}_{i}^{\dagger }({\rm{i}}{\partial }_{t}+{\lambda }_{i}(t)){h}_{i}+\displaystyle \sum _{i}{d}_{i}^{\dagger }({\rm{i}}{\partial }_{t}+{\lambda }_{i}(t)){d}_{i}\\ & & -\displaystyle \sum _{i}{\lambda }_{i}(t)-H.\end{array}\end{eqnarray}$
Another way to understand the relation between the Hamiltonian (4 ) and the Lagrangian (6 ) is as follows. According to (6 ), ${\pi }_{\lambda i}(t)=\tfrac{\delta L}{\delta {\dot{\lambda }}_{i}(t)}$, i.e. πλi(t) is the canonical conjugate field of λi(t). Therefore, according to the classical mechanics, the Lagrangian of the Hamiltonian (4 ) reads7 ) is exactly the same as (6 ).
$\begin{eqnarray}{L}_{\lambda }=\displaystyle \sum _{i}({\pi }_{\lambda i}{\dot{\lambda }}_{i}+{{\rm{\Pi }}}_{{{\rm{\Phi }}}_{i}}{\dot{{\rm{\Phi }}}}_{i})-{H}_{\lambda },\end{eqnarray}$
where ${{\rm{\Pi }}}_{{{\rm{\Phi }}}_{i}}$ are the canonical conjugate fields of Φi. The Lagrangian (2.2. Gauge symmetry
We now explain the reason to add the constraint ${\dot{\lambda }}_{i}(t)=0$ from the gauge symmetry point of view. In literature, instead of (6 ), the following Lagrangian is considered [9, 11]6 )10 ) where ξ is an arbitrary constant and integrating away πλi field, the path integral reads11 ) is not gauge invariant. In order to resolve this paradox, we recall the Faddeev–Popov quantization of the gauge theory. We insert 1 into the gauge invariant (9 ) to fix the redundant gauge degrees of freedom in terms of14 ) is gauge invariant. Comparing (14 ) and (11 ), they differ from a factor $\det \ ({\partial }_{t}^{2})$ after dropping N(ξ). At the present case, this determinant does not contain any fields and is a constant. This means that (11 ) is equivalent to (14 ). Therefore, up to a constant determinant, (11 ) is gauge invariant. However, for a non-Abelian gauge theory, the determinant in general is dependent on the gauge field and can not be dropped. This is why Faddeev–Popov ghost fields are introduced.
$\begin{eqnarray}\begin{array}{rcl}{L}_{\mathrm{GI}} & = & \displaystyle \sum _{i\sigma }{f}_{i\sigma }^{\dagger }({\rm{i}}{\partial }_{t}+{\lambda }_{i}(t)){f}_{i\sigma }+\displaystyle \sum _{i}({h}_{i}^{\dagger }({\rm{i}}{\partial }_{t}+{\lambda }_{i}(t)){h}_{i}\\ & & +\displaystyle \sum _{i}{d}_{i}^{\dagger }({\rm{i}}{\partial }_{t}+{\lambda }_{i}(t)){d}_{i}-\displaystyle \sum _{i}{\lambda }_{i}(t)-H.\end{array}\end{eqnarray}$
It was known that the electron operator ${c}_{i\sigma }^{\dagger }={f}_{i\sigma }^{\dagger }{h}_{i}+\sigma {f}_{i,-\sigma }{d}_{i}^{\dagger }$ is gauge invariant under $({h}_{i},{d}_{i},{f}_{i\sigma })\to {{\rm{e}}}^{-{\rm{i}}{\theta }_{i}}({h}_{i},{d}_{i},{f}_{i\sigma })$. LGI is invariant under this gauge transformation accompanied with ${\lambda }_{i}(t)\to {\lambda }_{i}(t)-{\dot{\theta }}_{i}$, i.e. λi(t) plays a role of a scalar gauge potential. There are redundant gauge degrees of freedom in the path integral $\begin{eqnarray}W^{\prime} =\int \displaystyle \prod _{i,t}{\rm{d}}{{\rm{\Phi }}}_{i}^{\dagger }(t){\rm{d}}{{\rm{\Phi }}}_{i}(t){\rm{d}}{\lambda }_{i}(t){{\rm{e}}}^{{\rm{i}}\int {\rm{d}}{{tL}}_{\mathrm{GI}}}.\end{eqnarray}$
One way to remove the redundant gauge degrees of freedom is replacing LGI by the Lagrangian ( $\begin{eqnarray}W=\int \displaystyle \prod _{i,t}{\rm{d}}{{\rm{\Phi }}}_{i}^{\dagger }(t){\rm{d}}{{\rm{\Phi }}}_{i}(t){\rm{d}}{\pi }_{{\lambda }_{i}}(t){\rm{d}}{\lambda }_{i}(t){{\rm{e}}}^{{\rm{i}}\int {\rm{d}}{{tL}}_{\lambda }}.\end{eqnarray}$
Making a transformation ${\dot{\lambda }}_{i}\to {\dot{\lambda }}_{i}+\xi {\pi }_{\lambda i}/2$ for equation ( $\begin{eqnarray}W\propto \int \displaystyle \prod _{i,t}{\rm{d}}{{\rm{\Phi }}}_{i}^{\dagger }(t){\rm{d}}{{\rm{\Phi }}}_{i}(t){\rm{d}}{\lambda }_{i}(t){{\rm{e}}}^{{\rm{i}}\int {\rm{d}}{{tL}}_{\mathrm{eff}}},\end{eqnarray}$
where $\begin{eqnarray}{L}_{\mathrm{eff}}={L}_{\mathrm{GI}}-\displaystyle \frac{1}{2\xi }\displaystyle \sum _{i}{\dot{\lambda }}_{i}^{2}(t).\end{eqnarray}$
This is a correct gauge fixing Lagrangian of the Abelian gauge theory but equation ( $\begin{eqnarray}1=\int \displaystyle \prod _{i,t}{\rm{d}}{\theta }_{i,t}\delta ({\dot{\lambda }}_{i}(t))\det \left(\displaystyle \frac{\delta {\dot{\lambda }}_{i}(t)}{\delta {\theta }_{j}(t^{\prime} )}\right),\end{eqnarray}$
and finally [43] $\begin{eqnarray}\begin{array}{rcl}1\cdot W^{\prime} & = & N(\xi )\int \displaystyle \prod _{i,t}{\rm{d}}{{\rm{\Phi }}}_{i}^{\dagger }(t){\rm{d}}{{\rm{\Phi }}}_{i}(t){\rm{d}}{\lambda }_{i}(t)\det \ ({\partial }_{t}^{2})\\ & & \times \exp \{{\rm{i}}\int {\rm{d}}{{tL}}_{\mathrm{eff}}\},\end{array}\end{eqnarray}$
where N(ξ) is an unimportant infinity constant. The path integral (2.3. BRST quantization
Historically, there is a standard approach to deal with the relativistic gauge theory with Dirac’s first-class constraints Tα(r, t) = 0 with $[{T}^{\alpha },{T}^{\beta }]={f}_{\gamma }^{\alpha \beta }{T}^{\gamma }$ for constants ${f}_{\gamma }^{\alpha \beta }$ [39, 40]. It is called the BRST quantization of a gauge theory which is the generalization of the Faddeev–Popov path integral quantization of a gauge theory. Fradkin and Vilkovisky pointed out that in such relativistic gauge theories, it is necessary to include the time derivative of all Lagrange multiplier fields to be the GFCs, i.e. ${\partial }_{t}{\lambda }_{\alpha }({\bf{r}},t)+{F}_{\alpha }({\pi }_{{\lambda }_{\alpha }})=0$ with an arbitrary function ${F}_{\alpha }({\pi }_{{\lambda }_{\alpha }})$ [39]. For Dirac’s first-class constrained systems, the BRST symmetry gives a criterion whether the GFC is correct. Applying their approach to the present case with an Abelian constraint Ti(t) = 0 and taking F(πλi) = ξπλi, the BRST invariant Lagrangian is given by11 ). Therefore, the BRST quantization of the gauge theory is exactly equivalent to the conventional path integral quantization. Notice that this BRST quantization may also be applied to non-Abelian gauge theory such as the SU(2) gauge theory of the slave boson [20]. In the quantization of the non-Abelian gauge theory, the determinant for the non-Abelian gauge theory will not be easily treated without introducing the Faddeev–Popov ghosts.
$\begin{eqnarray}{L}_{\mathrm{BRST}}={L}_{\mathrm{eff}}+\displaystyle \sum _{i}{\bar{u}}_{i}{\partial }_{t}^{2}{u}_{i},\end{eqnarray}$
where ui and ${\bar{u}}_{i}$ are the Faddeev–Popov ghost and anti-ghost fields, the fermionic fields obeying $\{{u}_{i},{\bar{u}}_{j}\}={\delta }_{{ij}}$. It is easy to check LBRST is invariant under the BRST transformations $\begin{eqnarray}{\delta }_{{\rm{B}}}{u}_{i}=0,{\delta }_{{\rm{B}}}{\bar{u}}_{i}=\epsilon {\dot{\lambda }}_{i}/\xi ,{\delta }_{{\rm{B}}}{\lambda }_{i}=-\epsilon {\dot{u}}_{i},\end{eqnarray}$
where ε is an anti-commuting constant with ε2 = 0 while the slave particles’ transformations are hi, di, fiσ varying a local phase ${{\rm{e}}}^{-{\rm{i}}\epsilon {u}_{i}(t)}=1-{\rm{i}}\epsilon {u}_{i}(t)$. It is easy to check that ${\delta }_{B}^{2}=0$. This is called the nilpotency of the BRST transformation. The BRST quantized path integral is given by $\begin{eqnarray*}{W}_{\mathrm{BRST}}=\int \displaystyle \prod _{i,t}{\rm{d}}{{\rm{\Phi }}}_{i}^{\dagger }(t){\rm{d}}{{\rm{\Phi }}}_{i}(t){\rm{d}}{\lambda }_{i}(t){\rm{d}}{\bar{u}}_{i}(t){\rm{d}}{u}_{i}(t){{\rm{e}}}^{{\rm{i}}\int {\rm{d}}{{tL}}_{\mathrm{BRST}}}.\end{eqnarray*}$
Integrating over the ghost fields, the path integral WBRST recovers the path integral (The benefits gained from the BRST quantization are that:
| 1. | (1) Because the BRST symmetry is a global symmetry with respect to the fermionic constant ε, one can define the conservation fermionic charge from Nother’s theorem $\begin{eqnarray}Q=\displaystyle \sum _{i}{u}_{i}({\dot{\lambda }}_{i}/\xi +{T}_{i})-\displaystyle \sum _{i}{\dot{u}}_{i}\displaystyle \frac{\delta }{\delta {\lambda }_{i}}.\end{eqnarray}$ All the physical states which are ghost-free obey $\begin{eqnarray}Q| \mathrm{Phys}\rangle =0.\end{eqnarray}$ This recovers Ti = 0 and ${\dot{\lambda }}_{i}(t)=0$ because ξ is an arbitrary constant. From the gauge theory point of view, ${\lambda }_{i}(t)=\bar{\lambda }$ is also a GFC, which removes the redundant degrees of freedom but the constraint Ti = 0 is relaxed to ⟨Ti⟩ = 0. This brings other unphysical degrees of freedom into the quantum state space. The gauge theory developed in [18–21] tried to solve this problem in a different way from ours. |
| 2. | (2) The nilpotency Q2 = 0 resembles the external differential operator d2 = 0 in the deRahm cohomology. The constraint ( |
| 3. | (3) Introducing the ghost fields greatly simplifies the quantization of the non-Abelian gauge theory which we will not involve in here. For the Abelian gauge theory considered in this paper, the ghost fields are decoupled to the gauge field and can be integrated away. Therefore, we will use the path integral ( $\begin{eqnarray}Z=\int \displaystyle \prod _{i,\tau }{\rm{d}}{{\rm{\Phi }}}_{i}^{\dagger }(\tau ){\rm{d}}{{\rm{\Phi }}}_{i}(\tau ){\rm{d}}{\lambda }_{i}(\tau ){{\rm{e}}}^{-{\int }_{0}^{\beta }{\rm{d}}\tau {L}_{\mathrm{eff}}},\end{eqnarray}$ where β = 1/T and ${\dot{\lambda }}_{i}(\tau )\equiv {\partial }_{\tau }\lambda (\tau )$. |
3. The Hubbard model at half-filling
To be concrete, we take the repulsive Hubbard model on a square lattice at half-filling as an example. The model Hamiltonian is given by
$\begin{eqnarray}H=-t\displaystyle \sum _{\langle {ij}\rangle ,\sigma }{c}_{i\sigma }^{\dagger }{c}_{j\sigma }+U\displaystyle \sum _{i}{n}_{i\uparrow }{n}_{i\downarrow },\end{eqnarray}$
where the hopping t is fixed in between nearest neighboring sites. U is the on-site Hubbard repulsion and ${n}_{i\sigma }={c}_{i\sigma }^{\dagger }{c}_{i\sigma }$. In the slave boson representation, the Hamiltonian reads $\begin{eqnarray}\begin{array}{rcl}{H}_{\lambda } & = & -t\displaystyle \sum _{\langle {ij}\rangle }[{\chi }_{{ij}}^{f}{\chi }_{{ji}}^{b}+{{\rm{\Delta }}}_{{ij}}^{f\dagger }{{\rm{\Delta }}}_{{ij}}^{b}+{\rm{h}}.{\rm{c}}.]+U\displaystyle \sum _{i}{d}_{i}^{\dagger }{d}_{i}\\ & & -\displaystyle \sum _{i}{\lambda }_{i}{T}_{i}-\mu \displaystyle \sum _{i}({d}_{i}^{\dagger }{d}_{i}-{h}_{i}^{\dagger }{h}_{i}),\end{array}\end{eqnarray}$
where ${\chi }_{{ij}}^{f}={\sum }_{\sigma }{f}_{i\sigma }^{\dagger }{f}_{j\sigma }$, ${\chi }_{{ij}}^{b}={h}_{i}^{\dagger }{h}_{j}-{d}_{i}^{\dagger }{d}_{j}$, ${{\rm{\Delta }}}_{{ij}}^{f}={\sum }_{\sigma }\sigma {f}_{i,-\sigma }{f}_{j\sigma }$, and ${{\rm{\Delta }}}_{{ij}}^{b}={d}_{i}{h}_{j}+{h}_{i}{d}_{j}$.As we have argued, to study the detailed properties the various phases, we need to do various MF approximations fluctuated by the spatial components of the gauge field which is not the task in this work. We only restrict on the fluctuation from λi(t) and examine the instability of the BCS MF states at half-filling for small U. The original Hubbard model at half-filling is metallic for small U, while the previous slave boson MF theory gave a SC phase [41, 42]. In this SC phase, charge and spin excitations are gapless and the slave bosons condense [41]. Neglecting the boson fluctuation of the condensate, the effective Lagrangian reads
$\begin{eqnarray}\begin{array}{rcl}{L}_{\mathrm{eff}}^{s} & = & \displaystyle \sum _{{\boldsymbol{k}}\sigma }{f}_{{\boldsymbol{k}}\sigma }^{\dagger }\omega {f}_{{\boldsymbol{k}}\sigma }-\displaystyle \frac{t{\chi }^{b}}{2}\displaystyle \sum _{{\boldsymbol{k}}\sigma }{k}^{2}{f}_{{\boldsymbol{k}}\sigma }^{\dagger }{f}_{{\boldsymbol{k}}\sigma }\\ & & -\displaystyle \sum _{i}\left(\displaystyle \frac{1}{2\xi }{\dot{\lambda }}_{i}^{2}-{\lambda }_{i}{n}_{{fi}}\right),\end{array}\end{eqnarray}$
where $t{\chi }^{b}=2t{\rho }_{h}={m}_{0}^{-1}$ with ρh the holon density and ${n}_{{fi}}={\sum }_{\sigma }{f}_{i\sigma }^{\dagger }{f}_{i\sigma }$. We take $\begin{eqnarray}{\lambda }_{i}=\bar{\lambda }+{{ga}}_{i},\end{eqnarray}$
where g is a constant which is arbitrary according to the constraint Ti = 0. The a-dependent part in ${L}_{\mathrm{eff}}^{s}$ reads $\begin{eqnarray}{L}_{a,f}=g\displaystyle \sum _{{\boldsymbol{kq}}\sigma }{a}_{{\boldsymbol{k}}}{f}_{{\boldsymbol{k}}+{\boldsymbol{q}}\sigma }^{\dagger }{f}_{{\bf{q}}\sigma }-\displaystyle \frac{{g}^{2}}{2\xi }\displaystyle \sum _{{\boldsymbol{k}}}{\dot{a}}_{{\boldsymbol{k}}}{\dot{a}}_{-{\boldsymbol{k}}},\,\end{eqnarray}$
where the first term is the interaction between the fluctuation field ${a}_{{\boldsymbol{k}}}={\sum }_{i}{a}_{i}{{\rm{e}}}^{-{\rm{i}}{\boldsymbol{k}}\cdot {{\boldsymbol{r}}}_{i}}$ and the spinon. The perturbation calculation for a small g finds that the renormalization constant ZF does not vanish at spinon Fermi surface and then the gapless spinon liquid is stable against the perturbation fluctuation.To see the instability, we rewrite La,f28 )
$\begin{eqnarray}\begin{array}{rcl}{L}_{a,f}(\omega ) & = & g\displaystyle \sum _{{\boldsymbol{k}}{\boldsymbol{q}}\nu \sigma }{a}_{{\boldsymbol{k}}}(\omega ){f}_{{\boldsymbol{k}}+{\boldsymbol{q}}\sigma }^{\dagger }(\omega +\nu ){f}_{{\bf{q}}\sigma }(\nu )\\ & & -{\omega }^{2}\displaystyle \frac{{g}^{2}}{2\xi }\displaystyle \sum _{{\boldsymbol{k}}}{a}_{{\boldsymbol{k}}}(\omega ){a}_{-{\boldsymbol{k}}}(-\omega ).\end{array}\end{eqnarray}$
Integrating away ak(ω), the effective interacting Lagrangian between the spinons reads $\begin{eqnarray}\begin{array}{rcl}{L}_{\mathrm{int}}^{s}\left(\omega \right) & = & \frac{\xi }{2{\omega }^{2}}\displaystyle \sum _{{\boldsymbol{k}}{\boldsymbol{q}}\nu \sigma }{f}_{-{\boldsymbol{k}}-{\boldsymbol{q}},\sigma }^{\dagger }\left(-\omega +\nu \right){f}_{-{\boldsymbol{k}},\sigma }\left(\nu \right)\\ & & \times \displaystyle \sum _{{\boldsymbol{q}}^{\prime} \nu ^{\prime} \sigma ^{\prime} }{f}_{{\boldsymbol{k}}+{\boldsymbol{q}}^{\prime} ,\sigma ^{\prime} }^{\dagger }\left(\omega +\nu ^{\prime} \right){f}_{{\boldsymbol{k}}\sigma ^{\prime} }\left(\nu ^{\prime} \right).\end{array}\end{eqnarray}$
The pairing Hamiltonian is then given by $\begin{eqnarray}\begin{array}{rcl}{H}_{\mathrm{pair}}\left(\omega \right) & = & -\frac{\xi }{2{\omega }^{2}}\displaystyle \sum _{{\boldsymbol{k}}{\boldsymbol{q}}\nu \sigma }{f}_{-{\boldsymbol{k}}-{\boldsymbol{q}},\sigma }^{\dagger }\left(-\omega -\nu \right){f}_{{\boldsymbol{k}}+{\boldsymbol{q}},-\sigma }^{\dagger }\left(\omega +\nu \right)\\ & & \times {f}_{-{\boldsymbol{k}},-\sigma }\left(-\nu \right){f}_{{\boldsymbol{k}},\sigma }\left(\nu \right).\end{array}\end{eqnarray}$
For the pairing parameter $\tilde{{\rm{\Delta }}}=\xi \langle {f}_{-{\boldsymbol{k}},\uparrow }{f}_{{\boldsymbol{k}},\downarrow }\rangle $, the MF energy is given by $\begin{eqnarray}{\omega }^{2}-\displaystyle \frac{4{\tilde{{\rm{\Delta }}}}^{2}}{{\omega }^{4}}-{\left(\displaystyle \frac{{{\boldsymbol{k}}}^{2}}{2{m}_{0}}-\bar{\lambda }\right)}^{2}=0.\end{eqnarray}$
When $\tilde{{\rm{\Delta }}}=0$, the spinon is gapless at spin Fermi surface. When $\tilde{{\rm{\Delta }}}\ne 0$, the spinon has a pairing gap is ${{\rm{\Delta }}}^{* }=\pm {\tilde{{\rm{\Delta }}}}^{1/3}$ and the spinon liquid is gapped. In this case, the spinon Fermi surface is not stable. Nor is the SC state of the Hubbard model at half-filling. It is metallic state because the charge excitation is still gapless. To see the paired spinon liquid has a lower ground state energy than that of the gapless spinon, we consider the ground state energy $\begin{eqnarray}{E}_{g}=\displaystyle \sum _{{\boldsymbol{k}}}\left(-\displaystyle \frac{2{\xi }_{{\boldsymbol{k}}}^{2}+4{\tilde{{\rm{\Delta }}}}^{2}/{\omega }^{2}}{2\omega }\right)+\cdots \equiv {E}_{g}^{* }+\cdots ,\end{eqnarray}$
where ${\xi }_{{\boldsymbol{k}}}=\tfrac{{{\boldsymbol{k}}}^{2}}{2{m}_{0}}-\bar{\lambda }$ and ⋯ is $\tilde{{\rm{\Delta }}}$-independent while ω can be solved from the dispersion relation ( $\begin{eqnarray}\begin{array}{rcl}{\omega }^{2} & = & \displaystyle \frac{1}{3}\left({\xi }_{{\boldsymbol{k}}}^{2}+\displaystyle \frac{{\xi }_{{\boldsymbol{k}}}^{4}}{\sqrt[3]{{\xi }_{{\boldsymbol{k}}}^{6}+54{\tilde{{\rm{\Delta }}}}^{2}+6\sqrt{3}| \tilde{{\rm{\Delta }}}| \sqrt{{\xi }_{{\boldsymbol{k}}}^{6}+27{\tilde{{\rm{\Delta }}}}^{2}}}}\right.\\ & & \left.+\sqrt[3]{{\xi }_{{\boldsymbol{k}}}^{6}+54{\tilde{{\rm{\Delta }}}}^{2}+6\sqrt{3}| \tilde{{\rm{\Delta }}}| \sqrt{{\xi }_{{\boldsymbol{k}}}^{6}+27{\tilde{{\rm{\Delta }}}}^{2}}}\right).\end{array}\end{eqnarray}$
In figure 1, we plot a typical curve of the $\tilde{{\rm{\Delta }}}$-dependent part of Eg, i.e. ${E}_{g}^{* }$ versus $\tilde{{\rm{\Delta }}}$. Here the lattice value of ξk is used and the summation of k runs over the first Brillouin zone. $\bar{\lambda }$ is converted to the spinon chemical potential μf. We see that the larger $\tilde{{\rm{\Delta }}}$ is, the lower ${E}_{g}^{* }$ is. On the other hand, in the MF theory, the gap $\tilde{{\rm{\Delta }}}$ is determined by the self-consistent equation $\begin{eqnarray}\tilde{{\rm{\Delta }}}=\displaystyle \frac{\xi }{{\omega }^{2}}\displaystyle \sum _{{\boldsymbol{k}}}\displaystyle \frac{\tilde{{\rm{\Delta }}}}{\sqrt{{\xi }_{{\boldsymbol{k}}}^{2}+{\left(\tfrac{2\tilde{{\rm{\Delta }}}}{{\omega }^{2}}\right)}^{2}}}.\end{eqnarray}$
In the large $\tilde{{\rm{\Delta }}}$ limit, ${\omega }^{2}\approx 2.314{\tilde{{\rm{\Delta }}}}^{2/3}$ and the self-consistent equation is approximated by $\begin{eqnarray}\displaystyle \sum _{{\boldsymbol{k}}}\displaystyle \frac{\xi }{2\tilde{{\rm{\Delta }}}}\approx 1,\end{eqnarray}$
i.e. the ξ-independent gap is ${\rm{\Delta }}=\tilde{{\rm{\Delta }}}/\xi \approx {N}_{e}/2$. Thus, the spinon is gapped. This means that the s-wave SC at half-filling for small U is unstable and the system is a conventional metal.
Figure 1. Numerical result of ${E}_{g}^{* }$ of $\tilde{{\rm{\Delta }}}$. We choose the dispersion of the normal state to be ${\xi }_{{\boldsymbol{k}}}=-t(\cos {k}_{x}+\cos {k}_{y})-{\mu }_{f}$, with t = 0.5 and μf = − 0.01. |
4. t–J model
We are going to the large U limit. It was known that the SC MF theory of the t–J model has a SC Tc ∼ 1000 K. Let us see if λi fluctuation can suppress it. The t–J model Hamiltonian on the square lattice is given by11 ) by ${L}_{\mathrm{eff}}^{{tJ}}$, one has the path integral for the t–J model. If we ignore the condensed holon fluctuation, the effective low-lying Lagrangian in the SC MF state is given by25 ). The holon is condensed and the SC gap function is ρhΔk. Similarly, after integral over a-field, there is an additional pairing term, Hpair in the MF Hamiltonian. The MF dispersion is then given by
$\begin{eqnarray}{H}_{t-J}=-t\displaystyle \sum _{\langle {ij}\rangle \sigma }{c}_{i\sigma }^{\dagger }{c}_{j\sigma }+J\displaystyle \sum _{\langle {ij}\rangle \sigma }\left({{\boldsymbol{S}}}_{i}\cdot {{\boldsymbol{S}}}_{j}-\frac{1}{4}{n}_{i}{n}_{j}\right),\end{eqnarray}$
where ${S}_{i}^{a}=\tfrac{1}{2}{\sum }_{\sigma ,\sigma ^{\prime} }{c}_{i\sigma }^{\dagger }{\sigma }_{\sigma \sigma ^{\prime} }^{a}{c}_{i\sigma ^{\prime} }$ are the spin operators and σa (a = x, y, z) are Pauli matrices. The constraint is that there is no double occupation on each lattice site. The slave boson decomposition is now ${c}_{i\sigma }^{\dagger }={f}_{i\sigma }^{\dagger }{h}_{i}$ and the local constraints are $\begin{eqnarray}{T}_{i}^{tJ}={h}_{i}^{\dagger }{h}_{i}+\displaystyle \sum _{\sigma }{f}_{i\sigma }^{\dagger }{f}_{i\sigma }-1=0.\end{eqnarray}$
The effective slave boson t–J Lagrangian then reads $\begin{eqnarray}\begin{array}{rcl}{L}_{\mathrm{eff}}^{tJ} & = & \frac{J}{4}\displaystyle \sum _{\langle {ij}\rangle }[| {\chi }_{{ij}}^{f}{| }^{2}+| {{\rm{\Delta }}}_{{ij}}^{f}{| }^{2}-({\chi }_{{ij}}^{f\dagger }\displaystyle \sum _{\sigma }{f}_{i\sigma }^{\dagger }{f}_{j\sigma }+{\rm{h}}.{\rm{c}}.)]\\ & & +\displaystyle \sum _{\langle {ij}\rangle }\frac{J}{4}[{{\rm{\Delta }}}_{{ij}}^{f}({f}_{i\uparrow }^{\dagger }{f}_{j\downarrow }^{\dagger }-{f}_{i\downarrow }^{\dagger }{f}_{j\uparrow }^{\dagger })+{\rm{h}}.{\rm{c}}.]\\ & & +\displaystyle \sum _{i}[{h}_{i}^{\dagger }({\rm{i}}{\partial }_{t}-\mu ){h}_{i}+\mu x]+\displaystyle \sum _{i\sigma }{f}_{i\sigma }^{\dagger }{\rm{i}}{\partial }_{t}{f}_{i\sigma }\\ & & -\displaystyle \sum _{i}\left(\frac{1}{2\xi }{\dot{\lambda }}_{i}^{2}-{\lambda }_{i}{T}_{i}^{{tJ}}\right)-t\displaystyle \sum _{\langle {ij}\rangle }{h}_{i}{h}_{j}^{\dagger }{f}_{i\sigma }^{\dagger }{f}_{j\sigma }.\end{array}\end{eqnarray}$
Replacing Leff in equation ( $\begin{eqnarray}\begin{array}{rcl}{\bar{L}}_{\mathrm{eff}}^{tJ-\mathrm{sc}} & = & \displaystyle \sum _{{\boldsymbol{k}}\sigma }{f}_{{\boldsymbol{k}}\sigma }^{\dagger }\omega {f}_{{\boldsymbol{k}}\sigma }-\displaystyle \sum _{{\boldsymbol{k}}\sigma }{\xi }_{{\boldsymbol{k}}}{f}_{{\boldsymbol{k}}\sigma }^{\dagger }{f}_{{\boldsymbol{k}}\sigma }\\ & & -\displaystyle \sum _{{\boldsymbol{k}}}({{\rm{\Delta }}}_{k}{f}_{{\boldsymbol{k}}\uparrow }^{\dagger }{f}_{-{\boldsymbol{k}}\downarrow }^{\dagger }+{\rm{h}}{\rm{.}}{\rm{c}}.)+{L}_{a,f},\end{array}\end{eqnarray}$
where ${\xi }_{{\boldsymbol{k}}}=-(J{\chi }^{f}/4+t{\rho }_{h})(\cos {k}_{x}+\cos {k}_{y})-{\mu }_{f}$ and ${{\rm{\Delta }}}_{{\boldsymbol{k}}}=J({{\rm{\Delta }}}_{x}\cos {k}_{x}+{{\rm{\Delta }}}_{y}\cos {k}_{y})$ with μf being the chemical potential of the spinon and ρh is the holon density. La,f is given by ( $\begin{eqnarray}{\omega }^{2}-{{\rm{\Delta }}}_{k}^{2}{\left(1+\displaystyle \frac{2\xi }{J{\omega }^{2}}\right)}^{2}-{\xi }_{{\boldsymbol{k}}}^{2}=0.\end{eqnarray}$
At the Fermi surface, the energy gap $\omega ({{\boldsymbol{k}}}_{{\boldsymbol{F}}})={{\rm{\Delta }}}_{{{\boldsymbol{k}}}_{{\boldsymbol{F}}}}^{* }$ is determined by $\begin{eqnarray}{\hat{{\rm{\Delta }}}}_{{\rm{F}}}^{6}-{\left({\hat{{\rm{\Delta }}}}_{{\rm{F}}}^{2}+{A}_{{\rm{F}}}\right)}^{2}=0,\end{eqnarray}$
where ${\hat{{\rm{\Delta }}}}_{{\rm{F}}}=\tfrac{{{\rm{\Delta }}}_{{{\boldsymbol{k}}}_{{\rm{F}}}}^{* }}{{{\rm{\Delta }}}_{{{\boldsymbol{k}}}_{{\rm{F}}}}}$ and ${A}_{{\rm{F}}}=\tfrac{2\xi }{J{{\rm{\Delta }}}_{{{\boldsymbol{k}}}_{{\rm{F}}}}^{2}}$.Because the MF approximation breaks the gauge symmetry, the calculation results are dependent on the gauge fixing parameter ξ. We treat ξ as a phenomenological parameter and so is AF(ξ). Obviously, Equation (38 ) has solutions in $0\leqslant {\hat{{\rm{\Delta }}}}_{{\rm{F}}}\leqslant 1$ for ξ ≠ 0 and ${\hat{{\rm{\Delta }}}}_{{\rm{F}}}=0$ or 1 only if ξ = 0. This means that the SC MF phase is stable but the gap parameter ${{\rm{\Delta }}}_{{k}_{{\rm{F}}}}^{* }$ is suppressed, e.g. when AF ∼ −10−2 ± 10−3, y ∼ 0.1. The renormalized pairing gap may be suppressed by one order from the MF gap ${{\rm{\Delta }}}_{{k}_{{\rm{F}}}}$. This may explain why ${T}_{c}\propto {\rho }_{{\rm{h}}}{{\rm{\Delta }}}_{{k}_{{\rm{F}}}}^{* }$ in cuprates is one order lower than the MF estimation.
We can directly solve the dispersion relation (37 ) and find thatAppendix ). The ground state energy of the BCS state is given byAppendix ).
$\begin{eqnarray}\omega ={\omega }_{{\boldsymbol{k}}}=\sqrt{{\xi }_{{\boldsymbol{k}}}^{2}+{{\rm{\Delta }}}_{{\boldsymbol{k}}}^{* 2}},\end{eqnarray}$
where ${{\rm{\Delta }}}_{k}^{* }$ is the renormalized gap (see $\begin{eqnarray}\begin{array}{rcl}{E}_{g}\left(\xi \right) & = & \frac{1}{2J}\displaystyle \sum _{{\boldsymbol{k}}}{{\rm{\Delta }}}_{k}^{* 2}\left(\xi \right)-\displaystyle \sum _{{\boldsymbol{k}}}{\omega }_{{\boldsymbol{k}}}\left(\xi \right)+\cdots \\ & \equiv & {E}_{g}^{* }+\cdots ,\,\end{array}\end{eqnarray}$
where ⋯ is ξ-independent terms. With the material data of cuprates, we plot the ξ-dependent part of the d-wave SC MF ground state energies ${E}_{g}^{* }$ varying as ξ in figure 2. We find that ξ ∼ −0.002 minimizes ${E}_{g}^{* }$. With this value of ξ, we have $\tfrac{{{\rm{\Delta }}}_{{\boldsymbol{k}}}^{* }}{{{\rm{\Delta }}}_{{\boldsymbol{k}}}}\sim 0.2$ for kx = 0, ky = π. This lowers the MF Tc a factor of one-fifth, i.e. from ${T}_{c}^{\mathrm{MF}}\sim 1000$ K to ${T}_{c}^{* }\sim 200$ K (for more details, see 
Figure 2. Numerical result of ${E}_{{\rm{g}}}^{* }$ of ξ. We choose the material data of cuprates [44], i.e. J ∼ 0.12 eV, χf ∼ 0.2 − 0.3, ρh ∼ 0.18 − 0.25, and μf ∼ −0.05 eV. The four curves in the diagram, from the top to bottom, correspond to Jχf/4 + tρh = 0.10 eV, 0.11 eV, 0.12 eV, and 0.13 eV, respectively. |
5. Conclusions and prospects
We have properly dealt with the local constraint conditions in the slave boson representation of the strongly correlated systems. We argued that as a gauge theory with Dirac’s first-class constraint, taking the GFC that removes the redundant gauge degrees of freedom must be consistent with the constraint. The BRST quantization is a consistent method to do that although the final path integral for the Abelian gauge theory is decoupled to the ghost fields. We have applied our theory to the Hubbard model at half-filling and found that the ground state of the system in small U is indeed a conventional metal. We showed that the MF s-wave SC state obtained by the slave boson representation in a previous study is not stable against the gauge fluctuation of the gauge field λi(t). For the strong coupling system, we studied the t–J model. We focused on the gauge fluctuation to the d-wave SC gap and found that it was substantially suppressed to a factor of one-fifth. For cuprates, this means that the MF SC Tc is lowered from 1000 to 200 K. As we have mentioned, the gauge fluctuation from the spatial components of the gauge field was not considered. It might further reduce Tc to be comparable to that of the cuprates materials.
Historically, the MF phase diagram of the t–J model was studied in the early days when the high Tc SC was found in cuprates. The gauge fluctuation to various MF states was studied. However, as we see here, the GFC might not be introduced properly because the spatial components of the gauge field also play a role of the Lagrange multiplier to the constraint on the currents. Additional constraints are also needed. This may endow the gauge field with dynamics. Then, the gauge invariant physical quantities can be calculated with perturbation theory. For instance, one can calculate the renormalized pairing gap by Dyson’s equation using the perturbation theory
$\begin{eqnarray}\displaystyle \frac{{{\rm{\Delta }}}_{{\boldsymbol{k}}}^{* }(\omega )}{{{\rm{\Delta }}}_{{\boldsymbol{k}}}}=\displaystyle \frac{{G}_{\uparrow \downarrow }^{-1}({\boldsymbol{k}},\omega )}{{G}_{0\uparrow \downarrow }^{-1}({\boldsymbol{k}})}=1-\displaystyle \frac{{{\rm{\Sigma }}}_{\uparrow \downarrow }({\boldsymbol{k}},\omega )}{{{\rm{\Delta }}}_{{\boldsymbol{k}}}},\end{eqnarray}$
where G↑↓ is the anomalous spinon Green’s function and ${G}_{0\uparrow \downarrow }^{-1}={{\rm{\Delta }}}_{k}$. Σ↑↓ is the anomalous self-energy of the spinon. This also implies Δk is suppressed to ${{\rm{\Delta }}}_{k}^{* }(\omega )$ by the fluctuation. However, to establish a complete gauge theory with all components of the gauge fields, many symmetry considerations must be taken care of. We leave them to the further works.Acknowledgments
This work is in memory of Professor Zhong-Yuan Zhu for his discussions with YY in the possible application of the BRST quantization to strongly correlated systems thirty years ago. The authors thank Professor Qian Niu for his insightful comments and resultful discussions. We are grateful to Jianhui Dai and Long Liang for useful discussions. This work is supported by NNSF of China with No. 12174067.
Appendix. Dispersion and renormalized pairing gap
The SC dispersion relation in the t–J model can be solved by equation (37 )A1 ) readsA4 ) is
$\begin{eqnarray}{\omega }^{2}-{{\rm{\Delta }}}_{{\boldsymbol{k}}}^{2}{\left(1+\displaystyle \frac{2\xi }{J{\omega }^{2}}\right)}^{2}-{\xi }_{{\boldsymbol{k}}}^{2}=0.\end{eqnarray}$
Defining $\begin{eqnarray*}{\hat{E}}_{{\boldsymbol{k}}}^{2}=\frac{{\omega }^{2}}{{{\rm{\Delta }}}_{{\boldsymbol{k}}}^{2}},\,\,{\hat{\xi }}_{{\boldsymbol{k}}}^{2}=\frac{{\xi }_{{\boldsymbol{k}}}^{2}}{{{\rm{\Delta }}}_{{\boldsymbol{k}}}^{2}},{A}_{k}=-\frac{2\xi }{J{{\rm{\Delta }}}_{{\boldsymbol{k}}}^{2}}.\end{eqnarray*}$
Equation ( $\begin{eqnarray}{\hat{E}}_{{\boldsymbol{k}}}^{2}-{\left(1-\displaystyle \frac{{A}_{{\boldsymbol{k}}}}{{\hat{E}}_{{\boldsymbol{k}}}^{2}}\right)}^{2}-{\hat{\xi }}_{{\boldsymbol{k}}}^{2}=0.\end{eqnarray}$
For ${E}_{{\boldsymbol{k}}}^{2}\gt 0$, we have $\begin{eqnarray}{\hat{E}}^{6}-({\hat{\xi }}_{{\boldsymbol{k}}}^{2}+1){\hat{E}}_{{\boldsymbol{k}}}^{4}+2{A}_{k}{\hat{E}}_{{\boldsymbol{k}}}^{2}-{A}_{{\boldsymbol{k}}}^{2}=0.\end{eqnarray}$
Define $y={E}_{{\boldsymbol{k}}}^{2}-({\hat{\xi }}_{{\boldsymbol{k}}}^{2}+1)/3$, the above equation becomes $\begin{eqnarray}{y}^{3}+{py}+q=0,\end{eqnarray}$
with $\begin{eqnarray*}\begin{array}{rcl}p & = & \displaystyle \frac{6{A}_{{\boldsymbol{k}}}-{\left({\hat{\xi }}_{{\boldsymbol{k}}}^{2}+1\right)}^{2}}{3},\\ q & = & \displaystyle \frac{-2{\left({\hat{\xi }}_{{\boldsymbol{k}}}^{2}+1\right)}^{3}+18({\hat{\xi }}_{{\boldsymbol{k}}}^{2}+1){A}_{{\boldsymbol{k}}}-27{A}_{{\boldsymbol{k}}}^{2}}{27}.\end{array}\end{eqnarray*}$
Then $\begin{eqnarray*}\begin{array}{rcl}{ \mathcal D } & = & \displaystyle \frac{{q}^{2}}{4}+\displaystyle \frac{{p}^{3}}{27}\\ & = & \displaystyle \frac{1}{108}{A}_{{\boldsymbol{k}}}^{2}\left(27{A}_{{\boldsymbol{k}}}^{2}-4{A}_{{\boldsymbol{k}}}(9{\xi }_{{\boldsymbol{k}}}^{2}+1)+4{\xi }_{{\boldsymbol{k}}}^{2}{\left({\xi }_{{\boldsymbol{k}}}^{2}+1\right)}^{2}\right).\end{array}\end{eqnarray*}$
It is easy to see that for a small Ak, up to ${ \mathcal O }({A}_{k}^{2})$, ${ \mathcal D }\gt 0$. Therefore, the solution of equation ( $\begin{eqnarray}\begin{array}{rcl}y & = & \sqrt[3]{-\displaystyle \frac{q}{2}+\sqrt{{\left(\displaystyle \frac{q}{2}\right)}^{2}+{\left(\displaystyle \frac{p}{3}\right)}^{3}}}\\ & & +\sqrt[3]{-\displaystyle \frac{q}{2}-\sqrt{{\left(\displaystyle \frac{q}{2}\right)}^{2}+{\left(\displaystyle \frac{p}{3}\right)}^{3}}}.\end{array}\end{eqnarray}$
Then $\begin{eqnarray}{\hat{E}}_{{\boldsymbol{k}}}^{2}=y+({\hat{\xi }}_{{\boldsymbol{k}}}^{2}+1)/3.\end{eqnarray}$
This gives the dispersion relation. Defining the renormalized gap by $\begin{eqnarray*}{{\rm{\Delta }}}_{{\boldsymbol{k}}}^{* 2}={E}_{{\boldsymbol{k}}}^{2}-{\xi }_{{\boldsymbol{k}}}^{2},\end{eqnarray*}$
we have $\begin{eqnarray}{{\rm{\Delta }}}_{{\boldsymbol{k}}}^{* 2}={{\rm{\Delta }}}_{{\boldsymbol{k}}}^{2}\left(y-\displaystyle \frac{1}{3}\left(2{\hat{\xi }}_{{\boldsymbol{k}}}^{2}-1\right)\right).\end{eqnarray}$
The ground state energy is given by $\begin{eqnarray}{E}_{g}=\displaystyle \frac{1}{2J}\displaystyle \sum _{{\boldsymbol{k}}}{{\rm{\Delta }}}_{{\boldsymbol{k}}}^{* 2}-\displaystyle \sum _{{\boldsymbol{k}}}{\omega }_{{\boldsymbol{k}}}+\cdots \equiv {E}_{g}^{* }+\cdots ,\end{eqnarray}$
where the ⋯ stands for the ξ-independent part. In the numerical calculation, we choose the parameters of cuprates, i.e. J ∼ 0.12 eV, χf ∼ 0.2 − 0.3, ρh ∼ 0.18 − 0.25, and Jχf/4 + tρh ∼ 0.1 − 0.13 eV [44]. The spinon chemical potential μf is determined by $\begin{eqnarray}1-{\rho }_{h}={\rho }_{f}=-\displaystyle \frac{1}{N}\displaystyle \frac{\partial G}{\partial {\mu }_{f}},\end{eqnarray}$
where ρf is the spinon density and G is the Gibbs free energy. To the zeroth order of ξ, it reduces to $\begin{eqnarray}{\rho }_{f}=\displaystyle \frac{1}{N}\displaystyle \sum _{{\boldsymbol{k}}}\left(1-\displaystyle \frac{{\xi }_{{\boldsymbol{k}}}}{{\omega }_{{\boldsymbol{k}}}}\tanh \left(\displaystyle \frac{{\omega }_{{\boldsymbol{k}}}}{2{k}_{{\rm{B}}}T}\right)\right),\end{eqnarray}$
where N is the number of total particles. By varying T from 200 to 1000 K, μf keeps ∼ −0.05 eV. From the numerical results (see figure 2 in the main text), ξ = −0.002 minimizes the ground state energy. At kx = 0, ky = π and Δk ∼ 0.17 eV, ${{\rm{\Delta }}}_{k}^{* }\sim 0.0318\,\mathrm{eV}$. Therefore, the critical temperature becomes ${T}_{c}\sim 1000\tfrac{{{\rm{\Delta }}}_{{\boldsymbol{k}}}^{* }}{{{\rm{\Delta }}}_{{\boldsymbol{k}}}}\sim 200$ K.