A Lagrangian constraint analysis of first order classical field theories with an application to gravity

We present a method that is optimized to explicitly obtain all the constraints and thereby count the propagating degrees of freedom in (almost all) manifestly first order classical field theories. Our proposal uses as its only inputs a Lagrangian density and the identification of the a priori independent field variables it depends on. This coordinate-dependent, purely Lagrangian approach is complementary to and in perfect agreement with the related vast literature. Besides, generally overlooked technical challenges and problems derived from an incomplete analysis are addressed in detail. The theoretical framework is minutely illustrated in the Maxwell, Proca and Palatini theories for all finite $d\geq 2$ spacetime dimensions. Our novel analysis of Palatini gravity constitutes a noteworthy set of results on its own. In particular, its computational simplicity is visible, as compared to previous Hamiltonian studies. We argue for the potential value of both the method and the given examples in the context of generalized Proca and their coupling to gravity. The possibilities of the method are not exhausted by this concrete proposal.


Introduction
It is hard to overemphasize the importance of field theory in high energy physics. Suffice it to recall that each and every of the fundamental interactions we are aware of as of yet -the gravitational, electromagnetic, strong and weak interactions-are described in terms of fields. Correspondingly, their dynamics are studied by means of field theory. Most often, this is done by writing a Lagrangian (or a Hamiltonian) density that is a real smooth function of the field components (and their conjugate momenta) and that is then subjected to the principle of stationary action. It is customary to encounter the situation where not all of the a priori independent quantities -field components and/or conjugate momenta-are conferred a dynamical evolution through the equations of motion. In such a case, the field theory is said to be singular or constrained. For instance, it is well-known that all gauge theories are singular.
In this work, we focus on singular classical field theories that are manifestly first order and analyze them employing exclusively the Lagrangian formalism. Non-singular theories are also in (trivial) reach. Throughout the paper, manifest first order shall stand for a Lagrangian that depends only on the field variables and their first derivatives. This implies the equations of motion are guaranteed to be second order at most. Within this framework, we present a systematic methodology that is optimized to determine the number of field components that do propagate, which we denominate physical/propagating modes/degrees of freedom. To do so, we explicitly obtain the constraints: specific functional relations among the field variables and their time derivatives that avoid the propagation of the remaining field components. Our approach is complementary to the similarly aimed procedures in [1][2][3] and is markedly distinct from, yet equivalent to, that in [4].
Apart from the intrinsic relevance of understanding and characterizing the constraint structure of those theories satisfying our postulates, an ulterior motivation for this investigation is to pave the way towards a consistent theory building principle. Indeed, theoretical physics is currently in need of new fundamental and effective field theories that are capable of accounting for experimental data; the strong CP problem, neutrino masses and the nature of the dark sector, to mention but a few of the most relevant examples. A recurrent and challenging obstacle in the development of well-defined field theories consists in guaranteeing the correct number of physical modes. In this context, most effort is devoted to avoiding the propagation of Ostrograski instabilities [5] -additional unphysical degrees of freedom, which we shall denote ghosts for short. The general problem is delineated in [6] and numerous realizations of this idea can be found, e.g. [7]. However, it is equally important to ensure the theory is not overconstrained, i.e. there are fewer than required propagating modes. Our subsequent prescription provides a solid footing to this (double) end and is presented in a comprehensible and ready to be used manner, with the goal of being useful to communities such as, but not limited to, theoretical cosmology and black hole physics. We describe how to convert the analytical tool here exposed into a constructive one, but the concrete realization of this idea is postponed to future investigations.
A specific materialization of the preceding general discussion (and the one we later on employ to ground our conversion proposal) is as follows. We recall that an earlier version of the method here augmented and refined already allowed for the development of the most general non-linear multi-vector field theory over four-dimensional flat spacetime: the Maxwell-Proca theory [8,9]. There, the inclusion of a dynamical gravitational field was beyond scope. The present work provides a sound footing for the study of singular field theories defined over curved backgrounds. Thus, it paves the way for the ghost-free coupling of Maxwell-Proca to gravity.
Bearing in mind the above future objective and in order to clarify the formal presentation of the method, we (re)analyze the simplest spin one and two theories by means of our proposed procedure: Maxwell, Proca and Einstein's gravity. While the former two are manifestly first order, the latter is not. Indeed, gravity, cast in the Einstein-Hilbert way, is a second order Lagrangian for the metric, up to a non-covariant boundary term. As such, it exceeds the domain of applicability of our approach. Favorably, this property can be circumvented taking advantage of the deluge of reformulations available for the theory. Among them, we single out the Palatini formalism -see [10] for a historical overview-, which considers the metric and the affine connection as a priori independent fields.
Our determination of the explicit constraints present in Palatini, while not yielding novel information about the theory, conforms a remarkable piece of work. Not only it is carried out minutely and can be readily seen to be computationally easier and shorter than the previously performed Hamiltonian studies, e.g. [11][12][13][14]. It also provides the basis for a consistent inclusion of matter fields. As such, we regard this comprehensive analysis as an intrinsically valuable result.
Organization of the paper. In the following section 2, we introduce the Lagrangian methodology we shall use throughout the paper. Our approach is complementary to the existing literature. In particular, it is equivalent to the recent proposal in [4], as argued and exemplified in section 5.1.
We proceed to employ it to analyze various well-known theories: Maxwell electromagnetism, together with the (hard) Proca action in section 3 and the Palatini formulation of gravity in section 4. Their study is cornerstone to understand the Maxwell-Proca theory [8,9] and paves the way to its consistent coupling to gravity. This is discussed in section 5.2.
We conclude in section 6, restating the instances when our method is most convenient and emphasizing two crucial aspects that are sometimes overlooked.

Conventions.
We work on a d-dimensional spacetime manifold M of the topology M ∼ = R × Σ. Namely, we assume M admits a foliation along a time-like direction. This is true for all (pseudo-)Riemannian manifolds. For simplicity, we consider Σ has no boundary. The dimension d is taken to be arbitrary but finite, with the lower bound d ≥ 2. Spacetime indices are denoted by the Greek letters (µ, ν, . . .) and raised/lowered with the metric g µν and its inverse g µν . We employ the standard short-hand notation ∂ µ := ∂ ∂x µ , where x µ := (x 0 , x 1 , . . . , x d−1 ) ≡ (x 0 , x i ), with i = 1, 2, . . . , d − 1, are spacetime local coordinates, naturally adapted to the foliation R × Σ. The dot stands for derivation with respect to time, so that for local functions f : M → R, we writeḟ := ∂ 0 f andf := ∂ 2 0 f . Brackets indicating symmetrization and antisymmetrization of indices are defined as T (µν) := (T µν + T νµ )/2 and T [µν] := (T µν − T νµ )/2, respectively. As is customary, summation over repeated indices should be understood at all times.

Exposition of the method
We begin by putting forward a coordinate-dependent, i.e. non-geometrical, Lagrangian approach to obtain all the constraints present in a manifestly first order classical field theory. Needless to say, there exists a vast literature on the topic: some standard references are [15]; but for its elegance and concision, we particularly recommend [16]. This section serves us to fix the notation used throughout the paper and provide a self-contained derivation of all our results. We stress that, although the method is not new per se, we are not aware of any reference where this material is comprehensively presented in a ready to be used manner and keeping the technicalities at a bare minimum, as we do here.
Our only assumptions shall be the principle of stationary action and finite reducibility. The first assumption is rather obviously a very mild one, but it is worth noting that this is not an essential requirement; for instance, see [17]. We will explain the second assumption shortly. For the time being, it suffices to note that, to our knowledge, the only known example of a classical field theory (of the kind here considered) not satisfying it is bosonic string field theory, both in its open [18] and closed [19] variants.
Given a Lagrangian density L within the above postulates, our analysis yields the constraint structure characterizing triplet t (N ) := (l, g, e). (2.1) Here, N is the number of a priori independent field variables in terms of which L is written. As such, N is equal to the dimension of the theory's configuration space, which we shall shortly introduce. The other numbers l, g and e are defined below. On shell, we obtain l: the total number of functionally independent Lagrangian constraints. Our analysis elaborates on the iterative algorithm presented in [1] and employed in appendix A of [8]. It is the suitable generalization to field theory of the coordinate-dependent method used in [2] for particle systems, which is in turn based on [20]. The non-trivial geometric extension to field theory of [2] was carried out in [3], where the discussion was extended to the treatment of off shell constraints as well. Thus, our discussion is complementary to all these references [1][2][3].
Off shell, we shall obtain g and e: the number of gauge identities and effective gauge parameters, respectively. Gauge identities are to be understood in the usual sense, as (differential) relations between certain functional variations of the action that identically vanish. By effective gauge parameters we mean the number of independent gauge parameters plus their successive time derivatives that explicitly appear in the gauge transformations. We determine g and e for theories where the gauge transformations are known a priori and provide suitable references that deal with the treatment of theories where the gauge transformations are unknown beforehand. Notice that knowledge of the gauge transformations for the field theory is not a necessary assumption, unlike the principle of stationary action and finite reducibility. However, this information considerably shortens the analysis and, being a feature of all the theories we shall explicitly consider, we have opted for only developing in detail such case.
Given the triplet t (N ) , the physical degrees of freedom n dof in the theory under study can be counted, employing the result derived in [2]: We will refer to (2.2) as the master formula, the way the authors of [2] themselves do. The remarkable feature about the previous counting is that it is purely Lagrangian, as opposed to the usually employed Hamiltonian formula attributed to Dirac. Here, (N 1 , N 2 ) denote the number of first and second class constraints, respectively. As a reminder, first (second) class constraints are those which do (not) have a weakly vanishing Poisson bracket with all of the constraints present in a given theory. Needless to say, the proven equivalence between the Lagrangian and Hamiltonian formulations of classical theories [1,21] is a most celebrated body of work. The two given prescriptions for the degree of freedom count in (2.2) and (2.3) are a particular materialization of this equivalence, which was further exploited in [2] to develop a one-to-one mapping between the Lagrangian parameters (l, g, e) and their Hamiltonian counterparts: stands for the number of so-called primary first class constraints, those first class constraints that hold true off shell. Using this information, the triplet t (N ) defined in (2.1) can be readily seen to admit the following equivalent Hamiltonian parametrization: An important comment is in order here. Our subsequently proposed Lagrangian approach to determine t (N ) does not guarantee n dof ∈ N ∪ {0}. This means that, even though all l, g and e in (2.2) are integers by definition, their sum need not be an even number. The reason is simple: we put forward an analytical tool, not a mechanism to detect (or even correct) ill posed theories. If, for some Lagrangian density L, a half-integer number of physical degrees of freedom is found upon correctly employing our prescription for t (N ) together with (2.2), then it must be concluded that the theory is unphysical. The (possibly non-trivial) modifications required on L for it to propagate an integer number of physical modes is a question beyond the scope of this manuscript 1 .
For the renowned examples in sections 3 and 4, we shall minutely determine the triplet t (N ) defined in (2.1) and then use (2.2) to explicitly count physical modes. Afterwards, we shall (partially) verify our results by comparing them to a representative subset of the Hamiltonian-based literature via (2.3) and (2.4). Additionally, the examples of section 3 shall be worked out in two different (but dynamically equivalent) Lagrangian formulations, based on distinct values N and N N N = N of the dimension of the configuration space. We will then see that, even though the constraint structure characterizing triplets don't coincide, the number of propagating modes n dof does match for both descriptions: t (N ) := (l, g, e) = t (N N N ) := (l l l, g g g, e e e), N − 1 2 (l + g + e) = n dof = N N N − 1 2 (l l l + g g g + e e e).
(2.6) This is because n dof is a physical observable, while (N, l, g, e) are not. Obviously, the same situation arises in the Hamiltonian picture as well, which we briefly illustrate at the end of section 4. In the following, we explain how to obtain the constraint structure characterizing triplet t (N ) in (2.1).

On shell Lagrangian constraints
Let C be the configuration space of a classical field theory. As usual, we take C to be a differentiable Banach manifold whose points are labeled by N real field variables Q A : We stress that A comprises all possible discrete indices that the real field variables have. For instance, if one considers Yang-Mills theory, A consists of both spacetime indices and color indices. If one wishes to entertain complex Yang-Mills, then the real and imaginary parts of each and every Yang-Mills field component must be counted separately in A. So, for SU (2) complex Yang-Mills theory in four spacetime dimensions, we would have that N = 2(4 · 3) = 24. Notice that Q A are real smooth functions of spacetime Q A = Q A (x µ ), but we will suppress this dependence all along, so as to alleviate notation. Thus, our notation matches that in [3] and leaves out the spacetime argument compared to the condensed notation introduced by DeWitt in [22] and extensively used in the literature, e.g. [23]. Then, T C is the tangent bundle of C, which is spanned by {Q A ,Q A }. We refer to (Q A ,Q A ,Q A ) as the generalized coordinates, velocities and accelerations of the theory, respectively. As already stated and common to most field theories, we assume that the dynamics are derivable from a principle of stationary action. In other words, the Euler-Lagrange equations E A ! = 0 for the field theory follow from the requirement that the action functional remains stationary under arbitrary functional variations δQ A = δQ A (x 0 , x i ) that vanish at times t 1 and t 2 on the spatial slice Σ: with δQ A (t 1 , x i ) = 0 = δQ A (t 2 , x i ). The above variational derivative is defined as where the latter equality is the on shell demand. This on shell requirement commences the iterative algorithm we shall employ to determine the Lagrangian constraints present in the theory. Here, L = L[Q A ] is the Lagrangian density. Observe that we have already restricted attention to manifestly first order field theories, i.e. we consider L depends only on Q A and its first derivatives ∂ µ Q A . The study of higher order field theories 2 -where L explicitly depends on ∂ n µ Q A , with n ≥ 2-lies beyond the scope of our 2 One may be tempted to evade the higher order character of a theory via the Ostrogradski prescription, i.e introducing additional generalized coordinates in a manner that results in a manifestly first order Lagrangian density. Such alteration of T C must be compensated through the inclusion of Lagrange multipliers that preserve the equivalence to the original setup. To do so consistently, one needs to either verify the so-called Ostrogradsky non-singularity condition or exploit alternative methods, as detailed in [24]. In view of these non-trivial subtleties, we restrict ourselves to the study of manifestly first order theories. present investigations. We omit the possible dependence of L on non-dynamical field variables, such as the spacetime metric in any special relativistic theory. The said dependence can be easily incorporated to our analysis, but it does not arise in the theories we discuss in this work. An important remark on notation follows. As introduced in (2.7), Q A is an ordered set of a priori independent field variables; it is neither a row nor a column vector. The same is true for E A in (2.10): this is the ordered set of Euler-Lagrange equations for the Q A field variables; not a vector. We have opted for a notation where the set indices are always assigned the same position when ascribed to a certain ordered set (for instance, upper position for the field variables Q A and lower position for the Euler-Lagrange equations E A ). The assignation is such that the Einstein summation convention employed throughout the paper is apparent. The only quantities that will show up in this section which have a definite character within matrix calculus are the following. The various Hessians, their Moore-Penrose pseudo-inverses and the Jacobians are all matrices. The null vectors of the Hessians are row vectors. Their transposed column vectors also show up. The row or column character of the ordered sets is then straightforwardly fixed according to dimensional analysis in all formulae.
As a practical starting point for our iterative method, it is convenient to recast the Euler-Lagrange equations (2.10) in the form where we have defined the so-called primary Hessian W AB := ∂Ȧ∂ḂL, as well as (2.12) To alleviate notation, we have introduced the following short-hands: , which we shall extensively employ henceforth. We focus on singular (or constrained) field theories next 3 . That is, we look at field theories described by a Lagrangian density whose primary Hessian has a vanishing determinant det(W AB ) = 0. This means that the rank of W AB (the number of linearly independent rows or columns) is not equal to its dimension N ; instead, it is reduced.
By definition it follows that, for singular Lagrangians, the N number of Euler-Lagrange equations in (2.11) can be split into two types. First, primary equations of motion: these are the R 1 := rank(W AB ) number of on shell second order differential equations that explicitly involve the generalized accelerations Q A . Second, primary Lagrangian constraints: these are the M 1 := dim(W AB ) − rank(W AB ) = N − R 1 number of on shell relations between the generalized coordinates Q A and their generalized velocitiesQ A . We stress an explicit dependence onQ A (Q A ) is not necessary for the primary Lagrangian constraints, they can be relations between the Q A 's (Q A 's) only. Consistency requires that these constraints are preserved under time evolution.
In the following, we obtain the said constraints and ensure the consistency of the field theory by means of an iterative algorithm. We refer to each iteration in the algorithm as a stage. In every stage, the above specified notions of equations of motion and Lagrangian constraints will arise. The algorithm closes when the preservation under time evolution of all Lagrangian constraints is guaranteed. Equivalently, when all n-th stage Lagrangian constraints are stable, for some finite integer n ≥ 2. An n-th stage Lagrangian constraint is said to be stable if its time derivative does not lead to a new (i.e. functionally independent) Lagrangian constraint in the subsequent (n + 1)-th stage. Below, we explain in detail the different manners in which the necessary stability of the functionally independent Lagrangian constraints may manifest itself.

Primary stage.
In order to determine the subset of M 1 number of primary Lagrangian constraints out of the set of all N number of Euler-Lagrange equations in (2.11), we first introduce a set of M 1 number of linearly independent null vectors γ I associated to the primary Hessian W AB : (2.14) We require that these form an orthonormal basis of the kernel of W AB , which amounts to imposing the normalization condition with T denoting the transpose operation. We stress that, even though in all the examples considered in sections 3 and 4 we have chosen null vectors that are constant, this is not a required feature for our formalism. Rather, this is just a possible choice in all the given examples that has been opted for due to its computational convenience. Only the normalization (2.15) is an essential requirement for the null vectors.
In full generality, the null vectors of all stages can have an explicit dependence on the field variables Q A and their first derivatives ∂ µ Q A . Then, the M 1 primary Lagrangian constraints are obtained by contracting the Euler-Lagrange equations E A in (2.11) with the above null vectors 4 . Namely, by performing the contraction with γ I : Notice that the last equality is a direct consequence of the on shell demand in (2.10) or equivalently in (2.11). Hence, the primary Lagrangian constraints are on shell constraints by definition. One can also see this through equivalence to the more familiar Hamiltonian analysis. It is common knowledge, e.g. [25], that primary Lagrangian constraints relate to secondary constraints in the Hamiltonian framework, which are on shell constraints by definition. The primary Lagrangian constraints in (2.16) need not be functionally independent from each other 5 . When they are, the field theory is said to be irreducible at the primary stage. Otherwise, the theory is reducible at the primary stage. Before we carry on, we must restrict attention to the functionally independent primary Lagrangian constraints ϕ I ! = 0, where I = 1, 2, . . . M 1 ≤ M 1 . Their number is given by M 1 = rank(J IΛ ), where the Jacobian matrix J IΛ is defined as This test can be easily related to the standard Hamiltonian framework: it is the pullback of the phase space regularity conditions in [26]. For the theories we are concerned with in this work, we verify M 1 = M 1 . Hence, all of the primary Lagrangian constraints in (2.16) must be considered in the following 6 . The vanishing of all the functionally independent primary Lagrangian constraints defines the so-called primary constraint surface T C 1 , which is a subspace of the moduli space T C 0 of the field theory: For brevity, we write Equalities that hold true in T C 1 (and not in the entire of the moduli space) shall be denoted ≈ 1 and referred to as primary weak equalities.
As previously noted, consistency requires us to not only enforce the primary Lagrangian constraints (2.19), but also to ensure that these are preserved under time evolution. Explicitly, E J :=φ J ! ≈ 1 0. This requirement starts the second iteration in the algorithm.

Secondary stage.
The freshly introduced demands E J ! ≈ 1 0 7 are known as the secondary Euler-Lagrange equations. In order to split them into secondary equations of motion and secondary Lagrangian constraints, it is convenient to write them as where we have defined We point out that, in obtaining this expressions, we have employed the on shell statement (2.11), so as to eliminate from (2.20) as much dependence on the generalized accelerationsQ A as possible 8 . Here, W IJ is the so-called secondary Hessian and the auxiliary matrix M AB is the Moore-Penrose pseudo-inverse (as detailed in [27]) of the primary Hessian. The latter is ensured to always exist and be unique. Its defining relations are 9 If M 1 < M1 and the functionally independent constraints are not straightforwardly identifiable, more work is required.
Indeed, there exists an iterative algorithm to extract the functionally independent subset of Lagrangian constraints from (2.16). This is explained in section IID of [3] and subsequently exemplified. When the said algorithm requires a(n) finite (infinite) number of iterations, we face a(n) finitely (infinitely) reducible theory. As already pointed out, the procedure here described requires, at the very least, the closure of the reducibility algorithm to proceed. Thus, infinitely reducible theories cannot be studied with the present formalism. We restate bosonic string field theory [18,19] is the only physically relevant example of an infinitely reducible theory we are aware of. 7 For clarity, we will use a notation where tilde quantities belong to the secondary stage and hat quantities pertain to the tertiary stage. This will be particularly helpful in section 4.1. 8 In the equivalent and more familiar Hamiltonian approach, this corresponds to solving as many generalized velocities as possible in terms of generalized coordinates and conjugate momenta:Q A =Q A (Q A , ΠA). 9 In [1], the first relation is referred to as completeness relation. There, both equations in (2.22) are further used to obtain the explicit form of the functionally independent secondary equations of motion. Unlike at the primary stage, functional independence is not guaranteed by construction. As in the first iteration earlier on, our interest lies in the form of the secondary Lagrangian constraints exclusively.
To gain some more intuition into M AB , we note that it constitutes a generalization of the standard matrix inverse. It is introduced so that W AB M BC and M AB W BC are orthogonal projections onto the image of W AB and M AB , respectively. For regular square matrices, the Moore-Penrose pseudo-inverse is equivalent to the standard matrix inverse: M = W −1 iff det(W ) = 0. If rank( W IJ ) = dim( W IJ ) = M 1 , no secondary Lagrangian constraints arise and thus the primary Lagrangian constraints are stable. In this case, we say that the consistency of the primary Lagrangian constraints (2.19) under time evolution is dynamically ensured, by a set of M 1 (necessarily functionally independent) secondary equations of motion E J = E J (Q A ). As a result, the total number of functionally independent Lagrangian constraints present in such field theories is l = M 1 . However, this is not what happens in the theories of our interest.
Generically, the rank of the secondary Hessian is smaller than its dimension. Consequently, M 2 := dim( W IJ ) − rank( W IJ ) of the equations in (2.20) are secondary Lagrangian constraints, whose consistency under time evolution must be ensured. This is done exactly as in the primary stage before. In other words, the analysis from equation (2.14) onwards is to be repeated.
In details, the M 2 number of linearly independent null vectors γ R of the secondary Hessian must be obtained: and chosen so that the normalization condition is satisfied. Then, these must be contracted with the secondary Euler-Lagrange equations in (2.20) to yield the secondary Lagrangian constraints in the theory, If the secondary Lagrangian constraints vanish when evaluated on the first constraint surface ϕ R ≈ 1 0, then the total number of functionally independent Lagrangian constraints is l = M 1 . Again, this is not what happens in (all of) the theories of our interest. As a consequence, we must proceed with the algorithm. First, we need to obtain the (subset of) ϕ R 's which are functionally independent among themselves when evaluated on the first constraint surface. Their number M 2 ≤ M 2 is given by and X Λ was introduced in (2.17). When M 2 = 0, we verify M 2 = M 2 for the theories we shall consider -so that they are irreducible theories at the secondary stage. Thus, all secondary Lagrangian constraints in (2.25) must be considered subsequently 10 .
The vanishing of the functionally independent secondary Lagrangian constraints defines the secondary constraint surface T C 2 ⊂ T C 1 ; which we write as ϕ R ! :≈ 2 0. Equalities holding true in T C 2 shall be denoted ≈ 2 10 When 0 < M 2 < M2, the iterative algorithm referenced in footnote 6 must be employed to extract the functionally independent secondary Lagrangian constraints from (2.25). and referred to as secondary weak equalities. It should be obvious that the secondary Lagrangian constraints are on shell constraints by definition.
Tertiary stage. Let W RS := ( γ R ) I (γ I ) A ∂Ȧ ϕ S be the tertiary Hessian. When the tertiary Hessian's rank does not match its dimension, the consistency under time evolution of M 3 := dim( W RS ) − rank( W RS ) number of the functionally independent secondary Lagrangian constraints is not (dynamically) guaranteed. Instead, it must be enforced through a third iteration of the just described procedure. We stress that it is essential to close the iterative algorithm in order to find the correct number l of functionally independent Lagrangian constraints.
For completeness, we provide the explicit expressions for all relevant quantities at some arbitrary stage of the algorithm in appendix A. These have not appeared in the literature, as far as we know.
Closure of the algorithm. In full generality and as already anticipated, our algorithm stops when all functionally independent Lagrangian constraints have been stabilized. This can happen in either of the following different manners: i Dynamical closure.
Firstly, it may happen when M n := dim(W (n) ) − rank(W (n) ) = 0 for some n-th stage Hessian W (n) , with n ≥ 2. This implies that no Lagrangian constraints arise at the n-th stage, since in this case W (n) has full rank and hence admits no null vector. Here, the consistency under time evolution of the previous stage's functionally independent Lagrangian constraints ϕ (n−1) is dynamically ensured, i.e. through the (necessarily functionally independent) n-th stage equations of motion. In other words, the functionally independent ϕ (n−1) 's are stable. This closure of the algorithm is exemplified in section 3.2.
ii Non-dynamical closure. Secondly, it may happen when M n > 0, but M n = 0, again with n ≥ 2. This implies that the n-th stage functionally independent Lagrangian constraints ϕ (n) 's do not define a new constraint surface, so that T C n ≡ T C n−1 . We differentiate two algebraically distinct scenarios: iia The ϕ (n) 's vanish identically in the (n − 1)-th constraint surface: ϕ (n) ∼ ∼ ∼ iib The ϕ (n) 's functionally depend on the (n − 1)-th stage functionally independent Lagrangian constraints. Schematically, In all the detailed cases, the total number of functionally independent Lagrangian constraints is given by where M a counts the number of functionally independent a-th stage Lagrangian constraints and n ≥ 2. We are not aware of any physically relevant example of a field theory where n is infinite.

Noteworthy considerations.
We restate that it is of utmost importance to close the iterative procedure in order to determine l. If the algorithm is not closed (only some or none of the constraints are stabilized), one can only give a lower bound on l. While this may be enough to ensure the absence of Ostrogradsky instabilities [5] in the field theory, it is insufficient to guarantee the propagation of a definite number of degrees of freedom. In such case, one can only infer an upper bound on n dof . This observation is further discussed and exemplified in section 6.
We also point out that, in general, the different stabilizations of the functionally independent Lagrangian constraints that we listed are all present in a given field theory. Namely, some functionally independent Lagrangian constraints in the theory are stabilized dynamically, while others are stabilized non-dynamically.

This is indeed what happens in our examples of sections 3.3 and 4.
Besides, we warn the readers against deceiving themselves regarding the ease of the exposed iterative algorithm. Even though our methodology is sound and rigorous and its logic is easy to follow, there can be no misapprehension as to the algebraic complexity of its implementation in concrete theories, most significantly those involving gravity. At last, we remark that the algorithm just exposed does not break covariance. Namely, if a field theory within our postulates is covariant, its study under the outlined iterative methodology will preserve this feature. Nonetheless, a suitable space and time decomposition of the a priori independent field variables and an evaluation of the Lagrangian constraints in the various constraint surfaces will generically break manifest covariance. This should not be confused with the loss of covariance.

Off shell gauge identities
We now obtain g and e, the two remaining numbers in the triplet t (N ) defined in (2.1) of our interest. To begin with, we notice that in the principle of stationary action (2.9), we have so far only considered that δS = 0 follows from the E A piece. However, δS = 0 may also follow from the δQ A piece. Subsequently, we briefly review the latter scenario: how the vanishing of δS may be a consequence of off shell identities stemming from a strict symmetry of the action. This kind of symmetry -gauge invariance-is only manifest through specific field variations δ θ Q A , in contrast to our previous consideration in section 2.1 of arbitrary δQ A 's. Correspondingly, we will differentiate between δ θ S and δS as well.
There are different methods to obtain the said off shell identities, but it is not our goal to provide an overview of them here. Our subsequent discussion summarizes and employs the approach put forward in [28] and later on adapted to exhibit manifest covariance in [14]. This adaptation makes it straightforward to apply [28] to any manifestly first order classical field theory, which is our framework.
Consider the field transformations Q A → Q A + δ θ Q A . Let the changes δ θ Q A be of the form where n ∈ N ∪ {0}, β is an (possibly collective) index that is to be summed over and the θ β 's and Ω β A 's are known as the gauge parameters and gauge generators of the transformation, respectively. The the θ β 's are real smooth functions of the spacetime coordinates x µ , while the Ω β A 's are defined in T C and as such are real smooth functions of (Q A ,Q A ). The former are unspecified, while the latter are to be determined. Introducing the above in (2.9) and operating, one finds that If, under the field variations (2.28) for some Ω β A 's, the action remains invariant δ θ S ≡ 0, then we have that holds true off shell (i.e. without making use of E A ! = 0). In such a case, (2.28) and (2.30) are known as the gauge transformations and gauge identities in the theory, respectively.
Given (2.28), g is equal to the number of different θ parameters there present. Equivalently, g is the number of linearly independent gauge identities (2.30). On the other hand, e is equal to the total number of distinct parameters plus their successive time derivatives (θ,θ,θ, . . .) that appear in (2.28). Obviously, e ≥ g.
The recursive construction of the gauge generators Ω β A has been a subject of vivid interest for decades. The approach in [29] is perhaps the most befitting to our own exposition, requiring only a suitable adaptation from particle systems to manifestly first order field theories that is devoid of conceptual subtleties. We shall not present the corresponding discussion here because, for the theories at hand, the explicit form of the gauge transformations is already known. This a priori knowledge allows us to effortlessly infer the generators Ω β A in all the subsequent examples. We stress that the determination of g and e is possible and has been made systematic in theories for which the gauge transformations are unknown from the onset. The calculations in such theories are more involved, but there is no theoretical obstacle that has to be overcome. To illustrate this point, the reader can consult [28] for the explicit derivation of the gauge generators in Yang-Mills theory and both the metric and Palatini formulations of General Relativity, by means of the formalism put forward in [29].
For the ease of the reader, we have schematically depicted the main line of reasoning behind this section 2 in figure 1.

Simple examples: vector field theories
This section is devoted to the study of some of the constraint structure characterizing triplets t (N ) that are possible for the theories (within the framework of section 2) describing the dynamics of a single vector field. Recall there are only two distinct types of vector fields that one can entertain classically: massless and massive. For simplicity, we will restrict to real Abelian vector fields and focus on their most elementary actions: Maxwell electromagnetism and the (hard) Proca theory, respectively. We shall consider two equivalent formulations of each of these theories, based on different numbers N and N N N = N of a priori independent field variables. Our forthcoming detailed analyses are based on the purely Lagrangian method described in the previous section 2 and thus serve to illustrate it.
Besides and as we shall explain in section 5.2, our forthcoming elementary calculations turn out to be enough to understand the complete set of manifestly first order (self-)interactions among an arbitrary δ θ S ≡ 0 (off shell) δS ! = 0 (on shell) Figure 1: Schematics of section 2. Here, (eqns., Lag. consts., f.i., num.) stand for equations, Lagrangian constraints, functionally independent and number, respectively. The computational challenge of the steps relating Lagrangian constraints to functionally independent Lagrangian constraints (represented with a double arrow), as well as the relevance of closing the iterative algorithm are further discussed in section 6. number of both Maxwell and (generalized) Proca [30] fields in four-dimensional flat spacetime [8,9]. This hints to the convenience of the proposed method, compared to other possible approaches; a point that shall be reinforced in the more elaborate examples of the next section 4 and discussed in the concluding section 6.
In the remaining of this section, we shall work on d-dimensional Minkowski spacetime, still for finite d ≥ 2. We will choose Cartesian coordinates with the mostly positive signature, so that g µν = η µν = diag(−1, 1, 1, . . . , 1). Subsequently, all spacetime indices shall be raised/lowered by η µν and its inverse η µν .

Maxwell electromagnetism
This renowned manifestly first order singular field theory describes an Abelian massless vector field and its linear interactions with sources in terms of N = d number of a priori independent field variables. As already stated, we take the Maxwell vector field (which we denote A µ ) to be real and consider the particularly simple case when there are no sources.

Lagrangian constraints.
The canonically normalized Lagrangian density of sourceless classical electromagnetism is The components of the Maxwell field constitute the generalized coordinates for this theory: Q A = {A µ }, so that A = 1, 2, . . . , d = dim(C) ≡ N , as already announced. As is well-known and can be easily calculated by means of (2.10), the Euler-Lagrange equations following from (3.1) are If we decompose the Maxwell field into its space and time components can be conveniently rewritten as where sum over repeated indices is to be understood and we have been careful to lower all indices with the flat metric η µν . It is then easy to see that the primary Hessian following from . and therefore manifestly possesses the symmetry dictated by its very definition: W AB = W BA . Further, its Moore-Penrose pseudo-inverse is given by Since the primary Hessian takes such an uncomplicated form, it readily follows that R 1 = 3 and thus M 1 = 1(= M 1 ) in this case. A convenient choice for the null vector of W AB amounts to (γ 1 ) A = δ 1 A . Then, the one and only primary Lagrangian constraint for the theory can be effortlessly calculated to take the explicit form This is the familiar Gauss law, telling us that, in the absence of sources, the electric field is divergenceless. Note that this is an on shell statement by construction. The Gauss law constraint straightforwardly yields a vanishing secondary Hessian W 11 ≡ 0, so that M 2 = M 1 = 1 and we choose ( γ 1 ) 1 = 1. With all this information, it is a matter of easy algebra to find the only secondary Lagrangian constraint: Therefore, M 2 = 0 and the end of the iterative algorithm is signalled according to the non-dynamical prescription in case iia. We have thus found that the total number of Lagrangian constraints for Maxwell electromagnetism is just l = M 1 + M 2 = 1.

Gauge identities.
Maxwell's theory enjoys an apparent U (1) gauge symmetry. Indeed, under the transformation A µ → A µ +∂ µ θ, the Lagrangian (3.1) remains invariant. Here, θ is the only gauge parameter, while (θ,θ) are the sole two effective gauge parameters present in the fields' transformation. Consequently, we have that g = 1 and e = 2. For completeness, we point out that the said transformation, when compared to (2.28) immediately allows us to read off the gauge generator of the symmetry. This is (Ω A ) ν = −δ Aν . When combined with the primary Euler-Lagrange equations (3.2) as indicated in (2.29), we can right away verify the off shell gauge identity we counted: = ∂ µ ∂ ν A µν ≡ 0.
Physical degrees of freedom. According to our prior analysis, which shows that the constraint structure of classical electromagnetism in its standard formulation with N = d is t (N ) M = (l = 1, g = 1, e = 2), (3.6) and making use of the master formula (2.2), we count n dof = d − 2 propagating modes. In d = 4, these correspond to the two polarizations of the photon. Exploiting the equalities in (2.4), we check that our counting corresponds to two first class constraints, one primary and one secondary. Therefore, our purely Lagrangian investigation is in perfect agreement with the standard literature, e.g. [31]. It also matches the Hamiltonian definition of the Maxwell field given in [8]: "a real Abelian vector field [. . .] associated with two first class constraints". This latter correspondence will play a role in section 5.2.

The (hard) Proca theory
We turn our attention to the Proca theory next, in the modern formulation of the original proposal in [32]. Namely, we focus on the (manifestly first order) field theory of a real Abelian vector field of mass m in the absence of any source described by N = d a priori independent field variables. The remark (hard) is to avoid ambiguity with respect to the Generalized Proca theory, discussed in section 5.2. We refer to the Proca field as B µ .

Lagrangian constraints.
The Lagrangian density of the said Proca theory is As in the Maxwell case earlier on, the components of the Proca field are the generalized coordinates: We thus see that A = 1, 2, . . . , d = dim(C) ≡ N here as well. The Euler-Lagrange equations following from (3.7) can be easily obtained as indicated in (2.10). The result is At this point, it is straightforward to see that the primary Hessian -and hence also its Moore-Penrose pseudo-inverse-is the same as for the Maxwell theory earlier on. This implies M 1 = 1(= M 1 ) and the associated null vector can again be chosen as (γ 1 ) A = δ 1 A . The primary Lagrangian constraint differs, though: The above once more leads to a vanishing secondary Hessian, so that M 2 = M 1 = 1 and ( γ 1 ) 1 = 1. The secondary Lagrangian constraint in this case takes the form Contrary to the Maxwell theory, (3.10) is obviously not a Lagrangian identity, so the algorithm is not closing here according to the prescription in case iia. Notice as well that ϕ 1 and ϕ 1 are functionally independent from each other, so that we are not in case iib of the general method either. Instead, we have M 2 = M 2 = 1 and we must move on to the tertiary stage.
It is easy to check that the tertiary Hessian following from (3.10) is W 11 = −m 2 . As such, its dimension and rank match (M 3 = 0 = M 3 ) and the algorithm closes according to the dynamical prescription in case i. Namely, the consistency of (3.10) under time evolution is ensured via a tertiary equation of motion and there are no tertiary constraints. As a result, we have obtained l = M 1 + M 2 + M 3 = 2 functionally independent Lagrangian constraints in the (hard) Proca theory.

Gauge identities.
The mass term for the Proca field explicitly breaks the U (1) gauge invariance of Maxwell electromagnetism. In our conventions, this means that there is no field transformation of the form (2.28) that leaves the action invariant. Therefore, there are no off shell identities associated to (3.7) and we have g = 0 = e.
Physical degrees of freedom. Using the (hard) Proca constraint structure for N = d t (N ) obtained before in the master formula (2.2), we count n dof = d − 1 degrees of freedom in the theory. By means of (2.4), it is immediate to certify that this corresponds to two second class constraints; as explicitly shown, for instance, in [33]. As with the Maxwell field before, we thus find agreement with the Proca field's definition given in [8]: "a real Abelian vector field [. . .] associated with two second class constraints". We will further comment on this connection in section 5.2 later on.

The Schwinger-Plebanski reformulation of Maxwell and Proca
In this section, we reanalyze the constraint structures of the above massless and massive vector field theories in a formulation with N N N = N = d a priori degrees of freedom. Specifically, we entertain the reformulation of sourceless classical electromagnetism originally proposed by Schwinger [34] and later on popularized by Plebanski [35] and employ it for the (hard) Proca theory simultaneously. In this setup, the real Abelian (covariant) vector field C µ -be it massless or massive-and its antisymmetric (contravariant) field strength F µν are regarded as independent at the onset: The aim of this section 3.3 is to determine the constraint structure characterizing triplets t P , so as to illustrate in a simple double-example the general claim in (2.6). Namely, these triplets differ from the previously determined ones t A clarifying remark follows. Classical electromagnetism as written in [34] is commonly called the manifestly first order formulation of electrodynamics. This refers to the order of its primary Euler-Lagrange equations, contrarily to our convention here, where the order refers to the Lagrangian density. For us, all examples in sections 3 and 4 are manifestly first order and as such can be investigated by means of the methodology in section 2. In view of this dissonance, we can already anticipate that there will be no primary equations of motion in our subsequent examples. The primary Euler-Lagrange equations, being first order, will not involve the generalized accelerationsQ A and so they will all be primary Lagrangian constraints. Further, this is possible iff the primary Hessian of the theories identically vanishes, as we shall see it does.

Lagrangian constraints.
Inspired by [34], we take the Lagrangian density Solving the latter for F µν and substituting the result into the former, we recover Maxwell's (3.2) or Proca's (3.8) equations of motion, depending on the value of m. Then, we say both formulations, in (3.13) and in (3.1) or (3.7) as pertinent, are dynamically equivalent, as foretold. We proceed to explicitly confirm our predictions. The primary Hessian following from (3.13) vanishes identically W AB ≡ 0, so R 1 = 0 and M 1 = N N N . We can choose its appropriate null vectors as (γ I ) A = δ I A . As a result, the primary Lagrangian constraints coincide with the primary Euler-Lagrange equations. These can be readily seen to be functionally independent among themselves. Consequently, the first constraint surface T C 1 coincides with the moduli space in this case. This set of circumstances can be summarized as or simply as M 1 = M 1 = N N N . Notice that W AB ≡ 0 immediately makes its Moore-Penrose pseudo-inverse vanish as well: M AB = 0. We encounter this same situation of a zero primary Hessian in both of the theories analyzed in section 4. We briefly depart from the application of the iterative algorithm in order to introduce an extremely useful notation that will be recurrent from now on. We wish to be able to refer to each kind of field variables in (3.12) individually. To this aim, we shall henceforth understand that the index A therein decomposes into two distinct sets of indices A ≡ A 1 A 2 , the first referring to the type of field variable and the second to the spacetime structure of each type of field variable. In this way, A 1 = 1, 2, . . . , 4 and we have Observe that we have employed the symbol [·] to visually split the A 1 index from the A 2 one. Back to the algorithm and putting into practice the above notation, we write the primary Lagrangian constraints as

as required by definition.
We go on to the secondary stage next. The secondary Hessian W IJ = ∂İ ϕ J , can be portrayed in our recently introduced notation as follows: where, for each entry of the secondary Hessian, we have placed the space-like tensorial indices of the field variables (primary Lagrangian constraints) labeled by I (J) to the left (right). A few explicit examples that should clarify our notation are 18 The only non-zero components in (3.18) are 20) which lead to a simple secondary Moore-Penrose pseudo-inverse M AB with non-zero elements j [ This corresponds to the transpose of (3.18). It is easy to see that R 2 := rank( W IJ ) = 2(d−1), which in turn implies that M 2 = (d 2 − 3d + 4)/2. We choose the suitably normalized linearly independent null vectors for (3.18) as ( γ R ) I = δ R I . The above results can be employed to determine the functionally independent secondary Lagrangian constraints ϕ R ! :≈ 2 0. First, we calculate ϕ R = ( γ R ) Iφ I and obtain where again the antisymmetry property [ ϕ 2 ] ij = −[ ϕ 2 ] ji required by definition is apparent. Evaluation on the first constraint surface then gives which respects the noted symmetry, as it must. In more detail, the evaluation has been carried out as follows. By setting to zero all ϕ's in (3.17), solving for (Ḟ i ,Ċ i ) and plugging the resulting expressions into (3.21). Next, we need to select only the functionally independent secondary constraints. It is obvious that the mass m plays a crucial role here, as could easily be anticipated in view of our results in the previous sections 3. 1 We turn to the time evolution of the functionally independent secondary constraints, i.e. we commence the tertiary stage. The tertiary Hessian can be succinctly expressed as where we have made use of the same notation as in (3.18) earlier on, so that Hence, the tertiary Hessian has full rank R 3 = M 2 and consequently M 3 = 0. Observe that this is true for both the m = 0 and the m > 0 cases. The functionally independent secondary constraints' consistency under time evolution is at this point dynamically ensured and the algorithm closes according to the prescription in case i. The total number of functionally independent Lagrangian constraints is We see that the mass m gives rise to one more functionally independent Lagrangian constraint, exactly as in the previous sections, where we found that l = 1 for electromagnetism, while l = 2 for the (hard) Proca theory.
Here, the remaining constraints that l l l counts are associated to the field strength F µν , as a result of having promoted it to a set of a priori independent field variables. Notice that there are d(d − 1) number of such supplementary constraints, two times the number of independent components in F ij . This duplicity makes it manifest that these fields are superfluous when describing the dynamics of the theory. In other words, no initial data is needed for them: F ij andḞ ij need not be specified at some initial time t 1 when solving their associated equations of motion. Yet another way to understand this is to map them to the Hamiltonian picture, where they correspond to second class constraints, as we shall shortly see.

Gauge identities.
Consider the following transformations of the field variables: Here, θ is an arbitrary parameter. It can be easily checked that, under the said transformations, the Lagrangian (3.13) remains invariant iff m = 0. Therefore, these are the very same gauge transformations of the massless theory that we noted in section 3.1, while the massive theory does not exhibit any kind of gauge symmetry. Straightforwardly, we count For completeness, we provide the gauge identity and generators for m = 0 next. Comparing (2.28) and (3.28), we can immediately read off the non-zero generators: Notice that here we have dropped the, in this case, single-valued β index from (2.28). Putting together (3.14) and (3.30) as indicated in (2.29), we readily confirm the gauge identity: Physical degrees of freedom.
We have now achieved our goal. Namely, we have shown that the constraint structure characterizing triplet for (3.13) is Substituting the quantities (3.32) into the master formula (2.2), we count

A comprehensive constraint analysis of Palatini theories
In the following, we apply the general framework presented in section 2 to the Palatini action. We split our calculations into the d > 2 and the d = 2 cases, as these are physically distinct theories. As we shall see, the former case is much more algebraically involved than the latter. However, compared to their equivalent Hamiltonian investigations, our Lagrangian approach shall prove much simpler in both instances. For concreteness, we specify our framework to be that of the metric-affine Palatini formulation of General Relativity, ordinarily ascribed to Palatini but firstly suggested by Einstein himself [10,37]. As such, we shall study a manifestly first order formulation of gravity based on N = d(d + 1) 2 /2 number of a priori independent degrees of freedom. Even though alternative manifestly first order formulations do exist, such as the tetradic-Palatini action (for example, see [38] and its recent canonical study [39]), inconvenient subtleties to our aims arise in those frameworks due to their geometric construction. For instance, unlike the metric, vielbeine are not required to be invertible. In such scenario, the strict equivalence between Palatini and Einstein's gravity is lost due to a singular vielbein and, in general, ends up in a dynamical manifestation of torsion [40]. Similar situations might arise in other manifestly first order formulations, like the Barbero-Holst action [41] or BF-like models [42], which happen to enclose the most celebrated Plebanski action. For a complete review on these topics, we refer the interested reader to [43].

Palatini in d > 2
The Palatini action in d > 2 is a well-known (re)formulation of the Einstein-Hilbert action, which is dynamically equivalent to it. This is explained shortly. Most significantly for us, Palatini is a manifestly first order formulation of General Relativity, which treats the spacetime metric g µν = g νµ and the affine connection Γ ρ µν = Γ ρ νµ as a priori independent variables. As such, and unlike Einstein-Hilbert, it readily allows for the application of the methodology introduced in section 2.

Lagrangian constraints.
The Palatini action is of the general form given in (2.8) and its Lagrangian density can be written as [44] Here, the independent variables h µν and G ρ µν are defined exclusively in terms of the spacetime metric and affine connection, respectively: and thus inherit their symmetry properties. The primary Euler-Lagrange equations for h µν and G λ µν following from (4.1) are Notice that these vanishings are on shell statements. Multiplying the second set of field equations by h µν and employing the identity h µν h νρ = δ µ ρ , one finds that Solving (4.4) implies that G ρ µν is fixed (on shell) to be a function of h µν and its first derivatives. The substitution of the resulting expression into (4.1) yields the second order formulation of General Relativity and we say d > 2 Palatini is dynamically equivalent to it.
It is natural and convenient to decompose the variables in (4.2) as follows: (4.5) The explicit form of the Lagrangian (4.1) in terms of the above variables is (4.6) We express the generalized coordinates of the Palatini Lagrangian in (4.6) as Henceforth, we shall employ the notation [·] introduced in section 3.3 for the collective index A above. In particular, see (3.16) and explanations around. This notation shall prove of utmost convenience. For instance, in this way, it is obvious that The primary Hessian following from (4.6) vanishes identically: W AB ≡ 0, as a result of having promoted the affine connection to a set of a priori independent field variables. This parallels the reformulations of classical electromagnetism and the (hard) Proca theory in section 3.3. In passing, we note that the primary Hessian is symmetric W AB = W BA , as it should by definition. Obviously, rank(W AB ) = 0 and we have M 1 = N = d(d + 1) 2 /2. This trivialization of the primary Hessian has a number of direct implications. First, it allows us to straightforwardly pick its suitably normalized null vectors to be (γ I ) A = δ I A . Second, it immediately makes its Moore-Penrose pseudo-inverse vanish as well: M AB = 0. Third, it becomes apparent that the primary Euler-Lagrange equations coincide with the primary Lagrangian constraints. All of these constraints turn out to manifestly be functionally independent from each other in this specific theory. In other words, the moduli space is the primary constraint surface in this case and we have M 1 = M 1 . Thus, exactly as in our examples of section 3.3 before, see (3.15). By means of the notation employed in (3.17), the explicit form of the ϕ I 's is where we have defined (4.14) In (4.13), the only non-zero components are Notice that the secondary Hessian is antisymmetric W IJ = − W JI , as it should by definition. It is easy to see that R 2 := rank( W IJ ) = d(d + 1), thus yielding M 2 = d(d 2 − 1)/2. This means that R 2 number of the functionally independent primary Lagrangian constraints are being dynamically stabilized at the secondary stage, while the remaining M 2 primary Lagrangian constraints are not stable: they lead to secondary Lagrangian constraints, which we proceed to determine.
To this aim, we first choose the suitably normalized linearly independent null vectors associated to the secondary Hessian as In our [·] notation, we have il , (4.17) where we have defined in terms of (4.12) as well as the following quantities: Observe that the appropriate symmetry [ ϕ 3 ] jk i = [ ϕ 3 ] kj i is manifest. To obtain the above, we have first computed ϕ R = ( γ R ) Iφ I . Then, we have evaluated the result in the first constraint surface. In practice, this means that we have substituted (Ġ,Ġ i , . . . , ∂ i h jk ) for their suitable weak expressions in terms of the generalized coordinates Q A , which follow from setting to zero (4.11).
To conclude the secondary stage, we calculate the Moore-Penrose pseudo-inverse of W IJ . It can be easily checked that this is M IJ = −( W IJ ) T .
Next, the consistency under time evolution of the above functionally independent secondary Lagrangian constraints is to be inspected at the tertiary stage. The first step is to calculate the tertiary Hessian W RS = ( γ R ) I ∂İ ϕ S . Employing the same conventions as in (4.13) before, we write where the non-zero components are The tertiary Lagrangian constraints are given by the requirement ϕ U = ( γ U ) R˙ ϕ R ! ≈ 2 0. In this case, the derivation with respect to time is particularly simple and coincides with the naively expected one, so that (4.24) 11 To determine the null vectors of WRS in (4.20), we considered the ansatz ( γU ) R = (a i , a i j , a i jk ), with a i jk = a i kj . Then, the equation ( γU ) R WRS = 0 results in a j i = ca l l δ j i and a j il = ca k kl δ j i , but does not impose any condition on a i . The first equation implies a i j = 0 for i = j. Setting j = l in the second equation yields a j ij = 0, which in turn implies a i jk = 0 for all i, j, k.

25
In our short-hand notation, we find it convenient to express these constraints as follows: where the operator O is defined as (4.26) Recall that ( τ 1 ) i and ( τ 2 ) ij are as introduced in (4.18). Following the procedure described under (4.19), the tertiary Lagrangian constraints in (4.25) can be evaluated on the first constraint surface T C 1 . After tedious algebraic manipulations, the weak tertiary Lagrangian constraints can be written exclusively in terms of the functionally independent secondary constraints as The above is a non-trivial result. Indeed, it becomes increasingly computationally challenging to evaluate Lagrangian constraints on constraint surfaces as one goes to higher stages. We elaborate on this topic and advice on how to handle the evaluations in section 6. Our results in (4.27) must be further evaluated on the second constraint surface T C 2 . Namely, in the subspace of T C 1 defined by the vanishing of (4.17). We thus see that which implies T C 3 ≡ T C 2 and there are no functionally independent tertiary constraints M 3 = 0. Consequently, the algorithm closes non-dynamically, according to case iib. We are finally able to obtain the result of interest from the analysis here presented. The number of functionally independent Lagrangian constraints for the Palatini theory in d > 2, when described in terms of N = d(d + 1) 2 /2 number of a priori independent field variables, is equal to Gauge identities. It is well-known (for instance, see [14]) that the Palatini action corresponding to the Lagrangian density (4.1) remains invariant under the following transformations of its independent variables: h µν → h µν + δ θ h µν and G ρ µν → G ρ µν + δ θ G ρ µν , with 30) where θ µ are the (unspecified) gauge parameters. Notice that the pertinent symmetries δ θ h µν = δ θ h νµ and δ θ G ρ µν = δ θ G ρ νµ are apparent in the precedent expressions. Of course, (4.30) is just the Palatini (re)formulation of the renowned diffeomorphism invariance of the Einstein-Hilbert action. This holds true off shell.
It is easy to see in (4.30) that the gauge parameters θ µ appear explicitly in all the gauge transformations ∀µ. Similarly, we note that the effective gauge parameters (θ µ ,θ µ ,θ µ ) are manifestly present in the gauge transformations ∀µ as well. By definition, it follows that g = d, e = 3d, (4.31) which are the off shell parameters we aimed to obtain in this short analysis. For completeness, we provide the gauge generators and confirm the gauge identities of d > 2 Palatini next. A direct comparison between (2.28) and (4.30), allows us to rewrite the latter as where we have introduced a bracket (·) to visually split the (in general collective) indices β and A for latter convenience. In view of these transformations, the gauge generators can easily be identified to be Combining (4.3) with the above as prescribed in (2.29) and working through, the gauge identities are obtained: Physical degrees of freedom.
Putting everything together, we can finally count the number of propagating modes present in the theory. Namely, employing (4.8), (4.29) and (4.31) in the master formula (2.2), we get When d = 4, we have that n dof = 2, corresponding to the two massless tensor's polarizations of the graviton. For d = 3, the widely known triviality is recovered, with no physical degrees of freedom being propagated. Our result is in perfect agreement with the counting performed in [13,14], where a purely Hamiltonian analysis was done. We have thus carried out another (non-trivial) explicit verification of the already noted equivalence between (2.2) and (2.3). This equivalence can be further verified as follows. It is explicitly shown in [13,14] that N 1 = 3d, N 2 = d(d − 1)(d + 2) and N (P) 1 = d for the d > 2 Palatini theory when (4.8) holds true. Substitution of these results in (2.4) readily confirms our own counting in (4.29) and (4.31). Besides, a direct comparison between the calculations in [13,14] and those presented in this section 4.1 unequivocally shows that our purely Lagrangian computation is an algebraically much simpler way to derive (4.36), from which the number of physical modes follows readily.
To sum up, we have derived the constraint structure characterizing triplet t (N ) Pa , with N = d(d + 1) 2 /2, of the Palatini theory in d > 2 dimensions in a purely Lagrangian approach and ratified its equivalence with a representative Hamiltonian analysis performed in the past. Mathematically, t (N ) (4.36) in the Lagrangian picture, while t (N ) in the Hamiltonian side -recall (2.5)-, both of which imply (4.35).

A special case: Palatini in d = 2
General Relativity, in its standard second order formulation, behaves drastically different in two dimensions. Specifically, it can be shown that Namely, the Einstein-Hilbert action is proportional to the Euler characteristic χ of the spacetime manifold M 2 , see e.g. [45]. The above implies that General Relativity is a topological theory in d = 2 and, accordingly, propagates no degrees of freedom; a fact that we shall explicitly verify in the following.
Turning to the Palatini Lagrangian in (4.1) for d = 2, we restate that this is not dynamically equivalent to two-dimensional Einstein's gravity (4.38). To see this, consider its corresponding Euler-Lagrange equations in (4.3). These are valid for d ≥ 2. However, recall that c := (d − 1) −1 , so that c = 1 in two dimensions. In this particular case, it is obvious that (4.4) cannot be solved as we said, i.e. G ρ µν ! = G ρ µν (h µν , ∂ ρ h µν ). As a result, the dynamical equivalence to Einstein's gravity is lost. A more general yet detailed argumentation can be found in [46].
Correspondingly, the dynamics of the two-dimensional Palatini action does not constitute a smooth limit of its higher dimensional counterpart. Namely, the Lagrangian (4.1) in d = 2 does not describe the evolution of the same family of fields as that very same Lagrangian in d > 2: these are two physically different theories. The easiest way to ratify this second inequivalence is to note that the counting of degrees of freedom in (4.35), when we set d = 2, yields a negative number of propagating modes, which is an unphysical result. Thus, a different constraint structure characterizing triplet is then to be expected. We proceed to determine this t (N =9) 2Pa next.

Lagrangian constraints.
As a starting point, we express the generalized coordinates of the Palatini theory in d = 2 as 28 in direct analogy to (4.7) earlier on. Next, we compute the first stage quantities associated to the d = 2 version of the Palatini Lagrangian (4.6). One can verify that the set of primary Lagrangian constraints thus obtained matches the consistent d = 2 evaluation of those for a generic dimension in (4.10) and (4.11). Comparatively, these two-dimensional constraints have a much simpler form, given by the vanishing of (4.41) The demand that the above be zero constitutes a set of nine scalar primary Lagrangian constraints (M 1 = 9), whose functional independence is rather obvious (M 1 = 9) -and can be ratified through the Jacobian test in (2.17). Therefore, such vanishing defines the primary constraint surface T C 1 of the theory, which coincides with the moduli space due to the primary Hessian being zero, as in the d > 2 case before. In other words, (4.10) holds true here as well.
The progress to the subsequent stage parallels that of the d > 2 case. The secondary Hessian is given by the skew-symmetric constant matrix The Hessian (4.42) has rank R 2 = 6 and so M 2 = 3. This means that six of the primary Lagrangian constraints are dynamically stabilized by the (functionally independent) secondary equations of motion. For the remaining three primary Lagrangian constraints, the algorithm must be pursued. We choose the suitably normalized linearly independent null vectors of (4.42) as ( γ R ) I = δ R+6 I , with R = 1, 2, 3. (4.43) Using (4.41) and (4.43), we obtain the three secondary Lagrangian constraints as the vanishing of Notice that the above are the total time derivatives of [ϕ 7 ] 1 , [ϕ 8 ] 1 1 and [ϕ 9 ] 11 1 in (4.41), respectively. It is easy to check that the secondary Lagrangian constraints are functionally dependent on the primary Lagrangian constraints. Specifically, Therefore, the secondary constraints vanish on T C 1 , and so T C 2 ≡ T C 1 . This in turn implies that there are no functionally independent secondary constraints: M 2 = 0. Here, the algorithm closes non-dynamically, as described in case iib. It follows that the total number of functionally independent Lagrangian constraints is l = M 1 + M 2 = 9. We note that this result does not correspond to setting d = 2 in (4.29).

Gauge identities.
Given the already pointed out inequivalence between the d > 2 and d = 2 Palatini theories, it is not too surprising that the gauge transformations (4.30) preserving the action (4.6) meet a non-smooth limit for d = 2. The argument is more subtle than that of the purely on shell inequivalence; for example, see [11]. We shall touch upon it shortly. It has been proven, e.g. [11,47], that the two-dimensional Palatini action is invariant under the field transformations h µν → h µν + δ θ h µν and G ρ µν → G ρ µν + δ θ G ρ µν , with Here, µν stands for the two-dimensional Levi-Civita symbol (we work with the convention 01 = 1) and θ µν = θ νµ , so there are three arbitrary gauge parameters that characterize the transformation. It is obvious that all the gauge parameters explicitly appear in the gauge transformations. Hence, g = 3. Their first time derivatives also show up, adding to a total number of effective gauge parameters e = 6. We point out that g does not match the value predicted in (4.31) for d = 2. There is a match for e, but this is purely coincidental. These numbers (g = 3, e = 6), in contrast to the naively expected ones (g = 2, e = 6) from the diffeomorphism transformations (4.31) in d > 2 Palatini, reflect the fact that d = 2 Palatini is associated to a comparatively larger symmetry group. Its connection to the d > 2 gauge group is not obvious, but finds its origin in the underlying two-dimensional geometry. Briefly recall the conformal flatness of two-dimensional spacetimes, i.e. g µν = Ω 2 η µν , g µν = Ω −2 η µν , with µ, ν = 0, 1, (4.48) for some conformal factor Ω = Ω(x µ ). Given this property, the variable h µν introduced in (4.2) simplifies to h µν := det(−g µν )g µν = det(−Ω 2 η µν )Ω −2 η µν = η µν . (4.49) Consequently, in the conformal frame, h µν is flat and det(h µν ) = −1, independent of det(g µν ). This latter equality can be expressed as an algebraic constraint: referred to as the metricity condition. We will soon get back to such condition. For a richer discussion on this topic, though, we refer the reader to [12]. In analogy to the higher dimensional case before, we provide the gauge generators and identities of d = 2 Palatini next. Comparing (2.28) to (4.47), we can conveniently rewrite the latter as with the gauge generators readily recognized as where the bar | notation delimits the symmetrized indices. Besides the the generators, the other element needed to determine the gauge identities are the primary Euler-Lagrange equations. These are given by the straightforward evaluation of (4.3) for d = 2. Let us refer to them as E (2) (h µν ) and E (2) (G ρ µν ) , respectively. Then, their merging together with (4.52) as indicated in (2.29) yields the gauge identities for d = 2 Palatini we were seeking, after some tedious yet elementary algebra: Observe that the manifest symmetry under the exchange α ↔ β makes the number of independent gauge identities coincide with our earlier counting: g = 3.
Physical degrees of freedom. Plugging (4.54) in the master formula (2.2), we confirm the well-known fact that there are no physical degrees of freedom propagated by the theory: n dof = 0. We restate that the above cannot be obtained by simply setting d = 2 in (4.35).
To wrap up this section, we check our results are in good agreement with some of the previously carried out Hamiltonian calculations. We begin our comparisons by looking into the approach closest to our own, the one in [11]. There, the quantities (h µν , G ρ µν ) were regarded as the N = 9 a priori independent field variables for d = 2 Palatini, exactly as we did here. Following the Dirac-Bergmann procedure, it was shown that N (P) 1 = 3, N 1 = 6 and N 2 = 6, which readily confirms our own independent findings. In [12], the metricity condition (4.50) was taken into account from the onset. As a result of incorporating this information in the form of additional terms preceded by two Lagrange multipliers in the Hamiltonian, their setup had N N N = 11 = N = 9 a priori independent field variables. It was there shown that, in such formulation, N N N but we find that n dof = 0 for both sets of numbers upon employing (2.3). As a last remark, we notice that our calculations in this section 4.2 are comparatively simpler than those in [11,12]. Namely, our approach is certainly to be preferred if the goal is to determine the constraint structure of the theory and thereby manifestly count its propagating modes.

31
The study of constrained systems was initiated in the thirties by Rosenfeld, in a sometimes overlooked work [48], nowadays acknowledged and revisited [49]. It was later greatly developed during the fifties [50] and has since been a very active field of theoretical research. As such, one may have the impression that the investigation of manifestly first order singular classical field theories must be an already closed subject. This is not true. There are ongoing advances in this fundamental topic, particularly within the Lagrangian picture. Besides the references already provided in section 2, the recent work [4] stands as a neat example. The methodology there put forward is equivalent to our own proposal, as we shall show in the next section 5.1.
To further reassure the reader of the topicality of our formalism, in section 5.2 we explain how our method lends itself to a conversion from an analytic machinery to a constructive one. Indeed, the Lagrangian building principle originally put forward in [8,9] finds in the contents of section 2 a solid footing for attempting the construction of novel theories. This argumentation is carried out in terms of a concrete application for clarity, but the general proposal is much broader. In particular, we explain that the less elaborated upon procedure in [8,9] was cornerstone for the development of the so-called Maxwell-Proca theory. This discussion justifies an interest in the calculations of section 3 well beyond a simple exemplification of the explicit usage of the proposed method. When gravity is to be involved, the examples in section 4 provide a useful possible basis.

On a recent equivalent Lagrangian approach
During the preparation of this manuscript, a novel Lagrangian approach to obtain the functionally independent Lagrangian constraints and count propagating modes in constrained systems (of the kind here considered) appeared [4]. The method therein is physically equivalent to that put forward in [2,3], which -as already mentioned-are complementary references to our own discussion in section 2. This equivalence can be easily verified, as both [4] and [2,3] provide a mapping between their proposed Lagrangian parameters and the usual numbers of different kinds of Hamiltonian constraints. We have checked this leads to a consistent mapping between their different Lagrangian parameters.
In our understanding, the method in [4] distinguishes itself because it introduces the notion of first and second class (functionally independent) Lagrangian constraints. In our language, these are easy to identify. They are the sum of the various functionally independent Lagrangian constraints arising at all prior stages whose algorithm finalizes non-dynamically (as in cases iia and iib) and dynamically (as in case i), respectively. This abstract definition is clarified in the following, by classifying the functionally independent Lagrangian constraints we found in all the given examples into first and second class Lagrangian constraints.
In the case of Maxwell electromagnetism, the primary Lagrangian constraint (3.4) we found is a first class Lagrangian constraint. This is because it leads to a secondary constraint (3.5) that is identically satisfied and so non-dynamically stabilized by means of the closure iia. In fact, this same example is worked out in [4] as well.
Next, consider the (hard) Proca theory. There, both the primary (3.9) and secondary (3.10) Lagrangian constraints we determined are second class Lagrangian constraints, since the algorithm closes dynamically at the next stage by means of case i. Such closure implies that the consistency under time evolution of the secondary constraint is determined through a tertiary equation of motion.
We move to Schwinger-Plebanski formulation of both electromagnetism and the (hard) Proca theory. In Turning to d = 2 Palatini, we see it is rather simple to reclassify the nine functionally independent Lagrangian constraints we obtained into first and second class. At the primary level, we notice that there are six velocity dependent Lagrangian constraints among the relations that follow from requiring the vanishing of (4.41). These are [ϕ a ] ! ≈ 1 0, where a = 1, 2, . . . , 6 and the tensorial indices outside the square brackets [·] have been omitted. Their stability is dynamically ensured (via the secondary equations of motion) and so these are second class Lagrangian constraints. The three remaining primary Lagrangian constraints, 7, 8, 9, are manifestly velocity independent. They show a trivial stability at the secondary stage, see (4.46). Accordingly, we identify these as three first class Lagrangian constraints.
At last, we reclassify the functionally independent Lagrangian constraints we found for d > 2 Palatini into first and second class Lagrangian constraints. Recall that we obtained M 1 = d(d + 1) 2 /2 functionally independent primary Lagrangian constraints, given by the vanishing of (4.11). Notice now that we can straightforwardly split these primary constraints into where the subscripts (n)v stand for (non-) velocity dependent constraints. The consistency under time evolution of the velocity dependent constraints is dynamically fixed at the secondary stage and so these are second class Lagrangian constraints. The remaining velocity independent constraints give rise to the M 2 = M 1(nv) number of functionally independent secondary Lagrangian constraints, which are equal to the vanishing of (4.17). Once more, it is trivial to differentiate between

Relation to the Maxwell-Proca theory and beyond
As we explicitly showed in section 3.1, in a purely Lagrangian formulation with as many a priori independent field variables as the dimension of the underlying flat spacetime, the constraint structure of the simplest theory for a single Maxwell field can be characterized by the triplet t in (3.11). Employing the results of [2], we also verified the corresponding Hamiltonian characterization of these two triplets. We thus checked that the Maxwell and (hard) Proca fields are associated with two first and second class constraints, respectively. Although usually Maxwell and Proca fields are defined in the latter Hamiltonian manner, in the following we take the former Lagrangian triplets as the vector fields' defining features. We stress both definitions are equivalent.
The manifestly first order completions of the Maxwell and (hard) Proca theories analyzed in sections 3.1 and 3.2 are non-linear electrodynamics (NLE) and the so-called generalized Proca (GP) or vector-Galileon theory 12 , respectively. NLE encompasses a large class of theories. The celebrated Born-Infeld theory [52] is part of it, but also the more recently proposed exponential [53] and logarithmic [54] electrodynamics, among others. Schematically, the Lagrangian density for NLE can be written as where L M is the Maxwell Lagrangian as introduced in (3.3) and f is a smooth real function. Notice that the above depends on the Maxwell field A µ exclusively through its field strength A µν -up to boundary terms. Indeed, it is well-known [35] that a more involved dependence is not possible, if the U (1) gauge symmetry is to be respected. This feature remains true even when coupling the Maxwell field to General Relativity [55]. Only a few fine-tuned terms that contract A µν with the Riemann tensor are possible in such case. It is not hard to convince oneself that the constraint structure of NLE is characterized by the triplet t (N ) M in (3.6). In other words, it has the same constraint structure as classical electromagnetism, in its standard formulation of section 3.1.
The GP theory was put forward in [30] and its complete Lagrangian was established in [56]. Again schematically, we may express it as (5.4) in d dimensions, where L P is the (hard) Proca Lagrangian in (3.7), g is a real smooth function and each T ν 1 ...νnρ 1 ...ρn is a certain smooth real object constructed out of the spacetime metric η µν , the d-dimensional Levi-Civita tensor µ 1 ...µ d and the Proca field B µ . Although GP has only been formulated for d = 4, its systematic construction allows for a straightforward inferring of (5.4). Here, the underlying key idea consists in supplementing the (hard) Proca Lagrangian with derivative self-interaction terms of the Proca field B µ .
This implies a non-local extension of the notion of mass for the vector field. As such, we regard GP as an effective classical field theory 13 . It can be readily inferred from the calculations in [8] that the constraint structure of GP is characterized by the triplet t (N ) P in (3.11). Namely, GP has the same constraint structure as the (hard) Proca theory, when the latter is formulated as in section 3.2.
Next, we consider a multi-field scenario, including n M number of Maxwell fields, as well as n P number of (generalized) Proca fields. In four-dimensional Minkowski spacetime, the Maxwell-Proca (MP) theory [8,9] is the complete set of manifestly first order (self-)interactions among an arbitrary number of real Abelian vector fields that propagates the correct number of degrees of freedom. These consistent interactions were derived by demanding that the constraint structure of each Maxwell and Proca field is characterized by the triplets t MP . Then, we say that the building principle of the theory is based on the requirement: ⊕ n P · t (N =d) P = (l = n M + 2n P , g = n M , e = 2n M ), (5.5) where in the last equality we have made use of (3.6) and (3.11).
At this point, it should be clear that our calculations of t in section 3, elementary as they are, can be used as a basis for the construction of non-trivial theories. Having a ready-to-be-used method optimized to obtain such triplets (i.e. the method explained in section 2 and graphically summarized in figure 1) is thus a powerful tool for the development of manifestly first order classical field theories where multiple fields of different spins (self-)interact.
For instance, an interesting open question is that of the consistent coupling of the MP theory to gravity. It is in principle possible to combine our calculations in all the previous sections to attempt this ambitious goal as follows. Let N 2 = (n M + n P )d + d(d + 1) 2 /2, with d ≥ 2 the dimension of the spacetime. A manifestly first order Lagrangian density L MP(2)Pa that describes the dynamics of n M number of Maxwell fields and n P number of (generalized) Proca fields in the presence of Einstein's gravity in terms of N 2 a priori independent field variables must be associated with a constraint structure characterizing triplet t where all the triplets on the right-hand side have already been calculated in this work; see (3.6), (3.11), (4.36) and (4.54). Substituting these results, we have that The conversion of any of the above necessary conditions into a Lagrangian density building principle is an algebraically involved exercise beyond the scope of our present investigations. We thus leave it for future works. A last remark is due. As we observed at the very end of section 2.1 and should be apparent from our calculations in section 4.1, it is in general a conceptually clear but algebraically non-trivial exercise to obtain the triplet t (N ) of a given Lagrangian density L within our framework. It is even more challenging to determine the (exhaustive) form of L from the necessary condition that it should be associated to a certain triplet t (N ) . The reason is that such inversion in the logic requires solving sets of coupled non-linear partial differential equations in most cases. Therefore, it is overwhelmingly convenient to use all freedom of choice available in order to simplify this task to the utmost. For instance, one is advised to choose constant null vectors for the Hessians at all stages, if possible. For the concrete research project here proposed, it may be the case that (5.6) is not the optimal starting point. It could happen that the equivalent demand with the right-hand side triplets as given in (3.32), (4.36) and (4.54), is a more befitting way to try to derive the set of consistent (self-)interactions of vector fields in a curved background. For the reasons given at the beginning of section 4, we believe that t (N ) (2)Pa is indeed a beneficial basis for the gravity piece above.

Conclusions
In the following, we summarize the results we have put forward in this manuscript. Then, we proceed to discuss their relevance and pertinence. At last, we comment on the increasing (in n) computational difficulty of evaluating Lagrangian constraints on constraint surfaces T C n and concretize the pathologies a theory may suffer from when the algorithm of section 2.1 is not verified to close.

Summary of results.
In section 2, we have collected and complemented results from the extensive literature on constrained systems and presented a self-contained and ready to be used method to determine all the constraints in a theory. By postulation, the theory is required to be described by a manifestly first order Lagrangian. We make the mild assumptions of the principle of stationary action and finite reducibility. When the theory is covariant, the iterative algorithm presented for the determination of the functionally independent Lagrangian constraints does not contravene this feature. Nonetheless, manifest covariance is generically lost in our approach. In sections 3 and 4, we have minutely exemplified the usage of our said procedure. In section 5, we have argued for the pertinence and contemporaneity of both the general formalism and the given examples. Indeed, an equivalent but different methodology has been put forward lately [4]. The examples of section 3 constitute the foundation of the also recent Maxwell-Proca theory [8,9] and those of section 4 can potentially form the basis for the consistent coupling of Maxwell-Proca to gravity.
Critical discussion of results. The procedure explained in section 2 presents two main appealing features. First, it is a coordinatedependent approach, as opposed to a geometrical one. It thus readily allows for its application, given a Lagrangian density satisfying the initial postulates, without having to work out any symplectic two-form.
With pragmatism in mind, section 2 has been written in a way that is (hopefully) accessible to a broad audience. Even though the method stands on a rigorous footing, the discussion has been made largely devoid of mathematical technicalities. Second, it is an intrinsically Lagrangian procedure, as opposed to a Hamiltonian or a hybrid one. The appeal of this characteristic resides in the fact that, in many areas of high energy theoretical physics, manifestly first order classical field theories are predominantly posed and studied in their Lagrangian formulation. This is the case for instance in cosmology, astrophysics, black hole physics and holographic condensed matter. In all these disciplines, GP, MP and allied theories, specially in the presence of gravity, have been convincingly argued to be of significant interest, e.g. [9,30,61]. As such, our proposed procedure avoids non-negligible obstacles that typically arise in the transformation from the Lagrangian to the Hamiltonian picture. Besides, as already noted in the end of sections 4.1 and 4.2, our Lagrangian approach is a computationally faster and simpler way to obtain the constraint structures of these theories, compared to representative Hamiltonian analyses. (The examples in section 3 are so effortless comparatively that they do not substantiate an analogous argumentation.) In more detail, implementing our algorithm in section 2.1 is considerably easier than carrying out a Hamiltonian counterpart algorithm based on the Dirac-Bergman [50] procedure. As the attentive reader will have already noticed in our explicit examples of section 4 and we shall address shortly, the most demanding step in our approach consists in evaluating the n-th stage Lagrangian constraints in the (n − 1)-th constraint surface, with n ≥ 1. An analogous evaluation is necessary within the Hamiltonian picture as well, where two additional hurdles arise. On the one hand, one must classify the Dirac constraints into first and second class. This entails calculating the Poisson brackets of all Dirac constraints, a generically challenging task in field theory because non-local algebras usually arise 14 , e.g. [13,14]. On the other hand, in the standard Hamiltonian transition from one stage to the next, novel constraints emerge and must be consistently included via Lagrange multipliers. Closure of the algorithm requires the determination of as many Lagrange multipliers as possible, which in turn implies the resolution of algebraic or even differential equations. Even in the comparatively benign algebraic scenario, finding a solution is an increasingly (in stage) laborious and non-trivial task that involves inverting field-dependent matrices with complicated spatial index structures.
For a suggestive utility of the examples in sections 3 and 4, the reader is referred to section 5.2. Recall that the proposal therein is illustrative of the general theory-construction idea outlined in the introduction section 1 and at the beginning of section 5.

Two final observations.
In the first of our observations, we bring to light a series of considerations that must be taken into account when applying our method. In particular, we wish to discuss the practical complications that field theories of the kind here considered commonly exhibit when their Lagrangian constraints are to be evaluated on the suitable constraint surface.
First, we debunk what naively may look like an ambiguity. Recall that any constraint surface T C n for some finite n ≥ 1 is defined by the weak vanishing of the functionally independent Lagrangian constraints at all prior stages: As a direct consequence of the above, one can determine a maximal set of functionally independent relations of the formQ Though it should be clear by now, we confirm the different role played by the generalized velocitiesQ A and the spacelike derivatives of the generalized coordinates ∂ i Q A . The former are independent coordinates on T C, while the latter are functionally related to the generalized coordinates Q A . This clarification becomes pertinent when evaluating the secondary Lagrangian constraints in T C 1 already. At this point (and in subsequent stages), derivatives of the form ∂ iQ A generically show up. In such expressions, one must first replace the primary weak expression for the generalized velocityQ A -if pertinent-and then apply the spatial derivative on it.
Having clarified this point, we notice that its consistent implementation leads to the following nested situation. Substitution ofQ A according to (6.2) in ∂ iQ A normally leads to the presence of terms of the form ∂ i Q B in (6.2). These are again prone to be evaluated in T C n and can in turn contribute terms depending on Q B 's in (6.2); etc. We emphasize that one must reach an expression where this nesting ceases to occur, before proceeding with the algorithm. Not doing so would imply a wrong evaluation of the Lagrangian constraints in T C n , may lead to a misidentification of the functionally independent Lagrangian constraints and will almost invariably yield wrong results at the following (n + 1)-th stage. In fact, a wrong evaluation will typically land the researcher in a physically inequivalent theory from the one he/she started with.
Additionally and normally, when evaluating some Lagrangian constraints in T C n , potentially contrived functions of the previous stages' functionally independent Lagrangian constraints also show up. To understand the difficulty their appearance implies, consider the tertiary Lagrangian constraints (4.25) we found for d > 2 Palatini. Their raw expressions, prior to any evaluation in a constraint surface, contain quantities f = f (Q A ,Q A , ∂ i Q A ) that vanish in T C 1 . However, recognizing such f 's as primary weak zeros is a challenging task. Specifically, where in the ⊃ relation we have omitted numerical factors and the ϕ's are as given in (4.11). In the expression for f ij , the first equality is non-trivial, while the subsequent primary weak equality is obvious. An analogous situation arises with other f 's that are based on both functionally independent primary (4.11) and secondary (4.17) Lagrangian constraints. A brute force resolution to identify all such f 's consists in putting forward the most general ansatz compatible with the tensorial character of each of the Lagrangian constraints one is trying to evaluate and comparing it to their explicit expressions. This is indeed how we laboriously arrived at (4.27). For the second and last observation, the reader should heed (2.2) and (2.27). We already stressed the importance of closing the iterative algorithm for obtaining the functionally independent Lagrangian constraints towards the end of section 2.1. Now, we are equipped to better grasp the implications of not doing so, mentioned in the introductory section 1. Most often, failure to close the algorithm will give rise to the propagation of unphysical modes. These are Ostrogradski instabilities [5], but we shall loosely refer to them as ghosts. Even after ensuring ghost-freedom, not closing the algorithm can lead to trouble: it may overconstrain the theory, so that fewer than the desired number of degrees of freedom are propagated.
Let us consider the MP theory [8,9] discussed in section 5.2 as a concrete framework to clarify the above two unwanted scenarios. For our present purposes, it will suffice to consider the case when there are no Maxwell fields n M = 0 and there are an arbitrary but finite number of Proca fields n P . Recall that, in the standard formulation, we already saw in section 3.2 that a Proca field is associated to l = M 1 +M 2 = 1+1 = 2 number of functionally independent Lagrangian constraints. Bear in mind that this is also true for a generalized Proca field.
We denote the natural generalization of the GP theory in (5.4) to a multi-field setup as L PP . L PP automatically leads to M 1 = n P number of functionally independent primary constraints. The consistency under time evolution of these constraints does not generically yield the M 2 = n P number of functionally independent secondary constraints one would naively expect. Only a fine-tuned subset of terms in L PP does, precisely the terms that are part of the MP theory. For all those terms, it was shown that no tertiary constraints arise (M 3 = 0) and the algorithm closes dynamically giving rise to l = 2n P . Therefore, the correct number of physical modes n dof = d − n P are present in the theory. (To obtain this result, notice that, since there are no gauge identities, g = 0 = e.) Notice that, if one studies only the primary stage for L PP , one will be deceived into thinking that the theory is valid, as it suitably extends the primary stage of GP. However, L PP has M 2 < n P in general and therefore l < 2n P and n dof > d − n P . The additional propagating modes are precisely the ghosts of the first scenario we warn against.
If one studies both the primary and secondary stages for L PP , then one can fine-tune the Lagrangian density so that M 2 = n P as desired. But these functionally independent secondary constraints in the fine-tuned theory are at this point not necessarily stable. Their consistency under time evolution could in principle lead to further functionally independent tertiary constraints, so that l > 2n P and n dof < d − n P . This would place us in the second unwelcome scenario. For the given example, it so happens that the fine-tuned L PP is associated to a full rank tertiary Hessian. Consequently, the functionally independent secondary constraints are dynamically stabilized without further fine-tunings of the theory. However, this cannot be assumed, it has to be checked, so as to ensure the theory is not overconstrained.
It is interesting to point out that in [62] our very same admonition against the overconstrained scenario is made, albeit in a different context. The authors look into second order field theories with no gauge symmetry and derive the necessary conditions for such Lagrangians to not propagate ghosts. They show that, in the presence of Lorentz symmetry, the existence of any number M 1 > 0 of functionally independent Lagrangian constraints automatically leads to the same number M 2 = M 1 of functionally independent Lagrangian constraints. They unequivocally recognize our second scenario: those M 2 are not necessarily stable, so one could be facing an overconstrained theory.

A Formulae at an arbitrary stage of the algorithm
In this appendix, we show the explicit expressions of all quantities involved in an arbitrary a-th stage of the iterative algorithm for irreducible theories presented in section 2.1. Needless to say, in the appropriate limit, the general expressions here given yield the primary and secondary stages' formulae there shown.
Let ϕ Aa ! :≈ a 0 be a set of M a number of functionally independent Lagrangian constraints in the a-th stage, with A a = 1, 2, . . . , M a . These constraints are relations between the generalized coordinates Q A and velocitiesQ A of the field theory under consideration. They define the so-called a-th constraint surface where T C 0 is the moduli space of the theory, defined in (2.18). In order to ensure the preservation of the said constraints under time evolution, we demand We refer to E Aa as the (a + 1)-th stage Euler-Lagrange equations. Next, we will explicitly write E Aa . But to do so, we must first define the following objects. Let W AaA b denote the (a + 1)-th stage Hessian. This is a square matrix of dimension M a that allows us to define M a+1 := dim(W AaA b ) − rank(W AaA b ). We refer to the M a+1 number of linearly independent null vectors associated to W AaA b as γ A a+1 . Explicitly, (γ A a+1 ) Aa W AaA b = 0, with A a+1 = 1, 2, . . . , M a+1 . We require them to fulfil the normalization condition so that they form a basis in the kernel of W AaA b . Here, T stands for the transpose operation. With the help of the above null vectors, the (a + 1)-th stage Hessian can be expressed in terms of the the functionally independent a-th stage Lagrangian constraints as follows: Finally, we introduce the auxiliary matrix M AaA b (the Moore-Penrose pseudo-inverse of W AaA b ), which always exists and is uniquely determined from the relations Using the above, the (a + 1)-th stage Euler-Lagrange equations in (A.2) can be written as where the expression (A.4) is to be employed for W AaA b and where we have (recursively) defined with α A as given in (2.12). To obtain the presented α A b , the previous a-th stage Euler-Lagrange equations are employed. These in turn depend on the (a − 1)-th stage Euler-Lagrange equations and so on. This is the origin of the noted recursion. Notice we therefore explicitly employ the primary Euler-Lagrange equations and so the expression (A.6) is an on shell statement. In order to reproduce the results in section 2.1 from the above discussion, the reader only needs to do the index replacements (A 0 ≡ A, B, . . .), (A 1 → I, J, . . .), (A 2 → R, S, . . .), etc., as well as take footnote 7 into account.