{'fraction_non_alphanumeric': 0.04371116348832226, 'fraction_numerical': 0.012238621611061394, 'mean_word_length': 3.9157372986369268, 'pattern_counts': {'":': 0, '<': 38, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 53, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Existence and Construction of Vessiot Connections

A rigorous formulation of Vessiot's vector field approach to the analysis of general systems of partial differential equations is provided. It is shown that this approach is equivalent to the formal theory of differential equations and that it can be carried through if, and only if, the given system is involutive. As a by-product, we provide a novel characterisation of transversal integral elements via the contact map.


Introduction
Constructing solutions for (systems of) partial differential equations is obviously difficult -in particular for non-linear systems.Élie Cartan [4] proposed to construct first infinitesimal solutions or integral elements. These are possible tangent spaces to (prolonged) solutions. Thus they always lead to a linearisation of the problem and their explicit construction requires essentially only straightforward linear algebra. In the Cartan-Kähler theory [3,12], differential equations are represented by exterior differential systems and integral elements consist of tangent vectors pointwise annihilated by differential forms.
Vessiot [28] proposed in the 1920s a dual approach which does not require the use of exterior differential systems. Instead of individual integral elements, it always considers distributions of them generated by vector fields and their Lie brackets replace the exterior derivatives of differential forms. This approach takes an intermediate position between the formal theory of differential equations [13,19,23] and the Cartan-Kähler theory of exterior differential systems. Thus it allows for the transfer of many techniques from the latter to the former one, although this point will not be studied here.
Vessiot's approach may be considered a generalisation of the Frobenius theorem. Indeed, if one applies his theory to a differential equation of finite type, then one obtains an involutive distribution such that its integral manifolds are in a one-to-one correspondence with the smooth solutions of the equation. For more general equations, Vessiot proposed to "cover" the equation with infinitely many involutive distributions such that any smooth solution corresponds to an integral manifold of at least one of them.
Vessiot's theory has not attracted much attention: presentations in a more modern language are contained in [5,24]; applications have mainly appeared in the context of the Darboux method for solving hyperbolic equations, see for example [27]. While a number of textbooks provide a very rigorous analysis of the Cartan-Kähler theory, the above mentioned references (including Vessiot's original work [28]) are somewhat lacking in this respect. In particular, the question of under what assumptions Vessiot's construction succeeds has been ignored.
The purpose of the present article is to close this gap and at the same time to relate the Vessiot theory with the key concepts of the formal theory like formal integrability and involution (we will not develop it as a dual form of the Cartan-Kähler theory, but from scratch within the formal theory). We will show that Vessiot's construction succeeds if, and only if, it is applied to an involutive system of differential equations. This result is of course not surprising, given the well-known fact that the formal theory and the Cartan-Kähler theory are equivalent. However, to our knowledge an explicit proof has never been given. As a by-product, we will provide a new characterisation of integral elements based on the contact map, making also the relations between the formal theory and the Cartan-Kähler theory more transparent. Furthermore, we simplify the construction of the integral distributions. Up to now, quadratic equations had to be considered (under some assumptions, their solution can be obtained via a sequence of linear systems. We will show how the natural geometry of the jet bundle hierarchy can be exploited for always obtaining a linear system of equations. This article contains the main results of the first author's doctoral thesis [6] (a short summary of the results has already appeared in [7]). It is organised as follows. The next three sections recall the needed elements from the formal theory of differential equations and provide our new characterisation of transversal integral elements. The following two sections introduce the key concepts of Vessiot's approach: the Vessiot distribution, integral distributions and flat Vessiot connections. Two further sections discuss existence theorems for the latter two based on a stepby-step approach already proposed by Vessiot. As the proofs of some results are fairly technical, the main text contains only an outline of the underlying ideas and full details are given in three appendices. The Einstein convention is sometimes used to indicate summation.

The contact structure
Before we outline the formal theory of partial differential equations, we briefly review its underlying geometry: the jet bundle and its contact structure. Many different ways exist to introduce these geometric constructions, see for example [9,18,20]. Furthermore, they are discussed in any book on the formal theory (see the references in the next section).
Let π : E → X be a smooth fibred manifold. We call coordinates x = (x i : 1 ≤ i ≤ n) of X independent variables and fibre coordinates u = (u α : 1 ≤ α ≤ m) in E dependent variables. Sections σ : X → E correspond locally to functions u = s(x). We will use throughout a "global" notation in order to avoid the introduction of many local neighbourhoods even though we mostly consider local sections.
Derivatives are written in the form u α µ = ∂ |µ| u α /∂x µ 1 1 · · · ∂x µn n where µ = (µ 1 , . . . , µ n ) is a multi-index. The set of derivatives u α µ up to order q is denoted by u (q) ; it defines a local coordinate system for the q-th order jet bundle J q π, which may be regarded as the space of truncated Taylor expansions of functions s.
The hierarchy of jet bundles J q π with q = 0, 1, 2, . . . possesses many natural fibrations which correspond to "forgetting" higher-order derivatives. For us particularly important are π q q−1 : J q π → J q−1 π and π q : J q π → X . To each section σ : X → E, locally defined by σ(x) = x, s(x) , we may associate its prolongation j q σ : X → J q π, a section of the fibration π q locally given by j q σ(x) = x, s(x), ∂ x s(x), ∂ xx s(x), . . . .
The geometry of the jet bundle J q π is to a large extent determined by its contact structure. It can be introduced in various ways. For our purposes, three different approaches are convenient. First, we adopt the contact codistribution C 0 q ⊆ T * (J q π); it consists of all one-forms such that their pull-back by a prolonged section vanishes. Locally, it is spanned by the contact forms Dually, we may consider the contact distribution C q ⊆ T (J q π) consisting of all vector fields annihilated by C 0 q . A straightforward calculation shows that it is generated by the contact fields Note that the latter fields, C µ α , span the vertical bundle V π q q−1 of the fibration π q q−1 . Thus the contact distribution can be split into C q = V π q q−1 ⊕ H. Here the complement H is an ndimensional transversal subbundle of T (J q π) and obviously not uniquely determined (though any local coordinate chart induces via the span of the vectors C (q) i one possible choice). Any such complement H may be considered the horizontal bundle of a connection on the fibred manifold π q : J q π → X (but not for the fibration π q q−1 ). Following Fackerell [5], we call any connection on π q the horizontal bundle of which consists of contact fields a Vessiot connection (in the literature the terminology Cartan connection is also common, see for example [15]).
For later use, we note the structure equations of the contact distribution. The only nonvanishing Lie brackets of the vector fields (1) are Note that this observation implies that the vertical bundle V π q q−1 is involutive.
commutes for any section σ. Because of its linearity over π q q−1 , we may also consider it a map Γ q : J q π → T * X ⊗ J q−1 π T (J q−1 π) with the local coordinate form Obviously, Γ q (ρ, ∂ x i ) = T ρ π q q−1 (C (q) i ) and hence (C q ) ρ = (T ρ π q q−1 ) −1 im Γ q (ρ) for any point ρ ∈ J q π. Note that all vectors in the image of Γ q (ρ) are transversal to the fibration π q−1 q−2 . One of the main applications of the contact structure is given by the following proposition (for a proof, see [6,Proposition 2.1.6] or [23,Proposition 2.2.7]). It characterises those sections of the jet bundle π q : J q π → X which are prolongations of sections of the underlying fibred manifold π : E → X . Proposition 1. A section γ : X → J q π is of the form γ = j q σ for a section σ : X → E if, and only if, im Γ q γ(x) = T γ(x) π q q−1 T γ(x) im γ for all points x ∈ X where γ is defined.
Thus for any section σ : X → E the equality im Γ q+1 j q+1 σ(x) = im T x (j q σ) holds and we may say that knowing the (q + 1)-jet j q+1 σ(x) of a section σ at some x ∈ X is equivalent to knowing its q-jet ρ = j q σ(x) at x as well as the tangent space T ρ (im j q σ) at this point. This observation will later be the key for the Vessiot theory.

The formal theory of differential equations
We are now going to outline the formal theory of partial differential equations to introduce the basic notation. Our presentation follows [23]; other general references are [13,14,19]. Definition 1. A differential equation of order q is a fibred submanifold R q ⊆ J q π locally described as the zero set of some smooth functions on J q π: Note that we do not distinguish between scalar equations and systems.
We denote by ι : R q ֒→ J q π the canonical inclusion map. Differentiating every equation in the local representation (4) leads to the prolonged equation R q+1 ⊆ J q+1 π defined by the equations Φ τ = 0 and D i Φ τ = 0 where the formal derivative D i is given by Iteration of this process gives the higher prolongations R q+r ⊆ J q+r π. A subsequent projection leads to R (1) q = π q+1 q (R q+1 ) ⊆ R q , which is a proper submanifold whenever integrability conditions appear. Definition 2. A differential equation R q is formally integrable if at any prolongation order r > 0 the equality R (1) q+r = R q+r holds. In local coordinates, the following definition coincides with the usual notion of a solution.

Definition 3.
A solution is a section σ : X → E such that its prolongation satisfies im j q σ ⊆ R q .
For formally integrable equations it is straightforward to construct order by order formal power series solutions. Otherwise it is hard to find solutions. A constitutive insight of Cartan was to introduce infinitesimal solutions or integral elements at a point ρ ∈ R q as subspaces U ρ ⊆ T ρ R q which are potentially part of the tangent space of a prolonged solution.
Definition 4. Let R q ⊆ J q π be a differential equation and ι : R q → J q π the canonical inclusion map. Let I[R q ] = ι * C 0 q diff be the differential ideal generated by the pull-back of the contact codistribution on R q (algebraically, I[R q ] is then spanned by a basis of ι * C 0 q and the exterior derivatives of the forms in this basis). A linear subspace U ρ ⊆ T ρ R q is an integral element at the point ρ ∈ R q , if all forms in (I[R q ]) ρ vanish on it.
The following result provides an alternative characterisation of transversal integral elements via the contact map. It requires that the projection π q+1 q : R q+1 → R q is surjective.
lies transversal to the fibration π q q−1 is an integral element at the point ρ ∈ R q if, and only if, a pointρ ∈ R q+1 exists on the prolonged equation R q+1 such that π q+1 q (ρ) = ρ and T ρ ι(U ρ ) ⊆ im Γ q+1 (ρ). Proof . Assume first that U ρ satisfies the given conditions. It follows immediately from the coordinate form of the contact map that then firstly T ρ ι(U ρ ) is transversal to π q q−1 and secondly that every one-form ω ∈ ι * C 0 q vanishes on U ρ , as im Γ q+1 (ρ) ⊂ (C q ) ρ . Thus there only remains to show that the same is true for the two-forms dω ∈ ι * (dC 0 q ). Choose a section γ : R q → R q+1 such that γ(ρ) =ρ and define a distribution D of rank n on R q by setting T ι(Dρ) = im Γ q+1 γ(ρ) for any pointρ ∈ R q . Obviously, by construction U ρ ⊆ D ρ . It follows now from the coordinate form (3) of the contact map that locally the distribution D is spanned by n vector fields X i such that ι * X i = C (q) i + γ α µ+1 i C µ α where the coefficients γ α ν are the highest-order components of the section γ. Thus the commutator of two such vector fields satisfies . The commutators on the right hand side vanish whenever µ i = 0 or µ j = 0, respectively. Otherwise we obtain −∂ u α µ−1 i and −∂ u α µ−1 j , respectively. But this fact implies that the two sums on the right hand side cancel each other and we find that ι * [X i , X j ] ∈ C q . Thus we find for any contact form ω ∈ C 0 q that Each summand in the last expression vanishes, as all appearing fields are contact fields. Hence any form ω ∈ ι * (dC 0 q ) vanishes on D and in particular on U ρ ⊆ D ρ . For the converse, note that any transversal integral element U ρ ⊆ T ρ R q is spanned by linear Hence the coefficients are of the form γ α µ,i = γ α µ+1 i and a section σ exists such that ρ = j q σ(x) and T ρ (im j q σ) is spanned by the vectors T ρ ι(v 1 ), . . . , T ρ ι(v n ). This observation implies that U ρ satisfies the given conditions. For many purposes the purely geometric notion of formal integrability is not sufficient, and one needs the stronger algebraic concept of involution. This concerns in particular the derivation of uniqueness results but also the numerical integration of overdetermined systems [21]. An intrinsic definition of involution is possible using the Spencer cohomology (see for example [22] and references therein for a discussion). We apply here a simpler approach requiring that one works in "good", more precisely: δ-regular, coordinates x. This assumption represents a mild restriction, as generic coordinates are δ-regular and it is possible to construct systematically "good" coordinates -see [11]. Furthermore, it will turn out that the use of δ-regular coordinates is essential for Vessiot's approach.
Thus, the symbol is the vertical part of the tangent space to R q . Locally, N q consists of all vertical vector fields m α=1 |µ|=q v α µ ∂ u α µ where the coefficients v α µ satisfy the following linear system of algebraic equations: The matrix of this system is called the symbol matrix M q . The prolonged symbols N q+r are the symbols of the prolonged equations R q+r with corresponding symbol matrices M q+r . The class of a multi-index µ = (µ 1 , . . . , µ n ), denoted cls µ, is the smallest k for which µ k is different from zero. The columns of the symbol matrix (6) are labelled by the v α µ . We order them as follows. Let α and β denote indices for the dependent coordinates, and let µ and ν denote multi-indices for marking derivatives. Derivatives of higher order are greater than derivatives of lower order: if |µ| < |ν|, then u α µ ≺ u β ν . If derivatives have the same order |µ| = |ν|, then we distinguish two cases: if the leftmost non-vanishing entry in µ − ν is positive, then u α µ ≺ u β ν ; and if µ = ν and α < β, then u α µ ≺ u β ν . This is a class-respecting order: if |µ| = |ν| and cls µ < cls ν, then u α µ ≺ u β ν . Any set of objects indexed with pairs (α, µ) can be ordered in an analogous way. This order of the multi-indices µ and ν is called the degree reverse lexicographic ranking, and we generalise it in such a way that it places more weight on the multi-indices µ and ν than on the numbers α and β of the dependent variables. This is called the term-over-position lift of the degree reverse lexicographic ranking. Now the columns within the symbol matrix are ordered descendingly according to the degree reverse lexicographic ranking for the multi-indices µ of the variables v α µ in equation (6) and labelled by the pairs (α, µ). (It follows that, if v α µ and v β ν are such that cls µ > cls ν, then the column corresponding to v α µ is left of the column corresponding to v β ν .) The rows are ordered in the same way with regard to the pairs (α, µ) of the variables u α µ which define the classes of the equations Φ τ (x, u (q) ) = 0. If two rows are labelled by the same pair (α, µ), it does not matter which one comes first.
We compute now a row echelon form of the symbol matrix. We denote the number of rows where the pivot is of class k by β (k) q , the indices of the symbol N q , and associate with each such row its multiplicative variables x 1 , . . . , x k . Prolonging each equation only with respect to its multiplicative variables yields independent equations of order q + 1, as each has a different leading term. Definition 6. If prolongation with respect to the non-multiplicative variables does not lead to additional independent equations of order q + 1, in other words if then the symbol N q is involutive. The differential equation R q is called involutive, if it is formally integrable and its symbol is involutive.
The criterion (7) is also known as Cartan's test, as it is analogous to a similar test in the Cartan-Kähler theory of exterior differential systems. We stress again that it is valid only in δ-regular coordinates (in fact, in other coordinate systems it will always fail).

The Cartan normal form
For notational simplicity, we will consider in our subsequent analysis almost exclusively firstorder equations R 1 ⊆ J 1 π. At least from a theoretical standpoint, this is not a restriction, as any higher-order differential equation R q can be transformed into an equivalent first-order one (see for example [23,Appendix A.3]). For these we now introduce a convenient local representation.
Definition 7. For a first-order differential equation R 1 the following local representation, a special kind of solved form, · · · · · · · · · · · · · · · · · · · · · · · · · · · is called its Cartan normal form. The equations of zeroth order, u α = φ α (x, u β ), are called algebraic. The functions φ α k are called the right sides of R 1 . (If, for some 1 ≤ k ≤ n, the number of equations is β (k) 1 = m, then the condition β (k) 1 < δ ≤ m is empty and no terms u δ k appear on the right sides of those equations.) Here, each equation is solved for a principal derivative of maximal class k in such a way that the corresponding right side of the equation may depend on an arbitrary subset of the independent variables, an arbitrary subset of the dependent variables u β with 1 ≤ β ≤ β 0 , those derivatives u γ j for all 1 ≤ γ ≤ m which are of a class j < k and those derivatives which are of the same class k but are not principal derivatives. Note that a principle derivative u α k may depend on another principle derivative u γ l as long as l < k. The equations are grouped according to their class in descending order.
Proof . For the proof, see [19] or [23,Section 9.4] and references therein. The strategy is to split the system into subsystems according to the classes of the equations in it (see below). The solution is constructed step by step; each step renders a normal system in fewer independent variables to which the Cauchy-Kovalevskaya theorem is applied. Finally, the condition that R 1 is involutive leads to further normal systems ensuring that the constructed functions are indeed solutions of the full system with respect to all independent variables.
Under some mild regularity assumptions the algebraic equations can always be solved locally. From now on, we will assume that any present algebraic equation has been explicitly solved, reducing thus the number of dependent variables. We simplify the Cartan normal form of a differential equation as given in Definition 7 into the reduced Cartan normal form. It arises by solving each equation for a derivative u α j , the principal derivative, and eliminating this derivative from all other equations. Again, the principal derivatives are chosen in such a manner that their classes are as great as possible. Now no principal derivative appears on a right side of an equation (whereas this was possible with the non-reduced Cartan normal form of Definition 7). All the remaining, non-principal, derivatives are called parametric. Ordering the obtained equations by their class, we again can decompose them into subsystems: Note that the values β (k) 1 are exactly those appearing in the Cartan test (7), as the symbol matrix of a differential equation in Cartan normal form is automatically triangular with the principal derivatives as pivots. Provided that δ-regular coordinates are chosen, it is possible to perform a closed form involution analysis for a differential equation R 1 in reduced Cartan normal form. We remark that an effective test of involution proceeds as follows (see for example [23,Remark 7.2.10]). Each equation in (10) is prolonged with respect to each of its non-multiplicative variables. The arising second-order equations are simplified modulo the original system and the prolongations with respect to the multiplicative variables. The symbol N 1 is involutive if, and only if, after the simplification none of the prolonged equations is of second-order any more. The equation R 1 is involutive if, and only if, all new equations even simplify to zero, as any remaining first-order equation would be an integrability condition.
In order to apply this test, we now prove two helpful lemmata. We introduce the set B := (α, i) ∈ AE × AE: u α i is a principal derivative , and for each (α, i) ∈ B we define Φ α i := u α i − φ α i . Using the contact fields (1), any prolongation of some Φ α i can be expressed in the following form. Lemma 1. Let the differential equation R 1 be represented in the reduced Cartan normal form given by equation (10). Then for any (α, i) ∈ B and 1 ≤ j ≤ n, we have Proof . By straightforward calculation; see [6,Lemma 2.5.5].
For j > i, the prolongation D j Φ α i is non-multiplicative, otherwise it is multiplicative. Now let j > i, so that equation (11) shows a non-multiplicative prolongation, and assume that we are using δ-regular coordinates. According to our test, the symbol N 1 is involutive if, and only if, it is possible to eliminate on the right hand side of (11) all second-order derivatives by adding multiplicative prolongations.
If the differential equation is not involutive, then the difference does not necessarily vanish but may yield an obstruction to involution for any (α, i) ∈ B and any i < j ≤ n. The next lemma gives all these obstructions to involution for a first-order system in reduced Cartan normal form.
Lemma 2. Assume that δ-regular coordinates are used. For an equation in Cartan normal form and indices i < j and α such that (α, i) ∈ B, we have the equality Proof . By a tedious but fairly straightforward computation, see [6, pages 48-55].
In line (12) we have collected all terms which are of lower than second order. Furthermore, none of the appearing second-order derivatives is of a form that it could be eliminated by adding some multiplicative prolongation. Hence, under our assumption of δ-regular coordinates, the symbol N 1 is involutive if, and only if, all the expressions in square brackets vanish. The differential equation R 1 is involutive if, and only if, in addition line (12) vanishes, as it represents an integrability condition. Thus Lemma 2 gives us all obstructions to involution for R 1 in explicit form. They will reappear in the proof of the existence theorem for integral distributions in Section 7.

The Vessiot distribution
By Proposition 1, the tangent spaces T ρ (im j q σ) of prolonged sections at points ρ ∈ J q π are always subspaces of the contact distribution C q | ρ . If the section σ is a solution of the differential equation R q , then by definition it furthermore satisfies im j q σ ⊆ R q , and therefore T (im j q σ) ⊆ T R q . Hence, the following construction suggests itself.
Note that the Vessiot distribution is not necessarily of constant rank along R q (just like the symbol N q ); for simplicity, we will almost always assume here that this is the case. This definition of the Vessiot distribution is not the one usually found in the literature. But the equivalence to the standard approach is an elementary exercise in computing with pull-backs.
For a differential equation given in explicitly solved form, the inclusion map ι : R q → J q π is available in closed form and can be used to calculate the pull-back of the contact forms. This has the advantage of keeping the calculations within a space of smaller dimension, namely the submanifold R q . Thereby regarding the differential equation as a manifold in its own right, we bar its coordinates to distinguish them from those of the jet bundle. Example 1. Consider the first-order system given by the representation R 1 : Then from the prolongations u xt − v x = 0 and u xt − w t = 0 follows the integrability condition w t = v x by elimination of the second-order derivative u xt . Thus, the first projection of R 1 is (in Cartan normal form) represented by It is not difficult to verify that the projected equation R 1 is represented by a system in solved form, it is natural to choose appropriate local coordinates for R (1) 1 , which we bar to distinguish them: x, t; u, v, w; v x , w x . The contact codistribution for J 1 π is generated by This space is annihilated by These seven one-forms ω 1 , . . . , ω 7 annihilate the Vessiot distribution V[R 1 ], which is spanned by the four vector fields In local coordinates on R 1 , these four vector fields becomē They satisfy ι * Xi = X i , as a simple calculation using the Jacobian matrix for T ι shows. The vector fieldsX i are annihilated by the pull-backs of the contact forms, 1 ). For a fully nonlinear differential equation R q , in particular an implicit one, this approach to compute its Vessiot distribution V[R q ] via a pull-back is in general not effectively feasible. However, applying directly our definition of V[R q ], it is easily possible even for such equations to determine effectively T ι(V[R q ]), in other words: to realize it as a subbundle of T (J q π)| Rq . The contact fields (1) form a basis for C q . It follows that for any vector fieldX ∈ V[R q ], coefficients If the differential equation R q is locally represented by Φ τ = 0, where 1 ≤ τ ≤ t, it follows from the tangency of the vector fields in V[R q ] that dΦ τ (ι * X ) = ι * X (Φ τ ) = 0 and thus the coefficient functions must satisfy the following system of linear equations: where 1 ≤ τ ≤ t. Note that this approach to determine the Vessiot distribution is closely related to prolonging the differential equation R q and requires essentially the same computations. Indeed, the formal derivative (5) can be written in the form and in the context of the order-by-order determination of formal power series solutions (see for example [23, Section 2.3]) these equations are considered an inhomogeneous system for the Taylor coefficients of order q + 1 depending on the lower order coefficients. Taking this point of view, we may call (21) the "projective" version of (22). In fact for n = 1, that is, for ordinary differential equations, this is even true in a rigorous sense.
Example 2. We consider the fully nonlinear first-order ordinary differential equation R 1 locally defined by (u ′ ) 2 + u 2 + x 2 = 1. The contact distribution C 1 is spanned by the two vector fields consists of all linear combinations of these two fields which are tangent to R 1 . Setting ω = xdx + udu + u ′ du ′ , we thus have to solve the linear equation ω(aX 1 + bX 2 ) = 0 in order to determine T ι(V[R q ]). Its solution requires a case distinction (which is typical for fully nonlinear differential equations). If u ′ = 0, then we find For u ′ = 0 and x = 0, the Vessiot distribution is spanned by the vertical contact field X 2 . Finally, for x = u ′ = 0 the rank of the Vessiot distribution jumps, as at these points the whole contact plane is contained in it.
The definition of the symbol N q and of the Vessiot distribution V[R q ], respectively, of a differential equation R q ⊆ J q π immediately imply the following generalisation of the above discussed splitting of the contact distribution C q = V π q q−1 ⊕ H (such a splitting of the Vessiot distribution is also discussed by Lychagin and Kruglikov [14,15] where the Vessiot distribution is called "Cartan distribution").
Proposition 4. For any differential equation R q , its symbol is contained in the Vessiot distribution: for some complement H (which is not unique).
Such a complement H is necessarily transversal to the fibration R q → X and thus leads naturally to connections: provided dim H = n, it may be considered the horizontal bundle of a connection of the fibred manifold R q → X . Definition 10. Any such connection is called a Vessiot connection for R q .
In general, the Vessiot distribution V[R q ] is not involutive (that is, closed under the Lie bracket; an exception are formally integrable equations of finite type [23, Remark 9.5.8]), but it may contain involutive subdistributions. If these are furthermore transversal (to the fibration R q → X ) and of dimension n, then they define a flat Vessiot connection. Lemma 3. If the section σ : X → E is a solution of the equation R q , then the tangent bundle is an n-dimensional transversal involutive subdistribution, then any integral manifold of U has locally the form im j q σ for a solution σ of R q .
Proof . Let σ be a local solution of the equation R q . Then it satisfies by Definition 3 im j q σ ⊆ R q and thus T (im j q σ) ⊆ T R q . Besides, by the definition of the contact distribution, for be an n-dimensional transversal involutive subdistribution. Then according to the Frobenius theorem, U has n-dimensional integral manifolds. By definition, T ι(V[R q ]) ⊆ C q | Rq ; this characterises prolonged sections. Hence, for any integral manifold of U there is a local section σ such that the integral manifold is of the form im j q σ. Furthermore, the integral manifold is a subset of R q . Thus it corresponds to a local solution of R q . This simple observation forms the basis of Vessiot's approach to the analysis of R q : he proposed to construct all flat Vessiot connections. Before we do this, we first show how integral elements appear in this program.
be a transversal subdistribution of the Vessiot distribution of constant rank k. Then the spaces U ρ are k-dimensional integral elements for all points ρ ∈ R q if, and only if, Proof . Let {ω 1 , . . . , ω r } be a basis of the codistribution ι * C 0 q . Then an algebraic basis of the ideal I[R q ] is {ω 1 , . . . , ω r , dω 1 , . . . , dω r }. Any vector field X ∈ U trivially satisfies ω i (X) = 0 by Proposition 3. For arbitrary fields X 1 , X 2 ∈ U, we have The first two summands on the right hand side vanish trivially and the remaining equation implies our claim.
We call a subdistribution U ⊆ V[R q ] satisfying the conditions of Proposition 5 an integral distribution for the differential equation R q . In the literature [24], the terminology "involution" is common for such distributions which, however, is confusing. Note that generally an integral distribution is not integrable; the name only reflects the fact that it consists of integral elements.
A general differential equations R q does not necessarily possess any Vessiot connection (not even a non-flat one). Their existence is linked to the absence of integrability conditions. More precisely, we obtain the following characterisation.
Proposition 6. Let R q be a differential equation. Then its Vessiot distribution possesses locally a direct decomposition with an n-dimensional complement H such that V[R q ] = N q ⊕ H if, and only if, there are no integrability conditions which arise as the prolongation of equations of lower order in the system.
Proof . Let R q be locally represented by such that the equations Φ τ (x, u (q) ) = 0 do not imply lower-order equations which are independent of the equations Ψ σ (x, u (q−1) ) = 0. Let u (q) denote the subset of all derivatives of order q only; then the Jacobi matrix ∂Φ τ (x, u (q) )/∂u (q) has maximal rank. If we proceed as in the last proof, then the ansatz (20) for the determination of the Vessiot distribution yields for the above representation the linear system Here, the matrix C µ α (Φ τ ) has maximal rank, too; thus the equations C And since no terms of order q are present in Recall that we consider the Vessiot distribution, and thus the linear system (24), only on R q . It follows that the subsystem C From the proof of this proposition follows that for the determination of the Vessiot distribution V[R q ], equations of order less than q in the local representation of R q can be ignored if there are no integrability conditions which arise from equations of lower order. It is the integrability conditions which arise as prolongations of lower order equations that hinder the construction of n-dimensional complements, while those which follow from the relations between cross derivatives do not influence this approach.
Remark 1. If one does not care about the distinction between different kinds of integrability conditions and simply requires that R q = R (1) q (meaning that no integrability conditions at all appear in the first prolongation of R q ), then one can provide a more geometric proof for the existence of an n-dimensional complement (of course, in contrast to Proposition 6, the converse is not true then).
The assumption R q = R (1) q implies that to every point ρ ∈ R q at least one pointρ ∈ R q+1 with π q+1 q (ρ) = ρ exists. We choose such a point and consider im Γ q+1 (ρ) ⊂ T ρ (J q π). By definition of the contact map Γ q+1 , this is an n-dimensional transversal subset of C q | ρ . Thus there only remains to show that it is also tangential to R q , as then we can define a complement by T ρ ι(H ρ ) = im Γ q+1 (ρ). But this tangency is a trivial consequence ofρ ∈ R q+1 ; using for example the local coordinates expression (3) for the contact map and a local representation Hence it is possible to construct for each point ρ ∈ R q a transversal complement H ρ such that There remains to show that these complements can be chosen so that they form a smooth distribution. Our assumption R q = R (1) q implies that the restricted projection π q+1 q : R q+1 → R q is a surjective submersion, that is, it defines a fibred manifold. Thus if we choose a section γ : R q → R q+1 and then always takeρ = γ(ρ), it follows immediately that the corresponding complements H ρ define a smooth distribution as required.
Example 3. Consider again the differential equation R 1 in Example 1. It is locally represented by the same equations as R (1) 1 , except that the integrability condition w t = v x is missing. The matrix of T ι for the system R 1 has eleven rows and eight columns-one column more than the symbol matrix for the system R (1) 1 has dimension 2 and is spanned by ∂ vx + ∂ wt and ∂ wx + ∂ vt . But the one-forms ω 1 , ω 2 and ω 3 (and their pull-backs are the same, and therefore the coordinate expressions for the vector fields X 1 and X 2 (and their representations 1 ) look alike (see Example 1 for their representations). The integrability condition w t = v x does not influence the results as it stems from the equality of the cross derivatives, u tx = v x and u xt = w t , and not from the prolongation of a lower order equation. Now consider for comparison the differential equation which is locally represented bỹ It arises from the system R The vector fields ∂ vx and ∂ wx generate the symbolÑ 1 ; any of its complements in V[R 1 ] is onedimensional and, as the dimension of the base space is two, none of them has the right dimension to be the horizontal space of a connection.
The reason for this is thatR 1 is not formally integrable, as the prolongation of the algebraic equation u = x leads to the integrability conditions u x = 1 and u t = 0. Projecting the prolonged equation gives Now the symbol vanishes, and so do the pull-backs of the contact forms: As the Lie brackets of ∂ x and ∂ t trivially vanish, TR Any n-dimensional complement H is obviously a transversal subdistribution of V[R q ], but not necessarily involutive. Conversely, any n-dimensional subdistribution H of V[R q ] is a possible choice as a complement. Choosing a "reference" complement H 0 with a basis (X i : 1 ≤ i ≤ n), a basis for any other complement H arises by adding some vertical fields to the vectors X i . We will follow this approach in the next section. For the remainder of this section we turn our attention to the choice of a convenient basis of V[R q ] that will facilitate our computations.
Let r := dim N q . As an intersection of two involutive distributions, the symbol N q is an involutive distribution, too. Hence, there exists a basis (Y k : 1 ≤ k ≤ r) for it whose Lie brackets vanish: where W is again an involutive distribution. It can be spanned by vector fields W 1 , W 2 , . . . , W s (where s = n k=1 β (k) q equals the number of principal derivatives) which are chosen such that we have [W a , W b ] = 0 for all 1 ≤ a, b ≤ s. In local coordinates, a particularly convenient choice for the fields Y k and W a exists. We first choose for any 1 ≤ k ≤ r a parametric derivative u α µ , that is (α, µ) / ∈ B, and set Y k := Y α µ := ι * (∂ u α µ ); then we choose for any 1 ≤ a ≤ s a principal derivative u α µ , that is (α, µ) ∈ B, and set W a := W α µ := ∂ u α µ . The reference complement H 0 is chosen as follows. Any basis of it must consist of n transversal contact fields. Since the fields C µ α are vertical, we can always use a basis (X 1 : 1 ≤ i ≤ n) of the formX with some coefficient functions ξ α iµ chosen such thatX i is tangential to R q . The fields C µ α also span the vertical bundle V π q q−1 and hence we may exploit the above decomposition for a further simplification of the basis. By subtracting from eachX i a suitable linear combination of the fields Y k spanning the symbol N q , we arrive at a basis (X i : As already mentioned above, the Vessiot distribution V[R q ] is not necessarily involutive. So it is not surprising that its structure equations are going to be important later. We may extend the above chosen basis (X i , Y k ) of V[R q ] to a basis of the derived Vessiot distribution, V ′ [R q ], by adding vector fields Z c , 1 ≤ c ≤ C : where, using (2), for each c we have some coefficients κ α cν such that Z c = κ α cν ∂ u α ν with |ν| = q − 1. By construction, the non-vanishing structure equations of V[R q ] take now the form for 1 ≤ i, j ≤ n and 1 ≤ k ≤ r, with smooth functions Θ c ij and Ξ c ik . (For the complete set of structure equations, we have to add [Y k , Y l ] = 0 for 1 ≤ k, l ≤ r.) Remark 2. Since the vector fields Z c , which appear on the right sides of the structure equations (26), may, for q = 1, span a proper subspace in ∂ u α : 1 ≤ α ≤ m , about the exact size of which we know nothing, we write them as linear combinations Z c =: κ α c ∂ u α . The structure equations then become Knowing the larger sets of coefficients Θ α ij , Ξ α ik , we can reconstruct the true structure coefficients Θ c ij , Ξ c ik by solving the overdetermined systems of linear equations This is always possible since the fields Z c are assumed to be part of a basis for the derived Vessiot distribution V ′ [R 1 ] and therefore linearly independent. Thus there exist some coefficient functions κ c α such that For our later proof of an existence theorem for integral distributions, we will have to analyse certain matrices with the coefficients Θ c ij and Ξ c ik as their entries. It turns out that this analysis becomes simpler, if we use the extended sets of coefficients Θ α ij and Ξ α ik instead.
For a first-order equation R 1 with Cartan normal form (8) satisfying the assumptions of Proposition 6 it is possible to perform this process explicitly. We choose as a reference complement H 0 the linear span of the vector fields for 1 ≤ i ≤ n. One verifies in a straightforward computation that (27)  Lemma 4. Let i < j, without loss of generality. Then we obtain for the extended set of structure coefficients Θ α ij in local coordinates on R 1 the following results: Proof . By a rather straightforward computation; see [6, pages 75 and 76].
We collect these coefficients into vectors Θ ij which have m rows each where the entries are ordered according to increasing α. It is useful to denote the symbol fields Y k = ι * (∂ Lemma 5. SetȲ k =:Ȳ β h , then the extended set of structure coefficients Ξ α ik in local coordinates on R 1 are Proof . Again by a rather straightforward computation; see [6, page 76].
We end this section with two technical remarks on how these coefficients are collected into matrices Ξ i . The examination of the ranks of these matrices Ξ i is basic for the proof of the existence Theorem 2 for integral distributions.
We collect the coefficients Ξ α ik into matrices Ξ i using i as the number of the matrix to which the entry Ξ α ik belongs, α as the row index of the entry and k as its column index. These matrices have m rows each, ordered according to increasing α; and they have r = dim N 1 columns each of which can be labelled by pairs (β, h) ∈ B or the symbol fieldsȲ k =Ȳ β h . More precisely, for 1 ≤ h ≤ n, we set Such a matrix with an upper index h collects all those Ξ α ik into a block where (α, i) ∈ B. For any 1 ≤ i ≤ n, we have m − β Furthermore, for 1 ≤ h ≤ n, we assemble the remaining terms Ξ α ik (which are those where (α, i) / ∈ B) in a matrix. As above, letȲ k =Ȳ β h , and denote any entry Ξ α ik accordingly, for the sake of introducing the following matrix, by h β Ξ α i . Now set columns. Let for any natural numbers a and b denote 0 a×b the a × b zero matrix. According to equation (29), we have Using the matrices [Ξ i ] h and [Ξ i ] h as blocks, we now build the matrix Taking into account what we have noted on its entries, this means A sketch of such a matrix Ξ i where the entries which may be different from zero are marked as shaded areas and −½ α (i) 1 as a diagonal line is given in Fig. 1.   Note that for a differential equation with constant coefficients all vectors Θ ij vanish and for a maximally over-determined equation there are no matrices Ξ i .
The unit block of α (i) 1 rows, ½ α (i) 1 , leads immediately to the estimate Example 5. For Example 1, we have n = 2 and therefore two matrices Ξ 1 and Ξ 2 . Their entries follow immediately from our results in Example 4. For i = 1, we have

Flat Vessiot connections
In this section, we develop an approach for constructing flat Vessiot connections which improves recent approaches [5,24,27]: we exploit the splitting of V[R q ] suggested by (23) to introduce convenient bases for integral distributions which yield structure equations that are particularly simple. Finally, we give necessary and sufficient conditions for Vessiot's approach to succeed.
Let R q locally be represented by the system Φ τ (x, u (q) ) = 0 where 1 ≤ τ ≤ t. Our goal is the construction of all n-dimensional transversal involutive subdistributions U within the Vessiot distribution V[R q ]. Taking for the Vessiot distribution the basis (X i , Y k ) where the vector fields Y k are the above mentioned basis of the symbol T ι(N q ) with vanishing Lie brackets and the vector fields X i are given in (25), we make for the basis (U i : 1 ≤ i ≤ n) of such a subdistribution U the ansatz . This ansatz follows naturally from our considerations so far, as the fields X i are transversal to the fibration over X and span a reference complement to the symbol and in N q all fields Y k are vertical.
Since the fields U i are in triangular form, the distribution U is involutive if, and only if, their Lie brackets vanish, and using the structure equations (26) this means: It follows from the definition of the fields Y k and Z c that they are linearly independent and so their coefficients must vanish for U to be involutive. Thus the Lie bracket (32) yields two sets of conditions for the coefficient functions ζ k i : a system of algebraic equations and a system of differential equations In the algebraic system (33) the true structure coefficients Θ c ij , Ξ c jk appear. For our subsequent analysis we follow Remark 2 and replace them by the extended set of coefficients Θ α ij , Ξ α jk . This corresponds to replacing (33) by an equivalent but larger linear system of equations which is, however, simpler to analyse.
Obviously, (33) is an inhomogeneous linear system to which any solution method for linear systems can be applied. We shall see in the next section that its structure (induced by the structure equations of the Vessiot distribution) allows us to decompose it into simpler subsystems which are considered step by step. Many papers on the Vessiot theory (see for example [5,24]) study at this stage a quadratic system which only by such a step-by-step approach can be reduced to a series of linear problems. The linearity of (33) is a simple consequence of our choice of a basis for V[R q ] which in turn exploits the natural splitting V[R q ] = N q ⊕ H.
Remark 5. The vector fields Y k lie in the Vessiot distribution V[R 1 ]. Thus, according to Proposition 5, the subdistribution U is an integral distribution if, and only if, the coefficients ζ i k satisfy the algebraic conditions (33). This observation permits us immediately to reduce the number of unknowns in our ansatz. Assume that we have values 1 ≤ i, j ≤ n and 1 ≤ α ≤ m such that both (α, i) and (α, j) are not contained in B, that is, u α i and u α j are both parametric derivatives (and thus obviously the second-order derivative u α ij , too). Then there exist two symbol fields Y k = ι * (∂ u α i ) and Y l = ι * (∂ u α j ). Now it follows from the coordinate form (3) of the contact map that the subdistribution U, spanned by the vector fields U h = X h + ζ k h Y k for 1 ≤ h ≤ n, can be an integral distribution if, and only if, for all 1 ≤ i < j ≤ n and 1 ≤ k, l ≤ r.
As the unknowns ζ k j may be understood as labels for the columns of the matrices Ξ h , this identification leads to a contraction of these matrices. The contracted matrices, denoted bŷ Ξ h , arise as follows: whenever ζ k j = ζ l i then the corresponding columns of Ξ h are added; their sum replaces the first of these columns, while the second column is dropped. Similarly, we introduce reduced vectorsζ h where the redundant components are left out. From now on, we always understand that in the equations above this reduction has been performed. Otherwise some results would not be correct (see Example 10 below).
For a more thorough outline of these technical details, see [6, Lemma 3.3.5 and Remark 3.3.6].
It is easy to show that the symbol of such an equation is always involutive and that the used coordinates are δ-regular. In order to simplify the notation, we use the following convention for the indices: 1 ≤ k, ℓ < n − r and n − r ≤ a, b ≤ n. Thus we may express our system in the concise form u a = φ a (x, u, u k ) and local coordinates on the submanifold R 1 are (x, u, u k ). The pull-back of the contact form ω = du − u i dx i generating the contact codistribution C 0 1 is ι * ω = du − u k dx k − φ a dx a and V[R 1 ] = X 1 , . . . ,X n ,Ȳ 1 , . . . ,Ȳ n−r−1 wherē The fieldsȲ k span the symbol N 1 and the fieldsX i our choice of a reference complement H 0 .
Now we make the above discussed ansatz U i = X i + ζ k i Y k for the generators of a transversal complement H. Modulo the Vessiot distribution V[R 1 ] we obtain for their Lie brackets The algebraic system (33) is now obtained by requiring that all the expressions in parentheses on the right hand sides vanish. Its solution is straightforward. The first subsystem (36a) implies the equalities ζ ℓ k = ζ k ℓ . This result was to be expected by the discussion in Remark 5: both u k and u ℓ are parametric derivatives for R 1 and thus we could have made this identification already in our ansatz for the complement. The second subsystem (36b) yields that ζ k a = ζ ℓ k Y ℓ (φ a ) + X k (φ a ). If we enter these results into the third subsystem (36c), then all unknowns ζ drop out and the solvability condition arises. Thus in this example the algebraic system (33) has a solution if, and only if, this condition is satisfied.
Comparing with the classical theory of such systems, one easily verifies that this solvability condition is equivalent to the vanishing of the Mayer or Jacobi bracket [u a − φ a , u b − φ b ] on the submanifold R 1 , which in turn is a necessary and sufficient condition for the formal integrability of the differential equation R 1 . Thus we may conclude that R 1 possesses n-dimensional integral distributions if, and only if, it is formally integrable (which in our case is also equivalent to R 1 being involutive).
Thus here involution can be decided solely on the basis of the algebraic system (33). We will show in the next section that this observation does not represent a special property of a very particular class of differential equations but a general feature of the Vessiot theory.

The existence theorem for integral distributions
Now the question arises, when the combined system (33), (34) has solutions? We begin by analysing the algebraic part (33). We use for this analysis a step-by-step approach originally proposed by Vessiot [28]. Our analysis will automatically reveal necessary and sufficient assumptions for it to succeed. As outlined in Remark 5, we replace ζ 2 byζ 2 since for the entries ζ . Thus we begin the construction of the integral distribution U by first choosing an arbitrary vector field U 1 and then aiming for another vector field U 2 such that [U 1 , U 2 ] ∈ T ι(V[R q ]). During the construction of the field U 2 we regard the components of the vector ζ 1 =ζ 1 as given parameters and the components ofζ 2 as the only unknowns of the system Since the components ofζ 1 are not considered unknowns, the system (37) must not lead to any restrictions for the coefficientsζ k 1 . Obviously, this is the case if, and only if, rankΞ 1 = rank (Ξ 1Ξ2 ). (38) Assuming that (38) holds, the system (37) is solvable if, and only if, it satisfies the augmented rank condition When we have succeeded in constructing U 2 , the next step is to seek a further vector field . Now the components of both vectorŝ ζ 1 andζ 2 are regarded as given, and the components ofζ 3 are regarded as the unknowns of the systemΞ Now this system is not to restrict the components of bothζ 1 andζ 2 any further. The interrelations between theζ i following from the condition (35) for the existence of integral distributions, given in Remark 5, are taken care of by contracting Ξ 3 intoΞ 3 . This implies that now the rank condition rank Ξ 1 has to be satisfied. If it is, then for 1 holds. Now we proceed by iteration. Given j − 1 vector fields U 1 , U 2 , . . . , U j−1 of the required form spanning an involutive subdistribution of T ι(V[R 1 ]), we construct the next vector field U j by solving the system Ξ 1ζj =Ξ jζ1 − Θ 1j , · · · · · · · · · · · · · · · · · · · · · (40) Again we consider only the components of the vectorζ j unknowns, and the system (40) must not imply any further restrictions on the components of the vectorsζ i for 1 ≤ i < j. The corresponding rank condition is Assuming that it holds, the equations (40) are solvable and yield solutions for the components ofζ j if, and only if, the system satisfies the augmented rank condition (42) Remark 6. We can simplify our calculations for a differential equation R q , represented by u α µ = φ α µ (x, u,û (q) ) where each equation is solved for a principal derivative (all pairwise different) u α µ with |µ| = q, and whereû (q) denotes the set of the remaining, thus parametric, derivatives of order q and less. (Example 1 discusses a system of that kind.) Now we can use the local coordinates on R q . For the generators of the symbol N q we may choose for (α, µ) ∈ B the vector fieldsȲ α µ = ∂ u α µ ; as a basis for the complement H ⊆ V[R q ], we can take the vector fields which satisfy equation (27). For the vector fields W a which appear in equation (25) we may choose the contact vector fields C µ α where (α, µ) ∈ B. The further procedure, solving first the structure equations, which now take the simple form (26), to find the generators of V ′ (R q ) and then solving equation (33) for the coefficient functions ζ j , is the same as in the case where the representation is not in the above solved form.
Example 7. Before we proceed with our theoretical analysis, let us demonstrate with a concrete differential equation that in general we cannot expect that the above outlined step-by-step construction of integral distributions works. In order to keep the size of the example reasonably small, we use the second-order equation 1 R 2 defined by the system u xx = αu and u yy = βu with two real constants α, β. Its symbol N 2 is not involutive. However, one easily proves that R 2 is formally integrable for arbitrary choices of the constants α, β. One readily computes that the Vessiot distribution V[R 2 ] is generated by the following three vector fields on R 2 : They yield as structure equations for V[R 2 ]: For the construction of a two-dimensional integral distribution U ⊂ V[R 2 ] we make as above the ansatz U i = X i + ζ 1 i Y 1 with two coefficients ζ 1 i . As we want to perform a step-by-step construction, we assume that we have chosen some fixed value for ζ 1 1 and try now to determine ζ 1 . Evaluation of the Lie bracket yields the equation A necessary condition for the vanishing of the right side is that ζ 1 1 = αu y . Hence one cannot choose this coefficient arbitrarily, as we assumed, but (43) determines both functions ζ 1 i uniquely. Note that the conditions on the coefficients ζ 1 i imposed by (43) are trivially solvable and thus R 2 possesses an integral distribution (which is even involutive). The problem is that it could not be constructed systematically with the above outlined step-by-step process.
The following theorem links the satisfaction of the rank conditions (41) and (42), and thus the solvability of the algebraic system (33) by the above described step-by-step process, with intrinsic properties of the differential equation R q and its symbol N q . It represents an existence theorem for integral distributions.
Theorem 2. Assume that δ-regular coordinates have been chosen for the local representation of the differential equation R q . Then the rank condition (41) is satisfied for all 1 ≤ j ≤ n if, and only if, the symbol N q is involutive. The augmented rank condition (42) holds for all 1 ≤ j ≤ n if, and only if, the differential equation R q is involutive.
The proof of Theorem 2 is given in Appendix A, as parts of it are rather technical. First, we transform the differential equation into an equivalent first-order system with a representation in the reduced Cartan normal form (10) and the same Cartan characters. This is always possible; see, for example, [23,Proposition A.3.1]. The basic idea of the proof is then that the rank condition (41) is equivalent to the vanishing of the obstructions to involution of the symbol in Lemma 2, and the augmented rank condition (42) is equivalent to the vanishing of the remaining integrability conditions (recall that Lemma 2 gives all obstructions to involution only if the used coordinates are δ-regular so that this assumption is necessary in Theorem 2). The concrete implementation of this idea requires some technical considerations concerning the transformation of the matrices (41) and (42) into row echelon form, working out their contractions and analysing the interrelation between these operations.
Example 8. We demonstrate the role of δ-regularity of the coordinates with the wave equation u xy = 0 in characteristic coordinates which are not δ-regular. A straightforward computation yields as generators for the Vessiot distribution the fields X 1 = C 2 , Y 1 = ∂ uxx and Y 2 = ∂ uyy . Following our approach, we make the ansatz Evaluation of the Lie bracket [U 1 , U 2 ] yields now the algebraic equations ζ 2 1 = 0 and ζ 1 2 = 0 and the differential equations U 1 (ζ 1 2 ) − U 2 (ζ 1 1 ) = 0 and U 1 (ζ 2 2 ) − U 2 (ζ 2 1 ) = 0. As in Example 7, we see that the step-by-step process is broken, since we obtain a condition on the coefficient ζ 2 1 = 0 which we should be able to consider a parameter here.
At this point, we have proven that integral distributions within the Vessiot distribution exist if, and only if, the algebraic conditions (33) are solvable, and that this is equivalent to the augmented rank condition (41) being satisfied. This in turn is the case precisely if the differential equation is involutive. Thus we have characterised the existence of Vessiot connections for Vessiot's [28] step-by-step approach.

The existence theorem for flat Vessiot connections
There remains to analyse the solvability, if we add the differential system (34). Its solvability is equivalent to the existence of flat Vessiot connections in that each flat Vessiot connection of R 1 corresponds to a solution of the combined system (33), (34). We first note that the set of differential conditions (34) alone is again an involutive system. Proposition 7. The differential conditions (34) alone represent an involutive differential equation of first order.
The proof, which is somewhat technical, is given in Appendix B. It is based on an evaluation of the Jacobi identity for the Lie brackets of the fields U i .
If the original equation R 1 is analytic, then the quasi-linear system (34) is analytic, too. Thus we may apply the Cartan-Kähler theorem (Theorem 1) to it which guarantees the existence of solutions. The problem is that the combined system (33), (34) is not necessarily involutive, as the prolongation of the algebraic equations (33) may lead to additional differential equations. Instead of analysing the effect of these integrability conditions, we proceed as follows. If we assume that R 1 is involutive, then we know from Theorem 2 that the algebraic equations (33) are solvable. Actually, the proof of the Theorem 2 produces an explicit row echelon form of the system matrix given on the right hand side of equation (41). Now we use the interrelations between the unknowns ζ (α,h) ℓ , where (α, h) ∈ B, to eliminate in the differential conditions (34) some of them by expressing them as linear combinations of the remaining ones. Let 1 ≤ i < j ≤ n; then for all (α, i) where β while for all (α, i) where β (j) 1 + 1 ≤ α ≤ m the correspondence (35), given in Remark 5, holds again (see [6,Corollary 3.3.23] for the detailed calculation). In this way, we plug the algebraic conditions into the differential conditions and thus get rid of them.
Example 9. For Example 1, where i = 1 and j = 2, we deduce for The rank conditions (38) and (39) are satisfied because v x − w t = 0, so that the matrix (Ξ 1 Ξ 2 − Θ 12 ) has a vanishing first row. (There is no contraction of matrices for j = 2.) There are no interrelations like those in equation (35) because there is no index value α for which both ι * (∂ u α 1 ) and ι * (∂ u α 2 ) would be symbol fields. The differential conditions (34) are − U 2 (ζ wx 1 ) = 0, and by substituting ζ vx 2 = ζ vx 1 and ζ wx 2 = ζ wx 1 , we can drop the algebraic conditions from the system. Since n = 2, there is no further procedure.
We can now prove the following existence theorem for flat Vessiot connections.
Theorem 3. Assume that δ-regular coordinates have been chosen for the local representation of the analytic differential equation R 1 . Then the combined system (33), (34) is solvable.
The idea of the proof is this: following the strategy we have just outlined, we eliminate some of the unknownsζ k i . Because of the simple structure of (34), it turns out that we must take a closer look only at those equations where the leading derivative is of one of the unknowns we eliminate. A somewhat lengthy but straightforward computation shows that these equations actually vanish. The remaining equations still form an involutive system. Thus we eventually arrive at an analytic involutive differential equation for the coefficient functionsζ k i which is solvable according to the Cartan-Kähler theorem. The details are worked out in Appendix C.
Example 10. Consider the first-order equation in the five independent variables x, y, z, s, t and the three dependent variables u, v, w. It is formally integrable, and its symbol is involutive with dim N 1 = 8. Thus R 1 is an involutive equation. For the matrices Ξ i , all of which are 3 × 8-matrices, we find For the first two steps in the construction of the fields U i , the rank conditions are trivially satisfied even for the non-contracted matrices. But not so in the third step where we have in the row echelon form of the arising 9×32-matrix in the 7th row zero entries throughout except in the 12th column (where we have −2) and in the 17th column (where we have 2). As a consequence, we obtain the equality ζ 4 1 = ζ 1 2 and the rank condition for this step does not hold. However, since both u x and u y are parametric derivatives and in our ordering Y 1 = ι * (∂ ux ) and Y 4 = ι * (∂ uy ), this equality is already taken into account in our reduced ansatz and for the matricesΞ i the rank condition is satisfied.
Note that the rank condition is first violated when the rank reaches the symbol dimension (which is 8). From then on, the rank of the left matrix in (41) stagnates at dim N 1 while the rank of the augmented matrix rises further. The entries breaking the rank condition differ by their sign, while their corresponding coefficients in Lemma 2 are collected into one sum and thus vanish. Here it shows that contracting the matrices, as explained in Remark 5, is necessary.

Conclusions
Vessiot's [28] original motivation for the introduction of his theory was to provide an alternative proof of the Cartan-Kähler theorem (the same holds for Stormark's presentation [24] of the theory). If one takes this point of view, then one may say that our proof of Theorem 3 is a "cheat", as it uses the Cartan-Kähler theorem instead of proving it. However, in our opinion, this point of view is the main reason why the corresponding proofs in [24,28] are so difficult to read (and actually incomplete as they neglect that the involutivity of the system is a necessary condition). The direct proof of the Cartan-Kähler theorem (for differential equations) given in textbooks on the formal theory like [19,23] is much simpler and more transparent; in particular, it makes clear where involution (as opposed to mere formal integrability) is needed. Thus, if the only goal consisted of proving such an existence and uniqueness theorem, then there would be no need to bother with Vessiot's theory.
We believe that the real value of Vessiot's theory lies in the fact that it provides via the distribution V[R q ] an additional geometric structure on an involutive differential equation R q which is very useful for the further analysis of the equation, that is, after its solvability has been established. One possible application is the investigation of certain forms of singular behaviour of solutions. In fact, the classical works of Arnold and collaborators (see [2] for an elementary introduction) on implicit ordinary differential equations are based on the Vessiot distribution (without using this terminology); for some further works also in the context of partial differential equations or numerical analysis see, for example, [16,25,26]. The basic idea here is that such behaviour mainly stems from the fact that at some points the considered involutive distributions cease to be transversal to the fibration R q → X so that integral manifolds can no longer be interpreted as prolonged solutions in the classical sense. However, it often makes sense to consider them a generalised form of (potentially multi-valued) solutions.
Another application concerns symmetry theory. In classical symmetry theory one always considers the differential equation R q as a submanifold of J q π and looks then for diffeomorphisms of J q π which are (i) compatible with the contact structure of J q π and (ii) leave R q invariant. It was only fairly late realized [1] that this approach yields only what is now called external symmetries and that in general further useful symmetries may exist. Using Vessiot's theory, we may represent a differential equation as the pair (R q , V[R q ]), that is, as a manifold together with a distribution on it -without any recourse to an ambient space (all relevant properties of J q π are captured in the Vessiot distribution V[R q ]). This point of view yields automatically all internal symmetries (and is equivalent to the symmetry theory of exterior differential systems).
We mentioned already in the Introduction that Vessiot's theory takes an intermediate position between the formal theory and the theory of exterior differential systems. It allows for the transfer of ideas from the latter to the former one without the need to rewrite a differential equation as an exterior system. For example, [23, Section 9.5] contains a formulation of prolongation structures without exterior forms using the Vessiot distribution instead 2 . It should also be possible to give an explicit proof of the equivalence of the formal theory and the Cartan-Kähler theory on the basis of the results presented here.
Finally, we comment on the practical side of the Vessiot theory. Throughout this article, we made a number of assumptions on the treated differential equation: we assumed that it is first rewritten as a first-order system, then that all present algebraic equations are explicitly solved and finally that each equation is solved for its principal derivative. While, from a theoretical point of view, these operations are always possible under fairly mild regularity assumptions, usually they cannot be performed effectively. However, the made assumptions are only needed in order to be able to prove results like Lemma 2 on the explicit form of the obstructions to involution. These expressions are already bad enough with our simplifying assumptions; for more general systems they would hardly be manageable.
A concrete differential equation need not be transformed to such a special normal form in order to apply the Vessiot theory. Even for fully implicit equations of arbitrary order, the determination of all integral distributions requires only linear algebra and is easily implemented in a computer algebra system. In fact, already Fackerell [5] reported of an implementation of Vessiot's theory. Within a recent project thesis [8], Globke implemented large parts of the theory in MuPAD. The problems and limitations are here exactly the same as in implementations of the Cartan-Kähler theory (see, for example, [10]).

A Proof of Theorem 2
In this section of the appendix, we prove Theorem 2, the existence theorem for integral distributions. It is in principle by straightforward matrix calculation and involves a tedious distinction of several cases and subcases: we have to compare the entries in the matrix on the right hand side of equation (42) for step j after turning it into row echelon form with the integrability conditions and the obstructions to involution as they are given in Lemma 2. For this purpose, we fix 1 ≤ i < j and consider the blockΞ i and the entries to its right in the complete matrix.
Since for a differential equation with an involutive symbol the obstructions to involution vanish and for an involutive differential equation the integrability conditions vanish, too, it follows that the augmented rank condition, stated in equation (42), is equivalent to the differential equation being involutive, which in turn is the case if, and only if, the algebraic conditions (33) are satisfied, which is necessary and sufficient for the existence of integral distributions within the Vessiot distribution.
If the order of the differential equation is q > 1, transform it into an equivalent first-order equation; for the details of this procedure, see [6, Subsection 2.5.1] or [23]. If R q is involutive, then so is R 1 (see [23] for a straightforward proof of this). We assume this first-order equation is represented in reduced Cartan normal form (10). We first prove the rank condition (41) for the homogeneous system. (The proof for the augmented rank condition (42) follows.) We proceed in two steps: We consider the complete matrix at the jth step, (It differs from the matrix on the right hand side of equation (41) in that it is not yet contracted.) The complete matrix at the jth step is built from (j − 1)j blocks: the stack of j − 1 matrices Ξ i , 1 ≤ i ≤ j − 1, on the left, with each of the Ξ i having another j − 1 blocks to its right, the ith of which being Ξ j and all the others being zero. For easier reference, let, for 1 ≤ i, k ≤ j − 1, be [i, k] the kth block right of a Ξ i . Then, for all 1 ≤ i, k ≤ j − 1, we have For convenience, we set [i, 0] := Ξ i . Let, for 1 ≤ g ≤ h ≤ n,  We first turn the complete matrix into row echelon form. Then we contract columns and consider the effect.
The matrices Ξ i are shown in equation (31) and have a block-structure of the form where A, R, S, T are adequate natural numbers or zero, and the blocks are the S × S unit matrix ½ S or the zero matrix 0 S×R . (Here the index r is just some index and not supposed to be the dimension of any symbol.) Then the substitution transforms the matrix (g ij ) into For the transformation of the complete matrix into row echelon form, we use this obvious method. Let 1 < j ≤ n be given. Choose 1 ≤ i < j. Now we eliminate all the non-trivial entries in the rows We turn to the block Non-trivial entries therein may appear according to equation (31) only for 1 ≤ h ≤ i in the blocks Then it follows that the substitution (47) leaves (a αr ) unchanged, turns (b as ) into zero as required and makes (c αt ) into (c αt + S s=1 b αs d st ). According to (48b) the entries For h = l, we have [h, h] = Ξ j . The structure of Ξ j , which is (31) with i replaced by j, implies that non-vanishing entries are possible in the blocks m This matrix has A = β for all 1 ≤ k ≤ j. Some of the β = m summands vanish because of the special entries (49). As a consequence, from β (j) 1 + 1 on, of all the summands C h γ (φ α i )C k δ (φ γ j ) in (51) at most one remains: none for k = j, exactly one for k = j, namely the one for γ = β (j) Using the Kronecker-delta, for all β These are the terms S s=1 b αs d st in the case h = l when 1 ≤ k ≤ j. Now for the subcase j + 1 ≤ k ≤ n: these blocks m again on account of the structure of Ξ j , which is (31) with i replaced by j, are zero (if they exist at all; they do not for k with α According to the structure of Ξ j , which is (31) with i replaced by j, the only entries c αt within the block which may not vanish are, for 1 ≤ k ≤ j, those of the form −C k δ (φ α j ). They make up the block So for any row index α and any column index δ = β (k) as the most general form of an entry. Since Since 1 ≤ α ≤ β 1. Let h > i. Then according to the structure of Ξ i , (31), all the C h γ (φ α i ) and C h δ (φ α i ) in (52) vanish. Since h = i, δ hi = 0, thus all of (52) vanishes.
2. Let h = i. We shall consider several subcases for h = i, and for h < i after that, which may be labelled by second-order derivatives u δ hk where h and k are combined as shown in It turns out that (52) is the common form of all the sums that appear in the squared brackets of Lemma 2 as the coefficients of second-order derivatives u δ hk . Not only can the case distinctions of the following argument, 2.(a), 2.(b) and 3.(a) to 3.(d2), be labelled by these u δ hk , but in fact a case labelled u δ hk is dealt with by using the fact that according to Lemma 2 the coefficient in square brackets of that same u δ hk vanishes for an involutive system. The correspondence of the cases and the terms u δ hk is given in the caption of Fig. 2.
The cases h = i with subcase k < i and h < i with subcase k < i need not be considered because the columns of the complete matrix are to be contracted when proving the rank condition. This contraction concerns those columns of the complete matrix which are labelled by the same second-order derivative, and each second-order derivative is used exactly once in the argument.
(a) For the first subcase of h = i let i ≤ k < j. Then (52) becomes This vanishes for a system with an involutive symbol according to Lemma 2,line (17).
(b) For the second subcase of h = i choose k = j. Then (52) becomes According to Lemma 2, line (19), this vanishes for a system with an involutive symbol.
3. Let h < i. We consider several subcases with regard to the relation between h < i, j and k.
(a) First choose k = j. Then (52) becomes According to Lemma 2, line (18), this vanishes for a system with an involutive symbol.
(c) For the third subcase choose k = i. Then (52) becomes where β (k) 1 + 1 ≤ δ ≤ m. Since we have h < k = i, for any such δ the cross derivative u δ ih = u δ hi labels two columns in the complete matrix: the one with label δ in where β (h) 1 + 1 ≤ δ ≤ m. In the current subcase, (54) becomes Since the columns of both (53) and (55) are labelled by the same second-order derivatives, we have to contract them, which means adding their new entries. For all β which, for a system with an involutive symbol, vanishes according to Lemma 2, line (15a), and for all β which, for a system with an involutive symbol, again vanishes according to Lemma 2, line (15b).
(d) For the fourth subcase choose k < i. Under this assumption, we have to distinguish two further subcases: k = h and k < h. The subcase k > h need not be considered since u δ kh = u δ hk and any second-order derivative is used only once to label a column in the contracted matrix.
(d2) At last consider h < k < i. Then k < j, and (52) becomes where β (k) 1 + 1 ≤ δ ≤ m. Since we have h < k < i, for any such δ the cross derivative u δ kh = u δ hk labels two columns in the complete matrix: the one with label δ in where β (h) 1 + 1 ≤ δ ≤ m. In the current subcase, (57) becomes where β (h) 1 + 1 ≤ δ ≤ m. Since the columns of both (56) and (58) are labelled by the same second-order derivatives, we have to contract them, which means adding their new entries. This yields 1 and vanishes according to Lemma 2, line (14a), and for β (k) 1 + 1 ≤ δ ≤ m and vanishes according to Lemma 2, line (14b).
What we have shown is that in the contracted matrix all the entries without some entry −1 in a negative unit block to their left may be eliminated by elementary row transformations if the equation has an involutive symbol. Thus under this assumption the rank condition (41) is satisfied. Now for the augmented rank condition, equation (42). To transform the augmented complete matrix into row echelon form, we use for each j and i the same procedure as for the transformation of the non-augmented complete matrix, except that now the matrices (c αt ) and (d st ) are augmented by one more column each as follows. Fix 1 < j ≤ n. Let 1 ≤ i < j. Then, according to the structure of the augmented complete matrix, for its transformation into row echelon form we have to eliminate for 1 ≤ h ≤ i the entries in the blocks [Ξ i ] h , as we did for the non-augmented complete matrix given in equation (45). We have to consider the effect of these transformations on the additional entries which make up the rightmost column in the augmented complete matrix. These are −Θ 1 ij , −Θ 2 ij , . . . , −Θ m ij , given in equation (28). Of these entries, only −Θ 1 ij , −Θ 2 ij , . . . , −Θ hj as its rightmost column, in accordance with the structure of the augmented complete matrix. Then the substitution (47) yields as the transformed entries c αt + S s=1 b αs d st the same as for the transformation of the non-augmented matrix except, of course, for the new last column. Now let denote t the index of this last column. Fix some row index 1 ≤ α ≤ β C h γ (φ α i )Θ γ hj .
The entries Θ γ hj can be taken from equation (28). Since β as the summand for each h to be added to c αt = −Θ α ij . Since 1 ≤ α ≤ β j (φ α i ). Therefore the entry c αt transforms into which is the integrability condition in line (12) of Lemma 2, except for the sign. The difference in sign comes from the fact that the augmented complete matrix describes the system of equations (40), where the entries from the row with index α in the matrix Ξ i are on the opposite side from the entries from the row with index α in the matrix Ξ j and the inhomogeneous term −Θ α ij , while in Lemma 2 for each i, j and α the corresponding equation is set to zero, if the differential equation is involutive.
This means that the augmented rank condition holds for all 1 ≤ i < j ≤ n if, and only if, the equation has an involutive symbol and is formally integrable.

B Proof of Proposition 7
We are going to show that the system (34) is involutive. Consider the basis (U i : 1 ≤ i ≤ n) of the distribution U given by U i = X i + ζ p i Y p with yet undetermined coefficient functions ζ p i ∈ F(R 1 ). For all 1 ≤ i ≤ n and 1 ≤ p ≤ r, the independent variables of the functions ζ p i are the coordinates on R 1 , which are x, u and all u α h such that (α, h) ∈ B. To apply a vector field U j = ∂ x j + · · · to a function ζ p i includes a derivation with respect to x j . We order the independent variables such that if j > i, then x j is greater than x i , and each x i is greater than all the variables u α and u α h where (α, h) ∈ B. For any equation H p ij within the system (34), the application of the vector field U j = ∂ x j + · · · to ζ p i yields ∂ζ p i /∂x j as the leader of that equation; therefore equation H p ij is of class j, and the equations of maximal class are H p in ; the equations of second highest class in the system are H p in−1 and so on. There are only equations H p ij of a class indicated by some index 2 ≤ j ≤ n.
From the Jacobi identity for vector fields U i , U j and U k where 1 ≤ i < j < k ≤ n, we have The structure equations (32) for the vector fields U h and the definitions of G c ij and H p ij in equations (33) and (34) imply that this is The combined system (33, 34) means that all G c ab = 0 and all H p ab = 0 which implies that U is involutive, which it is, being in triangular form, exactly if all [U a , U b ] = 0. This leaves only As part of a basis for V ′ [R q ], the vector fields Z c and Y p are linearly independent, which means their coefficients must vanish individually. So in particular U i (H p jk ) + U j (H p ki ) + U k (H p ij ) = 0. Under the assumption i < j < k, the term U k (H p ij ) contains derivations with respect to x k of U i (ζ p j ) and U j (ζ p i ). Thus, according to our order, this is a non-multiplicative prolongation, and the remaining terms are multiplicative prolongations. But since any non-multiplicative prolongation within the system (34) must be of such a form, it is a linear combination of multiplicative prolongations. Therefore, no integrability conditions arise from cross-derivatives (and none arise from a prolongation of lower order equations since all equations of the system are of first order).