Solution (source code)

= Solution

Take $L\ge2$, so neighbouring operators act at distinct sites. The identity
$$
S_n^xS_{n+1}^x+S_n^yS_{n+1}^y
=\tfrac12(S_n^+S_{n+1}^-+S_n^-S_{n+1}^+)
$$
expresses the <Hamiltonian operator> of this <quantum XY chain> in <spin ladder operators>. Substituting the <Schwinger boson representation> gives the requested quartic expression:
$$
\boxed{H=-\frac J2\sum_n
\left(a_n^\dagger b_n b_{n+1}^\dagger a_{n+1}
+b_n^\dagger a_n a_{n+1}^\dagger b_{n+1}\right).}
$$
Operators at different sites commute, allowing <normal ordering> without a correction here. For the first quartic term, the <inverse discrete Fourier transforms> give
$$
\begin{aligned}
\sum_n a_n^\dagger b_{n+1}^\dagger b_n a_{n+1}
&=\frac1{L^2}\sum_{n,p,q,r,t}
e^{i(p+q-r-t)n}e^{i(q-t)}A_p^\dagger B_q^\dagger B_r A_t\\
&=\frac1L\sum_{p,q,r,t}
\delta^{(L)}_{p+q,r+t}e^{i(q-t)}A_p^\dagger B_q^\dagger B_r A_t.
\end{aligned}
$$
Here $\delta^{(L)}$ is one when its two <momenta> agree modulo $2\pi$ and zero otherwise; it expresses lattice <momentum conservation>. The other quartic term has the same ordered operators with phase $e^{i(p-r)}$. On the support of this <Kronecker delta>, $p-r=t-q$ modulo $2\pi$, so the two phases sum to $2\cos(q-t)$. Thus the <Fourier vertex of a Schwinger boson chain> is
$$
\boxed{H=-\frac JL\sum_{p,q,r,t}
\delta^{(L)}_{p+q,r+t}\cos(q-t)\,
A_p^\dagger B_q^\dagger B_r A_t.}
$$
An equivalent three-sum expression sets $t=p+q-r$ modulo $2\pi$, so its coefficient is $\cos(r-p)$. The factor $1/L$ comes from four factors $L^{-1/2}$ and one site sum, while the factor of two from the two exchange terms cancels the original $1/2$. The formula is <Hermitian>: conjugation exchanges $p$ with $t$ and $q$ with $r$, and the coefficient is unchanged after imposing lattice <momentum conservation>. It also preserves each local total <boson> number in its real-space form, hence preserves the physical constrained subspace.