Gauge Field Tensor from Wilson Loop It is possible to introduce the gauge field in a QFT purely on geometric arguments. For simplicity, consider QED, only starting with fermions, and seeing how the gauge field naturally emerges. The observation is that the derivative of the Dirac field doesn't have a well-defined transformation, because:
$$n^\mu \partial_\mu \,\psi = \lim\limits_{\epsilon\rightarrow 0}\big[\psi(x+\epsilon n)-\psi(x)\big],$$
i.e. the derivative combines two fields at different spacetime points (having different transformation rules). We need to introduce a parallel transporter $U(y,x)$ that transforms as
$$U(y,x) \rightarrow e^{ig\alpha(y)}U(y,x) e^{-ig \alpha(x)},$$
such that we can adapt the defintion of the derivative into a covariant derivative, that transforms in a well-defined way:
$$n^\mu D_\mu \,\psi = \lim\limits_{\epsilon\rightarrow 0}\big[\psi(x+\epsilon n)- U(x+\epsilon n,x)\psi(x)\big].$$
From geometric arguments, it is straightforward to show that the parallel transporter is a Wilson line:
$$U(y,x) = \mathcal{P}\,e^{ig\int\limits_x^y dz^\mu\, A_\mu(z)},$$
which introduces a new field, namely the gauge field $A_\mu$. See e.g. Peskin & Schroeder Chapter  15 for more detail.
However.. Where the interaction term $\bar{\psi} A_\mu \psi$ emerged in a natural way, I totally don't see how the kinetic terms emerge. The standard way to proceed, is to consider a Wilson loop (a Wilson line on a closed path), and use Stokes' theorem:
$$\text{exp}\left\{ig\oint_\mathcal{C}dx^\mu\, A_\mu \right\} = \text{exp}\left\{ig\int_\Sigma dx^\mu \wedge  dx^\nu\,\left(\partial_\mu A_\nu - \partial_\nu A_\mu \right)\right\},$$
where of course $\partial_\mu A_\nu - \partial_\nu A_\mu \equiv F_{\mu\nu}$. In Peskin & Schroeder, they then consider a small rectangular loop, and see that in the limit $\epsilon\rightarrow 0$, $F_{\mu\nu}$ is invariant. But what's the point? I mean, the transformation law for $A_\mu$ is easily calculated from the definition of the Wilson loop:
$$A_\mu\rightarrow A_\mu+\partial_\mu \alpha,$$
making $F_{\mu\nu}$ invariant by definition:
$$F_{\mu\nu}\rightarrow \partial_\mu A_\nu - \partial_\nu A_\mu +\square \alpha-\square \alpha. $$
I would have liked to see a calculation, starting from a particular loop parameterisation, that naturally leads to the correct kinetic terms in the Lagrangian, as was the case for the interaction term. In other words
$$\text{exp}\left\{ig\oint_\mathcal{C}dx^\mu\, A_\mu \right\} \leadsto -\frac{1}{4}\left(F_{\mu\nu}\right)^2,$$
but I have no idea how to do it.
Or is the idea simply 'look I have found some quadratic derivative terms that are invariant, now let me fiddle a bit and put its square in $\mathcal{L}$'? If yes, then why do Peskin & Schroeder bother calculating a loop parameterisation (p484), if using Stokes' theorem would have been enough to find $F_{\mu\nu}$ somewhere?
EDIT
Based on the comments, a few clarifications. This is not a question about why we need kinetic terms or how they should like. I perfectly know how to construct the SM Lagrangian the standard way. You have a vector field $A_\mu$, you need kinetic terms so you use some Klein-Gordon like structure, adapt it a bit (because of gauge behaviour of $A_\mu$) and so on. Standard QFT.
But this is all manual, almost a bit trial and error, putting in terms because you know they work the way they should. It is like making the Lagrangian as a puzzle: just put in those pieces that fit. No problem with that, it works and is a widely accepted way of doing physics, but from a theoretical point of view, it is not so elegant.
But we can do it in a more elegant way.
If we start only from a Dirac field and the Dirac equation, then purely by mathematical and geometric arguments, the gauge field appears, as discussed above. The question is if we can also make the kinetic terms for the gauge field appear purely based on mathematical and geometric arguments. If you believe the QFT book by Peskin & Schroeder, you can, starting from a Wilson loop (as discussed above see p484-494). But then you end up with a factor $\epsilon^4$, the loop area squared, in front of the field tensor. You could restrict to a class of loops with area $S=1$, but there are two problems with that:


*

*In the calculation, the exponent and integral where expanded, using the fact that $\epsilon <\!\!< 1$. This contradicts our restriction $S=1$.

*This restriction invalidates the elegance somewhat, as now the loop area has to be fine-tuned. I expected to find something like 'in the limit $\epsilon\rightarrow 0$, the field tensor is what remains'.


So, is it not possible then to get the kinetic terms for the gauge field purely from geometric arguments?
If the answer is no, it doesn't make sense to brag about the natural emergence of the interaction term.
It is not elegant when the interaction term emerges naturally, but you have to hand-pick the kinetic terms. For one, how do you prove that the naturally emerging field (interaction terms) and the one you put in (kinetic terms) are the same field?
 A: There exists an extensive literature for discretization of the abelian and the non-abelian gauge theories, known as lattice QED and lattice QCD, respectively. Here we will only sketch the main idea.
Let us for simplicity use Euclidean signature $(+,+,+,+)$. A small Wilson-loop 
$$\tag{1} W~=~{\rm Tr}{\cal P}e^{ig\int_{\gamma}A}$$ 
lies approximately in a 2-plane. In 4 spacetime dimensions we have six 2-planes labelled by an antisymmetric double index $\mu\nu$, where $\mu,\nu=1,2,3,4$. 
The $F_{\mu\nu}^2$ term is proportional to the next-to-leading term in a small loop expansion of 
$$\tag{2} \prod_{\mu\nu}\frac{W_{\mu\nu}+{\rm c.c.}}{2}~=~1+{\cal O}(F_{\mu\nu}^2). $$
The ${\rm c.c.}$ (complex conjugate) is inserted to render the result real, and to remove the linear terms ${\cal O}(F_{\mu\nu})$ in the abelian case. [In the non-abelian case, the linear terms ${\cal O}(F_{\mu\nu})$ are also removed by tracelessness. See also Ref. 1.]
Various quantities, such as, the action, the fields, and the coupling constant are subject to rescalings and renormalization in order to reproduce the correct continuum theory. In particular when summing over all lattice points of spacetime we should divide with $a^4$, where $a$ is the lattice spacing.
References:


*

*M.E. Peskin & D.V. Schroeder, An Intro to QFT, p. 494.

*M. Caselle, Lattice Gauge Theories and the AdS/CFT Correspondence, arXiv:hep-th/0003119.
