I have a question about a quantum field theory (QFT) partition function in the following example. Consider a single-site Hubbard model with the Hamiltonian given by:
$$H = \varepsilon_I \left(a_\uparrow^\dagger a_\uparrow + a_\downarrow^\dagger a_\downarrow \right) + U a_\uparrow^\dagger a_\uparrow a_\downarrow^\dagger a_\downarrow.$$
Typically, if we want to write down the partition function, we sum the Boltzmann factors over the Fock space. In this case, we have four states: one empty, two single-filled, and one double-filled. The partition function is:
$$\mathcal{Z} = 1 + 2e^{-\beta\varepsilon_I} + e^{-\beta(2\varepsilon_I + U)}.$$
Now, I would like to obtain the same using QFT.
The first step is writing down the action:
$$ S\left[\bar\phi,\phi,\bar\psi,\psi\right] = \int_0^\beta d\tau \left[\bar\psi(\tau)\left(\partial_\tau+\varepsilon\right)\psi(\tau) + \bar\phi(\tau)\left(\partial_\tau+\varepsilon\right)\phi(\tau) + U\bar\phi(\tau)\phi(\tau)\bar\psi(\tau)\psi(\tau)\right]\,, $$ where $\psi$ and $\phi$ correspond to spin-up and spin-down states.
Defining the Fourier transforms of the fields using Matsubara frequencies gives $$\psi(\tau) = \frac{1}{\sqrt{\beta}}\sum_{i}\psi_ie^{-i\omega_i\tau}\,,\quad \bar\psi(\tau) = \frac{1}{\sqrt{\beta}}\sum_{j}\bar\psi_je^{i\omega_j\tau}\,,\quad \phi(\tau) = \frac{1}{\sqrt{\beta}}\sum_{k}\phi_ke^{-i\omega_k\tau}\,, \quad \bar\phi(\tau) = \frac{1}{\sqrt{\beta}}\sum_{l}\bar\phi_le^{i\omega_l\tau}\,.$$
This results in $$ S\left[\bar\phi,\phi,\bar\psi,\psi\right] = \sum_{j}\bar\psi_j\left(-i\omega_j+\varepsilon\right)\psi_j + \sum_{l}\bar\phi_l \left(-i\omega_l+\varepsilon\right)\phi_l +\frac{U}{\beta}\sum_{l,k,j,i} \bar\phi_l\phi_k\bar\psi_j\psi_i\delta_{\omega_l+\omega_j,\omega_i+\omega_k} \\ = \sum_{j}\bar\psi_j\left(-i\omega_j+\varepsilon\right)\psi_j + \sum_{l}\bar\phi_l \left(-i\omega_l+\varepsilon\right)\phi_l +\frac{U}{\beta}\sum_m\left[\sum_{l}\bar\phi_l\phi_{l-m}\right]\left[\sum_j\bar\psi_j\psi_{j+m}\right] \\ = \sum_{j}\bar\psi_j\left(-i\omega_j+\varepsilon\right)\psi_j + \sum_{l}\bar\phi_l \left(-i\omega_l+\varepsilon\right)\phi_l +\frac{U}{\beta}\sum_{m}\Phi_{-m}\Psi_m \nonumber \\ = \sum_{j,n}\bar\psi_j\left[\delta_{j,n}\left(-i\omega_j+\varepsilon\right)+\frac{U}{\beta}\Phi_{j-n}\right]\psi_{n} + \sum_{l}\bar\phi_l \left(-i\omega_l+\varepsilon\right)\phi_l \,. $$
The partition function is
$$ \mathcal{Z} = \int D(\bar\psi,\psi)D(\bar\phi,\phi)e^{-S} = \int D(\bar\phi,\phi) \det\left[-i \hat\omega + \varepsilon + \frac{U}{\beta}\hat\Phi\right]\exp\left[-\sum_l\bar\phi_l \left(-i\omega_l+\varepsilon\right)\phi_l\right]\,. $$
We can put the determinant inside the exponential:
$$ \mathcal{Z} = \int D(\bar\phi,\phi) \exp\left[-\sum_l\bar\phi_l \left(-i\omega_l+\varepsilon\right)\phi_l + \mathrm{tr} \ln \left(-i \hat\omega + \varepsilon + \frac{U}{\beta}\hat\Phi\right)\right]\,. $$
This expression looks similar to what one gets after using the Hubbard-Stratonovich transformation to get rid of a quartic field. The difference is that here, the $\phi$ field is fermionic, not bosonic.
I am not sure how to proceed from here to obtain the exact solution for the partition function. Generally, one would try to look for a saddle point and Taylor expand. However, in this case, there is a known solution.
Any suggestions are greatly appreciated.
