Lagrange multiplier in spin liquid mean-field theory (Paper by X.G. Wen) My question is about a step in this paper: PhysRevB.65.165113 (X.G. Wen) page 6
Or alternatively: PhysRevB.90.174417 page 3. All the papers concerning spin liquids and the projective symmetry group derive the mean field Hamiltonian in the same way. 
Background: The Hilbert space $(\mathbb C^2)^{\otimes d}$ of $d$ spin-$1/2$ particles on a lattice is embedded in a larger fermionic Fock space with creation operators $f_{i\uparrow}^\dagger, f_{i\downarrow}^\dagger$ for each site $i$. The Heisenberg model Hamiltonian 
\begin{align}
H = \sum_{ij} J_{ij} \, \vec S_i \cdot \vec S_j
\end{align}
is written in terms of the fermionic operators and a mean field decoupling is performed yielding
\begin{align}
H_{MF} = \sum_{ij} (\Psi_i^\dagger U_{ij} \Psi_j + \text{h.c.}) + \text{const.} , \qquad \text{where }\Psi_i = \begin{bmatrix}f_{i\uparrow}^\dagger &f_{i\downarrow} \end{bmatrix}^t 
\end{align}
This is the part that I understand. But then an additional Lagrange multiplier term is added to the Hamiltonian. The Fock space generated by $f_{i\alpha}$ is  bigger than $(\mathbb C^2)^{\otimes d}$. Only those states are "physical", i.e. correspond to states in $(\mathbb C^2)^{\otimes d}$, which satisfy a one-particle per site constraint:
\begin{align}
 \langle \psi|f_{i\uparrow}^\dagger f_{i\uparrow} + f_{i\downarrow}^\dagger f_{i\downarrow} |\psi \rangle = 1 
\end{align}
Wen says in his paper that such a constraint can be forced by adding the (site dependent?) Lagrange multiplier term
\begin{align}
+ \sum_{i} a_3 (f_{i\uparrow}^\dagger f_{i\uparrow} + f_{i\downarrow}^\dagger f_{i\downarrow} -1) + \left[(a_1 + i a_2) f_{i\uparrow}f_{i\downarrow} + \text{h.c.}\right]
\end{align}
to the Hamiltonian. $a_\mu$ are real numbers. The second paper furthermore tells that $a_\mu$ can be obtained by the condition 
\begin{align}
\partial E_g/ \partial a_\mu = 0
\end{align} 
where $E_g$ is the ground state energy of the mean field Hamiltonian. 
Question: How can such an additional term enforce a constraint on the eigenstates of the Hamiltonian? 
I see that the first term corresponds to the condition that each site is single occupied, and that the other two terms correspond to the conditions that no site is not and that no site is doubly occupied.
They look quite similar as Lagrange multiplier terms in multivariable calculus. I understand how Lagrange multipliers work there. But I don't understand how they work in quantum mechanics. Any help is highly appreciated.
Links to the two papers:


*

*https://journals.aps.org/prb/abstract/10.1103/PhysRevB.65.165113 

*https://journals.aps.org/prb/abstract/10.1103/PhysRevB.90.174417
arXiv: 


*

*https://arxiv.org/abs/cond-mat/0107071

*https://arxiv.org/abs/1407.4124
 A: 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, Wen'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, Wen's constraints are site-dependent. This is often rather inconvenient, particularly in numerical calculations. 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


*

*Start by making initial guesses for the values of the mean field parameters (e.g. $\chi$)

*Find values for the Lagrange multipliers such that the constraint is satisfied

*Calculate new expectation values ($\chi_{current}$)

*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)

*Return to step 2, and repeat until the parameters converge.

A: I have a different perspective that what Anyon said, which is interesting.
I regard this as a classic trick, going back at least to Polyakov if not earlier. I learned about it from E. Fradkin's book Field Theories in Condensed Matter Physics, which describes how you can use it to understand anti-ferromagnets, Mott insulators, and the like as gauge theories.
Consider the operators
$$Q(\theta_i) = \prod_i e^{i \theta (f_{i\uparrow}^\dagger f_{i\uparrow} + f_{i\downarrow}^\dagger f_{i\downarrow} - 1)},$$
which depend on a parameter $\theta_i \in \mathbb{R}/2\pi \mathbb{Z}$. The idea is that in a system with the local constraint
$$f_{i\uparrow}^\dagger f_{i\uparrow} + f_{i\downarrow}^\dagger f_{i\downarrow} - 1,$$
the $Q$ will generate gauge symmetries. Thus we can think about imposing the constraint as gauging a symmetry.
The idea is to use minimal coupling of the form
$$L = L_{\rm matter} + ja + L_{\rm gauge}$$
where $L_{\rm matter}$ is the Lagrangian of the unconstrained theory, $j$ is the current, $a$ is the gauge field, and $L_{\rm gauge}$ is some weakly-coupled gauge theory kinetic term. For the Hamiltonian we get the first that Wen writes (the other two terms are to preserve manifest $SU(2)$ symmetry, as Anyon says in the comments), where
$$j_i = f_{i\uparrow}^\dagger f_{i\uparrow} + f_{i\downarrow}^\dagger f_{i\downarrow} - 1,$$
and the $a_3$ is the time component of the gauge field. This component is important because it occurs without a canonical conjugate, and intuitively is free to fluctuate wildly. Thus, the energy spectrum will be bottomless unless the constraint is satisfied.
In field theory it's easier to see, since in the path integral, $a_3$ will occur only in a term
$$\exp{i \int a_3 (j - \nabla \cdot E)},$$
where $E_i = \partial_0 a_i$ is the electric field of $a$, which only depends on spatial components. Thus we see the Lagrangian is linear in $a_3$, and when we integrate it out, we obtain a delta function which imposes the constraint
$$j = \nabla \cdot E.$$
If we take the gauge field to be a background, we can set $E = 0$ and we obtain the constraint.
A neat trick you can now do is consider integrating out the matter to obtain an effective action for this gauge field and then letting it become dynamical. This often gives a dual formulation of the problem being studied in terms of a pure gauge theory. Every duality I know besides holography is of this form, actually, and it's the reason one expects to see degrees of freedom like gauge fields appearing in spin liquids.
