Well, there is of course the cheating way to use constraint fields

$$\mathcal{L} = \alpha_\nu \partial_\mu F^{\mu \nu} + \beta_\nu \partial_\mu *\!F^{\mu\nu}$$

The variation with respect to $\alpha$ and $\beta$ then give you Maxwell's equations, the variation with respect to F will give you equations for the constraint fields which are solved after solving the Maxwell equations (and we do not really care about them).

But the reason why this cannot be done without using auxiliary fields seems to me to be the more general fact that one cannot take a single bosonic field (that is, not a field and its antiparticle comrade) and give a Lagrangian generating a first-order differential equation. Technically, this has to do with how the gamma matrices are able to "double the index" on the gradient for fermions, $\partial_\mu \to \partial_\mu \gamma^\mu_{ab}$.

Actually, $F^{\mu\nu}, R^\mu_{\;\nu \kappa \lambda}$ are "the" proper covariant representations of fields of helicity $h=\pm 1, \pm 2$ and the fact that they do not have standalone Lagrangian principles can be seen as the reason why we need Gupta-Bleuler and similar business.