A mechanism for quantum-critical  Planckian metal phase in high-temperature cuprate superconductors

Yung-Yeh Chang
Institute of Physics, Academia Sinica, Taipei 11529, Taiwan
Physics Division, National Center for Theoretical Sciences, Taipei 10617, Taiwan
Department of Electrophysics, National Yang Ming Chiao Tung University, Hsinchu 30010, Taiwan
Khoe Van Nguyen
Department of Electrophysics, National Yang Ming Chiao Tung University, Hsinchu 30010, Taiwan
Kim Remund
Physics Division, National Center for Theoretical Sciences, Taipei 10617, Taiwan
Department of Electrophysics, National Yang Ming Chiao Tung University, Hsinchu 30010, Taiwan

Chung-Hou Chung
Physics Division, National Center for Theoretical Sciences, Taipei 10617, Taiwan
Department of Electrophysics, National Yang Ming Chiao Tung University, Hsinchu 30010, Taiwan
Center for Theoretical and  Computational Physics (CTCP), National Yang Ming Chiao Tung University, Hsinchu 30010, Taiwan

The mysterious metallic phase showing perfect $T$-linear resistivity and a universal scattering rate $1/\tau = \alpha_P k_B T /\hbar$ with a universal prefactor $\alpha_P \sim 1$ and logarithmic-in-temperature singular specific heat coefficient, so-called “Planckian metal phase" was observed in various overdoped high-$T_c$ cuprate superconductors over a finite range in doping. Here, we propose a microscopic mechanism for this exotic state based on quantum-critical bosonic charge Kondo fluctuations coupled to both spinon and a heavy conduction-electron Fermi surfaces within the heavy-fermion formulation of the slave-boson $t$-$J$ model. Using a controlled perturbative renormalization group (RG) analysis, we examine the competition between the pseudogap phase, characterized by Anderson's Resonating-Valence-Bond spin-liquid, and the Fermi-liquid state, characterized by the electron hoping (effective charge Kondo effect). We find a quantum-critical metallic phase with a universal Planckian $\hbar \omega/k_B T$ scaling in scattering rate near a localized-delocalized (pseudogap-to-Fermi liquid) charge  Kondo breakdown transition. Our results are in excellent agreement with the recent experimental observations on optical conductivity (without fine-tuning) in Nat. Commun. 14, 3033 (2023), universal doping-independent field-to-temperature scaling in magnetoresistance in Nature 595, 661 (2021), and the marginal Fermi-liquid spectral function observed in ARPES (Science 366, 1099 (2019)) as well as Hall coefficient in various overdoped cuprates in Nature 595, 661 (2021) and Annu. Rev. Condens. Matter Phys. 10, 409 (2019). Our mechanism offers a microscopic understanding of the quantum-critical Planckian metal phase observed in cuprates and its link to the pseudogap, d-wave superconducting, and Fermi liquid phases.

§ INTRODUCTION

(a) Schematic plots for the dispersive $f$-spinon (blue curve) and the weakly dispersive $\xi$ (orange curve) bands (corresponding to hole doping $\delta$), associated with their chemical potentials, $\mu_f$ and $\mu_\xi$. The green dashed line denotes the energy level (Lagrange multiplier) $\lambda$ for the slave boson. The shaded areas represent filling of the $f$- and $\xi$-bands. $a_0$ here denotes the lattice constant. (b) Fermi surfaces of the $f$-spinon band (blue) and of the $\xi$ band with different levels of hole doping (orange for hole doping $\delta$, green for $1+\delta$, and red in between $\delta$ and $1+\delta$). (c) Main plot: Schematic plot of strange-metal state. The blue dashed curves represent the Kondo-like hopping term. Upper left: Feynman diagram for the interaction vertex of the Kondo-like hopping term. Upper right: the RVB spin-singlet bond in the $f$-spinon band and the spinon-holon bound state (the $\xi$ field).

Over the recent three decades, metallic behavior that cannot be described within the Fermi liquid (FL) paradigm has commonly been observed in a wide variety of strongly correlated quantum materials. Yet, the emergence of such metals is poorly understood. This non-Fermi liquid (NFL) behavior often exists near a quantum phase transition, and shows “strange metal (SM)" phenomena with (quasi-)linear-in-temperature decreasing resistivity and divergent logarithmic-in-temperature specific heat coefficient as $T \to 0$.

Of particularly intriguing class of SM states is the “Planckian metal", observed in the normal state of unconventional superconductors, including cuprate superconductors <cit.>, iron pnictides and chalcogenides <cit.>, organic <cit.> and heavy-fermion compounds <cit.>, and twisted bilayer graphene <cit.>. It shows perfect $T$-linear scattering rate, reaching “Planckian dissipation limit" allowed by quantum mechanics, $1/\tau(T) = \alpha_P k_BT/\hbar$ with $\alpha_P \sim 1$. Since d-wave superconductivity (dSC) emerges out of this exotic state of matter with decreasing temperatures, it was argued that revealing the mystery of the Planckian metal state is the key to understanding the mechanism for high-$T_c$ superconductivity in cuprates <cit.>. Experimental observations in heavy-fermion superconductors suggest that $\alpha$ is a non-universal constant but depends on the strength of the Kondo correlations <cit.>. By contrast, in overdoped cuprates the Planckian perfect $T$-linear scattering rate with the same $\alpha$ persists from very high temperatures ($T \sim 300$K) to very low temperatures ($T \to 0$)  over a wide range  of samples near optimal doping, indicating  that the Planckian state in overdoped cuprates might be a universal feature. Similar universal behavior has also been observed in magic angle twisted bi-layer graphene <cit.>.

The Planckian behavior across the entire temperature range, in particular for the low temperature regime ($T\to 0$) in the  two mentioned systems is unlikely to be explained by phonons whose contributions are mainly at high temperatures. Meanwhile, perfect frequency-to-temperature ($\hbar \omega/k_B T$) scaling from optical conductivity measurement <cit.> and field-to-temperature ($B/T$ ) scaling in magnetoresistance as well as other quantum critical-like properties extending over a range from the critical doping to the end of the dSC dome in various hole-doped <cit.> and electron-doped  <cit.> cuprates strongly suggest that a Planckian state is a quantum critical “phase” <cit.>.  These observations lead to fundamental questions: What is the microscopic mechanism for this exotic phase of matter and its links to quantum criticality and the neighboring phases? The Planckain metal state lies in between the pseudogap and Fermi liquid phases with localized and itinerant characters of electrons, respectively. This points to an appealing scenario in which the Planckian metal phase may arise near a possible localized-delocalized quantum critical point due to competition between the pseudogap and Fermi liquid phases, and the d-wave superconductivity is  reached by condensing this quantum critical metal.

It is challenging to develop a controlled approach for this exotic state of matter. Recently, this state has been realized theoretically via controlled large-$N$ approach of the $t$-$J$ model with random hoping and exchange couplings, the “SYK model" <cit.> and in the Hubbard model by Quantum Monte Carlo <cit.>. Nevertheless, a microscopic understanding of this state in terms of critical spin and charge fluctuations within the well-established non-random Hubbard or t-J model framework has yet to be developed.

In this work, we address the microscopic mechanism for the Planckian metal phase in cuprates via different approaches from the previous ones <cit.> based on recently developed heavy-fermion perspective of the two-dimensional slave-boson t-J model, known to offer qualitative understanding of cuprate phase diagram,  in the formulation of the Kondo-Heisenberg lattice model <cit.>. Within our framework, the Planckian metal phase appears as a result of  local disordered bosonic charge fluctuations coupled to a fermionic spinon band and a heavy-electron band near a localized-deloclaized quantum critical point due to the competition between the  pseudogap and Fermi liquid phases.

In this approach, the hoping of electrons is expressed in terms of effective Kondo hybridization between a composite fermionic spinon-holon bound state representing conduction band and a charge-neutral gapless fermionic spinon band with a Fermi surface, while as the Heisenberg exchange coupling is described by resonating-valence-bond (RVB)  spin-liquid with both hoping of fermionic spinons and singlet pairing between them. The mean-field theory of this approach captures qualitatively the pseudogap, Fermi liquid, and d-wave superconducting phases. In particular, it advances the previous slave-boson t-J approach by capturing aspects of coherent quasi-particle excitation observed at nodal Fermi pocket in the pseudogap phase. Meanwhile, the spin-liquid state here is further stablized by coupling to charge Kondo hybridization, similar to the Kondo-stablized spin-liquid mechanism in heavy-fermion systems <cit.>. Here, by perturbative renormalization group analysis, we study the quantum phase transition of the model beyond mean-field and seek the possible emergence of Planckian metal state due to spin and charge fluctuations near criticality. Our approach is highly motivated by the striking similarity in strange metal phenomenology between cuprates <cit.> and heavy-fermion Kondo lattice systems <cit.>, in which the Fermi surface volume reconstructs over the entire strange metal region in both systems. This indicates a Kondo-breakdown-like physics observed in heavy-fermion systems may appear in cuprates where the SM state and Fermi surface reconstruction occur simultaneously near the Kondo breakdown QCP due to coupling of local charge fluctuations to the Fermi surface <cit.>. A related yet distinct heavy-fermion Kondo lattice approach was proposed in Ref. <cit.> to address the evolution from the pseudogap metal with small Fermi surfaces to the conventional Fermi liquid with a large Fermi surface.

Via the controlled RG analysis, a stable quantum critical Planckian strange metal (SM) phase over a finite range in doping with universal $T$-linear scattering rate ($\alpha \sim O(1)$ being a universal constant independent of microscopic couplings) is realized near a localized-delocalized Kondo breakdown transition. Therein, the local bosonic charge (effective Kondo) fluctuations coupled to composite fermionic conduction band and gapless fermionic spinons. The universal quantum critical $\omega /T$-scaling is found in this phase in dynamical scattering rate, in excellent agreement without fine-tuning with the optical conductivity measurement <cit.>, and the universal doping-independent field-to-temperature scaling in magnetoresistance over an extended doping range in Ref. <cit.> in the strange metal region of various overdoped cuprates. The marginal Fermi liquid single-electron spectral function and crossover in Fermi surface volume in this phase are both in good agreement with ARPES <cit.> and Hall measurements <cit.>, respectively. Our study indicates that this exotic phase is a quantum critical phase governed by the critical charge (effective Kondo) fluctuations at the localized-delocalized Kondo breakdown quantum critical point arising from the competition between the pseudogap (Cooper-pair formation) and Fermi liquid (electron hoping) phases. It provides an understanding of the strange metal state within the context of the global phase diagram of cuprate superconductors and reconciles the seemingly inconsistent scenarios of the strange metal state observed in cuprates.

The Feynman diagrams for the self-energy and vertex correction being used in the RG analysis. (a) The graphical representation of the bare propagator of various fields/operators.  Feynman diagrams of the (b)-(f) self-energy and (g) vertex correction for $H_J$.

§ RESULTS

§.§ Heavy-fermion formulation of the slave-boson $t$-$J$ model

We start from the Hamiltonian of  the slave-boson representation of $t$-$J$ model on a 2D lattice <cit.>, describing as $ H=H_{t}+H_{J}$ with $H_{t} =-t\sum_{\langle i,j \rangle , \sigma} c_{i\sigma}^{\dagger} c_{j\sigma}-\mu \sum_{i\sigma} c_{i\sigma}^{\dagger} c_{i\sigma}$ and $	H_J  = J_{H} \sum_{\langle i,j\rangle} \bm{S_{i}} \cdot\bm{S_{j}}.$ Here, $(t,\, \mu, \, J_H)$ denotes the (hopping strength, chemical potential, Heisenberg coupling), and $\langle i,j \rangle$ the nearest-neighboring sites, and the local spin operator $ \bm{S}_{i}=\frac{1}{2} \sum_{\sigma\sigma^\prime} c_{i\sigma}^{\dagger} \bm{\sigma}_{\sigma\sigma^\prime} c_{i\sigma^\prime} $.

Under the slave-boson representation $c_{i\sigma}^\dagger \to  f^\dagger_{i\sigma} b_i$ with $f_{i\sigma} \, (b_i)$ being fermionic charged-neutral spinon (bosonic spinless charged holon), we further factorize the $H_t$ and $H_J$ terms via Hubbard–Stratonovich transformation following Ref. <cit.>, i.e., $H_{t} \to t\sum_{\langle i,j\rangle,\sigma} \left[\left(f_{i\sigma}^{\dagger}b_{j}^{\dagger}+f_{j\sigma}^{\dagger}b_{i}^{\dagger}\right)\xi_{ij,\sigma}+H.c.\right]$ and $ H_{J} \to \sum_{\langle i,j\rangle\sigma}\left(-\chi_{ij}f_{i\sigma}^{\dagger}f_{j\sigma}+ \Delta_{ij}\tilde{\sigma}f_{i\sigma}^{\dagger}f_{j,-\sigma}^{\dagger}  +H.c.\right), ~ \tilde{\sigma} \equiv \text{sgn}(\sigma)$.

Note that we decompose the Heisenberg interaction into both the particle-hole and particle-particle forms. Here, $\xi_{ij,\sigma}$, $\chi_{ij}$, and $\Delta_{ij}$ represent the auxiliary Hubbard-Stratonovich fields living between sites $i$ and $j$, as a result, their degrees of freedom for the bond fields are twice as many as that for $b$ and $f$. The $\chi_{ij}$ and $\Delta_{ij}$ fields represent the spinon hopping and pairing bond fields, respectively, whose condensate play the role of the effective hopping and  pre-formed singlet d-wave Cooper pairing of spinons in the RVB spin liquid. The  model shows a U(1) gauge symmetry: $f_{i\sigma} \to f_{i\sigma} e^{i \theta_i}$, $b_i \to b_i e^{i \theta_i}$, $\chi_{ij} \to \chi_{ij} e^{i (\theta_i - \theta_j)}$ and $\xi_{ij,\sigma} \to \xi_{ij,\sigma} e^{i (\theta_i+\theta_j)}$. A finite and uniform mean-field $\chi=\langle \chi_{ij} \rangle$ is assumed. The $\xi_{ij,\sigma}$ field, previously proposed in Refs. <cit.>,  describes a gauge charge-2$e$ spinon-chargeon bound state; it further advances the original slave-boson approach by capturing the nodal Fermi pocket (arc) quasi-particle excitations in the pesudogap phase of cuprates. The hoping term $H_t$, after the decomposition, behaves as an effective Kondo coupling in the charge sector. If the slave boson $b$ gets Bose condensed, the $\xi$ fermion becomes a physical conduction electron. Terms containing slave boson occupation number $n_b^i$ are neglected in the Hubbard-Stratonovich transformation for $H_t$ and $H_J$. The constraints $b_{i}^{\dagger} b_{i} + \sum_{\sigma} f_{i\sigma}^{\dagger} f_{i\sigma} = 1$  is imposed <cit.>. Four distinct mean-field phases are realized in this model, depending on whether or not the two boson fields $b$ and $\Delta$ get Bose-condensed [see Fig. <ref>(b) and Ref. <cit.>]: the pseudogap phase, known as $Z_2$ fractionalized Fermi liquid or $Z_2$ FL$^\star$, (Landau Fermi liquid, FL) is reached when $\Delta \neq 0, \, \langle b \rangle =0 \, (\Delta =0,  \langle b\rangle \neq 0)$; while the U(1) FL$^*$ (d-wave superconducting) phase arises when $\langle b \rangle = \Delta =0 \, (\langle b \rangle \neq 0, \, \Delta \neq 0)$.  The $Z_2$ FL$^\star$ and U(1) FL$^\star$ phases realized here are examples of the previously proposed fractionalized Fermi liquid states in Ref. <cit.> in the context of heavy-fermion Kondo lattice systems with a small Fermi surface volume and charge-neutral fractionalized spinon excitations which carry a gauge charge characterizing the topological order. The spinons in these two fractionalized spin-liquid phases are deconfined due to the existence of a spinon Fermi surface, and are thus stable against U(1) gauge field fluctuations <cit.>. Meanwhile, these two fractionalized phases are energetically more stable compared to that in the earlier approach of slave-boson t-J model <cit.>  by the charge Kondo hybridization term, similar to the Kondo-stablized spin-liquid mechanism in heavy-fermion systems <cit.>.

We shall explore the phase diagram in U(1) FL$^*$ phase beyond mean-field by including fluctuations of both $t$- and $J$-terms. We will not include U(1) gauge fluctuations in our RG analysis since the spinons are deconfined and stable against them due to the presence of a spinon Fermi surface <cit.>.  The leading  effective action in units of the half-bandwidth of the f-spinon band $D\approx 4\chi \approx 4J_H$ beyond the mean-field level reads ($\hbar = k_B = 1$) <cit.>

\begin{align}
S=&-\sum_{k\sigma}f_{k\sigma}^{\dagger}\mathcal{G}^{-1}_f (k)f_{k\sigma}-\sum_{k}b_{k}^{\dagger}\mathcal{G}^{-1}_b (k)b_{k}  \nn
&+ \sum_{k\sigma a}\xi_{k\sigma}^{a\,\dagger} \mathcal{G}^{-1}_\xi(k) \xi_{k\sigma}^{a} + \sum_{k a}\varphi_{k}^{a\,\dagger}\mathcal{G}^{-1}_\varphi(k) \varphi_{k}^{a} \nn
& +\frac{g}{\sqrt{\beta N_{s}}}\sum_{kp\sigma a}\left(f_{k\sigma}^{\dagger}b_{p}^{\dagger}\xi_{k+p,\sigma}^{a}+H.c.\right)  \nn
& +\frac{J}{\sqrt{\beta N_{s}}} \sum_{kpa}\left[f_{k\uparrow}^{\dagger}f_{p \downarrow}^{\dagger}\varphi_{k+p}^{a}+H.c.\right]
\label{eq:action-S}
\end{align}

with $k = (\bm{k},\omega)$ and $p = (\bm{p},\nu)$, $\varphi$ being the fluctuating field for the d-wave RVB pairing order parameter $\Delta$. In the phases with $\langle b\rangle=0$, the bare hoping parameter $t$ is strongly suppressed by the disordered bosons, leading to an effective hoping $g \equiv 2t \sqrt{\delta}$. Here, $ J = 2J_H$ is the effective exchange, $\mathcal{G}_f (k)= (i\omega-\varepsilon_{\bm{k}})^{-1}$, $\mathcal{G}_\xi (k) = \zeta^{-1} \left( i\omega-\xi_{\bm{k}} \right)^{-1}$, $\mathcal{G}_b = (i\omega - \lambda)^{-1}$, and $\mathcal{G}_\varphi = (J/2)^{-1}$ denote the bare Green's functions (see Supplementary Notes 2). The momenta relevant for the bosonic $\varphi$ field are near the antinode. Our perturbative expansion in bare couplings $g$, $J$ is controlled since $g/D, \, J/D <1$ (for an  estimated $J_H/t \sim 0.3$).

To study the charge dynamics and transport properties, we go beyond static mean-field level of $\xi$ field by generating its dynamics and dispersion via second-order hoping process at fixed $g=g^*_0$ such that $\zeta^{-1} \equiv ( g^2 \rho_0/D)^{-1}$ appears as a prefactor in $\mathcal{G}_\xi$, with $\rho_0 = 1/D$ being the constant density of states at Fermi level for $f$ spinon band.

The $f$-spinon band is approximated by a linear-in-momentum dispersion, namely $\varepsilon_{\bm{k}}=h_{\bm{k}}-\mu_{f} \approx v |\bm{k}|$ with  $h_{\bm{k}}\equiv-2\chi\left(\cos k_{x}+\cos k_{y}\right)$ and $\mu_{f}=\mu-\lambda$ being the effective chemical potential for $f$ spinon ($\mu$ is the chemical potential for the original conduction electron, see Supplementary Notes 1). Here, the spinon band is fixed at half-filling, $\mu_f=0$.  The band structure for the $\xi$ fermion shows a hole-like dispersion, $\xi_{\bm{k}}=-\zeta\varepsilon_{\bm{k}}-\mu_\xi$. The Lagrange multiplier $\lambda >0$ is introduced to enforce the local constraint for slave boson. The slave boson ($b$ field) is effectively treated as a local boson with a flat band of energy $\lambda$ and with a negligible dispersive band (see Supplementary Notes 2). In the pseudogap and U(1) FL$^\star$ phases, the chemical potential of the $\xi$ band $\mu_\xi$ fixes hole doping $\delta$ for the system such that  $N_s^{-1} \sum_{ij,\sigma}  \langle \xi^\dagger_{ij, \sigma}\xi_{ij, \sigma} \rangle = \delta$. Note that Eq. (<ref>) is also applicable when d-wave preformed Cooper pair order, relevant for temperatures slightly above the superconducting transition temperature of cuprates, is replaced by the pair density wave state, considered as a hallmark of pseudogap phase at higher temperatures <cit.>.

(a) Plot of the RG flow for the Heisenberg and hopping terms. Paths I, II, III in (a) indicate the three types of phase diagram that can be accessed from the RG equations. (b): three possible distinct types of phase diagrams (I, II, III) corresponding to paths (I, II, III) from the RG flow diagram in (a). $T_N$, $T^*$, and $T_{FL}$ denote the temperature crossovers for the anti-ferromagnetic (AFM) order, pseudogap, and  FL states, respectively. In (b)-I, $\delta_c$ and $\delta_{c2}$ represent the critical points associated respectively with PSG-SM and SM-FL transitions once superconductivity is suppressed.

§.§ Renormalization group analysis

The perturbative RG analysis is applied to the effective action of the modified slave-boson model Eq. (<ref>) by considering the one-loop diagrams shown in Fig. <ref>. It is convenient to define the dimensionless bare coupling constants $g \to g/D$, $J \to J /D$, and absorb $1/\zeta$ in $\mathcal{G}_\xi$ into $g$ by a rescaling, $g \to \bar{g} = g/\sqrt{\zeta}, \, \xi \to \sqrt{\zeta} \xi $. Our RG analysis is perturbatively controlled as $\bar{g}<1$
(the bare value $g<g_0^*$ is set), and $J<1$.

Our RG approach can be further controlled by the $\epsilon$-expansion technique with a small parameter $\epsilon=d-z$ within the  convergence radius $|\epsilon|\leq 1$. We set the dynamical exponent $z=1$ due to the linearized dispersion of $f$, and spatial dimension $d=2$ here. The RG scaling equations of the running renormalized dimensionless couplings, $\bar{g}(\ell) \rho_0$ and $J(\ell) \rho_0$ read (see Supplementary Notes 3)

\begin{align}
\frac{d(\bar{g}(\ell)\rho_{0})}{d\ell}	 & =-\left(\frac{\epsilon}{2}\right)(\bar{g}(\ell)\rho_{0}) + \left(\bar{g}(\ell)\rho_{0}\right)^{3}, \nn
\frac{d(J(\ell)\rho_{0})}{d\ell}	& =-\left(\frac{\epsilon}{2}\right)(J(\ell)\rho_{0})+\left(J(\ell)\rho_{0}\right)\left(\bar{g}(\ell)\rho_{0}\right)^{2}\nn
&\qquad   +\frac{3}{2}\left(J(\ell)\rho_{0}\right)^{3},
\label{eq:beta-fun-shift}
\end{align}

where scaling parameter $\ell\equiv-\ln D>0$. The RG flow equations of Eqs. (<ref>) is shown in Fig. <ref>(a) where the critical fixed point occurs at $C_r = \left(J^*\rho_{0},\, \bar{g}^*\rho_{0} \right) = (\sqrt{\epsilon/2},\sqrt{\epsilon/3})$. For simplicity, the U(1) gauge-field fluctuations are not included in our RG analysis; nevertheless, its effect on transport and thermodynamic properties are addressed in Discussions (see also Supplementary Notes 5).

RG flow and phase diagram. As shown in Fig. <ref>(a), our one-loop RG flow  diagram supports a quantum critical point (QCP) at $C_r$, separating the four phases: the pseudogap PSG (FL) phase is reached when the Heisenberg coupling $J$ becomes relevant (irrelevant), while the hopping (effective Kondo) term $\bar{g}$ becomes irrelevant (relevant); the d-wave superconducting state (U(1) FL$^*$, SM) phase arises for both $\bar{g}$ and $J$ being relevant (irrelevant), respectively. Near $C_r$, the correlation length $\eta (\ell)$ shows a power-law divergence as:  $\eta (\ell)\sim \ell^{-\nu}$ with the exponent $\nu = 1$, indicating linear-in-doping crossover scales $T^*, \, T_{FL} \propto \eta^{-z}  \propto |\delta-\delta_c|$, corresponding to the PSG$\to$SM and SM $\to$ FL crossovers, respectively <cit.>. We shall see below that the systems with initial couplings flowing to the mean-field U(1) FL$^*$ fixed point show strange metal features with universal quantum critical Planckian scaling in electron scattering rate; therefore beyond mean-field under RG, it corresponds to the strange metal fixed point.

$\omega/T$ scaling of scattering rate and the effective mass enhancement. Theoretical fitting (black curve) for both (a) scattering rate $\hbar/\tau$ and (b) the effective mass enhancement  from our theory exhibits a universal scaling behavior $m^*/m - 1=g(x)$ as a function of $x\equiv \hbar\omega/k_B T$. Theoretical result of $m^*/m - 1$ (red dashed line) obtained from Eq. (S180) in Supplementary Notes 6B is added in (b) for comparison. The experimental data shown here is reproduced from Ref. <cit.>.

§.§ Physical observables near the QCP

§.§.§ Conductivity and scattering rate

The total conductivity $\sigma_{tot} = \rho^{-1}_{tot}$ can be computed via the rule by Ioffe and Larkin <cit.>:  $\sigma_{tot} = \sigma_\xi \sigma_f/ (4\sigma_\xi + \sigma_f)$. Because the slave boson is localized and has negligible dispersion, it does not contribute to the conductivity.

Since the effective mass of the $\xi$ band is much heavier than that for the $f$ spinons, $m_\xi/m_f \sim 1/\zeta \gg 1$,  resistivity is dominated by the $\xi$-band, leading to  $\rho_{tot} = 1/\sigma_{tot} \approx 1/\sigma_\xi= \rho_\xi$.

In the theory of local (momentum independent) self-energy, the vertex corrections are negligible, and the optical conductivity can be obtained directly from the self-energy, $\sigma(\omega) = \left( i\Phi(0)/\omega\right)\int^\infty_{-\infty} d\varepsilon \left[n_F (\varepsilon)-\right. \left. n_F (\varepsilon+\hbar \omega)\right] / [ \hbar \omega + \Sigma_c^* (\varepsilon)-\Sigma_c(\varepsilon+\hbar \omega)]$, where $\Phi (\varepsilon)$ is the transport function (see Ref. <cit.> for details). Therefore, electron transport time is equivalent to the relaxation time since $\sigma_\text{tot} \sim  \sigma_\xi$. The total scattering rate is hence dominated by the contribution from the $\xi$ field: $1/\tau_\text{tot} \sim  1/\tau_\xi$. In the SM state with $\langle b \rangle=\Delta=0$, the gauge-invariant (physical) electron operator can be constructed from $\xi_{ij}$ and $\varphi_{ij}$ as $c= \xi \sqrt{2} \varphi_{ij}^* \sqrt{J}$. The dynamical scattering (relaxation) rate at $T = 0$ is calculated via the imaginary part of the electron self-energy: $\tau(\omega) =-\hbar /2\Sigma_c^{\prime\prime}(\omega)$ and $\Sigma_c^{\prime\prime}(\omega)=(1/2) \Sigma^{\prime\prime}_{\xi}(\omega)$,  where $\Sigma^{\prime\prime}_{\xi}(\omega)$ at $T=0$ is obtained via  second-order RG renormalized perturbation close to the QCP with bare couplings being replaced by the renormalized ones i.e. $\bar{g}, \, \bar{g}^*_0\to \bar{g}(\ell)$, showing a universal local MFL self-energy insensitive to couplings <cit.>, including a constant and a linear-in-$|\omega|$ term, i.e.  $\Sigma_{\xi}^{\prime\prime}(\omega)= \alpha-\varsigma\left|\omega\right|$ with $\alpha \approx \frac{5}{2}\bar{g}^{2}\rho_{0} = \frac{5}{2} D$ and $ \varsigma=\frac{2}{\pi}$.

This local self-energy $\Sigma_\xi$ comes as a consequence of the local nature of the slave boson in our approach.

The constant $\alpha$ is generated by the fluctuating non-condensed local slave-bosons through self-energy of the $\xi$ field, reminiscent of electrons in metals being scattered by disordered random impurities.

§.§.§ Quantum critical phase with universal Planckian scattering rate

Close to the QCP in the scaling regime where conformal symmetry is present, the scattering rate at finite temperatures can be derived from the scattering rate $T = 0$,  $1/\tau(\omega,T = 0)$, by a conformal transformation <cit.>.  Strikingly, following path III of Fig. <ref>(a) over the  region where both $\bar{g}(\ell)$ and $J(\ell)$ are irrelevant, we discover a quantum critical Planckian strange metal “phase" with universal quantum-critical   $\omega/T$-scaling in the scattering rate that is independent of coupling constant (see Supplementary Notes 5):

\begin{align}
\frac{1}{\tau(x)}-\frac{1}{\tau_0} = \frac{4}{\pi} k_B T  f(x/2), ~ f(x/2)=\frac{x}{2}\coth \left(\frac{x}{4}\right)
\label{eq:scatt-rate-scaling}
\end{align}

with $x\equiv \hbar \omega/k_B T$ and $\tau_0 \equiv \tau (\omega=0,T=0)$. Surprisingly, as shown in Fig. <ref>(a), the universal Planckian scattering rate of Eq. (<ref>) is in excellent agreement with the recent optical conductivity measurement in Ref. <cit.> without fine-tuning. In the high-frequency, low-temperature limit $x \gg 1$, the scattering rate divided by $k_B T$ shows a universal scaling behavior, $\frac{\tau^{-1}(T}{k_B T} \approx (2/\pi)  x$. Conversely, in the DC-limit ($x \to 0$), the scattering rate manifests the Planckian scattering rate, revealing a universal feature that is insensitive to microscopic coupling constants: $1/\tau \approx \alpha_P k_B T/\hbar$ with $\alpha_P \approx 8/\pi$, in good agreement with DC-scattering rate estimated in various overdoped cuprates <cit.>. We find this universal feature of scattering rate that is insensitive to coupling constants in our theory is originated from the cancellation of  the same RG renormalized running coupling constant $g^2 (\ell)$ in $\Sigma_\xi$ and in the denominator of $\mathcal{G}_\xi$ (arising from the second-order hoping process that generates dynamics and energy dispersion of the $\xi$-band). Due to this cancellation, all initial (bare) values of $\bar{g}$ and $J$ which flow to the U(1) FL$^*$ fixed point at $\bar{g}=J=0$ form an extended “phase" by showing the same Planckian strange-metal behavior [see Eq. (<ref>)].  Our finding is in excellent agreement with  the doping-independent universal $B/T$ scaling in magnetoresistance observed over the extended strange metal phase of Tl2201  <cit.> (see Supplementary Notes 6).

Meanwhile, the electron mass renormalization  $m^*_e - m^*_e(0 )$ is well captured by  a universal scaling form: $m^*/m -1 = g(x)$ with the scaling function g(x) obtaining from the real part of the dynamical self-energy $\Sigma_{\xi}^\prime$ via Kronig-Kramers relation [see Fig. <ref>(b) and Supplementary Notes 7].

§.§.§ Single-particle spectral function

Our theory offers a mechanism for the phenomenological MFL ansatz <cit.>: The complete electron self-energy $\Sigma_\xi$, is reminiscent of the MFL behavior with a distinction that $\Sigma_\xi$ here is insensitive to coupling constant. The resulting single-particle spectral function $A(\omega, \bm{k}, T) \equiv -\pi^{-1}G^R_c (\omega, \bm{k}, T) = -\pi^{-1} \text{Re} \left[ \omega -\varepsilon_{\bm{k}}-\Sigma_c\right]^{-1} $ thus shows an excellent fit to the recent ARPES measurement for overdoped cuprates in the SM region <cit.> (see Fig. <ref> and Supplementary Notes 7).

§.§.§ Strange metal state: Signature of quantum critical point vs. quantum critical phase.

Our RG results offer further insights to the issue on the origin of the Planckian strange metal state.  The RG flow in Fig. <ref>(a) in general may lead to three distinct temperature-doping ($T\,,\delta$) phase diagrams by following the three paths I, II, and III depending on initial values of $t$ and $J$ with increasing doping, see Figs. <ref>(b)-I, <ref>(b)-II, <ref>(b)-III. Experimentally, the three paths can be tuned by magnetic fields: When the external magnetic field is weak, the couplings lie on path I in which a finite range of couplings on this path flow to the superconducting phase. This corresponds to Fig. <ref>(b)-I. At a critical magnetic field, couplings follow path II and pass through the QCP at $C_r$; consequently, the strange metal state at finite $T$ is controlled by a single QCP. The most interesting case, however, is path III, occurring at a larger field which completely suppresses superconductivity. For path III, there is  a finite range of doping where all initial values of $t$ and $J$  [blue solid line in path III of Fig. <ref>(a)] flow to the SM fixed point. In this case,  the Planckian strange metal behaviors at finite T persist over a finite doping range as $T \to 0$, corresponding to the quantum critical strange metal “phase". Experimental signatures of both quantum critical point <cit.> and quantum critical phase <cit.> were reported in different cuprates. These seemingly incoherent results can be coherently unified in our generic RG phase diagram: the former may follow a path close to path II, while the latter follows path III.

Spectral intensity for the SM phase: (a) Spectral weight as a function of frequency $\omega$ at the anti-nodal point, and (b) spectral weight as a function of momentum $k_y$ at a zero frequency. Experimental data OD86 (red) and OD81 (blue) are reproduced from Ref. <cit.>.  Greens curves in (a) and (b) corresponds to our theoretical fitting to OD81, with temperature being fixed at $T = 250$K. (c) Density plot for the spectral weight as a function of $\omega$ and $k_y$ (see supplementary Notes 7 for details.)

§.§.§ Thermodynamic property: divergence of specific heat coefficient at the QCP

We further computed the specific heat coefficient in the SM phase near the QCP showing logarithmic-in-$T$ MFL behavior: $\gamma \equiv C_V/T \sim -[A+ B(J)] \ln T$ with $A$ being a constant and $B(J)$ the power-law divergent prefactor arising from the antiferromagnetically short-ranged spin fluctuations (the $J$ term) being quantum critical: $B(J) \sim |J-J_c|^{-p}$ with $p = \nu/2 = 1/2$, in good agreement with the power-law divergent $\gamma$ coefficient $\gamma(T\to0) \sim |\delta-\delta_c|^{-m}$ with $m \sim 0.5$ observed in Eu-LSCO and Nd-LSCO near pseudogap critical point <cit.>. We find the dominating logarithmic-in-$T$ term $-B(J) \ln T$ is originated from the van Hove singularity of the half-filled $f$-spinon band on a 2D square lattice [see the spinon Fermi surface in Fig. <ref>(b) and Supplementary Notes 4], while the sub-leading constant term $-A \ln T$ is contributed from the self-energy $\Sigma_{c,f}$.

§.§.§ Fermi surface reconstruction across the SM region

Our heavy-fermion two-band formulation of the slave-boson t-J model  offers a natural explanation on the enlargement of Fermi surface (FS) volume from a smaller FS with hole density $\delta$ in pseudogap phase to a larger one with hole density $1+\delta$ in the FL phase $(\langle b \rangle =\sqrt{\delta}, ~ \xi_\sigma \sim  \sqrt{\delta} f_\sigma, ~ \sum_\sigma f^\dagger_{i\sigma} f_{i\sigma} = 1-\delta)$ across the strange metal region [Fig. <ref>(b)] <cit.>.

In the pseudogap and SM phases with non-condensed bosons, the  $\xi$-band shows a small hole Fermi surface (hole doping $\delta$). When bosons get condensed across the FL phase boundary, the Fermi surface volume undergoes a sudden jump from $\delta$ to $1+\delta$, in parallel to that occurred at Kondo-breakdown transition in heavy-fermion systems  <cit.>. This is due to the participation of the half-filled spinon band to the Fermi surface by Kondo hybridization between the $f$- and the $\xi$-bands; therefore the hole Fermi surface volume of a normal conduction electron band with a hole density $1+\delta$ (corresponding to a doped hole density $\delta$ relative to half-filling) is recovered, in accordance with the Luttinger theorem [see Fig. <ref>(b)]. This scenario offers a qualitative understanding of the rapid increase observed in the normal-state Hall number near the critical doping of YBCO and Nd-LSCO <cit.>. Nevertheless, a smooth crossover of the Fermi surface volume from $\delta$ to $1+\delta$ is possible in the SM phase via a mixture of two fluids: (i) normal metal with $T$-quadratic resistivity by partial boson condensation and a large Fermi surface volume ($1+\delta$), and (ii) strange metal with $T$-linear resistivity without boson condensation and with a small Fermi surface volume ($\delta$), as was observed in Ref. <cit.>.

§ DISCUSSIONS

The umklapp scattering is usually expected to be important for electrical resistivity by relaxing electron momentum (charge current). Here, we propose an alternative relaxation mechanism without umklapp scattering via local charge Kondo fluctuations where electrical charge current decays into charge-neutral spinon current and local charge fluctuations despite forward scattering <cit.>.

Due to the heavy $\xi$-band, we find that the sub-leading corrections to the scattering rate are negligible compared to the leading Planckian scattering rate mentioned above, within the temperature range of our interest: $1/\tau_{\xi,\rm{gauge}} \sim (T/\chi_d)^{4/3} \varepsilon_{F,\xi} \propto a  (T/D)^{4/3}$  with $a/D \sim O(\delta^2) \ll 1$ and  specific heat coefficient $\gamma(T) \sim (b/D) \cdot (T/D)^{-2/3}$ with $b \sim O(\delta^2) \ll 1$ contributed from U(1) gauge field <cit.>, see Supplementary Notes 5. Meanwhile, due to the local nature of the slave boson in our approach, the T-linear scattering rate contributed from U(1) gauge field found in Ref. <cit.> is absent here. The Planckain metal phase we find here can be further stablized in the generalized version of our model in large-N and multi-channel limit where fluctuations from gauge-field and higher order processes are suppressed.

The Planckian phase in transport and thermodynamical coefficients within our theoretical framework show distinct qualitative features: the former is insensitive to but the latter strongly depend on the microscopic couplings. This suggests possible breakup of low-energy excitations in this phase into charged particles (electrons) and the charge-neutral spinons; the former dominates in transport while the latter in thermodynamics, reminiscent of the Kondo breakdown QCP in heavy-fermion system <cit.>.

To reach the Planckian bound, it was argued in Ref. <cit.> that a strong inelastic scattering with local energy relaxation $\Delta E \sim k_B T$ is required so that the equilibration time $\tau_{eq} \sim \hbar/\Delta E \sim \hbar/k_B T$ reaches the Planckian time scale. The charge fluctuating effective Kondo term in our model offers a realization of such strong inelastic scattering channel with a local energy relaxation rate $1/\tau \sim k_B T/\hbar$, where electron loses its energy to local slave boson by coupling to a spinon band.

It was pointed out that the strong disorder coupled to the electron interactions is crucial to realize the Planckian strange metal state, as shown in the SKY model <cit.>. Within our two-band heavy-fermion slave-boson $t$-$J$ model, local fluctuating non-condensed slave boson plays the role of disorder embedded in electron interaction in the fluctuating Kondo hybridization. MFL and SM properties are thus  expected, reminiscent of the SKY models. Via the effective Kondo term, the fluctuating slave-boson also generates impurity-like scattering, leading to a residual scattering rate at $T = 0$ due to the breaking of the translational symmetry of the original Hamiltonian.

A formal link of our model to the SYK-like models deserves further investigations.

§ CONCLUSIONS

We provide a microscopic mechanism for the Planckian metal phase in the Kondo-Heisenberg formulation of the slave-boson t-J model based on the dynamical charge Kondo breakdown near the localized-to-delocalized phase transition in the form of critical charge (Kondo) fluctuations. The slave-boson $t$-$J$ model was first mapped onto an effective Kondo-Heisenberg model. Via perturbative RG analysis on this effective model,  we realized an extended quantum critical Planckian metal phase. Within the one-loop RG, we identified a quantum critical point separating pseudogap, $d$-wave superconducting, a normal Fermi liquid and a strange metal phases.  In particular, we find a stable strange metal phase in our RG phase diagram close to this critical point exhibiting $T$-linear scattering rate and $\omega/T$-quantum critical scaling and the Planckian behavior where $\alpha_P \sim O(1)$ that is independent of coupling constant, in excellent agreement with the optical conductivity and magnetoresistance measurements in various overdoped cuprates. Our theoretical predictions on the specific heat coefficient, effective mass enhancement, and single particle spectral function in the strange metal state agree well with the experimental observations and hence offers a microscopic mechanism for the marginal Fermi liquid ansatz. Our theoretical framework offer a promising route to reveal the mystery of quantum critical strange metal phase and its link to the generic phase diagram of high-$T_c$ cuprate superconductors.

§ METHODS

The slave-boson approach to the $t$-$J$ model. The model we consider here is the $t$-$J$ model $ H	=H_{t}+H_{J}$,

\begin{align}
H_{t}& =-t\sum_{\langle i,j\rangle,\sigma}c_{i\sigma}^{\dagger}c_{j\sigma}-\mu\sum_{i\sigma}c_{i\sigma}^{\dagger}c_{i\sigma}, \nn	H_J & = J_{H}\sum_{\langle i,j\rangle}\bm{S}_{i}\cdot\bm{S}_{j}.
\end{align}

$\bm{S}_{i}=\frac{1}{2}c_{i\sigma}^{\dagger}\bm{\sigma}_{\sigma\sigma^\prime}c_{i\sigma^\prime}$ represents the (dimensionless) spin operator at site $i$, $t$ is the hopping strength, and $\mu$ denotes the chemical potential of the $c$-band. Here, $\langle i,j\rangle$ means only the nearest-neighboring sites are considered. The Heisenberg coupling $J_{H}\sim t^{2}/U>0$. For low dopings, the doubly-occupied state is not energetically favorable. To project out the doubly-occupied state from the Hilbert space, we exploit the slave-boson technique, which splits a conduction electron as a product of a charge-neutral spinon $f_{i\sigma}$ and a spinless charged boson $b_i$, i.e., $c_{i\sigma}^{\dagger}\rightarrow f_{i\sigma}^{\dagger}b_{i}$ subject to the single-occupancy constraint for the $b_i$ and $f_{i\sigma}$ operators on each site <cit.>.

Near half-filling, the doubly-occupied state is not energetically favorable, thus will be ruled out from the Hilbert space. Therefore, a no double-occupancy constraint has to be enforced on every site, i.e. $\sum_{\sigma}c_{i\sigma}^{\dagger}c_{i\sigma}\leq 1$. In terms of the fermion operators, the inner product of two spin operators takes the form of

\begin{align}
\bm{S}_{i}\cdot\bm{S}_{j}=&-\frac{1}{4}\sum_{\sigma\sigma^\prime}\left(\tilde{\sigma}c_{i\sigma}^{\dagger}c_{j,-\sigma}^{\dagger}\right)\left(\tilde{\sigma^\prime}c_{j,-\sigma^\prime}c_{i\sigma^\prime}\right) \nn
&+\frac{1}{4}\sum_{\sigma \sigma^\prime}c_{i\sigma}^{\dagger}c_{j\sigma^\prime}^{\dagger}c_{j\sigma}c_{i\sigma^\prime},
\end{align}

where $\tilde{\sigma}\equiv \text{sgn}(\sigma)$.

Near half-filling, the doubly-occupied state is not energetically favorable, thus will be ruled out. In this circumstance, we may apply the slave-boson $c_{i\sigma}^{\dagger}\rightarrow f_{i\sigma}^{\dagger}b_{i}$ to project out the doubly occupied state from the Hilbert space, leading to the following single-occupancy constraint for the b and f operators on each site: $ b_{i}^{\dagger}b_{i}+\sum_{\sigma}f_{i\sigma}^{\dagger}f_{i\sigma}=1.$  This constrain can be incorporated into the whole Hamiltonian via applying the Lagrange multiplier,

\begin{align}
\sum_{i}\lambda_{i}\left(b_{i}^{\dagger}b_{i}+\sum_{\sigma}f_{i\sigma}^{\dagger}f_{i\sigma}-1\right).
\end{align}

In the slave-boson representation, the kinetic term $H_{t}$ and the chemical potential term $H_{\mu}$ can be expressed as ($t>0$ being the strength of hopping)

\begin{align}
H_{t}+H_{\mu}=-t\sum_{\langle i,j\rangle,\sigma}\left(f_{i\sigma}^{\dagger}b_{i}\right)b_{j}^{\dagger}f_{j\sigma}-\mu\sum_{i,\sigma}f_{i\sigma}^{\dagger}f_{i\sigma}.
\end{align}

The Heisenberg term $H_J$ can be expressed in the forms of particle-particle (p-p) and particle-hole (p-h) channels,

\begin{align}
& \text{p-p} \to \left(-\frac{J_{H}}{2}\right)\sum_{\langle i,j\rangle}\sum_{\alpha\beta\gamma\delta}\left(\varepsilon_{\alpha\beta}f_{i\alpha}^{\dagger}f_{j\beta}^{\dagger}\right)\left(\varepsilon^{\gamma\delta}f_{j\delta}f_{i\gamma}\right)+const., \nn
& \text{p-h} \to \left(-\frac{J_{H}}{2}\right)\sum_{\langle i,j\rangle}\sum_{\sigma,\sigma^\prime}\left(f_{i\sigma}^{\dagger}f_{j\sigma}\right)\left(f_{j\sigma^\prime}^{\dagger}f_{i\sigma^\prime}\right),
\end{align}

where $\varepsilon=i\sigma_{y}$ is an anti-symmetric tensor. In our mean-field analysis, we consider both of them.

§ ACKNOWLEDGMENTS

We acknowledge discussions with Cenke Xu, G. Kotliar, A. Georges, J. Schmalian, E. Berg, T. Devereaux, S. A. Kivelson, P. A. Lee, T. K. Lee, P. Coleman, N. Hussey, R. Greene, B. Michon, D. Van Der Marel. This work is supported by the National Science and Technology Council (NSTC) Grants 110-2112-M-A49-018-MY3, the National Center for Theoretical Sciences of Taiwan, Republic of China (C.- H.C.). Y.-Y.C. acknowledges the financial support from The 2023 Postdoctoral Scholar Program of Academia Sinica, Taiwan. C.-H. C. acknowledges the hospitality of Aspen Center for Physics, USA and Kavli Institute for Theoretical Physics, UCSB, USA where part of the work was done.

§ HUBBARD-STRATONOVICH (HS) TRANSFORMATION, MEAN-FIELD THEORY AND THE EFFECTIVE ACTION FOR THE U(1) FL$^*$ PHASE

The Hamiltonian of the Kondo-Heisenberg slave-boson approach to the $t$-$J$ model reads

\begin{align}
H=&-t\sum_{\langle i,j\rangle\sigma}\left(f_{i\sigma}^{\dagger}b_{i}\right)b_{j}^{\dagger}f_{j\sigma}-\mu\sum_{i\sigma}f_{i\sigma}^{\dagger}f_{i\sigma}+\sum_{i}\lambda_{i}\left(b_{i}^{\dagger}b_{i}+\sum_{\sigma}f_{i\sigma}^{\dagger}f_{i\sigma}-1\right) \nn
& -J_{H}\sum_{\langle i,j\rangle}\sum_{\sigma\sigma^\prime}\left(\tilde{\sigma}f_{i\sigma}^{\dagger}f_{j,-\sigma}^{\dagger}\right)\left(\tilde{\sigma^\prime}f_{j,-\sigma^\prime}f_{i\sigma^\prime}\right)-J_{H}\sum_{\langle i,j\rangle}\sum_{\sigma\sigma^\prime}\left(f_{i\sigma}^{\dagger}f_{j\sigma}\right)\left(f_{j\sigma^\prime}^{\dagger}f_{i\sigma^\prime}\right) ,
\label{eq:H-SB}
\end{align}

where $\sigma,\sigma^\prime$ refer to the spin-1/2 indices. Following Ref. <cit.>, we perform the HS transformation on the $H_{t}$ term. To do so, we first define the following operator

\begin{align}
\xi_{ij,\sigma}^{\dagger}=\left(f_{i\sigma}^{\dagger}b_{j}^{\dagger}+f_{j\sigma}^{\dagger}b_{i}^{\dagger}\right)/\sqrt{2},
\label{eq:def-HS-xi_ij}
\end{align}

and re-express the hopping term $H_{t}$ in terms of $\xi_{ij,\sigma}$. The $\xi_{ij,\sigma}$ field is defined on the middle of the bond connecting sites $i$ and $j$ such that

\begin{align}
\xi_{ij,\sigma}=\xi_{\sigma}\left(\frac{\bm{r}_{i}+\bm{r}_{j}}{2}\right)=\xi_{ji,\sigma}\np
\end{align}

On a two-dimensional square lattice, its Fourier component can be expressed as

\begin{align}
\xi_{\bm{k}\sigma}^{a}=\frac{1}{\sqrt{N_{s}}}\sum_{i}\xi_{i,i+\hat{e}_{a},\sigma}e^{-i\bm{k}\cdot(\bm{r}_{i}+\hat{e}_{a}/2)},
\end{align}

where $\hat{e}_{a}=\bm{r}_{j}-\bm{r}_{i}$ with $a = x,y$ being the label of the nearest neighbor sites for a 2D square lattice and $N_s$ the number of sites. It is straightforward to show that the $H_{t}$ term can be expressed as a quadratic form in $\xi_{ij,\sigma}$:

\begin{align}
\sum_{\langle i,j\rangle,\sigma}\xi_{ij,\sigma}^{\dagger}\xi_{ij,\sigma}	=\sum_{\langle i,j\rangle,\sigma}n_{f}^{i}n_{b}^{j}-H_{t}/t \Rightarrow  H_{t}=-t\sum_{\langle i,j\rangle,\sigma}\xi_{ij,\sigma}^{\dagger}\xi_{ij,\sigma}.
\end{align}

The $n_{f}^{i}n_{b}^{j}\sim(1-n_{b}^{i})n_{b}^{j}$ term shown above represents the holon-holon density interaction. In our theory, as we consider the low-doping regime, this interaction can be neglected. Upon performing the HS transformation, the hoping term takes the form of

\begin{align}
H_{t}	\to  -\frac{t}{\sqrt{2}}\sum_{\langle i,j\rangle,\sigma}\left[\left(f_{i\sigma}^{\dagger}b_{j}^{\dagger}+f_{j\sigma}^{\dagger}b_{i}^{\dagger}\right)\xi_{ij,\sigma}+H.c.\right]
+t\sum_{\langle i,j\rangle\sigma}\xi_{ij,\sigma}^{\dagger}\xi_{ij,\sigma},
\label{eq:Ht-HS-realspace}
\end{align}

which is analogous to the Kondo interaction of a Kondo lattice model. In the momentum space, the hoping term in the form of Eq. (<ref>) becomes

\begin{align}
H_{t} \to \frac{2t}{\sqrt{N_{s}}}\sum_{\bm{k}\bm{k}^{\prime}\sigma a}\left(f_{\bm{k}\sigma}^{\dagger}b_{\bm{k}^{\prime}}^{\dagger} y_{\bm{k}-\bm{k}^\prime}^{a}\xi_{\bm{k}^{\prime}+\bm{k},\sigma}^{a}+H.c.\right)+t\sum_{\bm{k}\sigma a}\xi_{\bm{k}\sigma}^{a\,\dagger}\xi_{\bm{k}\sigma}^{a}
\end{align}

\begin{align}
y_{\bm{k}-\bm{p}}^{a}&\equiv e^{i(\bm{k}-\bm{p})\cdot\hat{e}_{a}/2}.
\end{align}

To decompose the Heisenberg term ($H_J$), we introduce two HS fields, $\chi_{ij}$ and $\Delta_{ij}$, corresponding to the spinon bond field (particle-hole channel) and  spinon pairing field (particle-particle), respectively,

\begin{align}
H_{J}	 \to\left(-J_{H}\right)\sum_{\langle i,j\rangle}\left(\sum_{\sigma}\chi_{ij}f_{i\sigma}^{\dagger}f_{j\sigma}+H.c.-\left|\chi_{ij}\right|^{2}\right) +J_{H}\sum_{\langle i,j\rangle}\sum_{\sigma,\sigma^\prime}\left(\Delta_{ij}\varepsilon_{\sigma\sigma^\prime}f_{i\sigma}^{\dagger}f_{j\sigma^\prime}^{\dagger}+H.c.\right)  +J_{H}\sum_{\langle i,j\rangle}\left|\Delta_{ij}\right|^{2}.
\end{align}

In the following analysis, we assume a finite mean-field value of spinon bond field $\left\langle \chi_{ij}\right\rangle =\chi \neq 0$, which is uniformly distributed on the whole square lattice. The spinon pairing field $\Delta_{ij}$ can be introduced via $\Delta_{i,i+\hat{e}_{a}}\to\eta_{a}\Delta+J_{H}\varphi_{i,i+\hat{e}_{a}}$ where $\Delta \equiv \langle \Delta_{ij}\rangle$ denotes the mean-field value of the spinon pairing bond field and $\varphi_{i,i+\hat{e}_{a}}$ represents the fluctuation field of $\Delta$. Both the spinon bond field and the spinon  pairing field are defined on the middle of the links, similar to the $\xi$ field. Here, $\eta_{a}$ is a bond-dependent index which characterizes the underlying symmetry of lattice. We choose the $d$-wave pairing in our theory, corresponding to $\eta_{x}=1=-\eta_{y}$. The full mean-field Hamiltonian of $H_{J}$ reads

\begin{align}
H_{J}^{MF}	=\sum_{\bm{k},\sigma}h_{\bm{k}}f_{\bm{k}\sigma}^{\dagger}f_{\bm{k}\sigma}+\sum_{\bm{k},\sigma}\left(\Delta\tilde{\sigma}f_{\bm{k}\sigma}^{\dagger}f_{-\bm{k},-\sigma}^{\dagger}\sum_{a}\eta_{a}\cos k_{a}+H.c.\right) + \frac{2N_{s}\left|\chi\right|^{2}}{J_{H}}+\frac{2N_{s}\left|\Delta\right|^{2}}{J_{H}},
\end{align}

where $h_{\bm{k}}=-2J_{H}\chi\left(\cos k_{x}+\cos k_{y}\right)$ denotes the spinon dispersion.  The fluctuating part of $H_{J}$, denoted as $H_{J}^{fluc}$ is

\begin{align}
H_{J}^{fluc}\to J_{H}\sum_{\langle i,j\rangle,\alpha}\left(\varphi_{i,i+\hat{e}_{a}}\tilde{\alpha}f_{i\alpha}^{\dagger}f_{i+\hat{e}_{a},-\alpha}^{\dagger}+H.c.\right)+J_{H}\sum_{i,a}\left|\varphi_{i,i+\hat{e}_{a}}\right|^{2}.
\end{align}

The Fourier transform of the fluctuating  spinon bond field $\varphi_{\bm{k}}^{a}$ is defined as

\begin{align}
\varphi_{\bm{k}}^{a}=\frac{1}{\sqrt{N_{s}}}\sum_{i}\varphi_{i,i+\hat{e}_{a}}e^{-i\bm{k}\cdot(\bm{r}_{i}+\hat{e}_{a}/2)}.
\end{align}

The fluctuating part of $H_{J}$  can thus be expressed as

\begin{align}
H_{J}^{fluc}\to\frac{2J_{H}}{\sqrt{N_{s}}}\sum_{\bm{k}\bm{p}}\left[f_{\bm{k}\uparrow}^{\dagger}f_{\bm{p}\downarrow}^{\dagger}\sum_{a}\cos\left(\frac{k_{a}-p_{a}}{2}\right)\varphi_{\bm{k}+\bm{p}}^{a}+H.c.\right]+J_{H}\sum_{\bm{k}a} \left(\varphi_{\bm{k}}^{a}\right)^\dagger \varphi_{\bm{k}}^{a}.
\end{align}

The $H_{t}$ and $H_{J}$ shown above both contain a nonlinear (momentum-dependent) form factors, $y^{a}_{\bm{k-p}}$ and $\cos \left(\frac{k_a-p_a}{2}\right)$, respectively. To simplify our calculations, we first perform Taylor expansion of these momentum-dependent form factors as a polynomial in momentum, i.e. $f(k) = 1+ f^\prime(k=0) k^2 + (2!)^{-1}f^{\pp}(k=0)k^2 + O(k^3)$. By power-counting the scaling dimensions of these Taylor-expanded coupling constants, $t$ and J, we keep only the zeroth-order (constant) terms of them and neglect all the higher-order momentum terms. This approximation is justified since the higher-order momentum-dependent couplings are much more irrelevant than the constant terms under renormalization group (RG) transformation. We therefore use these approximate forms of $H_t$ and $H_J$ in our following analysis. They read

\begin{align}
H_{t}	\to \frac{g}{\sqrt{N_{s}}}\sum_{\bm{k}\bm{p}\sigma a}\left(f_{\bm{k}\sigma}^{\dagger}b_{\bm{p}}^{\dagger}\xi_{\bm{k}+\bm{p},\sigma}^{a}+H.c.\right)+\frac{g}{2}\sum_{\bm{k}\sigma a} \left(\xi_{\bm{k}\sigma}^{a}\right)^\dagger \xi_{\bm{k}\sigma}^{a},\nn
H_{J}^{fluc}	\to \frac{J}{\sqrt{N_{s}}}\sum_{\bm{k}\bm{p}a}\left(\varphi_{\bm{k}+\bm{p}}^{a}f_{\bm{k}\uparrow}^{\dagger}f_{\bm{p}\downarrow}^{\dagger}+H.c.\right)+\frac{J}{2}\sum_{\bm{k}a}\left(\varphi_{\bm{k}}^{a}\right)^\dagger \varphi_{\bm{k}}^{a},
\end{align}

with the newly defined couplings

\begin{align}
g\equiv2t,\quad J\equiv2J_{H}.
\end{align}

In the following, we will use the simplified $H_{J}$ and $H_{t}$ terms for further evaluations.

The electron hoping $t$, in general, depends on doping. To extract the doping dependence of hoping, we may rescale the whole Hamiltonian with the half bandwidth of the $f$-spinon, i.e.  i.e., $H \to H/D$ with $D=4 J_H  \chi \sim 4J_H$. In the absence of boson condensate,  $\langle b \rangle =0$, we expect the hoping parameter $t$ in $H_t$ is significantly suppressed. Below, we estimate the amount of the suppressed hoping term  upon perturbatively taking the expectation value of $H_t$, i.e. $\langle H_t \rangle$. Due to $\langle b \rangle =0$, the expectation value of $H_t$ in the unperturbed limit vanishes, i.e.  $\langle H_t \rangle = 0$; however, the second-order  contribution $\langle H_t^2 \rangle_0$ survives:

\begin{align}
\langle H_{t}^2\rangle_{0}&=\frac{2t^{2}}{D^{2}}\sum_{\langle i,j\rangle,\sigma}\sum_{\langle k,l\rangle,\sigma^\prime}\langle\xi_{ij,\sigma}^{\dagger}\left(b_{j}f_{i\sigma}+b_{i}f_{j\sigma}\right)\left(f_{k\sigma^\prime}^{\dagger}b_{l}^{\dagger}+f_{l\sigma^\prime}^{\dagger}b_{k}^{\dagger}\right)\xi_{kl,\sigma^\prime}\rangle_{0} \nn
= & \frac{2t^{2}}{D^{2}}\sum_{\langle i,j\rangle,\sigma}\sum_{\langle k,l\rangle,\sigma^\prime}\left\langle \xi_{ij,\sigma}^{\dagger}\xi_{kl,\sigma^\prime}\right\rangle \left[\begin{array}{c}
\langle b_{j}b_{l}^{\dagger}\rangle_{0}\langle f_{i\sigma}f_{k\sigma^\prime}^{\dagger}\rangle_{0}+\langle b_{j}b_{k}^{\dagger}\rangle_{0}\langle f_{i\sigma}f_{l\sigma^\prime}^{\dagger}\rangle_{0}\\
+\langle b_{i}b_{l}^{\dagger}\rangle_{0}\langle f_{j\sigma}f_{k\sigma^\prime}^{\dagger}\rangle_{0}+\langle b_{i}b_{k}^{\dagger}\rangle_{0}\langle f_{j\sigma}f_{l\sigma^\prime}^{\dagger}\rangle_{0}
\end{array}\right] \nn
= & \frac{2t^{2}}{D^{2}}\sum_{\langle i,j\rangle,\sigma}\sum_{\langle k,l\rangle,\sigma^\prime}\left\langle \xi_{ij,\sigma}^{\dagger}\xi_{kl,\sigma^\prime}\right\rangle \left[\begin{array}{c}
\delta_{jl}\delta_{ik}\delta_{\sigma\sigma^\prime}\left(1-\langle f_{i\sigma}f_{i\sigma}^{\dagger}\rangle_{0}\right)+\delta_{jk}\delta_{il}\delta_{\sigma\sigma^\prime}\left(1-\langle f_{i\sigma}f_{i\sigma}^{\dagger}\rangle_{0}\right)\\
+\delta_{il}\delta_{jk}\delta_{\sigma\sigma^\prime}\left(1-\langle f_{j\sigma}f_{j\sigma}^{\dagger}\rangle_{0}\right)+\delta_{jl}\delta_{ik}\delta_{\sigma\sigma^\prime}\left(1-\langle f_{j\sigma}f_{j\sigma}^{\dagger}\rangle_{0}\right)
\end{array}\right] \nn
= & \frac{4t^{2}}{D^{2}}\sum_{\langle i,j\rangle,\sigma}\left\langle \xi_{ij,\sigma}^{\dagger}\xi_{ij,\sigma}\right\rangle _{0}\left[\langle f_{i,-\sigma}f_{i,-\sigma}^{\dagger}\rangle_{0}+\langle f_{j,-\sigma}f_{j,-\sigma}^{\dagger}\rangle_{0}\right].
\end{align}

This plot illustrates the schematic of a $\xi$ field (pink eclipse) on a square lattice, with its position (blue cross) being defined in between a spinon field $f$ (black circle) and slave boson $b$ (black diamond). The lattice constant of the unit cell formed by the $\xi$ field is $1/\sqrt{2}$ times of the original square lattice.

To proceed $\left\langle H_t^2 \right\rangle_0, \,\left\langle \xi_{ij,\sigma}^\dagger \xi_{ij,\sigma}\right\rangle_0$ is treated in the homogeneous and uniform limit, $ \left\langle \xi_{ij,\sigma}^\dagger \xi_{ij,\sigma}\right\rangle_0 \longrightarrow \overline{\left\langle \xi_{ij,\sigma}^\dagger \xi_{ij,\sigma}\right\rangle_0}$ where  $\overline{\left\langle \xi_{ij,\sigma}^\dagger \xi_{ij,\sigma}\right\rangle_0}$ takes the average value, independent of spins and sites.

The value of $\overline{\left\langle \xi_{ij,\sigma}^{\dagger}\xi_{ij,\sigma}\right\rangle _{0}}$ can be obtained from the doping constraint. In the uniform limit, the doping constraint over a unit cell (consisting of two $ \xi$-bonds per spin flavor, see Supplementary Figure <ref>) reads

\begin{align}
\frac{2}{N_{s}}\sum_{\langle i,j\rangle\sigma}  \overline{\left\langle \xi_{ij,\sigma}^{\dagger}\xi_{ij,\sigma}\right\rangle _{0}}=\frac{2}{N_{s}}\times4N_{s}\times  \overline{\left\langle \xi_{ij,\sigma}^{\dagger}\xi_{ij,\sigma}\right\rangle _{0}}=\delta,
\end{align}

leading to $\overline{\left\langle \xi_{ij,\sigma}^{\dagger}\xi_{ij,\sigma}\right\rangle _{0}}=\frac{\delta}{8}$. The above result gives rise to

\begin{align}
\langle H_{t}^2\rangle_{0}&=\frac{4t^{2}}{D^{2}}\times\frac{\delta}{8}\sum_{\langle i,j\rangle,\sigma}\left[\langle f_{i,-\sigma}f_{i,-\sigma}^{\dagger}\rangle_{0}+\langle f_{j,-\sigma}f_{j,-\sigma}^{\dagger}\rangle_{0}\right].
\end{align}

The original hoping $t$ is thus modified by a prefactor proportional to $\sqrt{\delta}$:

\begin{align}
\frac{\langle H_{t}^2\rangle_{0}}{8N_{s}}=\frac{t^{2}\delta}{4D^{2}}\Rightarrow t\to t_{eff}=\frac{t\sqrt{\delta}}{2}.
\end{align}

The suppression of the hoping parameter by the non-condensed slave boson mentioned above is reminiscent of the reduction (by a factor of $\delta$ instead) in hoping term within the renormalized mean-field theory approach, i.e. $t \to t_{eff} = t \times g_t$ with $g_t \sim 2\delta/(1+\delta)^2$ (see, e.g., Ref. <cit.>).

The RG-relevant action of slave-boson approach to the t-J model for the U(1) FL$^*$ phase ($\Delta = 0$ and $\langle b \rangle = 0$ but $\chi \neq 0$) <cit.> is described as

\begin{align}
S=	& -\sum_{k\sigma}f_{k\sigma}^{\dagger}\left(i\omega-\varepsilon_{\bm{k}}\right)f_{k\sigma}^{\dagger}-\sum_{k}b_{k}^{\dagger}\left(i\omega-\lambda\right)b_{k} +\frac{g}{2}\sum_{k\sigma a}\left(\xi_{k\sigma}^{a}\right)^{\dagger}\xi_{k\sigma}^{a}  +\frac{J}{2}\sum_{ka}\left(\varphi_{k}^{a}\right)^{\dagger}\varphi_{k}^{a} \nn
& +\frac{g}{\sqrt{\beta N_{s}}}\sum_{kp\sigma a}\left(f_{k\sigma}^{\dagger}b_{p}^{\dagger}\xi_{k+p,\sigma}^{a}+H.c.\right) +\frac{J}{\sqrt{\beta N_{s}}}\sum_{kpa}\left[f_{k\uparrow}^{\dagger}f_{p\downarrow}^{\dagger}\varphi_{k+p}^{a}+H.c.\right]
\label{eq:action}
\end{align}

Here, the simplified notations $k \equiv (\omega_n, \, \bm{k})$ and $p \equiv (\nu_n, \, \bm{p})$ are used to represent both Matsubara frequency ($\omega_n, \, \nu_n$) and momentum ($\bm{k}, \, \bm{p}$) simultaneously. Based on the action $S$ of Eq. (<ref>) above, we define the  bare  (without RPA corrections) Green's functions of various fields as

\begin{align}
&\mathcal{G}_{f}(k)=\frac{1}{i\omega_{n}-\varepsilon_{\bm{k}}}, \quad
\mathcal{G}_{b}(k)=\frac{1}{i\omega_{n}-\lambda},\nn
& \mathcal{G}_{\xi}(k)=\frac{2}{-g},\quad\mathcal{G}_{\varphi}(k)=\frac{2}{-J}.
\end{align}

They are diagrammatically represented in Fig. 2(a) of the main text. Note that, in  <ref> below, we will use RPA to  perturbatively generate dispersion and dynamics of the $\xi$ field. Thus, the form of its Green's function will change accordingly.

§ DYNAMICS OF VARIOUS FIELDS

In this section, we perturbatively estimate the dispersion and dynamics of the slave boson $b$, the spinon-holon bond field $\xi$  and the spinon pairing field $\varphi$.

At the bare level, the $\xi$ fermion does not have dynamics. We generate its dynamics and band dispersion by   perturbatively calculating the leading non-trivial self-energy correction, $\Sigma_{\xi}$, via Random Phase Approximation (RPA) approach where the dynamics and dispersion of the $\xi$ fermion are generated by integrating out the higher energy modes in its self-energy <cit.>. We first evaluate $\Sigma_{\xi}$: The Green's function of the $\xi$ fermion to the second order in hoping term is given by

\begin{align}
\mathcal{G}_{\xi}^{(2)}(\bm{k},\omega_{n})	& =-\frac{g^{2}}{\beta}\sum_{p,q,\sigma}\sum_{p^{\prime},q^{\prime},\sigma^{\prime}}\langle\xi_{\bm{k}\sigma}(\omega_{n})\xi_{p+q,\sigma}^{\dagger}b_{q}f_{p\sigma}f_{p^{\prime}\sigma^{\prime}}^{\dagger}b_{q^{\prime}}^{\dagger}\xi_{p^{\prime}+q^{\prime},\sigma^{\prime}}\xi_{\bm{k}\sigma}^{\dagger}(\omega_{n})\rangle_{0} \nn
& =-\frac{g^{2}}{\beta}\sum_{p,q,\sigma}\sum_{p^{\prime},q^{\prime},\sigma^{\prime}}\left\langle \xi_{\bm{k}\sigma}(\omega_{n})\xi_{p+q,\sigma}^{\dagger}f_{p\sigma}f_{p^{\prime}\sigma^{\prime}}^{\dagger}\xi_{p^{\prime}+q^{\prime},\sigma^{\prime}}\xi_{\bm{k}\sigma}^{\dagger}(\omega_{n})\right\rangle _{0}\left\langle b_{q}b_{q^{\prime}}^{\dagger}\right\rangle _{0} \nn
& =-\frac{g^{2}}{\beta}\sum_{p}\mathcal{G_{\xi}}^{2}(\bm{k},\omega_{n})\mathcal{G}_{f}(\bm{p},p_{n})\mathcal{G}_{b}(\bm{k-p},\omega_{n}-q_{n}).
\end{align}

Thus, via the Dyson equation, $G_{\xi}^{-1}=\mathcal{G}_{\xi}^{-1}-\Sigma_{\xi}$, the full $G_{\xi}$ up to the second order approximation, $  G_{\xi}(\bm{k},\omega_{n})	 \approx\mathcal{G}_{\xi}(\bm{k},\omega_{n})+\mathcal{G}_{\xi}^{(2)}(\bm{k},\omega_{n}) $, leads to

\begin{align}
G_{\xi}^{-1}(\bm{k},\omega_{n})=\mathcal{G}_{\xi}^{-1}(\bm{k},\omega_{n})- \Sigma_{\xi}(\bm{k},i\omega_{n}),
\end{align}

where we obtain

\begin{align}
\Sigma_{\xi}(\bm{k},\omega_{n})=-\frac{g^{2}}{\beta}\sum_{\bm{p},p_{n}}\mathcal{G}_{f}(\bm{p},p_{n})\mathcal{G}_{b}(\bm{k}-\bm{p},\omega_{n}-p_{n}).
\end{align}

To acquire the dispersion of the $\xi$ fermion, we only focus on real part of $\Sigma_{\xi}$. We integrate out the higher energy modes of the conduction electron, leading to

\begin{align}
\Sigma_{\xi}(\bm{k},\omega_{n})	& =-\frac{g^{2}}{\beta}\sum_{\bm{p},p_{n}}\mathcal{G}_{f}(\bm{p},p_{n})\mathcal{G}_{b}(\bm{k}-\bm{p},\omega_{n}-p_{n}) \nn
% & =g^{2}\sum_{\bm{p}}\left[\frac{n_{F}(\varepsilon_{\bm{p}})}{-i\omega_{n}+\varepsilon_{\bm{p}}+\lambda}+\frac{n_{F}(i\omega_{n}-\lambda)}{i\omega_{n}-\lambda-\varepsilon_{\bm{p}}}\right] \nn
& =\sum_{\bm{p}}\left(\frac{g^{2}}{i\omega_{n}-\lambda-\varepsilon_{\bm{p}}}\right)+\sum_{\bm{p}}^{\sim}\left(\frac{g^{2}}{-i\omega_{n}+\varepsilon_{\bm{p}}+\lambda}\right),
\end{align}

where, in the above equation, $\tilde{\sum}_{\bm{p}}$ means summing over the relatively higher energy (quasi-momenta) states whose energy levels are slightly below the Fermi surface.

Letting $i\omega_{n}\to\omega+i\epsilon$ with $\epsilon>0$ being an infinitesimal number. The real part of first term is

\begin{align}
\sum_{\bm{p}}\left(\frac{g^{2}}{\omega-\lambda-\varepsilon_{\bm{p}}}\right)= & -g^{2}\rho_{0}\left(\int_{-D}^{-D+\varepsilon_{\bm{k}}}+\int_{D-\varepsilon_{\bm{k}}}^{D}\right)\frac{d\varepsilon}{-\omega+\lambda+\varepsilon} \nn
% =& \left(-g^{2}\rho_{0}\right)\left[\ln\left(\frac{-\omega+\lambda-D+\varepsilon_{\bm{k}}}{-\omega+\lambda-D}\right)+\ln\left(\frac{-\omega+\lambda+D}{-\omega+\lambda+D-\varepsilon_{\bm{k}}}\right)\right] \nn
\approx & \left(-g^{2}\rho_{0}\right)\left[\ln\left(1+\frac{\omega-\lambda-\varepsilon_{\bm{k}}}{D}\right)-\ln\left(1+\frac{-\omega+\lambda-\varepsilon_{\bm{k}}}{D}\right)\right] \nn
\approx &  \frac{-2g^{2}\rho_{0}}{D}\left(\omega-\lambda\right),
\end{align}

and the real part of the second term is given by

\begin{align}
\sum_{\bm{p}}^{\sim}\left(\frac{g^{2}}{-\omega+\varepsilon_{\bm{p}}+\lambda}\right)
= & g^{2}\rho_{0}\int_{-D}^{-D+\varepsilon_{\bm{k}}}\frac{d\varepsilon}{-\omega+\varepsilon+\lambda} \nn
=& g^{2}\rho_{0}\ln\left(\frac{-\omega+\lambda-D+\varepsilon_{\bm{k}}}{-\omega+\lambda-D}\right) \nn
% \approx &  g^{2}\rho_{0}\ln\left(1+\frac{\omega-\lambda-\varepsilon_{\bm{k}}}{D}\right) \nn
\approx & \frac{g^{2}\rho_{0}}{D}\left(\omega-\lambda-\varepsilon_{\bm{k}}\right).
\end{align}

The real part of $\Sigma_{\xi}(\omega+i\epsilon,\bm{k})$  then reads

\begin{align}
\Sigma_{\xi}(\omega+i\epsilon,\bm{k})\approx-\frac{g^{2}\rho_{0}}{D}\varepsilon_{\bm{k}}-\frac{g^{2}\rho_{0}}{D}\left(\omega-\lambda\right).
\end{align}

Next, we evaluate the imaginary part of $\Sigma_{\xi}$, denoted as $\Sigma_{\xi}^{\prime\prime}$. Using $(x\pm i\epsilon)^{-1}=\mathcal{P}(\frac{1}{x})\mp i\pi\delta(x)$ and assuming $\left|\omega-\lambda\right|\leq D$, we have

\begin{align}
\Sigma_{\xi}^{\prime\prime}	(\omega\pm i\epsilon,\bm{k})   & =   \sum_{\bm{p}}\left(\frac{g^{2}}{\omega-\lambda-\varepsilon_{\bm{p}}\pm i\epsilon}\right)^{\prime\prime}+\sum_{\bm{p}}^{\sim}\left(\frac{g^{2}}{-\omega+\varepsilon_{\bm{p}}+\lambda\pm i\epsilon}\right)^{\prime\prime} \nn
%=\mp\pi g^{2}\sum_{\bm{p}}\delta\left(\omega-\lambda-\varepsilon_{\bm{p}}\right)\pm\pi g^{2}\sum_{\bm{p}}^{\sim}\delta\left(-\omega+\varepsilon_{\bm{p}}+\lambda\right)
& =\left(\mp\pi g^{2}\rho_{0}\right)\int_{-D}^{D}d\varepsilon\delta\left(\omega-\lambda-\varepsilon\right)\pm\pi g^{2}\rho_{0}\int_{-D}^{0}\delta\left(-\omega+\varepsilon+\lambda\right) \nn
& =\begin{cases}
\mp\pi g^{2}\rho_{0} & \text{if \ensuremath{\lambda\leq\omega}\ensuremath{\leq D+\lambda}}\\
0 & \text{otherwise}
\end{cases}
\end{align}

Since $\lambda$ is positive, it implies that $\omega$ should also be positive. In the Matsubara frequency space, we have $\Sigma_{\xi}(\bm{k},i\omega_{n})$

\begin{align}
\Sigma_{\xi}(\bm{k},\omega_{n})
= &   \sum_{\bm{p}}\left(\frac{g^{2}}{i\omega_{n}-\lambda-\varepsilon_{\bm{p}}}\right)+\sum_{\bm{p}}^{\sim}\left(\frac{g^{2}}{-i\omega_{n}+\varepsilon_{\bm{p}}+\lambda}\right) \nn
= & g^{2}\rho_{0}\int_{-D}^{D}\frac{d\varepsilon}{i\omega_{n}-\lambda-\varepsilon}+g^{2}\rho_{0}\int_{-D}^{0}\frac{d\varepsilon}{-i\omega_{n}+\varepsilon+\lambda} \nn
= & g^{2}\rho_{0}\ln\left(\frac{i\omega_{n}-\lambda+D}{i\omega_{n}-\lambda-D}\right)+g^{2}\rho_{0}\ln\left(\frac{i\omega_{n}-\lambda}{i\omega_{n}-\lambda+D}\right) \nn
= & g^{2}\rho_{0}\ln\left(\frac{i\omega_{n}+D}{i\omega_{n}-D}\right)+g^{2}\rho_{0}\ln\left(\frac{i\omega_{n}-\lambda}{i\omega_{n}+D}\right) \nn
= & -i\pi g^{2}\rho_{0}\text{sgn}(\omega_{n})+\frac{1}{2}g^{2}\rho_{0}\ln\left(\frac{\omega_{n}^{2}+\lambda^{2}}{\omega_{n}^{2}+D^{2}}\right) \nn
& \qquad +ig^{2}\rho_{0}\left[\theta_{1}(\omega_{n})-\theta_{2}(\omega_{n})\right],
\end{align}

\begin{align}
\theta_{1}(\omega_{n})	=\tan^{-1}\left(\frac{\omega_{n}}{-\lambda}\right), \,\,
\theta_{2}(\omega_{n})	=\tan^{-1}\left(\frac{\omega_{n}}{D}\right).
\end{align}

Thus, we obtain the imaginary part of $\Sigma_{\xi}$

\begin{align}
\Sigma_{\xi}^{\prime\prime}(\bm{k},\omega_{n})
=& -\pi g^{2}\rho_{0}\text{sgn}(\omega_{n})+g^{2}\rho_{0}\left[\tan^{-1}\left(\frac{\omega_{n}}{-\lambda}\right)-\tan^{-1}\left(\frac{\omega_{n}}{D}\right)\right] \nn
=& -\pi g^{2}\rho_{0}\text{sgn}(\omega_{n})-g^{2}\rho_{0}\left[\frac{\omega_{n}}{\lambda}+\frac{\omega_{n}}{D}\right] \nn
\approx & -\pi g^{2}\rho_{0}\text{sgn}(\omega_{n}) \; ,
\end{align}

where, in the last line of the above equation, we have assumed $\lambda, \, D\gg \omega_n$. Combing $\Sigma_{\xi}^{\prime}$ and $\Sigma_{\xi}^{\prime\prime}$, the dressed retarded Green's function of the $\xi$ field up to the leading non-trivial self-energy correction is given by

\begin{align}
G_{\xi}^{R}(\omega+i\epsilon,\bm{k}) & 	=\frac{1}{-t-\Sigma_{\xi}^{\prime}(\omega+i\epsilon,\bm{k})-i\Sigma_{\xi}^{\prime\prime}(\omega+i\epsilon,\bm{k})} \nn
& \equiv\frac{1}{\omega_{\xi}-\xi_{\bm{k}}+i\Gamma_{\xi}},\quad\text{with }\text{\ensuremath{\lambda}\ensuremath{\leq\omega}\ensuremath{\leq D},}
\end{align}

where the quasi-particle width $\Gamma_{\xi}$, the generated  frequency $\omega_{\xi}$, and the dispersion $\xi_{\bm{k}}$ are given by

\begin{align}
&  \Gamma_{\xi}\equiv\pi g^{2}\rho_{0}, \nn
&  \omega_{\xi}\equiv\frac{g^{2}\rho_{0}}{D}\omega\equiv\zeta\omega, \nn
& \xi_{\bm{k}}= -\zeta\varepsilon_{\bm{k}}+\zeta\lambda+t
\end{align}

Here, $\zeta\equiv\frac{g^{2}\rho_{0}}{D}$, $\rho_{0}\sim1/D$ and $\varepsilon_{\bm{k}} = h_{\bm{k}} -\mu_f \approx v |\bm{k}| $ with $\mu_f \equiv \mu - \lambda $ being the effective chemical potential for $f$ spinon. Here, the $f$-spinon band is approximated by a linear-in-momentum dispersion near the effective chemical potential $\mu_f$. As we treat $g$ as a perturbation, we expect $\zeta=(g/D)^{2}\ll1$. In addition, as the $\xi$ fermion is dispersionless at the bare level, its bandwidth of the renormalized band structure falls in the range of  $-\zeta D \leq \zeta \varepsilon_{\bm{k}}\leq\zeta D$ for $\bm{k}$ to be satisfactory with $-D\leq\varepsilon_{\bm{k}}\leq D$. The value of $\zeta D=g^{2}/D\sim t\zeta\ll t$. This leads to $\xi_{\bm{k}}>0$.

We now generate and estimate  the dispersion of the slave boson $b$. At the bare level, the slave boson $b$ shows  only flat band with an atomic energy level $\lambda$. Here, we estimate the dispersion of the slave boson by computing its effective hoping through second-order perturbation in $H_t$. The effective hoping of slave boson can be obtained by contracting the $\xi$ and $f$ fermions of the following

\begin{align}
& g^2 \sum_{ i,j,\sigma}\sum_{ i^\prime,j^\prime,\sigma^\prime} \left\langle  f_{i\sigma}^{\dagger} \xi_{ij,\sigma}\xi^\dagger_{i^\prime j^\prime,\sigma^\prime}f_{i^\prime \sigma^\prime} \right\rangle b_{j}^{\dagger}b_{j^\prime}
\sim  g^2 \sum_{ i,j,\sigma}\sum_{ i^\prime,j^\prime,\sigma^\prime} \left\langle  f_{i\sigma}^{\dagger} f_{i^\prime \sigma^\prime} \right\rangle \left\langle \xi_{ij,\sigma}\xi^\dagger_{i^\prime j^\prime,\sigma^\prime} \right\rangle b_{j}^{\dagger}b_{j^\prime}.
\end{align}

Taking $\left\langle  f_{i\sigma}^{\dagger} f_{i^\prime \sigma^\prime} \right\rangle \sim \chi $ and $\left\langle \xi_{ij,\sigma}\xi^\dagger_{i^\prime j^\prime,\sigma^\prime} \right\rangle  \sim g^2/D^2 $, the dispersion of the slave boson field thus take the form $\varepsilon_b (\bm{k}) = \lambda + t_b (\cos k_x +\cos k_y)$ with the effective hoping  $t_b \propto g^4 \chi = t^4 \delta^2 \chi$ (in unit of D).  For small hole-doping $\delta$ and a weakly perturbative coupling $t$ where $t/D \approx 3/4$ with $D \approx 4J_H$, $t_b$ can be estimated as $t_b \sim t^4 \delta^2 \sim (3/4)^4 \delta^2 $ (in unit of D). We can fix the value of $\lambda$ such that $D>\lambda \gg |t_b|$, and therefore the effecitve hoping $t_b$ term of the generated slave boson band is negligible.

We now evaluate the leading-order self-energy correction to generate dynamics of the $\varphi$ field. The self-energy $\Sigma_{\varphi}$ proportional to $J^{2}$ is given by

\begin{align}
\Sigma_{\varphi}\left(\bm{k},\omega_{n}\right)	 =& \left(-\frac{J^{2}}{\beta}\right)\sum_{q_{n},\bm{q}}\mathcal{G}_{f}\left(\bm{k}-\bm{q},\omega_{n}-q_{n}\right)\mathcal{G}_{f}\left(\bm{q},q_{n}\right) \nn
= & \left(-\frac{J^{2}}{\beta}\right)\sum_{q_{n},\bm{q}}\frac{1}{iq_{n}-\varepsilon_{\bm{q}}}\frac{1}{i\omega_{n}-iq_{n}-\varepsilon_{\bm{k-q}}} \nn
%= & \left(\frac{J^{2}}{\beta}\right)\sum_{q_{n},\bm{q}}\frac{1}{iq_{n}-\varepsilon_{\bm{q}}}\frac{1}{-i\omega_{n}+iq_{n}+\varepsilon_{\bm{k-q}}} \nn
= & J^{2}\sum_{\bm{q}}\left[\frac{n_{F}(\varepsilon_{\bm{q}})}{-i\omega_{n}+\varepsilon_{\bm{q}}+\varepsilon_{\bm{k-q}}}+\frac{n_{F}(i\omega_{n}-\varepsilon_{\bm{k-q}})}{i\omega_{n}-\varepsilon_{\bm{q}}-\varepsilon_{\bm{k-q}}}\right] \nn
%=& J^{2}\sum_{\bm{q}}\left[\frac{n_{F}(\varepsilon_{\bm{q}})}{-i\omega_{n}+\varepsilon_{\bm{q}}+\varepsilon_{\bm{k-q}}}+\frac{1-n_{F}(\varepsilon_{\bm{k-q}})}{i\omega_{n}-\varepsilon_{\bm{q}}-\varepsilon_{\bm{k-q}}}\right] \nn
= & J^{2}\sum_{\bm{q}}\frac{1-n_{F}(\varepsilon_{\bm{k}-\bm{q}})-n_{F}(\varepsilon_{\bm{q}})}{i\omega_{n}-\varepsilon_{\bm{q}}-\varepsilon_{\bm{k-q}}}.
\end{align}

By setting the external momentum $\bm{k}=0$, we have

\begin{align}
\Sigma_{\varphi}\left(\omega_{n}\right)	=J^{2}\sum_{\bm{q}}\frac{1-2n_{F}(\varepsilon_{\bm{q}})}{i\omega_{n}-2\varepsilon_{\bm{q}}}=J^{2}\sum_{\bm{p}}\frac{1}{i\omega_{n}-2\varepsilon_{\bm{q}}}-2J^{2}\sum_{\bm{p}}^{\sim}\frac{1}{i\omega_{n}-2\varepsilon_{\bm{q}}}.
\end{align}

In the real frequency representation, the real part of $\Sigma_{\varphi}$ is

\begin{align}
\Sigma_{\varphi}^{\prime}\left(\bm{k},\,\omega + i\epsilon\right)	=& J^{2}\rho_{0}\left(\int_{D-\varepsilon_{\bm{k}}}^{D}+\int_{-D}^{-D+\varepsilon_{\bm{k}}}\right)\frac{d\varepsilon}{\omega-2\varepsilon}-2J^{2}\rho_{0}\int_{-D}^{-D+\varepsilon_{\bm{k}}}\frac{d\varepsilon}{\omega-2\varepsilon} \nn
= & \left(-\frac{J^{2}\rho_{0}}{2}\right)\left(\int_{2D-2\varepsilon_{\bm{k}}}^{2D}+\int_{-2D}^{-2D+2\varepsilon_{\bm{k}}}\right)\frac{du}{u-\omega}+\left(J^{2}\rho_{0}\right)\int_{-2D}^{-2D+2\varepsilon_{\bm{k}}}\frac{du}{u-\omega} \nn
= &\left(-\frac{J^{2}\rho_{0}}{2}\right)\ln\left|\frac{-\omega+2D}{-\omega+2D-2\varepsilon_{\bm{k}}}\right|+\left(-\frac{J^{2}\rho_{0}}{2}\right)\ln\left|\frac{\omega+2D-2\varepsilon_{\bm{k}}}{\omega+2D}\right|
+J^{2}\rho_{0}\ln\left|\frac{\omega+2D-2\varepsilon_{\bm{k}}}{\omega+2D}\right| \nn
= & \frac{J^{2}\rho_{0}}{2}\ln\left|\frac{-\omega+2D-2\varepsilon_{\bm{k}}}{-\omega+2D}\right|+\frac{J^{2}\rho_{0}}{2}\ln\left|\frac{\omega+2D-2\varepsilon_{\bm{k}}}{\omega+2D}\right| \nn
= & \frac{J^{2}\rho_{0}}{2}\ln\left[\frac{\left(2D-2\varepsilon_{\bm{k}}\right)^{2}-\omega^{2}}{(2D)^{2}-\omega^{2}}\right],
\end{align}

where the logarithmic term can be approximated as

\begin{align}
\ln\frac{\left(2D-2\varepsilon_{\bm{k}}\right)^{2}-\omega^{2}}{(2D)^{2}-\omega^{2}}	\approx & \ln\frac{\left(2D-2\varepsilon_{\bm{k}}\right)^{2}-\omega^{2}}{(2D)^{2}} \nn
\approx &  \ln\left[\frac{1}{\left(2D\right)^{2}}\left[(2D)^{2}-8D\varepsilon_{\bm{k}}-\omega^{2}\right]\right] \nn
= & \ln\left[1-\frac{2\varepsilon_{\bm{k}}}{D}-\frac{\omega^{2}}{4D^{2}}\right] \nn
= & -\frac{2\varepsilon_{\bm{k}}}{D}-\frac{\omega^{2}}{4D^{2}}.
\end{align}

The imaginary part of $\Sigma_{\varphi}$ vanishes in the wide-band and low-energy (frequency) limit, $D \gg \omega$:

\begin{align}
\Sigma_{\varphi}^{\prime\prime}\left(\bm{k},\omega+i\epsilon\right)	=-\pi J^{2}\rho_{0}\left(\int_{D-\varepsilon_{\mathbf{k}}}^{D}+\int_{-D}^{-D+\varepsilon_{\mathbf{k}}}\right)d\varepsilon\delta(\omega-2\varepsilon)+2\pi J^{2}\rho_{0}\int_{-D}^{-D+\varepsilon_{\mathbf{k}}}d\varepsilon\delta(\omega-2\varepsilon)
= 0.
\end{align}

Combining $\Sigma_{\varphi}^{\prime}$ and $\Sigma_{\varphi}^{\prime\prime}$, we have

\begin{align}
\Sigma_{\varphi}=\left(-\frac{J^{2}\rho_{0}}{2}\right)\left(\frac{2\varepsilon_{\bm{k}}}{D}+\frac{\omega^{2}}{4D^{2}}\right)\equiv-\varepsilon_{\varphi}(\bm{k})-\frac{\omega_{\varphi}^{2}}{D}.
\end{align}

\begin{align}
\omega_{\varphi}=\omega\sqrt{\frac{\kappa}{8}}\quad,\quad\varepsilon_{\varphi}(\bm{k})=\kappa\varepsilon_{\bm{k}}\quad,\text{with }\kappa\equiv\frac{J^{2}\rho_{0}}{D}=\frac{J^{2}}{D^{2}}.
\end{align}

The renormalized Green's function $G_{\varphi}$ is given by

\begin{align}
\end{align}

Since $J/t \ll 1$,  the generated frequency and dispersion which are of the order of $J^2/D^2$ can be thus neglected.

§ RENORMALIZATION GROUP (RG) ANALYSIS

In this section, we provide details of the RG analysis of our slave-boson approach to the Kondo-Heisenberg-like t-J model. The  RG procedure is primarily based on Ref. <cit.> .

§.§ Scaling dimension

The RG procedure starts from the scaling of the wave vector $\bm{k}$ and the frequency $\omega$, where

\begin{align}
k_i^{\prime}=e^{\ell}k_i,\quad\omega^{\prime}=e^{z \ell}\omega,
\end{align}

where $\ell$ is called the scaling factor, $i$ above denotes the component of $\bm{k}$, and $z$ represents the dynamical exponent. The bare scaling dimension for wave vector $k_i$ is chosen to be $[k_i]=1$. To preserve the Luttinger's theorem, the Fermi velocity has to be kept scale invariant. Since the dispersion of the $f$ spinon band is linear-in-momentum, we thus have

\begin{align}
\end{align}

The scaling dimension for the measure $d^{d}k$ for fermions with scale-invariant Fermi velocity is given by <cit.>

\begin{align}
\end{align}

The scale invariance for the action of the kinetic term of the $f$ spinon, $S_{0}^{f}$, implies that the bare scaling dimension of the  $f$ spinon operator is

\begin{align}
\end{align}

where $k \equiv (\omega,\, \bm{k})$. Since the bare slave boson $b_k$ is local, its scaling dimension can be thus obtained as

\begin{align}
\end{align}

The non-interacting part of the composite fermion field $\xi^a_{k\sigma}$ is given by

\begin{align}
S_{0}^{\xi}\sim\sum_{\sigma a}\int d\omega d^d k \left(\omega_{\xi}-\xi_{\bm{k}}\right)\left( \xi^a_{k\sigma}\right)^{\dagger}\xi^a_{k\sigma} \nc
\end{align}

where the renormalized frequency $\omega_{\xi}=\frac{g^{2}\rho_{0}}{D}\omega$ and renormalized dispersion $\xi_{\bm{k}}=-\frac{g^{2}\rho_{0}}{D}\varepsilon_{\bm{k}}+\text{const.}$. While performing the scaling procedure for $\xi_{\bm{k}}$ and $\omega_{\xi}$, we only scale $\varepsilon_{\bm{k}}$ and $\omega$. The prefactor $\frac{g^{2}\rho_{0}}{D}$ is not involved in the  scaling. This is equivalent to fixing $\frac{g^{2}\rho_{0}}{D}$ at its fixed point value, i.e. $\frac{g^{2}\rho_{0}}{D}\to \frac{(g^*)^{2}\rho_{0}}{D} $. We thus have the bare scaling dimension of the $\xi_{k\sigma}$,

\begin{align}
[\xi^a_{k\sigma}]=-\frac{3}{2} \np
\end{align}

In our theory, since the dispersion for the RVB pairing bond field $\varphi_{k}$ by perturbative correction is sub-leading,  we approximate the $\varphi^a_{k}$ field to be local. Hence, its bare scaling dimension is identical to that of the slave boson $b$, namely,

\begin{align}
[\varphi^a_k] = -\left(\frac{d+z}{2}\right).
\end{align}

The scaling dimension for various couplings can be computed in a similar approach. The action for $H_{t}$, denoted as $S_{t}$, takes the following form,

\begin{align}
S_{t}\sim \bar{g}\int d\omega d\nu\int d^{d}kd^{d}p\sum_{a\sigma}\left[f_{k\sigma}^{\dagger}b_{p}^{\dagger}\bar{\xi}_{k+p,\sigma}^{a}+H.c.\right].
\label{eq:S_t-action}
\end{align}

with $k=(\omega,\bm{k})$ and $p=(\nu,\bm{p})$. From Ref. <cit.>, for a special case $z=1$, the calculation of the scaling dimension for a boson-fermion coupling term can be performed via the usual way without invoking the patch scheme. The scaling dimension of $\bar{g}$ can be obtained similarly by simple power-counting, $ 0=[\bar{g}]+2z+1+d+[f_{k\sigma}]+[\xi_{k\sigma}]+[b_k]$. Consequently, we obtain

\begin{align}
[\bar{g}]	=\frac{4-d-3z}{2}=-\frac{1}{2} \np
\end{align}

Note that, in Eq. (<ref>), $\bm{k}$ denote the wave vector $\bm{k}$ for the $f$ spinon while $\bm{p}$ for the slave boson. Therefore, by Ref. <cit.>, we take $[d^dk] = 1$ and $[d^d p] = d$ in the above derivation for $[\bar{g}]	$. Following a similar approach, the bare scaling dimension of the Heisenberg coupling is given by

\begin{align}
[J]=\frac{4-d-3z}{2}=-\frac{1}{2} \np
\end{align}

§.§ Self-energy corrections

The self-energy correction to the bare Green's function of a given field begins with the following relation between the dressed (renormalized) Green's function $G$ of that field:

\begin{align}
G=\mathbb{G}\mathcal{G} \nc
\end{align}

with $\mathcal{G}$ being the bare Green's function and $\mathbb{G}$ being the corrections to $\mathcal{G}$ by perturbation expansion. From the Dyson equation, we can also show that $G^{-1}=\mathcal{G}^{-1}-\Sigma=\mathcal{G}^{-1}\left(1-\mathcal{G}\Sigma\right)$ with $\Sigma$ here being denoted as the self-energy correction, implying that $ \mathbb{G}=\frac{1}{1-\mathcal{G}\Sigma}$. For weak coupling theory, we can expand $\mathbb{G}$ as a power series,

\begin{align}
\mathbb{G}\approx1+\mathcal{G}\Sigma+O(\mathcal{G}\Sigma)^{2}.
\end{align}

Note that the self-energy correction $\Sigma_{\xi}^{g}$ for the $\xi$ field in the RG analysis requires an additional negative sign compared to $\Sigma_{\xi}^{g}$ from RPA. This difference arises from the minus sign  appearing in the $k$-dependent term in the dispersion $\xi_{\bm{k}}$ for the $\xi$ field computed from RPA. As a result, $\Sigma_{\xi}^{g}$ is changed to

\begin{align}
\Sigma_{\xi}^{g} (\omega_n)&=\left(\frac{\bar{g}^{2}}{N_{s}\beta}\right)\sum_{\bm{p},p_{m}}\frac{1}{ip_{m}-\varepsilon_{\bm{p}}}\cdot\frac{1}{i\omega_{n}-ip_{m}-\lambda} \nn
& = \left(-\frac{\bar{g}^{2}}{N_{s}}\right)\sum_{\bm{p}}\left[\frac{n_{F}(\varepsilon_{\bm{p}})}{-i\omega_{n}+\varepsilon_{\bm{p}}+\lambda}+\frac{n_{F}(i\omega_{n}-\lambda)}{i\omega_{n}-\lambda-\varepsilon_{\bm{p}}}\right].
\label{eq:Sigma_xi_RG}
\end{align}

Taking $-i\omega_{n}=-\omega-\lambda$, we have

\begin{align}
\Sigma_{\xi}^{g} (\omega)  = \left(-\bar{g}^{2}\rho_{0}\right)\int_{-D}^{0}\left(\frac{d\varepsilon}{-\omega+\varepsilon}\right)+\text{others}  = \bar{g}^{2}\rho_{0}\ln\frac{D}{\omega}+\text{others} \nc
\end{align}

where “others” in the above equation represents the terms with no logarithmic correction.

The self-energy of the $f$ spinon, $\Sigma_{f}^{g}$,  contributed from perturbation of $H_t$ is given by

\begin{align}
\Sigma_{f}^{g} (\bm{k}, \omega_n) &=\left(-\frac{\bar{g}^{2}}{N_{s}\beta}\right) \sum_{p}\overline{\mathcal{G}}_{\xi}\left(k+p\right)\mathcal{G}_{b}\left( p \right) \nn
& =\left(-\frac{\bar{g}^{2}}{N_{s}\beta}\right)\sum_{p_{m},\bm{p}}\frac{1}{i\omega_{n}+ip_{m}-\zeta^{-1}\xi_{\bm{k}+\bm{p}}}\cdot\frac{1}{ip_{m}-\lambda} \nn
& = \frac{\bar{g}^{2}}{N_{s}}\sum_{\bm{p}}\frac{n_{F}(\zeta^{-1}\xi_{\bm{p}})}{i\omega_{n}+\varepsilon_{\bm{p}}+\zeta^{-1}\varUpsilon+\lambda},
\end{align}

where $\varUpsilon\equiv-\zeta\lambda-t+\mu_{\xi}$. After substituting $i\omega_{n} + \lambda \to \omega$, we obtain

\begin{align}
\Sigma_{f}^{g} (\omega) = & \bar{g}^{2}\rho_{0}\int_{-\zeta^{-1}\varUpsilon}^{D}\frac{d\varepsilon}{\omega+\varepsilon} \nn
= & \bar{g}^{2}\rho_{0}\ln\left(\frac{D}{\omega-\zeta^{-1}\varUpsilon}\right).
\end{align}

From the result in RPA section, we know that the self-energy of the slave boson contributes no logarithmic divergence, $\Sigma^g_b \propto (\ln D)^0$.

Now, we evaluate the self-energy correction of various fields from the perturbation of the Heisenberg term, $H_J$. The self-energy of the f spinon is given by

\begin{align}
\Sigma_{f}^{J}(\bm{k},\omega_{n})=\frac{J^{2}}{N_{s}\beta}\sum_{q_{m},\bm{q}}\mathcal{G}_{f,\downarrow\downarrow}(\bm{q},q_{m})\mathcal{\bar{G}}_{\varphi}(\bm{k+q},\omega_{n}+q_{m}).
\end{align}

Using the $\varphi$ propagator

\begin{align}
\bar{\mathcal{G}}_{\varphi}(\omega_{n},\bm{k})=-1,
\end{align}

$\Sigma_{f}^{J}$ reads

\begin{align}
\Sigma_{f}^{J}(\bm{k},\omega_{n})	& =\left(-\frac{J^{2}}{N_{s}\beta}\right)\sum_{q_{m},\bm{q}}\frac{1}{ip_{m}-\varepsilon_{\bm{q}}}
=\left(-J^{2}\rho_{0}\right)\int_{-D}^{D}d\varepsilon n_{F}(\varepsilon) \nn
& =-J^{2}\rho_{0}\int_{-D}^{0}d\varepsilon
\end{align}

We find no logarithmic correction from $\Sigma_{f}^{J}$, i.e. $\Sigma_{f}^{J}  \sim\left(\ln D\right)^{0}$. Additionally, no imaginary part for $\Sigma_{f}^{J}$ exists.

We now evaluate the self-energy correction $\Sigma_{\varphi}$ for the fluctuating $\varphi$, which is given by

\begin{align}
\Sigma_{\varphi}^{J}\left(\bm{k},\omega_{n}\right)	& =\left(\frac{J^{2}}{N_{s}\beta}\right)\sum_{q_{n},\bm{q}}\mathcal{G}_{f}\left(\bm{k}-\bm{q},\omega_{n}-q_{n}\right)\mathcal{G}_{f}\left(\bm{q},q_{n}\right) \nn
&=\frac{J^{2}}{N_{s}\beta}\sum_{q_{n},\bm{q}}\frac{1}{iq_{n}-\varepsilon_{\bm{q}}}\frac{1}{i\omega_{n}-iq_{n}-\varepsilon_{\bm{k-q}}} \nn
& =\left(\frac{-J^{2}}{N_{s}}\right)\sum_{\bm{q}}\frac{1-n_{F}(\varepsilon_{\bm{k}-\bm{q}})-n_{F}(\varepsilon_{\bm{q}})}{i\omega_{n}-\varepsilon_{\bm{q}}-\varepsilon_{\bm{k-q}}}.
\end{align}

Since our purpose is the find to logarithmic correction due to the internal momentum integral, we thus can set the external momentum to zero. This leads to

\begin{align}
\Sigma_{\varphi}^{J}\left(\omega_{n}\right)=\left(\frac{-J^{2}}{N_{s}}\right)\sum_{\bm{q}}\frac{1-2n_{F}(\varepsilon_{\bm{q}})}{i\omega_{n}-2\varepsilon_{\bm{q}}}=\left(2J^{2}\rho_{0}\right)\int_{-D}^{0}\frac{d\varepsilon}{i\omega_{n}-2\varepsilon}.
\end{align}

Taking $i\omega_{n}\to\omega+i\epsilon$ and using $(x\pm i\epsilon)^{-1}=\mathcal{P}(\frac{1}{x})\mp i\pi\delta(x)$, the self-energy correction becomes

\begin{align}
\Sigma_{\varphi}^{J}\left(\omega\right)	& =\left(2J^{2}\rho_{0}\right)\int_{-D}^{0}\frac{d\varepsilon}{\omega-2\varepsilon+i\epsilon}=\left(-J^{2}\rho_{0}\right)\int_{-2D}^{0}\frac{d\varepsilon^{\prime}}{\varepsilon^{\prime}-\omega-i\epsilon} \nn
& =\left(-J^{2}\rho_{0}\right)\int_{-2D}^{0}d\varepsilon^{\prime}\left[\mathcal{P}(\frac{1}{\varepsilon^{\prime}-\omega})+i\pi\delta(\varepsilon^{\prime}-\omega)\right].
\end{align}

We find the real part for $\Sigma_{\varphi}^{J}$ is given by

\begin{align}
\left[  \Sigma_{\varphi}^{J}\left(\omega\right) \right]^\prime &	=\left(-J^{2}\rho_{0}\right)\int_{-2D}^{0}d\varepsilon^{\prime}\mathcal{P}(\frac{1}{\varepsilon^{\prime}-\omega}) \nn
& =\left(-J^{2}\rho_{0}\right)\ln\frac{\omega}{2D+\omega}\approx\left(-J^{2}\rho_{0}\right)\ln\frac{\omega}{2D} \nn
& \sim J^{2}\rho_{0}\ln\frac{D}{\omega},
\end{align}

which contributes a logarithmic correction that is relevant for further RG analysis. While the imaginary part can be evaluated as

\begin{align}
\left[\Sigma_{\varphi}^{J}\left(\omega\right) \right]^{\prime\prime}=\left(-J^{2}\pi\rho_{0}\right)\int_{-2D}^{0}d\varepsilon^{\prime}\delta(\varepsilon^{\prime}-\omega)=-J^{2}\pi\rho_{0}\Theta(-\omega).
\end{align}

§.§ Vertex corrections

Vertex correction is defined as the correction for the bare vertex function due to interactions. If we denote the full vertex function as $\Pi$ and the bare one as $\Pi^{0}$. The vertex correction, represented as $U$, is defined through

\begin{align}
\Pi = U \Pi^{0}.
\label{eq:def-vertex-correction}
\end{align}

While employing perturbation expansion, $U$ can be written as $U=1+U^{(1)}+U^{(2)}+U^{(3)}+\cdots$ with $U^{(n)}$ being the vertex correction from the $n$-th order perturbation.

It is straightforward to show that there is no vertex corrections up to one-loop order contributed from the hopping  term $H_{t}$. The only relevant one-loop vertex correction is coming from the Heisenberg term $H_J$. The vertex function of $H_J$, denoted as $\Pi_{J}$, is defined as

\begin{align}
\Pi_{J}\left(k,p\right)=\left\langle \varphi_{k+p}f_{k\uparrow}^{\dagger}f_{p\downarrow}^{\dagger}\right\rangle ,
\end{align}

where $k=(\bm{k},k_{n}),\, p=(\bm{p},p_{n})$. Expanding $\Pi_{J}$ above, it is easy to show that no zero-th order term is available. The first non-vanishing contribution is from the first-order expansion in $J$, given by

\begin{align}
\Pi_{J}^{(1)}\left(k,p\right)&=J\sum_{k^{\prime},p^{\prime}}\left\langle \varphi_{k+p}f_{k\uparrow}^{\dagger}f_{p\downarrow}^{\dagger}\varphi_{k^{\prime}+p^{\prime}}^{\dagger}f_{p^{\prime}\downarrow}f_{k^{\prime}\uparrow}\right\rangle _{0} \nn
& =  \left(-J\right)\mathcal{\bar{G}}_{\varphi}\left(k+p\right)\mathcal{G}_{f}^{\downarrow\downarrow}(p)\mathcal{G}_{f}^{\uparrow\uparrow}(k).
\label{eq:bare-vertex-Pi_J}
\end{align}

Thus, $\Pi_{J}^{(1)}\left(k,p\right)=\Pi_{J}^{0}$  is considered as the bare vertex function for the fluctuating Heisenberg term $H_J$. It is easy to demonstrate that $\Pi^{(2)}\propto J^{2}$ is zero. Thus, the lowest non-vanishing correction to the bare vertex function $\Pi_{J}^{(1)}$ comes from $\Pi^{(3)}\propto J^3$, which takes the following form

\begin{align}
\Pi_{J}^{(3)}(k,p)  = & \left(-\frac{3J^{3}}{3!\beta}\right)\sum_{k^{\prime},p^{\prime}}\sum_{k^{\prime\prime},p^{\prime\prime}}\sum_{k^{\prime\prime\prime},p^{\prime\prime\prime}}\left\langle \varphi_{k+p}f_{k\uparrow}^{\dagger}f_{p\downarrow}^{\dagger}\varphi_{k^{\prime}+p^{\prime}}^{\dagger}f_{p^{\prime}\downarrow}f_{k^{\prime}\uparrow}\varphi_{k^{\prime\prime}+p^{\prime\prime}}^{\dagger}f_{p^{\prime\prime}\downarrow}f_{k^{\prime\prime}\uparrow}\varphi_{k^{\prime\prime\prime}+p^{\prime\prime\prime}}f_{k^{\prime\prime\prime}\uparrow}^{\dagger}f_{p^{\prime\prime\prime}\downarrow}^{\dagger}\right\rangle_{0} \nn
= & \left(\frac{J^{3}}{\beta}\right)\sum_{p^{\prime}}\sum_{k^{\prime\prime}}\left\langle \varphi_{k+p}\varphi_{k^{\prime\prime}+p^{\prime}}\right\rangle _{0}\left\langle \varphi_{k+p^{\prime}}^{\dagger}\varphi_{k^{\prime\prime}+p}^{\dagger}\right\rangle _{0}\left\langle f_{k\uparrow}f_{k\uparrow}^{\dagger}\right\rangle _{0}\left\langle f_{p\downarrow}f_{p\downarrow}^{\dagger}\right\rangle _{0}\left\langle f_{k^{\prime\prime}\uparrow}f_{k^{\prime\prime}\uparrow}^{\dagger}\right\rangle _{0}\left\langle f_{p^{\prime}\downarrow}f_{p^{\prime}\downarrow}^{\dagger} \right\rangle_{0}.
\label{eq:Pi_J_3-1}
\end{align}

From Eqs. (<ref>), (<ref>) and (<ref>), the vertex correction of $ \Pi_{J}^{(1)}$ is related to the third-order perturbation,  denoted as $U_{J}^{(3)}(k,p)$, through $\Pi_{J}^{(3)}(k,p) = U_{J}^{(3)}(k,p)\Pi_{J}^{(1)}\left(k,p\right)$, with

\begin{align}
U_{J}^{(3)}(k,p)=\left(-\frac{J^{2}}{\beta}\right) &\sum_{p^{\prime}}\mathcal{\bar{G}}_{\varphi}\left(k+p^{\prime}\right)\mathcal{G}_{f}^{\uparrow\uparrow}(-k-p-p^{\prime}) \mathcal{G}_{f}^{\downarrow\downarrow}(p^{\prime}).
\end{align}

Using $\mathcal{\bar{G}}_{\varphi}=-1$, we have

\begin{align}
U_{J}^{(3)}(k,p)&=\frac{J^{2}}{\beta}\sum_{p^{\prime}}\mathcal{G}_{f}^{\uparrow\uparrow}(-k-p-p^{\prime})\mathcal{G}_{f}^{\downarrow\downarrow}(p^{\prime}) \nn
& = \frac{J^{2}}{\beta}\sum_{p^{\prime}}\mathcal{G}_{f}^{\uparrow\uparrow}(-p^{\prime})\mathcal{G}_{f}^{\downarrow\downarrow}(p^{\prime}) \nn
& =\frac{J^{2}}{\beta}\sum_{ip_{n}^{\prime},\bm{p}^{\prime}}\frac{1}{ip_{n}^{\prime}-\varepsilon_{\bm{p}^{\prime}}}\cdot\frac{1}{-ip_{n}^{\prime}-\varepsilon_{\bm{p}^{\prime}}} \nn
& =\left(-J^{2}\rho_{0}\right)\int_{-D}^{D}d\varepsilon\frac{n_{F}(\varepsilon)-n_{F}(-\varepsilon)}{2\varepsilon} \nn
& =\left(-J^{2}\rho_{0}\right)\ln\frac{D}{r}\Big|_{r\to0}+\text{others}.
\end{align}

We obtain a $\ln D$ term from the vertex correction of $H_J$.

In summary, the self-energy corrections of various fields from $H_{t}$ and $H_{J}$ are given by

\begin{align}
\Sigma_{\xi}^{g}(\omega) & \sim\bar{g}^{2}\rho_{0}\ln\frac{D}{\omega}, \nn
\Sigma_{b}^{g}(\omega)	& \sim(\ln D)^{0}, \nn
\Sigma_{f}^{g}(\omega)	& \sim\bar{g}^{2}\rho_{0}\ln\left(\frac{D}{\omega-\zeta^{-1}\varUpsilon}\right), \nn
\Sigma_{\varphi}^{J}(\omega) & \sim J^{2}\rho_{0}\ln\frac{D}{\omega}.
\end{align}

This leads to the following form of field renormalization $ \mathbb{G}$:

\begin{align}
\mathbb{G}_{f}	& \approx 1+ \left( \bar{g}^{2}\rho_{0}\right)\mathcal{G}_{f}\ln\left(\frac{D}{\omega-\zeta^{-1} \varUpsilon} \right)
\sim  1+ \left(\bar{g}^{2}\rho_{0}\right)\frac{1}{\omega}\ln\left(\frac{D}{\omega-\zeta^{-1}\varUpsilon}\right) \nc \nn
\mathbb{G}_{\xi} & \approx 1+ \bar{g}^{2}\rho_{0} \mathcal{\bar{G}}_{\xi}\ln\frac{D}{\omega} \sim 1+\bar{g}^{2}\rho_{0}\frac{1}{\omega}\ln\frac{D}{\omega}\nc \nn
\mathbb{G}_{\varphi} & \approx 1+ J^{2} \rho_{0} \mathcal{\bar{G}}_{\varphi}\ln\frac{D}{\omega} \sim 1+ J^{2}\rho_{0}\frac{1}{\omega}\ln\frac{D}{\omega} \nc \nn
\mathbb{G}_{b} & = 1 \np
\end{align}

Here, we have replaced the bare Green's function with $\mathcal{G}\sim1/\omega$. In addition, we find the perturbative expansion to the vertex function of $H_J$ and $H_t$ contributes the following vertex corrections:

\begin{align}
U_{J}(k,p)	& =1-J^{2}\rho_{0}\ln\frac{D}{r}\Big|_{r\to0}\nc \nn
U_{g}(k,p)	& =1 \np
\end{align}

§.§ The renormalization $z$-factors  and the RG $\beta$ functions

To compute the RG $\beta$ functions, we start by reducing the cutoff energy from $D\to D^{\prime}=yD$ with $0<y<1$. This generates the relations between the wave-function (Green's function) and vertex corrections between $D$ and $D^\prime$, which shows that <cit.>

\begin{align}
\mathbb{G}(\omega,D^{\prime},J_{i}^{\prime}) & =z_{\psi}\text{\ensuremath{\mathbb{G}}}(\omega,D,J_{i})\nc \nn U(\omega,D^{\prime},J_{i}^{\prime}) & =z_{U}^{-1}U(\omega,D,J_{i})
\label{eq:G-U-renormal}
\end{align}

with $z_\psi$ and $z_U$ being denoted as the renormalization factor for the Green's function for the $\psi$ field and vertex correction $U$, respectively. The renormalization factor for the $\xi$ field is given by

\begin{align}
z_{\xi}	& = \left(1+\bar{g}^{2}\rho_{0}\frac{1}{\omega}\ln\frac{D}{\omega}\right)^{-1}\left(1+\bar{g}^{2}\rho_{0}\frac{1}{\omega}\ln\frac{D^{\prime}}{\omega}\right) \nn
& =1-\frac{\bar{g}^{2}\rho_{0}}{\omega}\ln \ell \nc
\end{align}

where $\ln \ell \equiv\ln\frac{D}{D^{\prime}}$. As $\omega$ does not undergo a RG scaling [see Eq. (<ref>)], we can set, for convenience, $\omega=D\to1/\omega=1/D=\rho_{0}$. This yields

\begin{align}
z_{\xi}=1-\bar{g}^{2}\rho_{0}^{2}\ln \ell \np
\end{align}

Following a similar approach, we obtain the other two $z$ factors of fields,

\begin{align}
z_{f} & =1-\left(\bar{g}\rho_{0}\right)^{2}\ln \ell \nc \nn
z_{\varphi} & =1-\left(J\rho_{0}\right)^{2}\ln \ell \np
\end{align}

The $z$ factors for the couplings are

\begin{align}
& z_g =1 \nnc
&z_{J}=1+J^{2}\rho_{0}^{2}\ln \ell \np
\end{align}

§.§ The RG scaling equations ($\beta$ functions)

The renormalized couplings constants are defined as

\begin{align}
\bar{g}^{\prime}&=z_{f}^{-\frac{1}{2}}z_{b}^{-\frac{1}{2}}z_{\xi}^{-\frac{1} {2}}z_{g}\bar{g} \nnc
J^{\prime} & =z_{f}^{-1}z_{\varphi}^{-\frac{1}{2}}z_{J}J \np
\end{align}

Expanding $\bar{g}^\prime$ to the leading order in $\bar{g}\rho_0$, we obtain

\begin{align}
\bar{g}^{\prime}&	=\left(1-\left(\bar{g}\rho_{0}\right)^{2}\ln \ell\right)^{-\frac{1}{2}}\left(1-\bar{g}^{2}\rho_{0}^{2}\ln \ell\right)^{-\frac{1}{2}}\bar{g} \nn
& =\left(1+\frac{1}{2}\left(\bar{g}\rho_{0}\right)^{2}\ln \ell\right)\left(1+\frac{1}{2}\left(\bar{g}\rho_{0}\right)^{2}\ln \ell\right)\bar{g} \nn
& =\bar{g}+\bar{g}\left(\bar{g}\rho_{0}\right)^{2}\ln \ell \np
\end{align}

Similarly, $J^\prime$ shows

\begin{align}
J^{\prime}	& =\left[1-\left(\bar{g}\rho_{0}\right)^{2}\ln \ell\right]^{-1} \left[1-\left(J\rho_{0}\right)^{2}\ln \ell\right]^{-\frac{1}{2}}\left[1+J^{2}\rho_{0}^{2}\ln \ell \right]J \nn
& =\left[1+\left(\bar{g}\rho_{0}\right)^{2}\ln \ell \right] \left[1+\frac{1}{2}\left(J\rho_{0}\right)^{2}\ln \ell \right]\left(1+J^{2}\rho_{0}^{2}\ln \ell \right)J \nn
& =J+J\left[\left(\bar{g}\rho_{0}\right)^{2}+\frac{3}{2}\left(J\rho_{0}\right)^{2}\right]\ln \ell \np
\end{align}

Multiplying both sides of $J^{\prime}$ and $\bar{g}^{\prime}$ by $\rho_{0}$, we have

\begin{align}
\left(J\rho_{0}\right)^{\prime}	& =\bar{g}\rho_{0}+\left(\bar{g}\rho_{0}\right)^{3}\ln \ell \nnc
(\bar{g}\rho_{0})^{\prime}	& =J\rho_{0}+\left[\left(J\rho_{0}\right)\left(\bar{g}\rho_{0}\right)^{2}+\frac{3}{2}\left(J\rho_{0}\right)^{3}\right]\ln \ell \np
\end{align}

These lead to the following RG $\beta$ functions:

\begin{align}
\frac{d(\bar{g}\rho_{0})}{dl}	 & =-\left(\frac{d-z}{2}\right)(\bar{g}\rho_{0}) + \left(\bar{g}\rho_{0}\right)^{3} \nnc
\frac{d(J\rho_{0})}{dl}	& =-\left(\frac{d-z}{2}\right)(J\rho_{0})+\left(J\rho_{0}\right)\left(\bar{g}\rho_{0}\right)^{2}+\frac{3}{2}\left(J\rho_{0}\right)^{3} \np
\end{align}

The fixed points occur at  $\left(\left(\bar{g}\rho_{0}\right)^{2},\quad\left(J\rho_{0}\right)^{2}\right)=\left(0,\quad0\right),\quad Q:\left(0,\quad\frac{\epsilon}{3}\right),\quad P:\left(\frac{\epsilon}{2},\quad0\right)$
where $\epsilon \equiv d-z$. In our case, we set $d=2, \, z=1$ and thus  $\epsilon=1$. Note that the RG scheme is still reliable for $\epsilon=1$ though it makes the $\epsilon$ expansion divergent <cit.>.

§.§ Rescale the coupling constants

In this section, we will show how the coupling constants are  rescaled near the critical point $C_r$. Near the fixed point $P$, we fix $\bar{g}\rho_0$ at its fixed point value, namely $\bar{g}\rho_0 =\bar{g}^* \rho_0 = \sqrt{\epsilon/2} $. As $\bar{g}\rho_0$ is now fixed at a constant, it is equivalent of making the $\bar{g}\rho_0$  coupling marginal since it cannot flow anymore. Making $\bar{g}\rho_0$ marginal can be done by $H_t \sim \left( \bar{g}/e^{-\epsilon \ell/2}\right) e^{-\epsilon \ell/2} f^\dagger b^\dagger \xi  = \bar{g}^\prime f^{\dagger\,\prime} b^{\dagger \prime} \xi^\prime  $ with $\bar{g}^\prime = \bar{g} /e^{-\epsilon \ell/2} $ and $f^\prime = e^{-\epsilon \ell/2} f$ and $b^\prime = b$ and $\xi^\prime = \xi$. This leads to a shift in the bare scaling dimensions $[\bar{g}^\prime] = [\bar{g}]+\epsilon/2 = 0$ and $[f^\prime]=[f] +\epsilon/4$. We attribute the remaining exponential terms into $b$ and $\xi$ and change their scaling dimensions accordingly. In the meantime, due to the change of the scaling dimension of $f$, the scaling dimension of $J$ needs to be changed to $[J^\prime] = [J] - \epsilon/2$. To lighten the notation, here we will use the same symbols for the couplings after the rescaling as we did before the rescaling. Fixing $\bar{g}\rho_{0} = \left(\bar{g}\rho_{0}\right)^* = \sqrt{\epsilon/2}$  yields the following RG equation for $J\rho_0$ near $P$:

\begin{align}
\frac{d(J\rho_{0})}{dl}&=\left(-\frac{\epsilon}{2} -\frac{\epsilon}{2}\right)\left(J\rho_{0}\right)+\left(J\rho_{0}\right)\left[\left(\bar{g}\rho_{0}\right)^*\right]^{2}+\frac{3}{2}\left(J\rho_{0}\right)^{3} \nn
=& -\frac{\epsilon}{2} \left(J\rho_{0}\right)+\frac{3}{2}\left(J\rho_{0}\right)^{3},
\end{align}

which yields the same fixed point value for $ J\rho_0$, namely $\left[\left( J\rho_0\right)^*\right]^2  = \epsilon/3 $.

Following the similar approach, we evaluate the RG flow of $\bar{g}\rho_0$ near $Q$ by fixing $J\rho_0 =J^* \rho_0 = \sqrt{\epsilon/3}$. We make $J$ marginal by $[J^\prime] = [J] + \epsilon/2$ and $[\varphi^\prime] =[\varphi]-\epsilon/2 $. The RG equation for $\bar{g}$ does not change after this transformation. The RG flow diagram  of the rescaled coupling constants is shown in Fig. 3(a) of the main text.

§ SPECIFIC HEAT COEFFICIENT

Power-law divergence of the electronic specific heat coefficient $C_V/T$ for various samples around the optimal doping. Experimental data are reproduced from Ref. <cit.>. The blue line corresponds to the power-law fit for which we find a coefficient $b \approx 0.5$ in good agreement with our theoretical prediction [cf. Eq. (<ref>)].

In this section, we estimate the specific heat coefficient contributed solely from the half-filled f-spinon band. The entropy of a non-interacting many-fermion system with Hamiltonian $H_f=\sum_{\bm{k}\sigma}\varepsilon_{\bm{k}}f_{\bm{k}\sigma}^{\dagger}f_{\bm{k}\sigma}$ with  dispersion $\varepsilon_{\bm{k}}=-2 J_H \chi \left(\cos k_{x}+\cos k_{y}\right)$ is given by

\begin{align}
S_f(T)= & -2\sum_{\bm{k}}  \left[n_{F} (\varepsilon_{\bm{k}})\ln n_F (\varepsilon_{\bm{k}})
+\left(1-n_F \left(\varepsilon_{\bm{k}} \right)\right)\ln\left(1-n_{F} (\varepsilon_{\bm{k}})\right)\right] \nn
= &-2\int_{-D(T)}^{0}  d\varepsilon\rho(\varepsilon,T)\left[n_F(\varepsilon)\ln n_F(\varepsilon)
\end{align}

where $n_{F} (\varepsilon_{\bm{k}})=1/\left(e^{\varepsilon_{\bm{k}}/T}+1\right)$ and the temperature-dependent half-bandwidth $D(T) = 4 J_H \chi(T)$. The Boltzmann constant is set to be $k_{B}=1$, while the factor of $2$ in $S_f$ comes from the spins.

For the low-temperature regime, we can approximate $D(T) \approx D(T=0) \approx 4 J_H \chi(T = 0 ) = 4 J_H \chi$ with $\chi \equiv \chi(T=0)$, leading to $\rho(\varepsilon,T)\to\rho(\varepsilon,T=0)\equiv\rho_{0}(\varepsilon)$. The entropy becomes approximately

\begin{align}
S_f(T)\approx-2\int_{-D_{0}}^{0}d\varepsilon\rho_{0}(\varepsilon)\left[n_F(\varepsilon)\ln n_F(\varepsilon)
\end{align}

At half-filling, the density of states for the tight-binding dispersion $\varepsilon_{\bm{k}}$ on a  2D square lattice  exhibits a logarithmic divergence in the low-energy regime, i.e., $ \rho_{0}(\varepsilon) \sim - \left(\frac{1}{\Lambda} \right)\ln \left(\varepsilon / \varepsilon_0 \right)$ as $\varepsilon \to 0$, where the prefactor $\Lambda$, carrying the same unit of energy, is proportional to $J$, and $\varepsilon_0$ represents a constant energy scale. Consequently, the entropy in the low-energy regime is dictated by the logarithmically divergent density of states, given by

\begin{align}
S_f(T)	\sim &  \frac{2}{\Lambda}\int_{-D}^{0}d\varepsilon\ln \left(\varepsilon / \varepsilon_0 \right) \left[n_F(\varepsilon)\ln n_F(\varepsilon)   +\left(1-n_F(\varepsilon)\right)\ln\left(1-n_F(\varepsilon)\right)\right] \vspace{10pt} \nn
\approx & \frac{2}{\Lambda} \int_{-D/T}^{0}\left(Tdx\right)\ln\left( \frac{Tx}{\varepsilon_0}\right)\left[\frac{1}{e^{x}+1}\ln\left(\frac{1}{e^{x}+1}\right)+\left(1-\frac{1}{e^{x}+1}\right)\ln\left(1-\frac{1}{e^{x}+1}\right)\right] \vspace{10pt} \nn
= & \frac{2}{\Lambda} \int_{-\infty}^{0}\left(Tdx\right)\left[\ln \left( \frac{T}{\varepsilon_0}\right)+\ln x\right]F(x),
\end{align}

where in the last line of the above equation, we take the approximation $D/T\to\infty $  for $T \to 0$, and where we define

\begin{align}
F(x)\equiv & \frac{1}{e^{x}+1}\ln\left(\frac{1}{e^{x}+1}\right)  +\left(1-\frac{1}{e^{x}+1}\right)\ln\left(1-\frac{1}{e^{x}+1}\right).
\end{align}

For the low-temperature limit, the entropy can be compactly expressed as

\begin{align}
S_f(T)	=  \frac{\alpha^\prime T}{\Lambda}\ln \left( \frac{T}{\varepsilon_0}\right)+ \frac{\beta^\prime T}{\Lambda}
\label{eq:S-low-T}
\end{align}

with $\alpha^\prime, \, \beta^\prime$ being constant prefactors:

\begin{align}
\alpha^\prime \equiv 2\int_{-\infty}^{0}F(x)dx,\quad \beta^\prime \equiv2\int_{-\infty}^{0} dx F(x)\ln x.
\end{align}

The contribution of the non-interacting $f$ spinon to the specific heat coefficient at the low-temperature regime can be thus obtained by taking the temperature derivative of $S_f$ in Eq. (<ref>), $\frac{C^f_V}{T} = \frac{\partial S_f(T)}{\partial T}$, showing a logarithmic-in-T divergence

\begin{align}
\frac{C^f_V}{T} = \frac{\alpha^\prime}{\Lambda} \ln \left( \frac{\gamma^\prime T}{\varepsilon_0}\right)
\end{align}

with $\gamma^\prime = e^{\frac{\alpha^\prime + \beta^\prime}{\alpha^\prime}}$.

At low temperatures and near the QCP, the behavior of the specific heat coefficient $\frac{C^f_V}{T}$ has to be carefully estimated by considering the coupling constant renormalization. Here,  $\frac{C^f_V}{T}$ only depends on the Heisenberg coupling; thus, under the coupling constant renormalization $J \to J(\ell)$, we have $\frac{C^f_V}{T} \sim J^{-1}(\ell_0)$ at some fixed renormalization scale $\ell = \ell_0$. The renormalization of the Heisenberg coupling can be expressed as $J(\ell_0) = J e^{-\epsilon \ell_0 /2}$.  Using the relation of $e^{\ell_0} = \eta \sim |\delta - \delta_c|^{-\nu}$ with $\eta$ being the correlation length and $\nu = 1/z = 1$ the correlation length exponent in our theory,  $\frac{C^f_V}{T}$ near the QCP behaves as

\begin{align}
\frac{C^f_V}{T} \sim |\delta -\delta_c|^{-1/2}.
\label{eq:C_V^f-pwr}
\end{align}

The low-temperature behavior of  $\frac{C^f_V}{T} $ contributed from the half-filled $f$-spinon band exhibits a power-divergence as a function of $\delta -\delta_c$ with power-law exponent equal to $0.5$. A similar behavior has been experimentally observed in the cuprate superconductors <cit.>, see Supplementary Figure <ref>. Note that due to the divergence of the specific heat coefficient in Eq. (<ref>), the T-logarithmic marginal Fermi liquid behavior of the specific heat coefficient contributed from f-spinon band dominates over that contributed from the fluctuating charge (Kondo) term.

§ TRANSPORT PROPERTIES

§.§ Total electrical resistivity by Ioffe-Larkin's composition rule

Below, we use the Ioffe-Larkin composition rule <cit.> to compute the total resistivity (or conductivity) of the Kondo-Heisenberg approach to the slave boson t-J model [Eq. (1) of the main text]. To apply the Ioffe-Larkin composition rule in the calculation of the total resistivity, we need to analyze the gauge transformation and include both the  intrinsic U(1) gauge field and external gauge field in the form of vector potential of external magnetic field into this model. The model, described by the Lagrangian density below

\begin{align}
\mathcal{L}(\tau)	& =\sum_{i}b_{i}^{\dagger}\left(\partial_{\tau}+i\lambda_{i}\right)b_{i}+\sum_{i\sigma}f_{i\sigma}^{\dagger}\left(\partial_{\tau}+i\lambda_{i}-\mu\right)f_{i\sigma}+t\sum_{\langle i,j\rangle,\sigma}\left[\left(f_{i\sigma}^{\dagger}b_{j}^{\dagger}+f_{j\sigma}^{\dagger}b_{i}^{\dagger}\right)\xi_{ij,\sigma}+H.c.\right] \nn
& \qquad+\sum_{\langle i,j\rangle\sigma}\left(-\chi_{ij}f_{i\sigma}^{\dagger}f_{j\sigma}+H.c.\right)+\sum_{\langle i,j\rangle,\alpha}\left(\Delta_{ij}\tilde{\sigma}f_{i\sigma}^{\dagger}f_{j,-\sigma}^{\dagger}+H.c.\right)+{\rm others},
\end{align}

is invariant under the following gauge transformations

\begin{align}
b_{i}	&\to b_{i}e^{i\theta_{i}},\nn
f_{i\sigma}	& \to f_{i\sigma}e^{i\theta_{i}},\nn
\xi_{ij,\sigma}	& \to\xi_{ij,\sigma}e^{i(\theta_{i}+\theta_{j})},\nn
\Delta_{ij}	&\to \Delta_{ij}e^{i(\theta_{i}+\theta_{j})}, \nn
\chi_{ij}	&\to \chi_{ij}e^{i(-\theta_{i}+\theta_{j})}, \nn
\lambda_{i}	&\to \lambda_{i}-\partial_{\tau}\theta_{i}.
\end{align}

We fix $i\lambda_{i}$ and $\chi_{ij}$ at their saddle point solution such that  $i\lambda_{i}$ and $\chi_{ij}$ in our model are replaced by the mean-field values i.e. $i\lambda_{i}\to\lambda$ and $\langle\chi_{ij}\rangle \to \chi$. The fluctuations beyond the mean-field $\lambda$ and $\chi$ can be included in the effective Lagrangian by writing $i\lambda_{i} =\lambda+ia_{0}(\bm{r}_{i},\tau)$ and $\chi_{ij}=\chi e^{ia_{ij}}$, where $a_{0}(\bm{r}_{i},\tau)$ and $a_{ij}$ are the fluctuating fields for $\lambda$ and $\chi$ ($a_{ij}$ is the phase fluctuation for $\chi$). With the inclusion  of the fluctuations in the amplitude of  $\lambda$, the non-interacting parts for $f$ and $b$ can be rewritten as

\begin{align}
& \sum_{i}b_{i}^{\dagger}\left(\partial_{\tau}+i\lambda_{i}\right)b_{i}+\sum_{i\sigma}f_{i\sigma}^{\dagger}\left(\partial_{\tau}+i\lambda_{i}-\mu\right)f_{i\sigma} \nn
\longrightarrow &	\sum_{i}b_{i}^{\dagger}\left[\partial_{\tau}+\lambda+ia_{0}(\bm{r}_{i},\tau)\right]b_{i}+\sum_{i\sigma}f_{i\sigma}^{\dagger}\left[\partial_{\tau}+\lambda-\mu+ia_{0}(\bm{r}_{i},\tau)\right]f_{i\sigma}.
\end{align}

Note that in our further derivation of the conductivity (or the resistivity),  we will drop the non-interacting Lagrangian density of the $b$ field since its dynamics and dispersion are negligibly small. In addition, when the phase fluctuation $a_{ij}$ is included, the fermion hopping term becomes  $-\sum_{\langle i,j\rangle\sigma}\chi_{ij}f_{i\sigma}^{\dagger}f_{j\sigma}\to-\sum_{\langle i,j\rangle\sigma}\chi e^{ia_{ij}}f_{i\sigma}^{\dagger}f_{j\sigma}$, which, in the continuum limit, can be expressed approximately in the following familiar form,

\begin{align}
-\frac{1}{2m_{f}^{*}}\sum_{\sigma}\int d^{2}rf_{\sigma}^{\dagger}(\bm{r},\tau)\left[-i\nabla_{\bm{r}}+\bm{a}(\bm{r},\tau)\right]^{2}f_{\sigma}(\bm{r},\tau).
\end{align}

Within the saddle point solutions, the gauge transformations now becomes $a_{0}(\bm{r}_{i},\tau)\to a_{0}(\bm{r}_{i},\tau)-\partial_{\tau}\theta_{i}$ and $a_{ij}\to a_{ij}-\theta_{i}+\theta_{j}$. Thus, we observe that the fluctuating field $a_0(\bm{r},\tau)$ governs the gauge transformation for the scalar potential, while the phase field $a_{ij}$ governs the gauge transformation for the vector potential. The latter becomes apparent when we take the continuum limit of our model.

In the continuum limit, the phase field $a_{ij}$ can be approximate in the form of <cit.>

\begin{align}
a_{ij}(\tau)=(\bm{r}_{i}-\bm{r}_{j})\cdot\bm{a}\left(\frac{\bm{r}_{i}+\bm{r}_{j}}{2},\tau\right) \sim\int d\bm{r}\cdot\bm{a}(\bm{r},\tau),
\end{align}

which defines the gauge field $\bm{a}(\bm{r},\tau)$. The original form of the gauge transformation for $a_{ij}$ has now been modified to $\bm{a}(\bm{r},\tau) \to \bm{a}(\bm{r},\tau) -\nabla \theta (\bm{r},\tau)$.

The $\xi$ field should also couple to the gauge field since we introduce dynamics and dispersion to the $\xi$ field via Random Phase Approximation (RPA). This introduces an additional contribution to the conductivity via the Ioffe-Larkin composition rules, as described below. To demonstrate how the $\xi$ field couples with the gauge fields, we need to integrate out the $f$ spinon and the slave boson $b$, which generate terms of the form $\sim\xi_{ij,\sigma}^{\dagger}\xi_{jk,\sigma}\chi_{kl}\chi_{li}$ <cit.>. With a finite mean-field value for $\chi$, the $\xi$ field will couple to the phase field of $\chi_{ij}$, namely $a_{ij}$. In the continuum limit, this coupling of the $\xi$ field to the phase field is described by the following Lagrangian density  <cit.>,

\begin{align}
\sum_{\sigma}\int d^{2}r\xi_{\sigma}^{\dagger}(\bm{r},\tau)\left[\partial_{\tau}-2ia_{0}(\bm{r},\tau)\right]\xi_{\sigma}(\bm{r},\tau)-\frac{1}{2m_{\xi}^{*}}\sum_{\sigma}\int d^{2}r\xi_{\sigma}^{\dagger}(\bm{r},\tau)\left[-i\nabla_{\bm{r}}+2\bm{a}(\bm{r},\tau)\right]^{2}\xi_{\sigma}(\bm{r},\tau).
\end{align}

This term implies that $\xi$ is a charge-2$e$ fermion field. A similar term also exists for the pairing bond field $\Delta_{ij}$. However, we will neglect this contribution in our further analysis of the conductivity since the dynamics and dispersion of $\Delta_{ij}$ have been shown to be subdominant.

Next, we consider the effect of an external electromagnetic field coupling to the system. The coupling of the system to the electromagnetic field can be described by Peierls' substitution,

\begin{align}
tc_{i\sigma}^{\dagger}c_{j\sigma}\to te^{-ie\int_{i}^{j}d\bm{r}\cdot\bm{A}_{em}}c_{i\sigma}^{\dagger}c_{j\sigma}=t\sum_{\langle i,j\rangle,\sigma}\left[e^{-ie\int_{i}^{j}d\bm{r}\cdot\bm{A}_{em}}\left(f_{i\sigma}^{\dagger}b_{j}^{\dagger}+f_{j\sigma}^{\dagger}b_{i}^{\dagger}\right)\xi_{ij,\sigma}+H.c.\right].
\end{align}

We can couple the external vector potential $\bm{A}_{em}$ either to $f$, $b$, or $\xi$, but not to all simultaneously. Here, we choose $\bm{A}_{em}$ to couple to the $\xi$ field. The non-interacting part of the Lagrangian for the $\xi$ field shown above should be replaced by

\begin{align}
\sum_{\sigma}\int d^{2}r\xi_{\sigma}^{\dagger}(\bm{r},\tau) & \left[\partial_{\tau}-2ia_{0}(\bm{r},\tau)-iA_{em,0}(\bm{r},\tau)\right]  \xi_{\sigma}(\bm{r},\tau) \nn
&-\frac{1}{2m_{\xi}^{*}}\sum_{\sigma}\int d^{2}r\xi_{\sigma}^{\dagger}(\bm{r},\tau)\left[-i\nabla_{\bm{r}}+2\bm{a}(\bm{r},\tau)+\bm{A}_{em}(\bm{r},\tau)\right]^{2}\xi_{\sigma}(\bm{r},\tau).
\end{align}

The non-interacting effective Lagrangian, which includes all fields relevant for the following conductivity analysis coupled to the gauge fields, can be summarized as

\begin{align}
& \sum_{i\sigma}f_{i\sigma}^{\dagger}\left[\partial_{\tau}+\lambda-\mu+ia_{0}(\bm{r}_{i},\tau)\right]f_{i\sigma}-\frac{1}{2m_{f}^{*}}\sum_{\sigma}\int d^{2}rf_{\sigma}^{\dagger}(\bm{r},\tau)\left[-i\nabla_{\bm{r}}+\bm{a}(\bm{r},\tau)\right]^{2}f_{\sigma}(\bm{r},\tau)\nn
&\quad\qquad+\sum_{\sigma}\int d^{2}r\xi_{\sigma}^{\dagger}(\bm{r},\tau)\left[\partial_{\tau}-2ia_{0}(\bm{r},\tau)\right]\xi_{\sigma}(\bm{r},\tau)-\frac{1}{2m_{\xi}^{*}}\sum_{\sigma}\int d^{2}r\xi_{\sigma}^{\dagger}(\bm{r},\tau)\left[-i\nabla_{\bm{r}}+2\bm{a}(\bm{r},\tau)\right]^{2}\xi_{\sigma}(\bm{r},\tau).
\end{align}

Following the rules by Ioffe and Larkin <cit.>, we integrate out the $f$ and $\xi$ fermion fields to obtain the effective action of the gauge fields,

\begin{align}
S_{eff}	& \sim\sum_{\omega_{n},\bm{k}}\sum_{\mu=0,\perp}\Pi_{\mu}^{f}(\omega_{n},\bm{k})a_{\mu}(\omega_{n},\bm{k})a_{\mu}(-\omega_{n},-\bm{k}) \nn
&\qquad\qquad+\sum_{\omega_{n},\bm{k}}\sum_{\mu=0,\perp}\Pi_{\mu}^{\xi}(\omega_{n},\bm{k})\left[2a_{\mu}(\omega_{n},\bm{k})+A_{em,\mu}(\omega_{n},\bm{k})\right]\left[2a_{\mu}(-\omega_{n},-\bm{k})+A_{em,\mu}(-\omega_{n},-\bm{k})\right] \nn
&=\sum_{\omega_{n},\bm{k}}\sum_{\mu=0,\perp}\left[\Pi_{\mu}^{f}(\omega_{n},\bm{k})+4\Pi_{\mu}^{\xi}(\omega_{n},\bm{k})\right]\left|a_{\mu}(\omega_{n},\bm{k})\right|^{2}+\sum_{\omega_{n},\bm{k}}\sum_{\mu=0,\perp}\Pi_{\mu}^{\xi}(\omega_{n},\bm{k})\left|A_{em,\mu}(\omega_{n},\bm{k})\right|^{2} \nn
& \qquad\qquad+\sum_{\omega_{n},\bm{k}}\sum_{\mu=0,\perp}4\Pi_{\mu}^{\xi}(\omega_{n},\bm{k})A_{em,\mu}(\omega_{n},\bm{k})a_{\mu}(-\omega_{n},-\bm{k})
\end{align}

with $\Pi_{\alpha}^{f}\left(\Pi_{\alpha}^{\xi}\right)$ being the current-current correlation function for the $f \, (\xi)$ fermion. And $\mu \in 0,\, \perp$ denotes the scalar and transverse components. The transverse component is directly related to the conductivity. When further integrating out $a_{\mu}(\omega_{n},\bm{k})$, the resulting effective action for $\bm{A}_{em,\mu}$ describes the response of the system to the external electromagnetic field,

\begin{align}
S_{eff}\sim\sum_{\omega_{n},\bm{k}}\sum_{\mu=0,\perp}\Pi_{\mu}^{{\rm tot}}(\omega_{n},\bm{k})A_{em,\mu}(\omega_{n},\bm{k})A_{em,\mu}(-\omega_{n},-\bm{k}),
\end{align}

where $\Pi_{\mu}^{{\rm tot}}$ represents the correlation function of the total current, given by

\begin{align}
\Pi_{\mu}^{{\rm tot}} = \frac{\Pi_{\mu}^{\xi}\Pi_{\mu}^{f}}{\Pi_{\mu}^{f}+4\Pi_{\mu}^{\xi}}.
\end{align}

Since the transverse component of current-current correlation function is related to the conductivity through $\Pi_{\mu}^{f}(\omega_{n}, \bm{k}=0) = \sigma_{f}|\omega_{n}|$ (same for the $\xi$ field), the total conductivity is thus approximately  given by $\sigma_{{\rm tot}} = \frac{\sigma_{f}\sigma_{\xi}}{\sigma_{f} + 4\sigma_{\xi}}$, implying that the total resistivity is approximately given by

\begin{align}
\rho_{{\rm tot}} \approx 4\rho_{f}+\rho_{\xi}.
\end{align}

From Drude's model, the electrical resistivity is proportional to the effective mass, we thus obtain $\rho_{{\rm tot}}\approx\rho_{\xi}$ since the ratio of the effective masses for the $f$ spinon and $\xi$ field is

\begin{align}
m_{\xi}^{*}/m_{f}^{*}\sim\frac{\frac{1}{\zeta\chi}}{\frac{1}{\chi}} = \frac{1}{\zeta} \propto \frac{(D/t)^2}{\delta}  \gg 1 \nc
\end{align}

for doping $0.2<\delta <0.3$ and $t/D \approx 3/4$.

The resistivity from the $\xi$ fermion can be further spitted into two contributions: one from the scattering of the $\xi$ fermions by $H_t$, denoted as $\rho_{\xi,g}$, and the other from scattering of the $\xi$ fermions with gauge fields, denoted as $\rho_{\xi,\rm{gauge}}$. Calculation of $\rho_{\xi,g}$ is provided in the section below. The electron scattering rate $1/\tau_{\xi, \rm{gauge}}$ contributed from the gauge field, and the corresponding conductivity $\sigma_{\xi, \rm{gauge}}$ and resistivity $\rho_{\xi, \rm{gauge}}$ related by $\rho_{\xi, \rm{gauge}} = 1/\sigma_{\xi, \rm{gauge}} \propto 1/\tau_{\xi, \rm{gauge}}$ within a similar slave-boson approach has been estimated in Ref. <cit.>:

\begin{align}
\tau^{-1}_{\xi,\rm{gauge}} &\sim \left(\frac{T}{\chi_d}\right)^{4/3} \varepsilon_{F,\xi}
\sim  \left( \frac{T}{J} \right)^{4/3}  \frac{t^2 \delta^2}{D} ,
\end{align}

where $\chi_d \sim (m_f^*)^{-1} \sim J$ and the Fermi energy for the $\xi$ band $\varepsilon_{F,\xi} \propto \delta \times (g^2/D^2) \times D \sim t^2 \delta^2 / D$. Note that the estimate in Ref. <cit.> is applicable here for the $\xi$-band since the coupling of the gauge-field to the fermion band  ($\xi$-band here) and the gauge-field propagator are of the same form as that in our model. The scattering rate $1/\tau_{\xi, \rm{gauge}}$ receives a sub-leading $T$-super-linear power-law correction $1/\tau_{\xi, \rm{gauge}} \sim \bar{a} \times (T/D)^{4/3}$  with a negligible prefactor  $\bar{a}/D \sim (4t^2/D^2) \times \delta^2 \gg 1$. For the relevant temperature range $0<T<300$ K, it is apparent that $1/\tau_{\xi,\rm{gauge}}$ is negligible compared to the leading Planckian scattering rate $1/\tau = \alpha_P k_B T/\hbar$ with $\alpha_P \sim O(1) $ from the hopping (effective Kondo) term:

\begin{align}
\frac{\tau_{\xi, \rm{gauge}}^{-1}}{\tau_{\xi,g}^{-1}}\sim  \left( \frac{T}{D}\right)^{1/3} \left(\frac{t}{D}\right)^2  \delta^2 \ll 1
\end{align}

for $0<T<300$ K, $D \sim  1$ eV, $\delta = 0.2 \sim 0.3$ and  $1/\tau_{\xi, g}  \sim  T $ (set $k_B = \hbar = 1$ here).

§.§ Evaluation of $\Sigma_{f}$ and its scattering rate

To evaluate the scattering rate for the $f$ spinon, we start by calculating its self-energy contributed from the perturbation of the hopping term, $H_t$, denoted as $\Sigma_{f}$:

\begin{align}
\Sigma_{f}^{g}(\bm{k},\omega_{n})&=\left(-\frac{g^{2}}{N_{s}\beta}\right)\sum_{p_{m}}\sum_{\bm{p}}\mathcal{G}_{\xi}\left(p\right)\mathcal{G}_{b}\left(p-k\right) \nn
&  =-\frac{g^{2}}{N_{s}\beta}\sum_{p_m}\sum_{\bm{p}}\frac{1}{ip_{\xi}-\xi_{\bm{p}}-i\Sigma_{\xi}^{\prime\prime}(ip_{m})}
\times \frac{1}{ip_{m}-i\omega_{n}-\lambda}.
\label{eq:Sigma_g_f-scattering-rate}
\end{align}

Note that, due to the negative sign in front of the (k-dependent) dispersion $\xi_{\bm{k}}$, the imaginary part of the self-energy $\Sigma_{\xi}^{\prime\prime}$ for the $\xi$ field in the above equation will require an additional minus sign, which differs from $\Sigma_{\xi}$ evaluated via RPA. Similar situation also occurs for the calculation of $\Sigma_{\xi}$ in the RG section. Thus, we take $\Sigma_{\xi}^{\prime\prime}(ip_{m})=\pi g^{2}\rho_{0}\text{sgn}(p_{m})=\pi D\zeta\text{sgn}(p_{m})$. Due to the sign function in $\Sigma_{\xi}^{\prime\prime}$, Eq. (<ref>) can be expressed as  $\Sigma^g_{f}(\bm{k},i\omega_{n})=\Sigma_{f}^{g,>} (\bm{k},i\omega_{n}) +\Sigma_{f}^{g,<} (\bm{k},i\omega_{n})$ where $\Sigma_{f}^{g,\gtrless}$ corresponds to the $p_{m} \gtrless 0$ branch in the Matsubara sum:

\begin{align}
\Sigma_{f}^{g,\gtrless} (\bm{k},\omega_{n}) 	& \equiv-\frac{g^{2}}{\beta}\sum_{p_{m}\gtrless0}\sum_{\bm{p}}\frac{1}{ip_{\xi}-\xi_{\bm{p}}\mp i\pi g^{2}\rho_{0}}\cdot\frac{1}{ip_{m}-i\omega_{n}-\lambda} \nn
& =-\frac{g^{2}}{\zeta\beta}\sum_{p_{m}\gtrless0}\sum_{\bm{p}}\frac{1}{ip_{m}+\left(\varepsilon_{\bm{p}}-\bar{\mu}_{\xi}\right)\mp iD\pi}\cdot\frac{1}{ip_{m}-i\omega_{n}-\lambda}.
\label{eq:Sigma_f-RL}
\end{align}

To analytically evaluate Eq. (<ref>), we need to analyze the following contour integral, defined as

\begin{align}
\mathcal{S}_{F}^{\gtrless}	& =-\frac{g^{2}}{\zeta}\oint_{C_{>}}\frac{dz}{2\pi i}\,\frac{n_{F}(z)}{z+\left(\varepsilon_{\bm{p}}-\bar{\mu}_{\xi}\right)\mp iD\pi}  \times \frac{1}{z-i\omega_{n}-\lambda},
\label{eq:S_F-RL}
\end{align}

where $C_{>}$ encircles the upper hemisphere of the $z$ plane counterclockwise where $p_{m}>0$. There are two approaches to evaluate this contour integral either by residue theorem and line integral along $C_{>}$. On one hand, $\mathcal{S}_{F}^{>}$ can be  evaluated by using the residue theorem, namely

\begin{align}
\mathcal{S}_{F}^{>}	= &   \frac{g^{2}}{\zeta}\Bigg[\frac{1}{\beta}\sum_{p_{m}>0}\frac{1}{ip_{m}+\left(\varepsilon_{\bm{p}}-\bar{\mu}_{\xi}\right)-iD\pi}\cdot\frac{1}{ip_{m}-i\omega_{n}-\lambda}-\frac{n_{F}(i\omega_{n}+\lambda)\Theta(\omega_{n})}{i\omega_{n}+\left(\varepsilon_{\bm{p}}-\bar{\mu}_{\xi}\right)+\lambda-iD\pi} \nn
&\qquad\qquad\qquad+\frac{n_{F}(-\varepsilon_{\bm{p}}+\bar{\mu}_{\xi}+iD\pi)}{i\omega_{n}+\varepsilon_{\bm{p}}-\bar{\mu}_{\xi}+\lambda-iD\pi}\Bigg] \nn
= & \frac{g^{2}}{\zeta}\Bigg[\frac{1}{\beta}\sum_{p_{m}>0}\frac{1}{ip_{m}+\left(\varepsilon_{\bm{p}}-\bar{\mu}_{\xi}\right)-iD\pi}\cdot\frac{1}{ip_{m}-i\omega_{n}-\lambda}+\frac{n_{B}(\lambda)}{i\omega_{n}+\left(\varepsilon_{\bm{p}}-\bar{\mu}_{\xi}\right)+\lambda-iD\pi} \nn
&\qquad\qquad\qquad+\frac{n_{F}(-\varepsilon_{\bm{p}}+\bar{\mu}_{\xi}+iD\pi)}{i\omega_{n}+\varepsilon_{\bm{p}}-\bar{\mu}_{\xi}+\lambda-iD\pi}\Bigg] \nn
= & \frac{g^{2}}{\zeta}\Bigg[\frac{1}{\beta}\sum_{p_{m}>0}\frac{1}{ip_{m}+\left(\varepsilon_{\bm{p}}-\bar{\mu}_{\xi}\right)-iD\pi}\cdot\frac{1}{ip_{m}-i\omega_{n}-\lambda}+\frac{n_{F}(-\varepsilon_{\bm{p}}+\bar{\mu}_{\xi}+iD\pi)}{i\omega_{n}+\varepsilon_{\bm{p}}-\bar{\mu}_{\xi}+\lambda-iD\pi}\Bigg].
\end{align}

On the other hand, $\mathcal{S}_{F}^{>}$ can be alternatively analyzed via directly doing the line integral along  $C_{>}$ on the complex plane: first perform the integral along the real axis and then along the upper semicircle counterclockwise. The integrand  takes the approximate form $\sim\frac{1}{e^{\beta z}+1}\frac{1}{z^{2}}$. For the integral over $-\infty\leq|z|\leq 0$, integral along $[-\infty,0]$ gives $\mathcal{S}_{F}^{>}$ to be zero since the integrand is approximated as $\frac{1}{e^{\beta z}+1}\frac{1}{z^{2}}\to\frac{1}{z^{2}}$ when $T\to0$. While integrating over $0\leq|z|\leq\infty$, the integrand is approximately given by $\frac{1}{e^{\beta z}+1}\frac{1}{z^{2}}\to\frac{1}{e^{\beta z}}\frac{1}{z^{2}}$. It decays even faster than the previous one and thus the integral over this range is also zero. The second part of this integral is also zero as we take the radius of the semicircle to infinity. From the second approach, we have $\mathcal{S}_{F}^{>}=0$.

Combing these, we have

\begin{align}
\left(-\frac{g^{2}}{\zeta\beta}\right) \sum_{p_{m}>0}\frac{1}{ip_{m}+\left(\varepsilon_{\bm{p}}-\bar{\mu}_{\xi}\right)-iD\pi}\times \frac{1}{ip_{m}-i\omega_{n}-\lambda} 	= \left(\frac{g^{2}}{\zeta}\right)\frac{n_{F}(-\varepsilon_{\bm{p}}+\bar{\mu}_{\xi}+iD\pi)}{i\omega_{n}+\varepsilon_{\bm{p}}-\bar{\mu}_{\xi}+\lambda-iD\pi}.
\end{align}

Thus, $\Sigma_{f}^{g,>}(\bm{k},\omega_{n}>0)$ reads

\begin{align}
\Sigma_{f}^{g,>}(\bm{k},\omega_{n}>0)
& =-\frac{g^{2}}{\zeta\beta}\sum_{p_{m}>0}\sum_{\bm{p}}\frac{1}{ip_{m}+\left(\varepsilon_{\bm{p}}-\bar{\mu}_{\xi}\right)+iD\pi}
\frac{1}{ip_{m}-i\omega_{n}-\lambda} \nn
& =\frac{g^{2}}{\zeta}\sum_{\bm{p}}\frac{n_{F}(-\varepsilon_{\bm{p}}+\bar{\mu}_{\xi}+iD\pi)}{i\omega_{n}+\varepsilon_{\bm{p}}-\bar{\mu}_{\xi}+\lambda-iD\pi}\nn
& =\frac{g^{2}\rho_{0}}{\zeta}\int_{-D+\bar{\mu}_{\xi}}^{D+\bar{\mu}_{\xi}}d\varepsilon^{\prime}\frac{n_{F}(\varepsilon^{\prime}+iD\pi)}{i\omega_{n}-\varepsilon^{\prime}+\lambda-iD\pi},
\end{align}

where we have make the following change  $-\varepsilon+\bar{\mu}_{\xi}=\varepsilon^{\prime}$. By absorbing $\lambda$ into $\omega$ and performing the analytic continuation, namely $i\omega_{n}\to\omega-\lambda+i\epsilon$ with $\epsilon \to 0^+ > 0 $, we arrive $(-D+\bar{\mu}_{\xi}<0)$

\begin{align}
\Sigma_{f}^{g,>}(\omega>0)	=\frac{g^{2}\rho_{0}}{\zeta}\int_{-D+\bar{\mu}_{\xi}}^{D+\bar{\mu}_{\xi}}d\varepsilon^{\prime}\frac{n_{F}(\varepsilon^{\prime}+iD\pi)}{\omega-\varepsilon^{\prime}-iD\pi}.
\end{align}

For $T\to0$, the Fermi function exhibits the following asymptotic result,

\begin{align}
n_{F}(\varepsilon^{\prime}+iD\pi)	=\begin{cases}
1, & \text{for }-D^{\prime}\leq\varepsilon\leq0\\
0, & \text{for }-D+\bar{\mu}_{\xi}\leq\varepsilon\leq D+\bar{\mu}_{\xi}.
\end{cases}
\end{align}

This gives rise to

\begin{align}
\left[\Sigma_{f}^{g,>}(\omega>0)\right]^{\prime\prime}  = & \frac{g^{2}D\pi\rho_{0}}{\zeta}\int_{-D+\bar{\mu}_{\xi}}^{0}\frac{d\varepsilon^{\prime}}{\left(\omega-\varepsilon^{\prime}\right)^{2}+(D\pi)^{2}} \nn
= & \left(-\frac{g^{2}\rho_{0}}{\zeta}\right)\left[\tan^{-1}\left(\frac{\omega}{\pi D}\right)-\tan^{-1}\left(\frac{D-\bar{\mu}_{\xi}+\omega}{\pi D}\right)\right] \nn
\approx & \left(-\frac{1}{\rho_{0}}\right)\left[\frac{\omega}{\pi D}-\tan^{-1}\left(\frac{1}{\pi}\right)\right],
\end{align}

where the assumption $D\gg\omega,\, \bar{\mu}_{\xi} \, \lambda$ has been made in the last line of the above equation. Finally, we reach

\begin{align}
\left[\Sigma_{f}^{g,>}(\omega>0)\right]^{\prime\prime}=\alpha-\varsigma\omega.
\end{align}

with $\alpha=D\tan^{-1}\left(\frac{1}{\pi}\right)\approx D/3$ (where $\tan^{-1}\left(\frac{1}{\pi}\right) \approx 1/3$) and $\varsigma=\frac{1}{\pi}$.
Next, we turn to the $p_{m}<0$ branch of $\Sigma_f$. Following a similar approach, we need to analyze the following contour integral by the residue theorem,

\begin{align}
\mathcal{S}_{F}^{<}&=-\frac{g^{2}}{\zeta}\oint_{C_{<}}\frac{dz}{2\pi i}\,\frac{n_{F}(z)}{z+\left(\varepsilon_{\bm{p}}-\bar{\mu}_{\xi}\right)+iD\pi}\cdot\frac{1}{z-i\omega_{n}-\lambda},
\end{align}

where $C_{<}$ encircles the lower hemisphere where $p_{m}<0$ of the complex $z$-plane counterclockwise. Below, we skip all the immediate steps and directly jump to the results. The final result of $\Sigma_{f}^{g,<}$ can be simply obtained by replacing $iD\pi$ with $-iD\pi$ in $\Sigma_{f}^{g,>}$. We finally have

\begin{align}
\Sigma_{f}^{g,<}(\omega_{n}<0)&=\frac{g^{2}\rho_{0}}{\zeta}\int_{-D+\bar{\mu}_{\xi}}^{D+\bar{\mu}_{\xi}}d\varepsilon^{\prime}\frac{n_{F}(\varepsilon^{\prime}-iD\pi)}{i\omega_{n}-\varepsilon^{\prime}+\lambda+iD\pi},
\end{align}

leading to

\begin{align}
\left[\Sigma_{f}^{<\prime\prime}(\omega<0)\right]^{\prime\prime}&=\left(-\frac{g^{2}D\pi\rho_{0}}{\zeta}\right)\int_{-D+\bar{\mu}_{\xi}}^{0}\frac{d\varepsilon^{\prime}}{\left(\omega-\varepsilon^{\prime}\right)^{2}+(D\pi)^{2}}.
\end{align}

We therefore conclude

\begin{align}
\left[\Sigma_{f}^{g,<}(\omega<0)\right]^{\prime\prime} =-\left(\alpha-\varsigma\omega\right).
\end{align}

In summary, combining the results of $\Sigma_{f}^{g,\gtrless}$ above, the imaginary part of $\Sigma_f^g$ reads

\begin{align}
\left[\Sigma_{f}^{g}(\omega)\right]^{\prime\prime}=-\text{sgn(\ensuremath{\omega})}\left(\alpha-\varsigma\omega\right).
\label{eq:sigma_f-img-scatteringrate}
\end{align}

Due to the constant term in $\left[\Sigma_{f}^{g}(\omega)\right]^{\prime\prime}$ of Eq. (<ref>), we note that, although our Kondo-Heisenberg-like lattice model exhibits translational symmetry, various fields that are local at the bare level still exhibit local characteristics even after perturbation corrections. This implies that the translational symmetry of this model is still broken, analogous to the single-impurity problem.

§.§ Scattering rate for the spinon-holon bound field $\xi$

In this section, we calculate the resistivity contributing from the $\xi$ field: The self-energy of the $\xi$ field is given by

\begin{align}
\Sigma_{\xi}(ik_{n})	& =\left(\frac{\bar{g}^{2}}{N_{s}\beta}\right)\sum_{p}G_{f}(p)\mathcal{G}_{b}(k-p) \nn
& =\left(\frac{\bar{g}^{2}}{\beta}\right)\sum_{ip_{m}}\left(\frac{1}{N_{s}}\sum_{\bm{p}}G_{f}(p)\right)\mathcal{G}_{b}(k-p).
\label{eq:Sigma_xi}
\end{align}

Similarly, we also need to introduce an additional minus sign into $\Sigma_{\xi}$ here, as we did during the calculation of $\Sigma_{\xi}$ in the RG analysis, as depicted in Eq. (<ref>).  Here, $\Sigma_{\xi}$ can be justified to be local by iteration due to the local $b$ boson.

To evaluate $\Sigma_{\xi}$, we first need to evaluate the dressed Green's function of the $f$ spinon, defined as $G_{f}(ip_{m}) = N^{-1}_s \sum_{\bm{p}} G_f (p)$. This can be explicitly computed as

\begin{align}
G_{f}(p_{m})	 = & \frac{1}{N_{s}}\sum_{\bm{p}}\frac{1}{ip_{m}-\varepsilon_{\bm{p}}-\Sigma_{f}(p_{m})} \nn
= & \left(-\rho_{0}\right)\ln\frac{D-ip_{m}+\Sigma_{f}(p_{m})}{-D-ip_{m}+\Sigma_{f}(p_{m})} \nn
= & \left(-\rho_{0}\right)\ln\frac{D+\Sigma_{f}^{\prime}(p_{m})-ip_{m}+i\Sigma_{f}^{\prime\prime}(p_{m})}{-D+\Sigma_{f}^{\prime}(p_{m})-ip_{m}+i\Sigma_{f}^{\prime\prime}(p_{m})}\nn
= & \left(-\rho_{0}\right)\ln\left[\frac{\sqrt{\left(1+\Sigma_{f}^{\prime}/D\right)^{2}+\left(p_{m}-\Sigma_{f}^{\prime\prime}\right)^{2}/D^{2}}}{\sqrt{\left(1-\Sigma_{f}^{\prime}/D\right)^{2}+\left(p_{m}-\Sigma_{f}^{\prime\prime}\right)^{2}/D^{2}}}e^{i(\theta_{1}-\theta_{2})}\right] \nn
= & \left(-\rho_{0}\right)\ln\frac{\sqrt{\left(1+\Sigma_{f}^{\prime}/D\right)^{2}+\left(p_{m}-\Sigma_{f}^{\prime\prime}\right)^{2}/D^{2}}}{\sqrt{\left(1-\Sigma_{f}^{\prime}/D\right)^{2}+\left(p_{m}-\Sigma_{f}^{\prime\prime}\right)^{2}/D^{2}}} -i\rho_{0}(\theta_{1}-\theta_{2})
\end{align}

\begin{align}
&\theta_{1}=\tan^{-1}\left(\frac{-(p_{m}-\Sigma_{f}^{\prime\prime})}{D+\Sigma_{f}^{\prime}(p_{m})}\right), \nn
\end{align}

where $\Sigma_{f}=\Sigma_{f}^{J}+\Sigma_{f}^{g}$ can be justified to be local. In the wide-band limit where $D\gg|p_{m}|,\,\Sigma_{f}^{\prime},\,\Sigma_{f}^{\prime\prime}$,  we have

\begin{align}
& \theta_{1}	\approx\tan^{-1}\left(\frac{-(p_{m}-\Sigma_{f}^{\prime\prime})}{D}\right)=\begin{cases}
2\pi, & p_{m}-\Sigma_{f}^{\prime\prime}>0\\
0 & p_{m}-\Sigma_{f}^{\prime\prime}<0
\end{cases}, \nn
\quad
&  \theta_{2}	\approx\tan^{-1}\left(\frac{-(p_{m}-\Sigma_{f}^{\prime\prime})}{-D}\right)=\begin{cases}
\pi, & p_{m}-\Sigma_{f}^{\prime\prime}>0,\\
\pi & p_{m}-\Sigma_{f}^{\prime\prime}<0.
\end{cases}
\end{align}

Thus, in the wide-band limit, we have $\theta_{1}-\theta_{2}=\pi\text{sgn}(p_{m}-\Sigma_{f}^{\prime\prime})$. Now we deal with the leading order approximation of $G_{f}$ in terms of small $\frac{(p_{m}-\Sigma_{f}^{\prime\prime})}{D}$.