Maxwell's equations can be derived from a Lagrangian formulation using the Lagrangian term (modulo some constants) $$\mathcal L=-\frac{1}{4}F^{\mu\nu}F_{\mu\nu}-A_\mu J^\mu.$$ Focusing on the free term for the moment, I've seen mentioned (though I can't find a source right now) that the $F^2$ term is the only possible gauge-invariant Lagrangian for electromagnetism.
Is this the case? How can we prove it?
Not sure if this is what you're looking for, but here is a exhaustive search for all gauge invariant objects in electromagnetism: Say $f[A_{\mu},\partial_{\mu} A_{\nu},\partial_{\mu}\partial_{\nu} A_{\rho},...]$ is a gauge invariant object. Here $f$ could itself have indices which have been suppressed. Under a gauge transform:
$$A'_{\mu}(x)=A_{\mu}(x)+\partial_{\mu}\theta$$ Gauge invariance means that: $$f[A'_{\mu},\partial_{\mu} A'_{\nu},\partial_{\mu}\partial_{\nu} A'_{\rho},...]=f[A_{\mu},\partial_{\mu} A_{\nu},\partial_{\mu}\partial_{\nu} A_{\rho},...]$$ One can of course expand the left hand side to first order in $\theta$ $$f[A'_{\mu},\partial_{\mu} A'_{\nu},\partial_{\mu}\partial_{\nu} A'_{\rho},...]=f[A_{\mu},\partial_{\mu} A_{\nu},\partial_{\mu}\partial_{\nu} A_{\rho},...]+\frac{\delta f}{\delta A_{\mu}}\partial_\mu\theta+\frac{\delta f}{\delta \partial _\mu A_{\nu}}\partial_\mu\partial_\nu \theta+\frac{\delta f}{\delta \partial _\mu \partial _\nu A_{\rho}}\partial_\mu\partial_\nu\partial_\rho\theta+....$$ Since $\theta$ is an arbitrary function of $x$, each of the offending terms should separately vanish. $$\frac{\delta f}{\delta A_{\mu}}=0$$ For the next term: $$\frac{\delta f}{\delta \partial _\mu A_{\nu}}\partial_\mu\partial_\nu \theta=0$$ This means $\frac{\delta f}{\delta \partial _\mu A_{\nu}}$ is antisymmetric in $\mu, \nu$. i.e. $f=f[F_{\mu,\nu},....]$ For the next term, $$\frac{\delta f}{\delta \partial _\mu \partial _\nu A_{\rho}}\partial_\mu\partial_\nu\partial_\rho\theta=0$$ This means $\frac{\delta f}{\delta \partial _\mu \partial _\nu A_{\rho}}$ is antisymmetric in $\nu, \rho$(it can't be antisymmetric in any other pair of indices). i.e. $f=f[F_{\mu,\nu},\partial_{\rho}F_{\mu \nu}....]$ So to first order in the gauge parameter, $f$ must be a function of the field strength and its derivatives. Since we know that these quantities are gauge invariant to all orders in the gauge parameter, this is an exact result.
So the upshot is: any function of $F_{\mu \nu}$ and its derivatives is gauge invariant. This is the only way to form a gauge invariant object.
A sensible quantum theory requires we exclude higher derivatives of $F_{\mu \nu}$ and exclude higher powers of $F_{\mu \nu}$(higher than 2). Poincare invariance requires us to soak up all the indices, therefore requiring at least two powers of $F_{\mu \nu}$. This leaves us with only two terms in the Lagrangian: $$\mathcal{L} = c_1F_{\mu \nu}F^{\mu \nu}+c_2 \epsilon_{\mu \nu \rho \sigma}F^{\mu \nu}F^{\rho \sigma}$$
It is up to a trivial constant the only gauge invariant Lagrangian that reproduces the Maxwell equations. Since the Maxwell equations are linear in $F^{\mu\nu}$ only quadratic functions of $F^{\mu\nu}$ are possible. There are two possibilities, namely $F_{\mu\nu}F^{\mu\nu}$ and $\epsilon_{\mu\nu\rho\sigma} F^{\mu\nu} F^{\rho\sigma}$. The Maxwell equations can be recovered if the first but not the second term is present in the Lagrangian.
