Deriving Lagrangian density for electromagnetic field In considering the (special) relativistic EM field, I understand that assuming a Lagrangian density of the form
$$\mathcal{L} =-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}  + \frac{1}{c}j_\mu A^\mu$$
and following the Euler-Lagrange equations recovers Maxwell's equations. 
Does there exist a first-principles derivation of this Lagrangian? A reference or explanation would be greatly appreciated!
 A: I'm almost 100% sure the Lagrangian is an assumption of the theory.  It cannot be derived.  I don't have any references for this claim.  I just know that from every course I've been taught and every book I've read, the Lagrangian (assuming it is being used at all) is where you start.  It is the "first principle" in this case.
A: You can find the answer from the book "Differential Geometry and Lie Groups for Physicists" by Marian Fecko.
In geometrical language, an action of a field $F\in\Omega^{p}(M)$ on some n-dim Riemannian manifold $(M,g)$ should be understand as an 'inner-product'
$$\int_{M}F\wedge\ast_{g}F,$$
where $p<n$, and $\ast_{g}$ is the hodge star operator, i.e. 
$$\ast_{g}:\Omega^{p}(M)\rightarrow\Omega^{n-p}(M)$$
so that the action is diffeomorphism invariant and its density is a scalar. 
For example, the action of a free scalar field takes this form:
$$S[\phi]=\int_{M}d\phi\wedge\ast_{g}d\phi+\frac{m}{2}\int_{M}\phi\wedge\ast_{g}\phi=\int_{M}\sqrt{|g|}d^{n}x\left\{\partial_{\mu}\phi\partial^{\mu}\phi+\frac{m\phi^{2}}{2}\right\}.$$
When the 'worldsheet' manifold is 1-dim, only possible fields are 1-forms. One can consider an action takes the following form 
$$\int_{\mathbb{R}}A_{\mu}\frac{dx^{\mu}}{ds}ds$$
where $A=A_{s}ds=A_{\mu}\frac{dx^{\mu}}{ds}ds$ is a 1-form on the worldline, whose hodge star dual field is not defined.
A: Ultimately the reasoning must be that (as you stated) it must be constructed so the Euler-Lagrange equations are Maxwell's equations. So in a sense you have to guess the Lagrangian that produces this as is done here for example.
However you can get some guidance from the fact that we need to construct a Lagrangian for a massless non self interacting field. So we need a gauge and lorentz invariant combination of the 4-vector potential which only has a kinetic term (quadratic in derivatives of the fields). You are then not left with many options apart from $F^{\mu\nu}F_{\mu\nu}$. The source term is then trivial to add in if needed.
A: You can use the symmetries of E&M to show that there's essentially only one reasonable candidate that needs to be checked:

The action we seek should be Lorentz invariant, gauge invariant, parity and time-reversal invariant, and no more than second order in derivatives. The only candidate is [the Maxwell action]. [Srednicki QFT pg. 334.]

