Noether currents for the BRST tranformation of Yang-Mills fields The Lagrangian of the Yang-Mills fields is given by
$$
  \mathcal{L}=-\frac{1}{4}(F^a_{\mu\nu})^2+\bar{\psi}(i\gamma^{\mu} 
  D_{\mu}-m)\psi-\frac{1}{2\xi}(\partial\cdot A^a)^2+
  \bar{c}^a(\partial\cdot D^{ab})c^b
$$
where the metric is $(-,+,+,+)$, and the conventions are the following: 
$$
  [D_{\mu},D_{\nu}]=-igF_{\mu\nu},\quad 
  D_{\mu}=\partial_{\mu}-igA^a_{\mu}t^a, \quad 
  D^{ab}_{\mu}=\delta^{ab}\partial_{\mu}-gf^{abc}A^c_{\mu}
$$
Let $\epsilon$ be an infinitesimal anticummuting parameter, and consider the BRST transformation: 
$$
 \delta\psi=ig\epsilon c^at^a\psi,\quad
 \delta A^a_{\mu}=\epsilon D^{ab}_{\mu}c^b,\quad
 \delta c^a=-\frac{1}{2}g\epsilon f^{abc}c^bc^c,\quad
 \delta\bar{c}^a=\frac{1}{\xi}\epsilon\partial^{\mu}A^a_{\mu}
$$
I have calculated the corresponding Noether current as 
$$
  j_{BRST}^{\mu}=-g\bar{\psi}\gamma^{\mu}c^at^a\psi-F^{a\mu\nu}D^{ab}_{\nu}c^b-
  \frac{1}{\xi}(\partial\cdot A^a)D^{ab\mu}c^b+
  \frac{1}{2}gf^{abc}(\partial^{\mu}\bar{c}^a)c^bc^c
$$
I am not sure whether the result is correct or not, so I would like to check that $\partial_{\mu}j^{\mu}_{BRST}=0$. Even though I have used the equation of motion 
$$
  \partial_{\mu}F^{a\mu\nu}=-g\bar{\psi}\gamma^{\nu}t^a\psi-
  gf^{abc}A^b_{\mu}F^{c\mu\nu}-\frac{1}{\xi}\partial^{\nu}
  (\partial\cdot A^a)-gf^{abc}(\partial^{\nu}\bar{c}^b)c^c
$$
$$
  (i\gamma^{\mu}D_{\mu}-m)\psi=0,\quad \partial^{\mu}D^{ab}_{\mu}c^b=0
$$
and spent about four hours, I still cannot get it right. Could someone help me check this? Thanks a lot.
 A: The crucial point is the anticommutation of $\epsilon$ with fermion fields ($\psi,\bar\psi,\bar{c}^a,c^a$).

First, we will rewrite $\mathcal{L}$ as (for simplicity, we define $B^a \equiv \xi^{-1}\partial^\mu A_\mu^a$)
$$\tag{1} \mathcal{L}=-\frac{1}{4}(F^a_{\mu\nu})^2+\bar{\psi}(i\gamma^{\mu} 
  D_{\mu}-m)\psi-\frac{\xi}{2}B^aB^a-
  \partial^{\mu}\bar{c}^a D_{\mu}^{ab}c^b $$
which differs from the original $\mathcal{L}$ by a total derivative,
$$\tag{2}  \partial^{\mu}(\bar{c}^a D_{\mu}^{ab}c^b) $$
so $\delta\mathcal{L}$ no longer equals $0$, but
$$\tag{3}  \delta\mathcal{L} = -\delta \partial^{\mu}(\bar{c}^a D_{\mu}^{ab}c^b) = \partial^{\mu}( -\delta\bar{c}^a D_{\mu}^{ab}c^b) = \partial^{\mu}\left(-\epsilon B^a D_{\mu}^{ab}c^b\right) \equiv \partial^{\mu} K_{\mu} $$
We will use Jacobi identity occasionally,
$$\tag{4} f^{abd}f^{dce} + f^{bcd}f^{dae} + f^{cad}f^{dbe} = 0 $$

We will use right derivative in the following calculations. So the Noether current is defined as
$$\tag{5} \epsilon j^{\mu} \equiv \frac{\partial\mathcal{L}}{\partial (\partial_{\mu}\psi)}\delta\psi + \frac{\partial\mathcal{L}}{\partial (\partial_{\mu}\bar\psi)}\delta\bar\psi + \frac{\partial\mathcal{L}}{\partial (\partial_{\mu}c^a)}\delta c^a + \frac{\partial\mathcal{L}}{\partial (\partial_{\mu}\bar{c}^a)}\delta\bar{c}^a + \frac{\partial\mathcal{L}}{\partial (\partial_{\mu}A_{\nu}^a)}\delta A_{\nu}^a - K^{\mu} $$
Now we will calculate the individual parts of the current, moving $\epsilon$ to the left of each expression.  We'll get an extra minus sign if $\epsilon$ passes a fermion field,
$$\begin{aligned}
\frac{\partial\mathcal{L}}{\partial (\partial_{\mu}\psi)}\delta\psi &= (\bar\psi i \gamma^{\mu}) (ig \epsilon c^a t^a \psi) = \epsilon\, g \bar\psi \gamma^{\mu} c^a t^a \psi \\
\frac{\partial\mathcal{L}}{\partial (\partial_{\mu}\bar\psi)}\delta\bar\psi &= 0 \\
\frac{\partial\mathcal{L}}{\partial (\partial_{\mu}c^a)}\delta c^a &= (-\partial^{\mu}\bar{c}^a)\left(-\frac{1}{2}g\epsilon f^{abc}c^bc^c\right) = -\frac{1}{2} \epsilon\, g f^{abc} (\partial^{\mu}\bar{c}^a) c^bc^c \\
\frac{\partial\mathcal{L}}{\partial (\partial_{\mu}\bar{c}^a)}\delta\bar{c}^a &= (g^{\mu\nu}D_\nu^{ab} c^b) \left(\epsilon B^a\right) = - \epsilon  (g^{\mu\nu}D_\nu^{ab} c^b) B^a = K^{\mu}\\
\frac{\partial\mathcal{L}}{\partial (\partial_{\mu}A_{\nu}^a)}\delta A_{\nu}^a &= \left(-F^{a\mu\nu} - g^{\mu\nu} B^a\right)(\epsilon D_\nu^{ab} c^b) = \epsilon \left(-F^{a\mu\nu} - g^{\mu\nu}B^a\right)(D_\nu^{ab} c^b)\\
\end{aligned} \tag{6}$$
Inserting results from $(6)$ and $(3)$ into $(5)$ gives
$$\tag{7} j^\mu = \left(-F^{a\mu\nu} - g^{\mu\nu}B^a\right)D_\nu^{ab} c^b -\frac{1}{2} g f^{abc} (\partial^{\mu}\bar{c}^a) c^bc^c + g \bar\psi \gamma^{\mu} c^a t^a \psi $$
It's easy to derive the equations of motion,
$$\begin{aligned}
0= \frac{\partial\mathcal{L}}{\partial A_\nu^a} - \partial_\mu \frac{\partial\mathcal{L}}{\partial(\partial_\mu A_\nu^a)} &= D_\mu^{ab}F^{b\mu\nu} + \partial^\nu B^a + g \bar\psi \gamma^\nu t^a \psi + g f^{abc}(\partial^\nu \bar{c}^b) c^c \\
0= \frac{\partial\mathcal{L}}{\partial\bar\psi} - \partial_\mu \frac{\partial\mathcal{L}}{\partial(\partial_\mu \bar\psi)} &= -i \gamma^\mu\partial_\mu\psi - g A_{\mu}^a \gamma^\mu t^a \psi + m \psi \\
0= \frac{\partial\mathcal{L}}{\partial\psi} - \partial_\mu \frac{\partial\mathcal{L}}{\partial(\partial_\mu \psi)} &= -i \partial_\mu \bar\psi \gamma^\mu + g A_\mu^a \bar\psi \gamma^\mu t^a - m \bar\psi \\
0= \frac{\partial\mathcal{L}}{\partial c^a} - \partial_\mu \frac{\partial\mathcal{L}}{\partial(\partial_\mu c^a)} &= D_\mu^{ab} \partial^\mu \bar{c}^b \\
0= \frac{\partial\mathcal{L}}{\partial\bar{c}^a} - \partial_\mu \frac{\partial\mathcal{L}}{\partial(\partial_\mu \bar{c}^a)} &= -\partial^\mu D^{ab}_\mu c^b \\
\end{aligned} \tag{8a-e}$$
Now we will check the validation of $\partial_\mu j^\mu = 0 $,
$$\begin{aligned}
\partial_\mu \Bigl[ g \bar\psi \gamma^\mu c^a t^a \psi \Bigr] &\stackrel{(8b,8c)}{=} -\Bigl[ g \bar\psi \gamma^\nu t^a \psi \Bigr] D_\mu^{ad} c^d \\
\partial_\mu \left[ - \frac{1}{2}g f^{abc} (\partial^\mu \bar{c}^a) c^b c^c \right] &\stackrel{(4,8d)}{=} -\Bigl[ g f^{abc}(\partial^\nu \bar{c}^b) c^c \Bigr] D_\mu^{ad} c^d \\
\partial_\mu \Bigl[ (-F^{a\mu\nu} -B^a g^{\mu\nu}) D_\nu^{ab} c^b  \Bigr] &\stackrel{(8e)}{=} -(\partial_\mu F^{a\mu\nu} + \partial^\nu B^a) D_\nu^{ad}c^d - F^{a\mu\nu} \partial_\mu D_\nu^{ad}c^d \\
&= -(D_\mu^{ab} F^{b\mu\nu} + \partial^\nu B^a) D_\nu^{ad}c^d \\
&\quad - (gf^{abc}A_\mu^c F^{b\mu\nu}) D_\nu^{ad} c^d - F^{a\mu\nu} \partial_\mu D_\nu^{ad}c^d \\
&\stackrel{(4)}{=} -(D_\mu^{ab} F^{b\mu\nu} + \partial^\nu B^a) D_\nu^{ad}c^d \\
\end{aligned} \tag{9a-c}$$
We'll get
$$ \partial_\mu j^\mu \stackrel{(9a+9b+9c,8a)}{=} -\left(  \frac{\partial\mathcal{L}}{\partial A_\nu^a} - \partial_\mu \frac{\partial\mathcal{L}}{\partial(\partial_\mu A_\nu^a)} \right) D_\nu^{ad}c^d \tag{10}$$
This is indeed the off-shell Noether identity, and is zero when the equations of motion are fulfilled.
