{'fraction_non_alphanumeric': 0.048910547532785634, 'fraction_numerical': 0.013956217788437005, 'mean_word_length': 4.187268395976707, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 61, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Lagrangian Constraints and Differential Thomas Decomposition

In this paper we show how to compute algorithmically the full set of algebraically independent constraints for singular mechanical and field-theoretical models with polynomial Lagrangians. If a model under consideration is not singular as a whole but has domains of dynamical (field) variables where its Lagrangian becomes singular, then our approach allows to detect such domains and compute the relevant constraints. In doing so, we assume that the Lagrangian of a model is a differential polynomial and apply the differential Thomas decomposition algorithm to the Euler-Lagrange equations.


Introduction
Models with singular Lagrangians play a fundamental role in modern elementary particle physics and quantum field theory (cf. [1,2]). Singularity of such models is caused by local symmetries of their Lagrangians. Gauge symmetry is the most important type of local symmetries and it is imperative for all physical theories of fundamental interactions. The local symmetry transformations of differential equations describing dynamical systems or field-theoretical models relate their solutions satisfying the same initial (Cauchy) data. For dynamical systems with only one independent variable the initial data include (generalized) coordinates and velocities whereas for field-theoretical models they include the field variables, their spatial and the firstorder temporal derivatives ('velocities'). The presence of local symmetries in a singular model implies that its general solution satisfying the initial data depends on arbitrary functions.
A distinctive feature of singular Lagrangian models is that their dynamics is governed by the Euler-Lagrange equations which have differential consequences in the form of internal (hidden) constraints for the initial data. This is in contrast to regular constrained dynamics whose constraints are external with respect to the Euler-Lagrange equations.
Given a model Lagrangian, it is very important to verify whether it is singular, and if so to compute internal constraints that are hidden in the Euler-Lagrange equations. Knowledge of such constraints is necessary in constrained dynamics for the local symmetry analysis [1] and for quantization [2] of the singular model. As it was observed first in [3,4] in the framework of the Dirac formalism for investigation and quantization of constrained dynamical systems and gauge field theories, the Dirac procedure of the constraint calculation is nothing else than completion of the Hamiltonian equation system to involution. In this respect, the Hamiltonian constraints are integrability conditions [5]. The algorithmic aspects of computation and separation of Hamiltonian constraints into the first and second classes were studied in [6]. Computation of independent Lagrangian constraints as integrability conditions for the Euler-Lagrange equations is essential for well-posedness of initial value problems [5,7,8]. Such constraints are detected in the course of completion of the Euler-Lagrange equations to involution.
In the present paper we consider Lagrangian models whose Lagrangians (mechanics) and Lagrangian densities (field theory) are differential polynomials. Under this condition we show that the Thomas decomposition [9,10], being a characteristic one (cf. [11]) for the radical differential ideal generated by the polynomials in Euler-Lagrange equations, provides an algorithmic tool for the computation of Lagrangian constraints. Unlike the traditional linear algebra based methodology used in theoretical and mathematical physics [1] for computing linearly independent Lagrangian constraints, our approach is not only fully algorithmic but also takes into account rank dependence of the Hessian matrix on the dynamical (field) variables and outputs the complete set of algebraically independent constraints. The Thomas decomposition splits the system of Euler-Lagrange equations into a set of involutive subsystems and automatically determines the complete set of constraints for each output subsystem which has a certain rank deficiency of the Hessian. Moreover, each output subsystem admits a formally well-posed Cauchy problem [5,7,12] in terms of the output constraints. In contrast to other types of characteristic decompositions of (radical) differential ideals, a Thomas decomposition consists of differential systems whose solution sets are pairwise disjoint. This is accomplished by the use of inequations (cf. Example 4.6), allowing an irredundant description of the cases exhibiting different sets of constraints.
This paper is organized as follows. Section 2 contains a short description of our research objects. In Section 3 we outline the standard approach [1] to the computation of Lagrangian constraints. Our approach to this problem is based on the general algorithmic technique of differential Thomas decomposition [9,10] outlined in Section 4. In Section 5 we adopt this technique to the application to constrained dynamics and give a rigorous algebraic description of this adaptation. In Section 6 we illustrate its application by a number of examples including both field-theoretic and dynamical system models. Some concluding remarks are given in Section 7.

Preliminaries
In the framework of their field-theoretical description most of the fundamental laws of physics can be understood via the action [1, 2] where L is a Lagrangian density depending on the field variables and their first-order partial derivatives. In terms of its density L the Lagrangian L is expressed as The principle of least action (Hamilton's variational principle) implies the (field) Euler-Lagrange equations Here and in what follows we use Einstein's summation convention, i.e., summation over all possible values of the repeated indices occurring in a single term is assumed.
For dynamical systems the action and the Euler-Lagrange equations read [1, 2, 8] where q a are (generalized) coordinates and q a t ≡ dq a /dt are their velocities . A Lagrangian is regular if its Hessian matrix or Hessian defined as is invertible, and singular otherwise. Let K be a differential field of characteristic zero. By definition, K is a field containing the rational numbers Q with commuting derivations δ 0 , δ 1 , . . . , δ n , i.e., additive maps δ i : K → K satisfying the Leibniz rule Moreover, let R n = K{ϕ 1 , . . . , ϕ m } be the partial differential polynomial ring with differential indeterminates ϕ 1 , . . . , ϕ m , endowed with the set ∆ = {∂ t , ∂ x 1 , ∂ x 2 , . . . , ∂ x n } of commuting derivations. In other words, R n is the polynomial ring with infinitely many indeterminates, where the action of ∂ t , ∂ x 1 , ∂ x 2 , . . . , ∂ x n on an indeterminate increases the corresponding exponent i, j 1 , . . . , j n by one and on K coincides with the action of δ 0 , δ 1 , . . . , δ n . We identify ϕ i with ∂ 0 t ∂ 0 x 1 . . . ∂ 0 x n ϕ i and also indicate differentiation by subscripts; e.g., ϕ 1 t,t,x 1 = ∂ 2 t ∂ x 1 ϕ 1 . We shall denote by R 0 the ordinary differential polynomial ring with differential indeterminates ϕ 1 , . . . , ϕ m and ∆ = {∂ t }.
In our study of the Euler-Lagrange equation systems (1) and (2) we shall restrict our consideration to field-theoretical and mechanical models such that L ∈ R n and L ∈ R 0 , respectively.
It should be noted that this restriction is consistent with most of fundamental physical field theories [1, 2] whose (singular) Lagrangian densities are differential polynomials in the field variables and their first-order partial derivatives over a field of constants.

Linear algebra based approach
In terms of the Hessian (3) the set of Euler-Lagrange equations in (1) and (2) can be written as where the differential polynomials H i, j and P i contain derivations in t of order at most 1.
If the Hessian H is singular, then there are differential consequences of equations in (5) that do not contain the second order derivatives in t. These consequences, called constraints, are of principal importance in the study of singular mechanical and field-theoretical models (cf.
[1]- [8]). Given a Lagrangian, the standard approach used in physics and mechanics to compute Lagrangian constraints (cf. [1]) can be formulated as the following procedure:

Computation of constraints based on linear algebra
Step 1 Compute the Hessian H in accordance to (3). Then derive the set E of Euler-Lagrange equations (5) of cardinality m := |E| and let C := { } .
Step 2 Compute the rank r of the Hessian taking into account equations in E.
Step 3 If r = m, then go to Step 6. Otherwise, go to the next step.
Step 4 Compute a basis V of the nullspace of H, set up and enlarge the equation set Step 5 Set m := r and go to Step 2.
Step 6 Return C.
The elements in C that have differential order at most 1 are called Lagrangian (cf., for example, [1, 13]), and the set of Lagrangian constraints is formed by a largest functionally independent subset.
The above procedure suffers from the following algorithmic drawbacks.
(a) The entries of the Hessian are ordinary or partial differential polynomials, i.e., we have H i, j ∈ R 0 or H i, j ∈ R n (cf. Section 2). If Gaussian elimination is applied in Steps 2 and 4, the computation is performed in the field of fractions of R 0 or R n in general. Division by a non-zero differential polynomial which vanishes on the solution set of the given system has to be prevented. Such differential polynomials are not detected by the above procedure (cf. also (c) below). 4 On the other hand, the rings R 0 and R n are not principal ideal domains. For this reason the Hessian H does not admit a Smith normal form in general, which would provide an algorithmic way to compute rank(H) in Step 2. In the case of a multivariate polynomial matrix its Smith normal form can be defined and computed only in very special cases (cf. [14] and the references therein) that are not related to the Hessian of the general form (3). Similar remarks apply to the computation of a basis of the nullspace in Step 4 (cf., for example, [15] and its bibliography). (b) Generally, and we illustrate this fact by examples in Section 6, the rank of the Hessian may vary from one solution subspace of the Euler-Lagrange equations (1) or (2) to another one. Accordingly, the set of Lagrangian constraints may depend on a solution domain. (c) Computation of rank(H) in Step 2 assumes simplification (reduction) of the Hessian modulo (solution of) the system of Euler-Lagrange equations and its extension, respectively. Algorithmically, such a simplification can be performed by means of a characteristic decomposition [11] of the radical differential ideal generated by the left hand sides of the Euler-Lagrange equations. Thus, the procedure has to be extended with such a kind of decomposition. (d) The output set C of constraints has to be further processed to extract the set of Lagrangian constraints.
Below we show that differential Thomas decomposition [9,10] applied to the Euler-Lagrange equations (5) provides a fully algorithmic way to compute an algebraically independent set of Lagrangian constraints that takes into account the dependence of the rank of the Hessian on the solution domain.

Thomas decomposition and simple systems
We sketch the idea of Thomas decomposition of systems of polynomial ordinary or partial differential equations (ODEs or PDEs) and its computation. This method traces back to work by J. M. Thomas in the 1930s [16]. For more details, we also refer to [9,10]. Let be a system of polynomial ordinary or partial differential equations and inequations, i.e., the left hand sides p 1 , . . . , p s , q 1 , . . . , q h are elements of the differential polynomial ring R 0 or R n over the differential field K of characteristic zero. A fundamental problem is to determine all power series solutions of (6). Writing the unknown functions ϕ 1 , . . . , ϕ m (of either t or of t, x 1 , . . . , x n ) as power series with indeterminate coefficients, we are led to substitute this ansatz into (6) and compare coefficients. However, this naive procedure does not take into account conditions on the Taylor coefficients of ϕ 1 , . . . , ϕ m which arise as non-trivial consequences by differentiating the equations in (6), taking linear combinations of their left hand sides, etc. Therefore, system (6) needs to be transformed into an equivalent system which incorporates all integrability conditions and which is in this sense formally integrable (cf. [5]). For nonlinear systems of PDEs certain case distinctions are necessary in general, resulting in a finite family of formally integrable PDE systems (e.g., a differential Thomas decomposition) such that the union of their solution sets equals the solution set of the original PDE system. Differential Thomas decomposition has various applications, allowing, e.g., to decide membership to the radical differential ideal which is generated by the given PDE system or to solve certain differential elimination problems [10].
In order to identify the most significant variable in each non-constant (differential) polynomial, a total ordering ≻ on the set of (partial) derivatives of all orders of all unknown functions ϕ 1 , . . . , ϕ m is fixed beforehand. It is assumed to respect differentiation and that differentiation increases the significance of terms. The ordering ≻ is referred to as ranking.
We distinguish two stages of the computation of a differential Thomas decomposition, which may be intertwined in practice: its algebraic part and its differential part. The algebraic part of Thomas' algorithm eliminates variables (representing unknown functions and their derivatives) in equations and inequations by applying polynomial division and ensures square-freeness of all left hand sides. For these reductions the corresponding left hand sides are considered as univariate polynomials in their most significant variables, which are commonly referred to as their leaders, with coefficients that are polynomials in lower ranked variables.
More precisely, if two distinct equations p i = 0 and p j = 0 have the same leader v and , then there exist (differential) polynomials c 1 , c 2 , and r such that and r is zero or its leader is ranked lower than v or its leader is v and deg v (r) < deg v (p j ). This pseudo-division eliminates the highest power of v in p i , where c 1 may be chosen as a suitable power of the initial init(p j ) of p j , i.e., the coefficient of the highest power of v in p j . Replacing the equation p i = 0 with the equation r = 0 in the system requires that c 1 does not vanish on the solution set of the system, if the solution set of the system is to be maintained by this operation. For this reason, Thomas' algorithm possibly splits the system into two systems which are subsequently treated in the same way as the original system, which is discarded. These two systems are defined by inserting the inequation init(p j ) 0 and the equation init(p j ) = 0, respectively. For each two distinct inequations q i 0 and q j 0 with the same leader the least common multiple (lcm) r of q i and q j is computed and the two inequations are replaced with r 0. This process involves pseudo-divisions as above and also distinguishes cases of vanishing or non-vanishing initials.
For each pair p i = 0, q j 0 of an equation and an inequation with the same leader the greatest common divisor r of p i and q j is computed, where polynomials are considered again as univariate polynomials in their leaders and it is required that initials of pseudo-divisors do not vanish on the solution set of the system. If p i divides q j , then the system is inconsistent and is discarded. If r ∈ K \ {0}, then the inequation q j 0 is removed from the system. Otherwise, the equation p i = 0 is replaced with the equation p i /r = 0.
Finally, further reductions and splittings of systems may be necessary to ensure squarefreeness of all left hand sides (as univariate polynomials), i.e., non-vanishing of their discriminants. Since K is of characteristic zero, square-freeness of a univariate polynomial is characterized by the condition that the polynomial and its derivative have no root in common. Adhering to the recursive representation of polynomials, a computation of the greatest common divisor of a left hand side and its derivative with case distinctions as above allows to determine the square-free part.
The above procedure is a variant of Euclid's algorithm for multivariate polynomials, terminating after finitely many steps. If the fact that the variables occurring in p i and q j represent unknown functions and their derivatives is neglected, then the result is a finite collection of algebraic systems with the following property. Definition 4.1. An algebraic system S as in (6) is said to be simple (with respect to ≻), if the following three conditions are satisfied, where S ≺v is the subsystem of S consisting of those equations and inequations whose leader is ranked lower than the variable v.  . . , p s , q 1 , . . . , q h in a simple algebraic system are algebraically independent over K, i.e., the zero polynomial is the only polynomial with s+h indeterminates and coefficients in K which vanishes under substitution of p 1 , . . . , p s , q 1 , . . . , q h .
The differential part of Thomas' algorithm works towards formal integrability. For each two distinct equations p i = 0 and p j = 0 whose leaders ∂ k ϕ a : x n ϕ a are derivatives of the same unknown function ϕ a , the algorithm considers the derivatives of p i and p j which both have leader lcm(∂ k , ∂ l ) ϕ a . In an appropriate linear combination of these two derivatives the highest power of this leader is eliminated. Such a linear combination is referred to as a cross-derivative of p i and p j . The algorithm checks whether pseudodivision of this cross-derivative modulo (derivatives of) the equations of the system yields the zero polynomial. Non-zero remainders resulting from this reduction are added as new equations to the system, and the augmented system is treated by the algebraic part of Thomas' algorithm again. A method which traces back to M. Janet [17] allows an efficient organization of all crossderivatives. By distinguishing admissible and non-admissible derivations for each equation, each cross-derivative is considered only once.
The differential part performs pseudo-divisions as in (7), which however involve derivatives of equations of the system in general. If the equation p i = 0 is replaced with the equation r = 0, then we require that c 1 does not vanish on the solution set of the system. If a proper derivative of p j is subtracted, c 1 may be chosen as a suitable power of the initial of that derivative, which is the separant of p j , namely, the partial derivative of p j with respect to its leader. Since the discriminant of a differential polynomial (with respect to its leader) is essentially the resultant of the polynomial and its separant, the algebraic part of Thomas' algorithm ensures that separants do not vanish on the solution set of the system.
The result of Thomas' algorithm is a finite collection of differential systems with the following property. Definition 4.3. A differential system S as in (6) is said to be simple (with respect to ≻), if the following three conditions are satisfied.
(a) The system S is simple as an algebraic system (with respect to ≻, cf. Def. 4.1). (b) The cross-derivative of each pair of distinct equations whose leaders involve the same unknown function reduces to zero modulo the equations of the system and their derivatives.

7
(c) The left hand side q j of every inequation is reduced modulo the left hand sides of the equations, in the sense that no pseudo-division of q j modulo any p i is possible.

Remark 4.4.
Note that each splitting of a system S into two systems S 1 and S 2 performed by Thomas' algorithm is arranged so as to ensure that the solution set of S is the disjoint union of the solution sets of S 1 and S 2 .
Definition 4.5. Let S be a differential system as in (6). A finite collection S 1 , . . . , S k of simple differential systems whose solution sets form a partition of the solution set of S is called a differential Thomas decomposition of S (with respect to ≻).
A Thomas decomposition is not uniquely determined and its computation depends on the choice of the ranking ≻ and many other choices. Among the possible improvements of the above outline of Thomas' algorithm is the option to factorize the left hand sides of equations and to split a system according to a non-trivial factorization.
The algorithmic details of the Thomas decomposition method have been worked out in [9] (cf. also [10]). Implementations have been developed by T. Bächler and M. Lange-Hegermann as Maple packages AlgebraicThomas and DifferentialThomas and are freely available [18].
We illustrate the Thomas decomposition technique on a simple example. Example 4.6. Let us consider the PDE given by the following Hessian determinant: The Thomas decomposition method is applied over the differential polynomial ring K{u} with one differential indeterminate u and commuting derivations ∂ x , ∂ y , ∂ z , where K is a differential field of characteristic zero, e.g., K = Q so that the restrictions of ∂ x , ∂ y , ∂ z to K are all zero. We choose the degree-reverse lexicographical ranking satisfying ∂ x u ≻ ∂ y u ≻ ∂ z u, i.e., Then the initial of the PDE (8) is u z,z u y,y − u 2 y,z . We distinguish the cases whether this initial vanishes or not. In the first case we obtain the following simple differential system under the assumption that the initial u z,z of the initial of (8) does not vanish (no integrability condition needs to be checked for the single equation): Leaders of differential polynomials are underlined where not obvious. If u z,z vanishes in the first case, a simplification yields the second simple differential system y,z u x,x + 2 u y,z u x,z u x,y − u y,y u 2 x,z = 0 , u y,z 0 , u z,z = 0 8 after pseudo-reduction of a single cross-derivative to zero modulo the system. If the initial of (8) vanishes, we again distinguish the cases whether u z,z vanishes or not. If not, the pseudo-division u z,z 2 u y,z u x,z u y,z − u z,z u 2 x,y − u y,y u 2 x,z + u 2 x,z u z,z u y,y − u 2 y,z = u z,z u x,y − u y,z u x,z 2 allows to replace the PDE (8) with the expression in the bracket on the right hand side (ensuring square-freeness). We obtain the simple differential system after confirming that the integrability condition for the two equations is satisfied by applying pseudo-reduction modulo the system. Finally, if the initial of the PDE (8) vanishes, then u z,z = 0 implies u y,z = 0, and (8) reduces to u y,y u 2 x,z = 0. Distinguishing the cases whether u y,y vanishes or not, we obtain the simple differential systems Thus, we have the following tree of case distinctions (suppressing branches leading to inconsistent systems): u z,z u y,y − u 2 y,z 0 u z,z u y,y − u 2 y,z = 0 u z,z 0 u z,z = 0 u z,z 0 u z,z = 0 u y,y 0 u y,y = 0 As opposed to the Thomas decomposition (T 1 )-(T 5 ) with pairwise disjoint solution sets, e.g., the package DifferentialAlgebra in Maple 17 computes, using the Rosenfeld-Gröbner algorithm, the following characteristic decomposition for (8): The solution sets of these differential systems are not disjoint. The third and fourth systems have solutions f 1 (x) + f 2 (x) y + c z in common, where f 1 and f 2 are arbitrary analytic functions and c is an arbitrary constant.

Remark 4.7.
Due to its disjointness the Thomas decomposition partitions the solution space of the Euler-Lagrange equations (1) or (2). In doing so, it may occur that the field-theoretic or mechanical model, which is globally non-singular, is singular in some domains of the solution space. In those domains there are hidden Lagrangian constraints that are revealed by the Thomas decomposition (cf. Section 6 for examples). In terms of these local constraints one can pose an initial value (Cauchy) problem providing existence and uniqueness of local analytic solutions [5,7,12]. The local analyticity provides a smooth dependence of the solution on the initial data and thereby makes the initial-value problem well-posed (cf. [19]). It should be noted that the well-posedness of the initial value problem has been proven for an involutive simple differential system [7] related to the domain under consideration. For another type of characteristic decomposition which is not disjoint a proof of the well-posedness for Cauchy problem is not known in the literature.

Remark 4.8.
In some applications one might be interested only in one branch of the tree of case distinctions. For instance, if the differential ideal which is generated by the left hand sides of the given system of PDEs is prime, the resulting Thomas decomposition contains a simple system which is most generic in a precise sense (cf. [10,Subsection 2.2.3]). Clearly, if only the generic simple system is required, the computation can be restricted to that branch.

Adaptation of Thomas decomposition to constrained dynamics
Hereafter, we shall write the underlying ring of differential polynomials in (4) without subscript, since it will be clear from the context what ring is meant, partial or ordinary. We fix a ranking ≻ on R which satisfies The leader of any non-constant differential polynomial p in R, i.e., the differential indeterminate or derivative ∂ j t ∂ k 1 x 1 · · · ∂ k n x n ϕ a which is maximal with respect to ≻ among those occurring in p, is denoted by ld(p).
Remark 5.1. Condition (9) implies that for all i, j, k 1 , . . . , k n , l 1 , . . . , l n ∈ Z ≥0 and for all a, b ∈ {1, . . . , m} we have x 1 · · · ∂ l n x n ϕ b . The claim follows from applying ∂ k t to the left hand side and the right hand side of Definition 5.2. Let ψ ∈ {ϕ 1 , . . . , ϕ m } be minimal with respect to ≻. A non-constant differential polynomial p ∈ R is said to be a ∂ 2 t ψ-condition if holds. (2), their differential consequence p = 0 is called generalized Lagrangian constraint if p is a ∂ 2 t ψ-condition. If a constraint holds for all solutions of (1) or (2) it is called global, and local otherwise. A non-singular Lagrangian system with local constraints is locally singular.

Remark 5.4.
For singular models of mechanics, generalized Lagrangian constraints are just Lagrangian ones. For singular (or locally singular) field-theoretical models, Lagrangian constraints, as we underlined in Section 3, being of the first differential order, form a subset of the set of generalized Lagrangian (or local Lagrangian) constraints.
We consider the subring of R as a differential ring with derivations ∂ x 1 , . . . , ∂ x n and the ranking induced by ≻.
Remark 5.5. Every ∂ 2 t ψ-condition is an element of R ord(t)<2 . Remark 5.6. Let I be a differential ideal of R, i.e., an ideal of R which is closed under all derivations of R. The set of ∂ 2 t ψ-conditions in I is closed under taking R ord(t)<2 -linear combinations and applying an arbitrary number of the derivations ∂ x 1 , . . . , ∂ x n (cf. Remark 5.1). In other words, the set of ∂ 2 t ψ-conditions in I is a differential ideal of R ord(t)<2 . We denote by S = the set of left hand sides of equations in a differential system S . The saturation of a differential ideal E of R by a differential polynomial q ∈ R is defined as the differential ideal E : q ∞ := { p ∈ R | q k p ∈ E for some k ∈ Z ≥0 }. Proposition 5.7. Let S be a simple differential system over R with respect to a ranking satisfying (9) and let ψ ∈ {ϕ 1 , . . . , ϕ m } be minimal with respect to ≻. Let E be the differential ideal of R which is generated by S = and let q be the product of the initials and separants of all elements of S = . Moreover, let S ord(t)<2 := { p ∈ S = | ld(p) ≺ ∂ 2 t ψ } , E ord(t)<2 the differential ideal of R ord(t)<2 which is generated by S ord(t)<2 , and q ord(t)<2 the product of the initials and separants of all elements of S ord(t)<2 . Then the set of ∂ 2 t ψ-conditions in I := E : q ∞ is the radical differential ideal E ord(t)<2 : (q ord(t)<2 ) ∞ of R ord(t)<2 . 11 Proof. Let p be a ∂ 2 t ψ-condition in I. Since S is a simple differential system, repeated differential pseudo-reduction of p modulo S = yields zero (cf. [10,Prop. 2.2.50]). Since we have ld(p) ≺ ∂ 2 t ψ, only elements of S ord(t)<2 are chosen as pseudo-divisors. Their initials and separants are elements of R ord(t)<2 and the coefficients which are used for the pseudo-reductions are in R ord(t)<2 as well. This shows that the set of ∂ 2 t ψ-conditions in I is contained in E ord(t)<2 : (q ord(t)<2 ) ∞ . Conversely, every element of E ord(t)<2 : (q ord(t)<2 ) ∞ is clearly a ∂ 2 t ψ-condition in I.
We recall that the radical of a differential ideal I of R is defined as the differential ideal be a (not necessarily simple) system of polynomial partial differential equations and inequations, where p 1 , . . . , p s , q 1 , . . . , q h ∈ R, s, h ∈ Z ≥0 . We denote by E the differential ideal of R which is generated by p 1 , . . . , p s and by T the set of be a differential Thomas decomposition of (10) with respect to a ranking satisfying (9) and let ψ ∈ {ϕ 1 , . . . , ϕ m } be minimal with respect to ≻. For i = 1, . . . , r, let where E (i) ord(t)<2 is the differential ideal of R ord(t)<2 which is generated by and q (i) ord(t)<2 is the product of the initials and separants of all elements of S (i) ord(t)<2 . Then we have where E (i) is the differential ideal of R which is generated by (S (i) ) = and q (i) is the product of the initials and separants of all elements of (S (i) ) = , where i ∈ {1, . . . , r} (cf. [10, Prop. 2.2.72]). From Proposition 5.7 we conclude that, for each i ∈ {1, . . . , r}, the set of ∂ 2 t ψ-conditions in E (i) : (q (i) ) ∞ is equal to T (i) .
In the examples of the next section it is a very simple task to trace the origin of the ∂ 2 t ψconditions arising in the differential Thomas decomposition.

Examples
In this section we apply the Thomas decomposition to four models. Two of these models are field-theoretic and the other two are dynamical systems.
Example 6.1. The first model is the (1+1)-dimensional version of the chiral Schwinger model [20,21]. In the framework of Hamiltonian formalism the last model was studied, for example, in [2,22]. We consider the Lagrangian density in the form ( [22], Equation 19 ) Here, e, a are parameters, t, x are the independent variables and ϕ 1 = A 0 , ϕ 2 = A 1 , ϕ 3 = φ are the dependent variables. Let ≻ be the ranking satisfying (9) and such that

The Euler-Lagrange equations read
The Hessian (3) of this field-theoretic model is H = diag(1, 0, 1) and, hence, the model is singular. The Euler-Lagrange equations (11) are linear. In this case Thomas' algorithm applied to the system performs its completion to involution and outputs the Janet involutive form of (11) (cf. [9,10]). Denote the first and second equation in (11) by P and Q, respectively. Then the first step of the completion procedure consists in detection of the (underlined) leaders, calculation of the cross-derivative ∂ x P − ∂ t Q = 0 , and then elimination (reduction) of terms in the left hand side modulo equations in (11), as explained in Section 4. As a result of the completion, Thomas' algorithm yields The first equation in (12) (marked in boldface) is a Lagrangian constraint.

Example 6.2.
To show a field-theoretic situation when the differential Thomas decomposition computes local (generalized) Lagrangian constraints (cf. Definition 5.3) for a (globally) nonsingular (but locally singular) model, consider the following Lagrangian density: The independent variables are t, x and the dependent variables are ϕ 1 = u, ϕ 2 = v, ϕ 3 = w. We choose the ranking ≻ satisfying (9) such that The density (13) generates the Euler-Lagrange equations Let A, B, C denote the first, second and third equation, respectively. Considering the (underlined) leaders of these equations, Thomas' algorithm first replaces C with C − A. This corresponds to a row reduction of the Hessian: The leaders of the Euler-Lagrange equations are then u t,t , v t,t , w t,t . Since these differential polynomials involve pairwise distinct dependent variables, formal integrability of the system requires no further check (i.e., technically speaking, all derivations are admissible in the sense of Janet division). Only the possible vanishing of the initial of the second equation has be taken into account (compare also with the entry at position (2, 2) of the Hessian H (1) ). The field model (14) is not singular. For this reason the generic case leads to a simple system not containing any constraints (cf. (T 1 ) below). We conclude, using Theorem 5.8, that among the consequences of the given PDE system (14) there exist no (generalized) Lagrangian constraints.
There are precisely two possibilities for the initial of the second equation to vanish: w x = 0 or u = 0. All local generalized Lagrangian constraints for these more special consequences of the given PDE system (14) arise from the imposed equations or their integrability conditions [5].
In the first case (i.e., w x = 0), the integrability condition for the non-zero equations in (14) . We obtain three simple systems distinguishing the cases whether the initial 2 u v t + 1 of the last equation vanishes (cf. (T 4 ) below), or has vanishing initial u (cf. (T 3 ) below), or has nonvanishing initial u (cf. (T 2 ) below).
In the second case (i.e., u = 0, w x 0), the integrability condition for We obtain two simple systems distinguishing the cases whether v t vanishes or not (cf. (T 5 ) and (T 6 ) below, respectively).
Summarizing we obtain the following tree of case distinctions (14) where (T 1 ), . . . , (T 6 ) are the simple differential systems of the following Thomas decomposition of (14) and where for every simple subsystem we mark its local generalized Lagrangian constraints in boldface: Since the differential Thomas decomposition is disjoint (cf. Section 4), the simple differential subsystems (T 1 ), . . . , (T 6 ) (cf. Definition 4.3) partition the solution space of the Euler-Lagrange equations (14), and reveal local Lagrangian constraints.
Let A and B denote the first and second equation, respectively. According to the (underlined) leaders of these equations, Thomas' algorithm first performs a pseudo-reduction of B modulo A; more precisely, B is replaced with the remainder This computation corresponds to a row reduction of the Hessian: We obtain a zero row (indeed, the determinant of H (1) vanishes), which is reflected by the fact that all terms involving differentiation order 2 in (16) cancel. Hence, the model (15) is singular. The solution set of the system does not change when B is replaced with the remainder in (16) if the coefficient q 2 in (16) does not vanish on the solution set. For this reason, Thomas' 16 algorithm splits the original system (15) into one system incorporating the condition q 2 0, where B is replaced with q 2 1 − q 2 2 = 0, and a complementary system containing the new condition q 2 = 0. By taking the factorization q 2 1 − q 2 2 = (q 1 − q 2 )(q 1 + q 2 ) into account, the first system is split again into two complementary systems (T 1 ) and (T 2 ). For both systems, Thomas' algorithm reduces A modulo q 1 − q 2 = 0 or q 1 + q 2 = 0, respectively. Both remainders are divided by the non-vanishing factor 2 q 2 . Note that for ODE systems no formal integrability check is necessary, and the initials of the resulting left hand sides in (T 1 ) and (T 2 ) do not vanish. By reducing A and B modulo q 2 = 0 in the remaining case, we obtain the third simple differential system (T 3 ) of the Thomas decomposition. The Lagrangian constraints we find are q 1 − q 2 = 0 and q 1 + q 2 = 0 in (T 1 ) and (T 2 ), respectively (in boldface below).
Note that the local Lagrangian constraints in the simple systems (T 1 ) and (T 2 ) can be combined in a single global constraint q 2 1 − q 2 2 = 0 . ). Its Lagrangian is given by The independent variable is t, the dependent variables are y 1 = q 1 , y 2 = q 2 , and k, λ, and µ are non-zero parameters. We choose the ranking ≻ such that q 2 ≺ q 1 ≺ (q 2 ) t ≺ (q 1 ) t ≺ (q 2 ) t,t ≺ (q 1 ) t,t ≺ . . . Remark 6.5. We note that in our examples Thomas' algorithm is applied to a system of nonlinear differential equations which are linear in their leaders, i.e., in their highest derivatives with respect to the chosen ranking ≻. Therefore, one of the first steps in a computation of a differential Thomas decomposition for the given PDE system corresponds to the row reduction of the corresponding Hessian, which is part of the approach based on linear algebra recalled in Section 3.

Conclusion
In this paper we have shown that differential Thomas decomposition applied to the Euler-Lagrange equations of field-theoretical (1) or mechanical (2) models with polynomial Lagrangians is a proper algorithmic tool, implemented in Maple [18], for computing algebraically independent Lagrangian constraints or generalized Lagrangian constraints, respectively. In doing so, the decomposition being disjoint generates a partition of the solution space of the Euler-Lagrange equations. In the case of a singular or locally singular model the complete set of the corresponding local and algebraically independent constraints is computed for every element in the partition.
One can also apply another characteristic decomposition of the radical differential ideal generated by the Euler-Lagrange equations, for example, that based on the Rosenfeld-Gröbner algorithm and implemented in the diffalg library 1 of Maple (cf. [11] and its bibliography). Maple has one more built-in splitting algorithm [24] for nonlinear systems of differential equations (command rifsimp, a part of the package DEtools), that can also be used for the detection of Lagrangian constraints. However, diffalg and rifsimp do not yield disjoint decompositions in general, and the Lagrangian constraints that are inherent to different output subsystems may interfere.
The ranking (9) is a Riquier ranking [25,26] (cf. also [5]). Thus, by the Riquier existence theorems [27,28], a simple differential subsystem with nonempty equation set provides the existence and uniqueness of formal power series solution satisfying certain initial (Cauchy) data (cf. [5,12]). If the simple system under consideration is singular, and hence contains (generalized) Lagrangian constraints, then their presence is to be taken into account by the initial data, and the second-order derivatives in time ('accelerations') are uniquely defined. Differential Thomas decomposition can also be applied to detect and compute hidden constraints in differential-algebraic (and also in partial differential-algebraic) equations (DAEs). The presence of hidden algebraic constraints, i.e., those algebraic constraints that are not explicitly given in the system, is the main obstacle in numerical solving of DAEs (cf., for example, [29]).

Acknowledgements
The contribution of the first author (V.P.G.) was partially supported by grant No.13-01-00668 from the Russian Foundation for Basic Research.