condensed matter - Physical Meaning of the Gutzwiller Constraints

I have a doubt on the constraints for the expecation values obtained by Bünemann et all.

First i want to introduce my notation

To analytically solve a tight-binding model, $$\begin{equation} \hat{H}= \sum_{i,j,\sigma,\sigma'} t_{i,j}^{\sigma,\sigma'} \hat{c}^{\dagger}_{i;\sigma} \hat{c}_{j;\sigma'} + \sum_{i} \hat{H}_{i,at} \label{eq:Hubbard}\tag{1} \end{equation}$$ we can use a variational method made by Gutzwiller 1 , 2. He introduces this wave function, in the case of a single band system $$\begin{equation} \left|\varPsi_G\right>=\prod_{i}\left[1 + \left( g-1\right) \hat{D}_i\right]\left|\varPhi_0\right> \label{eq:GW1b_fin} \tag{2} \end{equation}$$ where $g$ is a real number between $0$ and $1$ that plays the role of a variational parameter, $\left|\varPhi_0\right>$ is a Slater determinat on which we can apply Wick's theorem and $\hat{D}_i = \hat{n}_{i,\uparrow}\hat{n}_{i,\downarrow}$ is the double occupation operator.

Here Bunemann et all. extend $\eqref{eq:GW1b_fin}$ to a multiband system, and it becomes $$\begin{equation} \left|\varPsi_G\right>=\prod_{i}\left[1 + \sum_{\Gamma} \left( \lambda_{i,\Gamma}-1\right) \hat{m}_{i,\Gamma}\right]\left|\varPhi_0\right> \label{eq:P_GW_multiband} \end{equation}$$ where $$\begin{equation} \hat{m}_{i,\Gamma} = \left|\Gamma\right>_{i} \left<\Gamma\right|_{i} = \prod_{\sigma \in \Gamma} \hat{n}_{i,\sigma} \prod_{\sigma \in \overline{\Gamma}} \left(1-\hat{n}_{i,\sigma}\right) \end{equation}$$ while $\lambda_{i,\Gamma}$ plays now the role of variational parameters. $\Gamma$ is the atomic eignestate of the atomic Hamiltonian.

They realize that in the limit of infinite coordination number, average values on the Gutzwiller wavefunction can be computed exactly if the following constraints are satised $$\begin{align} \left<\varPsi_{0}\left|\hat{P}^{\dagger}_{i}\hat{P}_{i}\right|\varPsi_{0}\right>&=1 \label{a}\tag{3}\\ \left<\varPsi_{0}\left|\hat{P}^{\dagger}_{i}\hat{P}_{i}\hat{c}^{\dagger}_{i,\alpha}\hat{c}_{i,\beta}\right|\varPsi_{0}\right>&=\left<\varPsi_{0}\left|\hat{c}^{\dagger}_{i,\alpha}\hat{c}_{i,\beta}\right|\varPsi_{0}\right> \label{b}\tag{4} \end{align}$$

It is not clear the physical meaning of these constraints, given also Here by Metzner.

Mathematically, applying Wick's theorem to $\eqref{b}$

$$\begin{align*} \left<\varPsi_{0}\left|\hat{P}^{\dagger}_{i}\hat{P}_{i}\hat{c}^{\dagger}_{i,\alpha}\hat{c}_{i,\beta}\right|\varPsi_{0}\right> &= \left<\varPsi_{0}\left|\hat{P}^{\dagger}_{i}\hat{P}_{i}\right|\varPsi_{0}\right> \left<\varPsi_{0}\left|\hat{c}^{\dagger}_{i,\alpha}\hat{c}_{i,\beta}\right|\varPsi_{0}\right> \\ +& {\left<\varPsi_{0}\left|\hat{P}^{\dagger}_{i}\hat{P}_{i}\hat{c}^{\dagger}_{i,\alpha}\hat{c}_{i,\beta}\right|\varPsi_{0}\right>_{contractions}} \end{align*}$$ assuming that $\varPsi_{0}$ is normalized, we have that the sum of all Wick's contractions of $\hat{c}^{\dagger}_{i,\alpha}\hat{c}_{i,\beta}$ with $\hat{P}^{\dagger}_{i}\hat{P}_{i}$ vanishes