Solution (source code)

= Solution

Let $\varepsilon=\hbar\omega$, use $M$ time slices and $q=e^{-\beta\varepsilon/M}$. The exact normal-ordered coherent Gaussian has cyclic matrix $A=1-qP$, where $P$ shifts one slice. Its eigenvalues are $1-qe^{2\pi in/M}$ and
$$
\det A=1-q^M=1-e^{-\beta\varepsilon},\qquad Z_{\rm normal}=\frac1{1-e^{-\beta\varepsilon}}.
$$
Restoring the constant energy gives
$$
\boxed{Z=e^{-\beta\varepsilon/2}Z_{\rm normal}=\frac1{2\sinh(\beta\varepsilon/2)},}
$$
consistent with the coordinate path integral. Alternatively, pairing positive and negative <Matsubara frequencies> symmetrically yields a determinant proportional to $x\prod_{n\ge1}[1+(x/2\pi n)^2]=2\sinh(x/2)$. This symmetric prescription already incorporates the zero-point convention and must not receive the vacuum factor twice. A naive unregulated frequency product can hide this difference.