Decoupling of ghost fields in axial-gauge QCD After quantizing QCD using the Faddeev-Popov "prescription", we end up with the original QCD Lagrangian plus the gauge-fixing term,
\begin{equation}
-\frac{1}{2\alpha}(n\cdot A)^2,
\end{equation}
and the ghost fields action
\begin{equation}
S_\mathrm{g}(\phi,\bar{\phi},A)=
\int\bar{\phi}(x)\bigl([n\cdot A(x),\phi(x)]+n\cdot\mathrm{d}\phi(x)\bigr)\,\mathrm{d}x.
\end{equation}
It is usually said that, using the axial gauge, the ghost fields decouple from the gauge field.
As long as $A$ appears in the ghost fields action $S_\mathrm{g}(\phi,\bar{\phi},A)$, a ghost-gluon vertex is created, so ghosts don't go away.
In $S_\mathrm{g}(\phi,\bar{\phi},A)$, $A$ appears in the product $n\cdot A$: I thought that the gauge condition $n\cdot A=0$ would help to eliminate this term, effectively removing $A$ from $S_\mathrm{g}(\phi,\bar{\phi},A)$.
But wouldn't this mean that the gauge fixing term is zero, too? Surely it cannot be, or we would be back at the beginning of the whole gauge-fixing procedure.
Also, the way the Faddeev-Popov prescription is usually presented in the literature, in order to "create" the gauge-fixing term, it requires a modification of the gauge condition $n\cdot A=0$ to $n\cdot A-\nu=0$ where $\nu$ is some $\mathrm{su}(N)$-valued function (just like $A$), then an integration on $\nu$ using a Gaussian weight, which in the end becomes the gauge-fixing term.
But then $n\cdot A$ isn't zero, so the relative term in the ghost action shouldn't even cancel, if I'm guessing correctly.
Exactly then how can I prove that the ghost fields really decouple?
 A: I believe that the answer to your question lies in thinking about the whole path integral, rather than just the Lagrangian of QCD. The purpose of using the Faddeev-Popov prescription is not necessarily to modify the form of the Lagrangian itself, but to ensure that the path integral is not summing over field configurations which are "equivalent", i.e. field configurations which differ only by a gauge transformation. Indeed, you are right that for the choice of gauge $n\cdot A =0$ the QCD part of the Lagrangian is identical to its form before we gauge fixed, and this will manifest in particular in the propagator for the gauge field, which will simply become
$$ G^{\mu \nu}_{ab}(p) = \delta_{ab} \frac{g^{\mu \nu}}{p^2+i\varepsilon}$$
i.e. precisely what you might have guessed for the propagator had you not considered gauge fixing at all. The difference, however, is that now your path integral measure has changed. Whereas without gauge fixing one would have
$$Z_{\text{unfixed}} = \int \mathcal{D}\bar\psi \mathcal{D}\psi \mathcal{D}A \exp(iS_{\text{QCD}}[\psi, \bar\psi, A])$$what one actually wants is to "gauge fix" the measure, which yields the result you have presented, namely
$$Z_{\text{QCD}} = \int \mathcal{D}\bar\psi \mathcal{D}\psi \mathcal{D}\bar\phi\mathcal{D}\phi \mathcal{D}A \exp\left(iS_{\text{QCD}}[\psi, \bar\psi, A] + iS_{\text{ghost}}[\phi, \bar\phi, A] - \frac{i}{2\alpha} (n \cdot A)^2\right)$$
The point here is that by picking the gauge $n \cdot A = 0$, you remove the gauge-ghost coupling, effectively turning $S_{\text{ghost}}[\phi, \bar\phi, A]$ into just $S_{\text{ghost}}[\phi, \bar\phi]$ and removing the new term which depends on $\alpha$, but you still have a different path integral measure which is properly gauge invariant. 
As to your question on decoupling of the ghost fields, I'm not entirely sure what you mean by this, as in a general gauge the ghosts do not decouple from the gauge fields. Indeed, you will have a gauge-ghost vertex which introduces ghost propagators as internal legs in your Feynman diagrams. The relationship between the diagrams containing ghosts and those that do not will be characterized by the Ward identities. In your gauge choice $n\cdot A = 0$, the ghost terms in these Ward identities will drop out, but the gauge invariance is still preserved by what remains in the Ward identity. 
As a reference to all of this, Peskin and Schroeder presents a fairly detailed discussion of QCD and quantization of non-Abelian gauge theories in general. That may be worth checking out, if you haven't already.
A: In the path integral with a $R_{\xi}$-gauge-fixing term ${\cal L}_{GF}=-\frac{\chi^2}{2\xi}$, the axial gauge-fixing condition $\chi=n\cdot A\approx 0$ is only imposed in a quantum average sense. In general the gauge-fixing condition may be violated by quantum fluctuations, except in the Landau gauge $\xi=0^+$, where such quantum fluctuations are exponentially suppressed (in the Wick-rotated Euclidean path integral). Therefore, only in the Landau gauge $\xi=0^+$, we may remove $n\cdot A$ from the Faddeev-Popov (FP) term. In this case the FP ghosts decouple from the gluon-field, cf. OP's question.
A: If you want, a more concrete way to show that the ghosts decouple is in the diagrammatic expansion. The gauge-fixed Lagrangian is 
$$\mathcal{L} = - \frac14 (F_{\mu\nu}^a)^2 - \frac{1}{2 \alpha} (n^\mu A_\mu)^2 - \bar{c} n^\mu D_\mu c.$$
Inverting the quadratic part of the Lagrangian to find the propagator for $A_\mu$, in the same way as usual, gives
$$i \Delta^{\mu\nu}_{ab}(k) = \frac{i \delta_{ab}}{k^2 + i \epsilon} \left(\eta^{\mu\nu} - \frac{k^\mu n^\nu + n^\mu k^\nu}{k \cdot n} - \frac{\alpha k^2 - n^2}{(k \cdot n)^2} k^\mu k^\nu \right).$$
Note that the interaction vertex between the ghosts and gauge boson contains the factor $n_\mu A^\mu$, and hence is proportional to
$$n_\mu \Delta^{\mu\nu}_{ab}(k) = - \frac{\alpha}{k \cdot n} \delta_{ab} k^\nu.$$
Taking the limit $\alpha \to 0$ where the gauge condition is exactly imposed, all diagrams where a ghost attaches to a gluon vanish, so the ghosts decouple.
