Solution (source code)

= Solution

The pair products are $PQ=N_aN_b$ and $QP=(N_a+1)(N_b+1)$. Hence
$$
K_1^2+K_2^2=\frac12(PQ+QP)
=N_aN_b+\frac12(N_a+N_b+1),
$$
whereas $K_3=(N_a+N_b+1)/2$. Subtraction gives
$$
\boxed{C=\frac{(N_a-N_b)^2-1}{4}}.
$$
The number difference commutes with pair creation, pair annihilation and $K_3$, since each pair operation changes both occupations equally. Its <boson-pair Casimir operator> $C$ therefore commutes with every generator. In a fixed-difference sector $N_a-N_b=d$, the Casimir <eigenvalue> is $(d^2-1)/4$, including $-1/4$ for equal occupations.