The Faddeev-Popov Lagrangian This is a non-abelian continuation of this QED question.
The Lagrangian for a non-abelian gauge theory with gauge group $G$, and with fermion fields and ghost fields included is given by
$$
\mathcal{L}=\overline{\psi}(i\gamma ^\mu D_\mu -m)\psi+\sum _{k=1}^{\dim (G)}\left[ \frac{-1}{4}F_{\mu\nu}^kF^{\mu \nu k}+\frac{1}{2\xi}(\partial _\mu A^{\mu k})^2+\sum _{i=1}^{\dim (G)}\overline{c}^k(-\partial _\mu D^{\mu ki})c^i\right] .
$$
How does the third term (the one that contains $\xi$) come into the picture?  The first terms is the standard fermion Lagrangian interacting with a gauge field, the second term is the standard term for gauge bosons, and the fourth arises because of the introduction of the ghost fields . . . but what about the third?
 A: I would think that almost every textbook that includes the Faddeev-Popov Lagrangian also explains it but you may also find explanations on the web, e.g. one by Faddeev himself in Scholarpedia:

http://www.scholarpedia.org/article/Faddeev-Popov_ghosts

The FP ghosts are needed to restore unitarity at one-loop level and are the key new players of the "modern covariant (BRST) quantization" of theories with gauge symmetries. Their existence may be most easily explained in Feynman's path integral approach to gauge theories. At the end, we need to "gauge fix" the gauge symmetry i.e. to choose a particular representative of physically equivalent field configurations in order to avoid (infinite-fold) multiple counting. This means that we're effectively inserting a delta-functional to the path integral.
However, $\delta(kx)$ isn't the same thing as $\delta(x)$: it is $|k|$ times smaller. Similarly, for a multi-dimensional delta-function or delta-functional, the ratio is given by a Jacobian (determinant of the matrix of derivatives). The only legitimate delta-functional would be one that imposes a particular (trivial) gauge transformation. However, the gauge-fixing conditions want to make other choices such as $A_3=0$ and a corresponding Jacobian has to be inserted to convert this delta-functional to the right one. The Jacobian is a determinant that may be expressed as a path integral over new fermionic fields.
One may also motivate the need for FP ghosts by discussing the BRST quantization based on $Q$, a nilpotent BRST charge obeying $Q^2=0$, a useful tool to describe physical states in all theories with gauge symmetries. Physical states are cohomologies of $Q$. This template of the answer automatically eliminates both the states violating the Gauss' constraint (and its generalizations) as well as states that are "pure gauge" and the FP ghosts are needed to define such a $Q$.
The term $(\partial\cdot A)^2/2\xi$ is a particular gauge-fixing term that  eliminates the gauge redundancy and imposes a gauge-fixing condition, in this case it "softly" imposes the Lorenz (not Lorentz!) condition. We may imagine that aside from the equations of motion, one also imposes an additional constraint, the Lorenz condition, on the gauge field. But aside from this term which could be replaced by another one if we decided for a different gauge-fixing condition, one still needs to include the FP terms, at least in non-Abelian gauge theories.
A: The third term is the gauge fixing term, you can think of $1/2\xi$ as a Lagrange multiplier. The EOM for $\xi$ implement the gauge fixing. That is the intuitive picture, the quantization of gauge systems can be treated at various levels of sophistication, luckily Yang-Mills theory is a relatively easy case.
A: The third term is a gauge fixing term. At this point in time, the Wikipedia article on gauge fixing has has a nice section on it, here: https://en.wikipedia.org/wiki/Gauge_fixing#R.CE.BE_gauges These express the so-called $R_ξ$ gauges - with $\xi=1$ called the Feynman–'t Hooft gauge, and the $\lim\xi\to0$ being the Landau gauge, the limit taken after completing calculations.
