its pairing Hamiltonian as $$\begin{aligned}
H_{pairing}=-\frac{3}{4}J_{H}\left(\sum_{k}\gamma_{k}f_{k\uparrow}^{\dag}f_{-k\downarrow}^{\dag}\right)\left(\sum_{k'}\gamma_{k'}f_{-k'\downarrow}f_{k'\uparrow}\right).\nonumber\end{aligned}$$ Then, follow the spirit of classic BCS mean-field theory,[@Xiang] we introduce pairing order parameter $$\begin{aligned}
\Delta=-\frac{3}{4}\sum_{k}\gamma_{k}\langle f_{k\uparrow}^{\dag}f_{-k\downarrow}^{\dag}\rangle=-\frac{3}{4}\sum_{k}\gamma_{k}\langle f_{-k\downarrow}f_{k\uparrow}\rangle,\nonumber\end{aligned}$$ and the corresponding mean-field Hamiltonian reads $$\begin{aligned}
H_{pairing}=J_{H}\Delta\sum_{k}\gamma_{k}(f_{k\uparrow}^{\dag}f_{-k\downarrow}^{\dag}
+f_{-k\downarrow}f_{k\uparrow})+\frac{4J_{H}\Delta^{2}}{3}.\label{eq5}\end{aligned}$$ Note that such pairing Hamiltonian is identical to ones in Refs. , where the RVB-like approach is utilized,[@Anderson2004] i.e. one assumes a spin-singlet pairing order of auxiliary fermion in real-space as $$\Delta_{ij}=\langle f_{i\uparrow}^{\dag}f_{j\downarrow}^{\dag}-f_{i\downarrow}^{\dag}f_{j\uparrow}^{\dag}\rangle  \nonumber$$

Obviously, although our method is more like the conventional antiferromagnetic spin-fluctuation-exchange theory, it is equivalent to the RVB-pairing theory at least in the present mean-field level, which emphasizes the core role of local magnetic interaction.[@Hu2012] In some sense, the mean-field pairing model misses the effect of original pairing interaction, but if we include the fluctuation of pairing interaction, one may identify which microscopic interaction is responsible for pairing.

Now, adding this pairing Hamiltonian into the normal state model Eq. \[eq1\], one can get the desirable model for heavy fermion superconductivity. In practice, the pairing between different heavy fermion quasi-particle bands is usually neglected due to the mismatch of Fermi surface as done in Ref. . Thus, the pairing Hamiltonian is simplified as $$\begin{aligned}
H_{pairing}&&=J_{H}\Delta\sum_{k}\gamma_{k}(\beta_{k}^{2}(A_{k\uparrow}^{\dag}A_{-k\downarrow}^{\dag})
+\alpha_{k}^{2}(B_{k\uparrow}^{\dag}B_{-k\downarrow}^{\dag})\nonumber\\
&&+\alpha_{k}\beta_{k}(A_{k\uparrow}^{\dag}B_{-k\downarrow}^{\dag}+B_{k\uparrow}^{\dag}A_{-k\downarrow}^{\dag})+\mathrm{H}.\mathrm{c}.)+J_{H}\frac{4\Delta^{2}}{3}\nonumber\\
&&\simeq\sum_{k}[J_{H}\Delta_{k}^{A}(A_{k\uparrow}^{\dag}A_{-k\downarrow}^{\dag}+A_{-k\downarrow}A_{k\uparrow})\nonumber\\
&&+J_{H}\Delta_{k}^{B}(B_{k\uparrow}^{\dag}B_{-k\downarrow}^{\dag}+B_{-k\downarrow}B_{k\uparrow})]+J_{H}\frac{4\Delta^{2}}{3}.\label{eq6}\end{aligned}$$

Here we have defined two 