\equiv c^\dagger_{i\alpha} c^{\;}_{i\alpha}$ is the electron number operator.

The quasi-particle density matrix is given by the Boltzmann distribution $\rho^{\rm qp} = e^{-\beta \mathcal{H}^{\rm qp}} / Z^{\rm qp}$, where $\mathcal{H}^{\rm qp}$ is the renormalized tight-binding Hamiltonian $$\begin{aligned}
    \label{eq:H_renorm}
    \mathcal{H}^{\rm qp} &=&  \sum_{\langle ij \rangle, \alpha} \left( \mathcal{R}^{\;}_{i,\alpha} \mathcal{R}^*_{j,\alpha}\, t_{ij} c^\dagger_{i,\alpha} c^{\;}_{j, \alpha} + \mbox{h.c.} \right) \nonumber \\
     & & \qquad - \sum_{i, \alpha} \mu_{i, \alpha}\, c^\dagger_{i, \alpha} c^{\;}_{i, \alpha}, \end{aligned}$$ where $t_{ij}$ depends on the displacements through Eq. (\[eq:t\_ij\]), and the renormalization factors are [@lanata12] $$\begin{aligned}
    \label{eq:R}
    \mathcal{R}_{i, \alpha} = \frac{ {\rm Tr}\left(\Phi^\dagger_i c^\dagger_{i, \alpha} \Phi^{\;}_i c^{\;}_{i, \alpha} \right) }{  \sqrt{n_{i, \alpha} (1 - n_{i, \alpha})} }.\end{aligned}$$ and $n_{i,\alpha} = {\rm Tr}(\Phi^\dagger_i \, \hat{n}_{i, \alpha} \Phi^{\;}_i)$ is the quasi-particle number at the $i$-th atom.

The minimization of $\mathcal{F}$ with respect to $\Phi_i$, subject to constraint ${\rm Tr}(\Phi^\dagger_i \Phi^{\;}_i) = 1$ can be recast into a nonlinear eigenvalue problem $\mathcal{H}^{\rm SB}_i \Phi_i = \epsilon_i \Phi_i$ for each atom [@lanata12]. Since $\Phi_i$ is diagonal in the case of single-band Hubbard model, we use the same notation $\Phi_i = (e_i, p_{i, \uparrow}, p_{i, \downarrow}, d_i)$ to denote the wavefunction or the coherent state of the slave bosons. The embedded Hamiltonian has the form $$\begin{aligned}
    