Jane, $\partial_\tau$ is clearly a derivative with respect to a bosonic time $\tau$, so it commutes with everything else (except for functions of $\tau$ itself, with which it has a nonzero commutator), rather than anticommutes. Only if both objects have a fermionic character (if both of them are Grassmann-odd), they anticommute with one another (or they have an anticommutator that can be evaluated).
There is no sign error in the formulae, however. You asked a good question: how do you get from
$$-\partial_\tau \bar\psi \psi$$
to
$$+\bar\psi\partial_\tau \psi$$
One needs a bit of patience to answer this question. Note that in the two expressions, a different variable is differentiated. In the first one, it's $\bar\psi$ that is differentiated; in the second one, it's $\psi$.
You can't just move derivatives around. Even for bosonic functions, $uv'$ is something else than $u'v$, isn't it?
So the two expressions are not "obviously equal", not even up to a sign, and to convert one to the other, you must carefully integrate by parts. Note that
$$\partial_\tau (\bar\psi \psi) = \partial_\tau\bar\psi \psi + \bar\psi\partial_\tau\psi. $$
This "Leibniz rule" proceeded just like for the derivative of products of bosonic factors because I had to bubble $\partial_\tau$ through the $\psi$'s, and $\partial_\tau$ is a bosonic object. If I were writing down a Leibniz rule for a Grassmannian derivative, I would have to change the sign everytime the derivative would bubble through a Grassmann-odd factor.
But here we deal with bosonic $\tau$-derivatives so the Leibniz rule is just like it has always been. So it implies that up to a total derivative - namely the left-hand side $\partial_\tau (\bar\psi \psi)$ that integrates to zero over the periodic Euclidean time - the two terms on the right hand side are opposite to one another. That's where the minus sign came from.
I) The Euclidean path integral with $\hbar=1$ reads
$$ Z~=~\int_{DBC} \!{\cal D}x ~e^{-S},\tag{1}$$
with Dirichlet boundary conditions (DBC)
$$ x(0)~=~0~=~x(T).\tag{2}$$
We expand the real periodic variable $x\in\mathbb{R}$ in Fourier series$^1$
$$ \begin{align} x(t) ~=~&\frac{a_0}{2}+\sum_{n\in\mathbb{N}} \left\{a_n \cos(\omega_n t) + b_n \sin(\omega_n t)\right\}\cr
~=~& \sum_{n\in\mathbb{Z}} c_n e^{i\omega_n t}, \cr
\omega_n~:=~&\frac{2\pi n}{T},\qquad a_n,b_n~\in~ \mathbb{R},\cr
c_n~\in~& \mathbb{C}, \qquad c_n^{\ast}~=~c_{-n}. \end{align} \tag{3} $$
The DBC (2) becomes
$$ \sum_{n\in\mathbb{Z}} c_n~=~0\qquad\Leftrightarrow\qquad c_0 ~=~ -2{\rm Re}\sum_{n\in\mathbb{N}}c_n .\tag{4}$$
The action for a free non-relativistic point particle with mass $m=1$ reads:
$$ S~=~\frac{1}{2} \int_0^T \!dt~ \dot{x}^2~=~T\sum_{n\in\mathbb{N}}\omega_n^2 |c_n|^2 . \tag{5} $$
II) We know that the proper normalization of the path integral (1) is
$$ Z~=~\frac{1}{\sqrt{2 \pi T}}.\tag{6} $$
This can e.g. be deduced (without introducing fudge factors!) from the (semi)group property of Feynman path integrals, cf. this Phys.SE post and links therein. Up till now we have basically just restated what OP wrote in his question.
III) Now we would like to repeat the same calculation using Fourier series, i.e. work with the Matsubara frequencies. In this answer, we will not explore the (semi)group property, but just do a quick and dirty calculation using various fudge factors, and see what we get. Since this is homework, the explanation will be somewhat brief.
To make heuristic sense of the path integral (1), we will use the following zeta function regularization rules:
$$ \prod_{n\in\mathbb{N}} a ~=~\frac{1}{\sqrt{a}} \quad\text{and}\quad \prod_{n\in\mathbb{N}} n ~=~\sqrt{2\pi}, \tag{7}$$
stemming from the zeta function values
$$ \zeta(0)~=~-\frac{1}{2}
\quad\text{and}\quad
\zeta^{\prime}(0)~=~-\ln\sqrt{2\pi} , \tag{8} $$
respectively. Now let the path integral measure be
$$ \begin{align} {\cal D}x~:=~&\delta\left(B\sum_{n\in\mathbb{Z}} c_n\right) A\mathrm{d} c_0 \prod_{n\in\mathbb{N}} A^2 \mathrm{d}^2c_n \cr
~\stackrel{(7)}{=}~&\frac{1}{B}\delta\left(\sum_{n\in\mathbb{Z}} c_n\right) \mathrm{d} c_0 \prod_{n\in\mathbb{N}} \mathrm{d}^2c_n , \end{align}\tag{9} $$
where $A$, $B$ are fudge factors. Interestingly, eq. (9) is independent of the $A$-fudge factor! After performing the delta function integration and the Gaussian integrals, we find
$$ \begin{align} Z~=~&\frac{1}{B} \prod_{n\in\mathbb{N}} \frac{\pi}{T\omega_n^2}\cr
~=~&\frac{1}{B} \prod_{n\in\mathbb{N}} \frac{T}{4\pi n^2}\cr
~\stackrel{(7)}{=}~&\frac{1}{B\sqrt{\pi T}}. \end{align}\tag{10}$$
Apparently we should chose the fudge factor $B=\sqrt{2}$ in order to achieve the correct normalization (6).
--
$^1$ Note that the sine (cosine) modes (3) correspond trivially (non-trivially) to the even (odd) modes in my Phys.SE answer here, respectively.
Best Answer
Yes, for this particular Matsubara sum you mentioned, taking different limits will lead to two results differed by 1. This is because the summation in consideration does not converge, which can be seen from the following integral (the continuous limit of the sum) considering the ultraviolet divergence (the large $\omega$ behavior) $$\frac{1}{2\pi}\int_{-\infty}^{\infty}\mathrm{d}\omega\frac{1}{i\omega-\epsilon}\sim\frac{1}{2\pi i}\int_{-\infty}^{\infty}\mathrm{d}\omega\frac{1}{\omega}\rightarrow\infty.$$ Asking for the result of an essentially divergent sum will not even lead to a definite answer. Introducing the factor $e^{i\omega_n\tau}$ is to control the convergence of the sum, but the result will then depend on how one chooses to control the convergence, i.e. $\tau\rightarrow0^+$ (controlling the left-half complex plane) or $0^-$ (controlling the right-half complex plane).
The two results will be different by a shift of 1. Even this difference "1" has a physical meaning, reflecting the ambiguity in defining Fermion particle number. Because the summation $$\frac{1}{\beta}\sum_n\frac{1}{i\omega_n-\epsilon}=G(\tau=0)=-\mathcal{T}_{\tau}\langle\psi(\tau=0)\bar{\psi}\rangle$$ corresponds to the Green's function at imaginary time 0, which physically means to count the number of Fermions on the energy level $\epsilon$. Here $\mathcal{T}_{\tau}$ stands for the time ordering operator. However $\tau=0$ can be either understood as $\tau\rightarrow0^+$ or $\tau\rightarrow0^-$ (corresponding to the two ways to control the convergence), which becomes very critical here because the time ordering operator will then order the two cases differently: $$-\mathcal{T}_{\tau}\langle\psi(\tau=0^+)\bar{\psi}\rangle = -\langle\psi\bar{\psi}\rangle,$$ $$-\mathcal{T}_{\tau}\langle\psi(\tau=0^-)\bar{\psi}\rangle = \langle\bar{\psi}\psi\rangle.$$ According to the anticommutation relation of the Fermion operators, we have $\bar{\psi}\psi=1-\psi\bar{\psi}$, that is why it is expected that the two results will differ by a constant 1. In fact $\bar{\psi}\psi$ counts the number of particles, while $\psi\bar{\psi}$ counts the number of holes. It becomes a problem whether to define the Fermion number (with respect to the vacuum state) as the particle number or as the negative of the hole number. This depends on our definition of the vacuum state: whether we treat the vacuum as no particle state or as a state filled with particles (like the concept of Fermi sea). We know both choices are acceptable, just as we can either define electron as a particle in the empty space or as a hole in the anti-electron Fermi sea. Because of this ambiguity, one needs to specify the choice when the $\tau\rightarrow0$ limit is taken. In conclusion, both ways of taking the limit are legitimate, and the difference in the results is just a matter of the definition of Fermion vacuum state.
Similar discrepancy also exists in the Bosonic case, which is again related to the definition of Boson vacuum state. In general, such discrepancy roots in the divergence of the Matsubara summation in consideration. If the summation itself converges to a definite result, then the result will be unique no matter which limit is chosen to take. For example, if one tries to calculate $$\frac{1}{\beta}\sum_n\frac{1}{(i\omega_n-\epsilon)^2},$$ the result will always be $\beta/4 \mathrm{sech}^2 (\beta\epsilon/2)$ (for Fermionic case) no matter $\tau\rightarrow0^+$ or $0^-$.