Gauge covariant derivative, - how do I get the field? Suppose I wish to create a gauge covariant derivative from
$$
\psi(x)\to e^{ia(x)}\psi(x)
$$
I first note that the usual derivate is not covariant:
$$
\partial_x(e^{ia(x)}\psi(x))=i(\partial _xa(x))e^{ia(x)}\psi(x)+e^{ia(x)}\partial_x\psi(x)
$$
To restaure covariance, I define $D_x$ such that $D_x(e^{ia(x)}\psi(x))= e^{ia(x)}\partial_x \psi(x)$
$$
D_x:=\partial_x-i\partial_x a(x)
$$
Here https://en.wikipedia.org/wiki/Gauge_covariant_derivative#Quantum_electrodynamics, it states that the correct answer is :
$$
D_u:=\partial_u-ieA_u-i\partial_u a
$$
I do not understand where does the term $ieA_u$ come from? I thought the gauge derivative was used to obtain the field?
 A: 
To restore covariance, I define $D_x$ such that $D_x(e^{ia(x)}\psi(x))= e^{ia(x)}\partial_x \psi(x)$
$$D_x:=\partial_x-i\partial_x a(x)$$

This isn't right. The idea is that we want $D\psi(x)\mapsto e^{ia(x)}D\psi(x)$ under change of gauge. Your expression isn't right for two reasons:

*

*We don't want to replace $D$ with $\partial$ on the right-hand side.  What if we perform two gauge transformations back to back?  If after the first one we only have a partial derivative acting on $\psi$, then we're out of luck.

*We can't bake a particular gauge function into $D$, because the gauge functions can be arbitrary.

The solution is to replace $\partial$ with $D:= \partial -iA$, where $A$ is an auxiliary field (which we will presumably figure out how to calculate later).  When we perform a gauge transformation $\psi(x) \mapsto e^{ia(x)}\psi(x)$, we simultaneously transform $A \mapsto A+\partial a(x)$.  If we do so, then obtain
$$D\psi(x):= (\partial-iA)\psi(x) \mapsto (\partial -i[A+\partial a(x)])e^{ia(x)}\psi(x) = e^{i\alpha(x)}D\psi(x)$$
as desired.

I thought the gauge derivative was used to obtain the field?

Our desire to have a gauge covariant derivative is why we need the field $A$, which has the specific transformation property $A\mapsto A+\partial a(x)$ under changes of gauge.  Because $A$ has an observable effect on the universe, presumably there must be a way to compute it.  A reasonable idea would be to build a gauge-invariant scalar out of it and use it as a Lagrangian density; the simplest such scalar is
$$\mathcal L_A = -\frac{1}{4}F^2= -\frac{1}{4}F_{\mu\nu}F^{\mu\nu}$$
where $F_{\mu\nu} := \partial_\mu A_\nu - \partial_\nu A_\mu$.  The resulting Euler-Lagrange equations for $A$ are precisely the Maxwell equations (expressed in terms of the 4-potential), where the coupling to $\psi$ occurs through $D_x$ (e.g. $\overline{\psi} \gamma^\mu D_\mu \psi$ for a spinor field).
