The nature of the Maxwell equations (including their Lagrangian aspects) can be clarified by viewing these equations as part of a larger family of classical field theories in physics. For simplicity I will focus on just four theories here:

- Vacuum classical Yang-Mills theory;
- Abelian vacuum classical Yang-Mills theory with gauge group $U(1)$, aka the vacuum Maxwell equations;
- The vacuum Einstein equations;
- The vacuum linearised Einstein equations.

However one could vastly generalise the discussion if desired, by adding a source matter term to any of these equations, or by coupling with other field theories to obtain such classical field theories as Maxwell-Klein-Gordon, Maxwell-Chern-Simons, Einstein-Maxwell, Yang-Mills-Higgs, Einstein-Vlasov, nonlinear sigma fields (wave maps), etc.. One can also similarly discuss quantum field theories but I will not attempt to do so here.

These equations all describe the dynamics of one or more underlying geometric objects, expressible in local coordinates by tensors (which I will call the *primary fields*):

- For Yang-Mills the primary field is the connection one-form $A_\alpha$ (whose components take values in a Lie algebra ${\mathfrak g}$).
- For Maxwell the primary field is the vector potential $A_\alpha$ (whose components take values in ${\bf R}$).
- For Einstein the primary field is the metric $g_{\alpha \beta}$.
- For linearised Einstein the primary field is the metric perturbation $h_{\alpha \beta}$.

For coupled systems one could have more than one primary field, but here we focus on simpler models in which there is just one primary field in each theory. From this primary field one can then create some secondary fields:

- For Yang-Mills one can form the curvature tensor $F_{\alpha \beta} = \partial_\alpha A_\beta - \partial_\beta A_\alpha + [A_\alpha, A_\beta]$ (normalisation conventions can vary from author to author here).
- For Maxwell one can form the electromagnetic field $F_{\alpha \beta} = \partial_\alpha A_\beta - \partial_\beta A_\alpha$.
- For Einstein one can form the Christoffel symbols $\Gamma_{\alpha \beta}^\gamma = \frac{1}{2} g^{\gamma \delta} (\partial_\beta g_{\gamma \alpha} + \partial_\alpha g_{\gamma \beta} - \partial_\gamma g_{\alpha \beta})$ and the Riemann curvature tensor $R_{\alpha \beta \gamma}^\delta = \partial_\beta \Gamma_{\alpha \gamma}^\delta - \partial_\alpha \Gamma_{\beta \gamma}^\delta + O(\Gamma^2)$ (where I omit the precise nature of the quadratic term $O(\Gamma^2)$ for simplicity).
- For linearised Einstein one can form the linearised Christoffel symbol $\Gamma_{\alpha \beta}^\gamma = \frac{1}{2} \eta^{\gamma \delta} (\partial_\beta h_{\gamma \alpha} + \partial_\alpha h_{\gamma \beta} - \partial_\gamma h_{\alpha \beta})$ and the linearised Riemann curvature tensor $R_{\alpha \beta \gamma}^\delta = \partial_\beta \Gamma_{\alpha \gamma}^\delta - \partial_\alpha \Gamma_{\beta \gamma}^\delta$.

These secondary fields are not unconstrained, but by their definition obey various Bianchi-type identities (as well as some symmetry and antisymmetry identities not mentioned here):

- For Yang-Mills one has the Bianchi identity $D_\alpha F_{\beta \gamma} + D_\beta F_{\gamma \alpha} + D_\gamma F_{\alpha \beta} = 0$, where $D_\alpha = \partial_\alpha + [A_\alpha,]$ is the covariant derivative.
- For Maxwell one has the Gauss-Faraday law $\partial_\alpha F_{\beta \gamma} + \partial_\beta F_{\gamma \alpha} + \partial_\gamma F_{\alpha \beta} = 0$.
- For Einstein one has the second Bianchi identity $\nabla_\alpha R_{\beta \gamma \delta}^\sigma + \nabla_\beta R_{\gamma \alpha \delta}^\sigma + \nabla_\gamma R_{\alpha \beta \delta}^\sigma = 0$.
- For linearised Einstein one has the linearised second Bianchi identity $\partial_\alpha R_{\beta \gamma \delta}^\sigma + \partial_\beta R_{\gamma \alpha \delta}^\sigma + \partial_\gamma R_{\alpha \beta \delta}^\sigma = 0$.

These identities are best thought of as geometric identities rather than true laws of motion, despite the fact that some components of these equations would involve a time derivative when viewed in standard coordinate systems. (All of these Bianchi identities can be derived mathematically from the Jacobi identity $[\nabla_\alpha,[\nabla_\beta, \nabla_\gamma]] + [\nabla_\beta,[\nabla_\gamma,\nabla_\alpha]] + [\nabla_\gamma,[\nabla_\alpha,\nabla_\beta]]=0$ applied to the covariant derivative $\nabla_\alpha$ that is naturally associated to the field theory.) In the case of Maxwell and linearised Einstein, the Bianchi-type identities can be expressed purely in terms of secondary fields, but in general we see that the Bianchi identities involve both primary and secondary fields.

To each of these field theories one can associate a Lagrangian which is typically expressed in terms of both primary and secondary fields:

- For Yang-Mills the functional is $\frac{1}{2} \int F^{\alpha \beta} \cdot F_{\alpha \beta}\ d\eta$ (using an invariant inner product on ${\mathfrak g}$, and with $d\eta$ denoting the Minkowski volume form).
- For Maxwell the functional is $\frac{1}{2} \int F^{\alpha \beta} F_{\alpha \beta}\ d\eta$.
- For Einstein the functional is the Einstein-Hilbert action $\frac{1}{2} \int R\ dg$ ($R = g^{\alpha \beta} R_{\alpha \gamma \beta}^\gamma$ being the scalar curvature).
- For linearised Einstein the functional is the quadratic component of the Einstein-Hilbert action applied to an infinitesimal perturbation $g = \eta + \varepsilon h$; it has a moderately complicated form, see e.g., equations (6), (7) of this paper of Menon.

One can express these Lagrangians purely in terms of the primary field by expanding every appearance of a secondary field in terms of primary fields. In the case of Yang-Mills and Maxwell one can also express the functional purely in terms of secondary fields, but this is a fluke that is not expected to occur for more general field equations.

By varying these functionals with respect to the primary fields (or with respect to the secondary fields using the Bianchi-type identities as constraints that produce Lagrange multipliers, in the cases where those Bianchi identities can be expressed purely in terms of secondary fields), one obtains the equations of motion as Euler-Lagrange equations:

- For Yang-Mills the equations of motion are $D^\alpha F_{\alpha \beta} = 0$.
- For Maxwell the equations of motion are the Gauss-Ampere law $\partial^\alpha F_{\alpha \beta} = 0$.
- For Einstein the equations of motion are $R_{\alpha \beta \gamma}^\beta - \frac{1}{2} R g_{\alpha \gamma} = 0$.
- For linearised Einstein the equations of motion are $R_{\alpha \beta \gamma}^\beta - \frac{1}{2} R \eta_{\alpha \gamma} = 0$ (with $R = \eta^{\alpha \beta} R_{\alpha \gamma \beta}^\gamma$ being the linearised scalar curvature).

In the latter three cases one can express the equations of motion purely in terms of secondary fields, but again this is something of a fluke that should not be expected in general (particularly when working with nonlinear theories). Note that in all four cases the insertion of matter into the theory would add a matter term to the Lagrangian and thus a source term to the equations of motion, but would not affect the geometric Bianchi-type identities, which highlights a key difference between the equations of motion and the Bianchi-type identities.

In all of these cases, there is also a gauge invariance that prevents one from observing the primary fields directly.

- For Yang-Mills, one can replace the connection $A_\alpha$ by the gauge-equivalent connection $U A_\alpha U^{-1} - (\partial_\alpha U) U^{-1}$ without affecting the physical observables.
- For Maxwell, one can similarly replace the vector potential $A_\alpha$ by the gauge-equivalent potential $A_\alpha + \partial_\alpha \lambda$ without affecting the physical observables.
- For Einstein, one can apply a diffeomorphic change of coordinates to spacetime and to the metric $g^{\alpha \beta}$ without affecting the physical observables (the general principle of relativity).
- For linearised Einstein, one can replace the metric perturbation $h_{\alpha \beta}$ by $h_{\alpha \beta} + \partial_\alpha X_\beta + \partial_\beta X_\alpha$ for any vector field $X^\alpha$ (which can be viewed as a linearised (or infinitesimal) diffeomorphism) without affecting the physical observables.

Because of this gauge invariance, it is initially tempting in all four cases to try to ignore the primary field as much as possible and try to formulate the physics in terms of the secondary fields (which tend to be much more directly physically measurable, being less affected by gauge transforms). However, despite the fact that in some cases the primary field can indeed be expunged from the equations of motion and/or the Lagrangian, there is no *a priori* reason why the geometric Bianchi-type identities should then be subsumed into the Euler-Lagrange equations (though one can always do so artificially, as explained in my other answer). The fact that the Bianchi identities and the equations of motion look so similar in the case of Maxwell's equations (when specialised to 3+1 spacetime dimensions) is a misleading red herring; once one sees how things work for all the other field equations one sees that these two systems of equations should really be treated as being of a very different nature.

4more comments