Solution (source code)

= Solution

Define the Grassmann spin density $\mathbf j_n=\bar\psi_n\boldsymbol\sigma\psi_n$. The <Pauli matrices> and Grassmann reordering give
$$
\mathbf j_n^2=-6D_n,\qquad 3UD_n=-\frac U2\mathbf j_n^2.
$$
For example $j_z^2=-2D_n$; each of the other two components gives the same result. For $U>0$, complete the square independently at each site and time slice:
$$
 e^{\Delta\tau U\mathbf j^2/2}=\left(\frac{\Delta\tau}{2\pi U}\right)^{3/2}\int_{\mathbb R^3}d^3M\,\exp\left[-\Delta\tau\left(\frac{\mathbf M^2}{2U}-\mathbf M\cdot\mathbf j\right)\right].
$$
This normalized <spin-vector Hubbard–Stratonovich decoupling> yields
$$
\boxed{S=S_0+\int_0^\beta d\tau\sum_n\left[\frac{\mathbf M_n^2}{2U}-\sum_{\alpha\beta}\bar\psi_{n\alpha}\mathbf M_n\cdot\boldsymbol\sigma_{\alpha\beta}\psi_{n\beta}\right].}
$$
The auxiliary <bosonic> <quantum field> is real and periodic, while the fermionic fields remain antiperiodic. All three vector components are required to reproduce the original $3U$ interaction; replacing this exact transformation by a single-component identity would change its coefficient.