Edit: It turns out that problem encountered in this question is not limited to BdG Hamiltonians.
I am having trouble in using the coherent state path integral approach to calculate the free energy. For example, consider a Bogoliubov-de Gennes (BdG) Hamiltonian:
$$ \mathsf{H}[b^\dagger,b] = \Delta(b^\dagger_1 b^\dagger_2 + b_2 b_1) + \lambda (b^\dagger_1 b_1 + b^\dagger_2 b_2) \quad (\lambda > \Delta > 0) $$
involving two bosons $\{b_i\}_{i=1}^2$. After diagonalizing the Hamiltonian with a Bogoliubov transformation, I can show that the free energy (up to constants independent of parameters $\Delta, \lambda$ in $\mathsf{H}$) is
$$ F = \frac{2}{\beta} \ln(1 - e^{-\beta E}) + E - \lambda, \quad E \equiv \sqrt{\lambda^2 - \Delta^2} $$
Now I describe my attempt to use the coherent state path integral. Since $\mathsf{H}$ is already normal-ordered ($b^\dagger$'s in front of $b$'s), writing down the coherent state path integral representation of the partition function is straightforward:
$$ \begin{align*} Z &= \int D[\bar{b},b] e^{-S[\bar{b},b]} \\ S[\bar{b},b] &= \int_0^\beta d\tau \, \Big[ \bar{b}_1 \partial_\tau b_1 + \bar{b}_2 \partial_\tau b_2 + \mathsf{H}[\bar{b},b] \Big] \end{align*} $$
Note that in this representation $b(\tau), \bar{b}(\tau)$ are commuting complex numbers. The integration measure $D[\bar{b},b]$ is (up to a normalization constant)
$$ D[\bar{b},b] = \prod_{0 \le \tau \le \beta} \prod_{a=1}^2 d\bar{b}_a(\tau) db_a(\tau) $$
Define the bosonic "spinor"
$$ \Phi = \begin{bmatrix} b_1 \\ \bar{b}_2 \end{bmatrix} \ \Rightarrow \ \bar{\Phi} = (\bar{b}_1, b_2) $$
I write $S[\bar{b},b]$ as
$$ \require{cancel} \begin{align*} S[\bar{b},b] &= \int_0^\beta d\tau \, \Big[ \bar{b}_1 \partial_\tau b_1 - b_2 \partial_\tau \bar{b}_2 + {\partial_\tau(\bar{b}_2 b_2)} \\ &\qquad + \Delta(\bar{b}_2 \bar{b}_1 + b_1 b_2) + \lambda (\bar{b}_1 b_1 + \bar{b}_2 b_2) \Big] \\ &= \int_0^\beta d\tau \, \bar{\Phi}(\tau) G^{-1}(\tau) \Phi(\tau) \end{align*} $$
Here I dropped total $\tau$-derivatives in $S$, and introduced the Green's function matrix
$$ G^{-1}(\tau) = \begin{bmatrix} \partial_\tau + \lambda & \Delta \\ \Delta & -\partial_\tau + \lambda \end{bmatrix} $$
The free energy is found by integrating over $b$. To deal with $\partial_\tau$, I Fourier transform to frequency space:
$$ b_a(\tau) = \frac{1}{\sqrt{\beta}} \sum_\omega b_{a\omega} e^{i\omega\tau}, \quad \bar{b}_a(\tau) = \frac{1}{\sqrt{\beta}} \sum_\omega \bar{b}_{a\omega} e^{-i\omega\tau} $$
where $\omega$ sums over boson Matsubara frequencies $\{2\pi n/\beta \ |\ n \in \mathbb{Z}\}$. This transformation is unitary, so I can change the integration measure to
$$ D[\bar{b},b] = \prod_{\omega} \prod_{a=1}^2 d\bar{b}_{a\omega} db_{a\omega} $$
The Fourier transform of $\Phi$ is
$$ \Phi = \begin{bmatrix} b_1 \\ \bar{b}_2 \end{bmatrix} = \frac{1}{\sqrt{\beta}} \sum_\omega \Phi_\omega e^{i\omega \tau}, \quad \Phi_\omega = \begin{bmatrix} b_{1\omega} \\ \bar{b}_{2,-\omega} \end{bmatrix} $$
The action becomes
$$ \begin{align*} S[\bar{b},b] &= \frac{1}{\beta} \int_0^\beta d\tau \sum_{\omega,\omega'} \bar{\Phi}_\omega e^{-i\omega\tau} G^{-1}(\tau) e^{i\omega'\tau} \Phi_{\omega'} \\ &\equiv \sum_{\omega,\omega'} \bar{\Phi}_\omega G^{-1}_{\omega \omega'} \Phi_{\omega'} \end{align*} $$
where the frequency-space Green's function is
$$ \begin{align*} G^{-1}_{\omega \omega'} &\equiv \frac{1}{\beta} \int_0^\beta d\tau \, e^{-i\omega\tau} G^{-1}(\tau) e^{i\omega'\tau} \\ &= \frac{1}{\beta} \int_0^\beta d\tau \, e^{-i\omega\tau} \begin{bmatrix} i\omega' + \lambda & \Delta \\ \Delta & -i\omega' + \lambda \end{bmatrix} e^{i\omega'\tau} \\ &= \begin{bmatrix} i\omega + \lambda & \Delta \\ \Delta & -i\omega + \lambda \end{bmatrix} \delta_{\omega \omega'} \equiv \mathcal{G}^{-1}_\omega \delta_{\omega \omega'} \end{align*} $$
Then the partition function is
$$ \begin{align*} Z &= \int D[\bar{b},b] \exp\bigg\{ - \sum_\omega \bar{\Phi}_\omega \mathcal{G}^{-1}_\omega \Phi_\omega \bigg\} \\ &= \prod_\omega \int d\bar{\Phi}_\omega \, d\Phi_\omega \exp[ - \bar{\Phi}_\omega \mathcal{G}^{-1}_\omega \Phi_\omega ] \end{align*} $$
where $d\Phi_\omega = db_{1\omega} d\bar{b}_{2,-\omega}$ and $d\bar{\Phi}_\omega = d\bar{b}_{1\omega} db_{2,-\omega}$. The integration is Gaussian, and yields
$$ Z = \text{const} \times \prod_\omega (\det \mathcal{G}^{-1}_\omega)^{-1} $$
Here and thereafter $(\text{const})$ is a constant number independent of parameters $\mu, \Delta$ in the Hamiltonian. In this formula the "const" comes from the normalization of the integration measure. Then the free energy is (define $E = \sqrt{\lambda^2 - \Delta^2}$)
$$ \begin{align*} F &= -\frac{1}{\beta} \ln Z = \frac{1}{\beta} \sum_\omega \ln \det \mathcal{G}^{-1}_\omega + \text{const} \\ &= \frac{1}{\beta} \sum_\omega \ln(\omega^2 + E^2) + \text{const} \\ &= \frac{1}{\beta} \sum_\omega \ln[\beta^2(\omega^2 + E^2)] + \text{const} \\ &= \frac{2}{\beta} \ln(1 - e^{-\beta E}) + E + \text{const} \end{align*} $$
The last equality can be found in, say, eq. (25.34) in Quantum Field Theory for the Gifted Amateurs. But the constant is independent of $\lambda$, contradicting the results from Bogoliubov transformation. This shift by $\lambda$ may look harmless, but it will cause problems if $\mathsf{H}$ is obtained by a mean field approximation, and $\Delta, \lambda$ are to be determined by minimizing the free energy. So it is important to get the $\Delta, \lambda$ dependence right in $F$.
What is the mistake I made in the derivation?