Lagrange multipliers in quantum systems are usually (always?) implemented on the level of expectation values - often specifically ground state expectation values. That's also the case for these two papers. So, given an extended Hamiltonian $$ H'= H_\mathrm{MF}+H_\mathrm{constraints}$$ where $H_\mathrm{constraints}$ contains the Lagrange-multiplier-enforced constraints, we want to find a ground state of $|\psi\rangle$ of $H'$ that satisfies the constraints. For reasonable constraints this will be possible by finding appropriate values for the Lagrange multipliers. Note that $|\psi\rangle=|\psi_0\rangle+|\delta\psi\rangle$, where $|\psi_0\rangle$ would be the ground state of $H_\mathrm{MF}$. That is, ground state expectation values $\langle\psi|\hat{O}|\psi\rangle$ of some operator $\hat{O}$ will depend on the Lagrange multipliers $a_i$, and we can tune the $a_i$ until we find expectation values compatible with the desired constraints.
Of course, we could use the same approach for excited states, but then we i) would expect different values for the Lagrange multipliers, and ii) might require additional constraints to make sure the state is orthogonal to lower-energy states. In other words, <PERSON>'s Lagrange multiplier terms are enough to enforce the constraints on the ground state, but you shouldn't expect them to be enough to enforce the constraint on the full Hilbert space. Compare the situation with the most familiar example of a Lagrange multiplier in stat. mech. or condensed matter systems - the chemical potential $$\mu \left( \int |\Psi \left(\mathbf{r}\right)|^2 d\mathbf{r} - N \right). $$ It can be used to enforce a fixed particle number $N=<\hat{N}>$ in a thermal state, but excited states are allowed to have other occupations.
Additional comments:
As you note, <PERSON>'s constraints are site-dependent. This is often rather inconvenient, particularly in numerical calculations. 976
Hence, often when searching for a translationally invariant mean field solution, the site-dependent local constraints can be relaxed to global ones, i.e. you might go from $\sum_i \lambda_i \left( f_{i\sigma}^\dagger f_{i\sigma} - 1\right)$ to $\lambda \left(\sum_i f_{i\sigma}^\dagger f_{i\sigma} - N\right)$, which now enforces half-filling on average. This is commonly done in the $t-J$ model theory for high-$T_c$ superconductivity, for example, see e.g. this review. Ideally, using the relaxed constraint can allow an easier path to a solution that satisfies the local constraints.
It might also be instructive to look at how one might solve such a constrained mean-field problem numerically and self-consistently. One algorithm would be
1. Start by making initial guesses for the values of the mean field parameters (e.g. $\chi$)
2. Find values for the Lagrange multipliers such that the constraint is satisfied
3. Calculate new expectation values ($\chi_{current}$)
4. Update guesses for the mean field parameters (e.g. by conservative mixing $\chi_{new}=(1-x)\chi_{old} + x\chi_{current}$ where $x$ is small)
5. Return to step 2, and repeat until the parameters converge. 915