The Variational Bi-Complex for Systems of Semi-Linear Hyperbolic PDEs in Three Variables

This paper extends, to a class of systems of semi-linear hyperbolic second order PDEs in three variables, the geometric study of a single nonlinear hyperbolic PDE in the plane as presented in [Anderson I.M., Kamran N., Duke Math. J. 87 (1997), 265-319]. The constrained variational bi-complex is introduced and used to define form-valued conservation laws. A method for generating conservation laws from solutions to the adjoint of the linearized system associated to a system of PDEs is given. Finally, Darboux integrability for a system of three equations is discussed and a method for generating infinitely many conservation laws for such systems is described.


Introduction
This paper belongs to the field known broadly as the geometric study of partial differential equations, which seeks to understand differential equations through the study of properties which remain invariant under particular groups of transformations. The subject has its roots in the foundational works of Lie, Darboux, Cartan and others. It was Cartan who recast partial differential equations geometrically as exterior differential systems. In doing so, the solutions to partial differential equations were realized as integral manifolds of corresponding exterior differential systems. More recently, the field of geometric PDEs has undergone a number of important developments. These include efforts to obtain explicit solutions and solution algorithms to specific classes of PDEs [3,35], the investigation of links between PDEs and the geometry of the submanifolds which their level sets define [22,23,31], as well as the study and computation of invariants such as conservation laws [4,14,34]. It is this last area, conservation laws, with which this paper is concerned.
The geometric approach to conservation laws has an extensive history and literature, which we will not attempt to describe in any detail, referring the reader instead to [24] for a comprehensive account of the subject. We will however provide a brief, non-exhaustive overview of some of the themes and contributions that are relevant to this paper. Our approach to conservation laws finds its origin in the study of the cohomology determined by the C-spectral sequence introduced by Vinogradov in [32] and [33]. In addition to being used to characterize conservation laws in terms of cohomology, this construction has facilitated work in other important aspects of the study of differential equations including the inverse problem of the calculus of variations (see for example [5,15]) and Euler-Lagrange operators as studied by Tulczyjew in [30]. Tsujishita and Duzhin built upon the work of Vinogradov to study conservation laws of the BBM equation [17], and Tsujishita expanded the applications of the C-spectral sequence to include topics such as the study of characteristic classes and Gel'fand-Fuks cohomology [28] and conservation laws of the Klein-Gordon equation [27].
In [1], Anderson gives a comprehensive treatment of the variational bi-complex, which emerged out of the works on the C-spectral sequence mentioned above, and the horizontal cohomology of which will serve as the natural framework for our present study of conservation laws. It should be pointed out that the variational bi-complex also lends itself to the study of various other topics such as the equivariant version of the inverse problem of the calculus of variations, Riemannian structures, and the method of Darboux integrability for a scalar second order PDE. A large body of work has been amassed in the study of conservation laws by utilizing the variational bi-complex. Of particular importance for the present discussion, we note that Anderson and Kamran performed an extensive study of the conservation laws of hyperbolic scalar second order PDEs in the plane in [4].
In the contemporaneous work of Bryant and Griffiths [6,7], conservation laws were studied from the distinct, yet related, perspective of the characteristic cohomology of exterior differential systems. In particular, local invariants of an exterior differential system are shown to govern the properties of the system's characteristic cohomology. This approach carries the study of the variational bi-complex and the C-spectral sequence into the realm of exterior differential systems, where the independent and dependent variables of a system of partial differential equations are treated equally.
To demonstrate the connection between this approach and that of Vinogradov, we refer the reader to Vinogradov's "two line theorem" [32] which indicates that a system of PDEs of Cauchy-Kovalevskaya type in n independent variables will have trivial horizontal cohomology of horizontal degree ≤ n − 2 in the associated variational bi-complex. This result can be recovered from a fundamental theorem of [6] regarding characteristic cohomology (see [21,Section 10.4]). It can also be generalized by using Vinogradov's spectral sequence [29] as well as characteristic cohomology techniques [6]. The literature on characteristic cohomology is too extensive to discuss in detail, but we will note that the notion of a hyperbolic exterior differential system was introduced and studied in [8] and [9]. Additional examples of work in this field include that of Clelland [14] and Wang [34], each of whom studied conservation laws using the approach of characteristic cohomology, the former studying second order parabolic PDEs in one dependent and three independent variables, and the latter considering third order scalar evolution equations.
As mentioned previously, our study of conservation laws will take place in the setting of the variational bi-complex. In this framework, differential forms on the jet bundle of infinite order J ∞ (E) of a fibered manifold π : E → M are bi-graded and the (unconstrained) variational bi-complex is defined using the exterior derivative split into horizontal and vertical components, denoted by d H and d V respectively. The constrained variational bi-complex is associated to a given partial differential equation or system of equations, R, by taking the pullback of the unconstrained variational bi-complex to the infinite prolongation of the equation manifold, R ∞ . Then the study of equivalence classes of conservation laws of a given PDE (or system of PDEs) corresponds to the study of the horizontal cohomology of this constrained bi-complex, denoted by H r,s (R ∞ , d H ) where (r, s) is referred to as the bi-degree of the conservation law. In this context, a classical conservation law will have the form where d H ω = 0 and the M i are functions defined on R ∞ . We call ω a trivial conservation law if there exists an (r − 1, 0) form γ on R ∞ such that d H γ = ω. In other words, the classical conservation laws determine the cohomology classes of H r,0 R ∞ , d H where two conservation laws lie in the same cohomology class if they differ by a trivial conservation law. The notion of a higher-dimensional, or contact form-valued, conservation law will then arise when ω is a representative of a cohomology class for s ≥ 1. In this case the terms M i would be type (0, s) contact forms. Analogous to the classical case, ω ∈ H r,s R ∞ , d H is said to be trivial if there exists a form γ of bi-degree (r − 1, s) such that d H γ = ω. To illustrate these ideas, we provide the following simple example in two independent variables: Example 1.1. Consider the Liouville equation An example of a classical conservation law for this equation is given by the cohomology class u xx − 1 2 u 2 x dx ∈ H 1,0 R ∞ , d H , while an example of a contact form valued conservation law would be [(θ xx − u x θ x ) ∧ dx] ∈ H 1,1 R ∞ , d H , where θ I = du I − p j=1 u I,j dx j .
We will now restrict our attention to the particular class of PDEs with which this paper is concerned. Specifically, we undertake the study of involutive systems of three nonlinear, hyperbolic equations of the following form: Note that the expression f ij depends on all three independent variables, but only depends on the derivatives of u with respect to x i and x j . Systems of partial differential equations within this class appear in several interesting contexts. Linear systems of this type arise in the parametrization of Cartan submanifolds in Euclidean space as carried out by Kamran and Tenenblat [22]. This work built upon that of Chern, who in [12,13] gave a generalization of the Laplace transformation to a class of n-dimensional submanifolds in projective space, which he termed Cartan submanifolds as they had previously been studied by Cartan in [11]. The submanifolds in Chern's study admit a parametrization by a conjugate net, which in Euclidean space implies that the functions giving the parametrization satisfy an overdetermined system of second order PDEs, These fit into the class of systems studied in [22]. Other examples of applications of systems of the form (1.1) include the study of semi-Hamiltonian systems of hydrodynamic type (see [21,26], and [16]), and (1.1) are also a special case of the nonlinear Darboux-Manakov-Zakharov systems studied by Vassiliou in [31]. We proceed to give a synopsis of how this paper is organized. An overview of the essential background material needed to set the stage for the subsequent sections is carried out in Section 2 which includes the theory of jet bundles and the variational bi-complex. In Section 3 we develop the Laplace adapted coframe which provides the framework for our study and describe a method for generating (2, s) conservation laws for systems (1.1). Darboux integrability is discussed in Section 4, followed by concluding remarks and directions for future research in Section 5. The systems (1.1) are described in the context of Cartan's structural classification [10] of involutive systems of three equations in one dependent and three independent variables in Appendix A and the full structure equations and Lie bracket congruences for the Laplace adapted coframe are given in Appendix B.
In Section 3 the construction of the Laplace adapted coframe utilizes the methods presented in [4], where a version of the Laplace transformation that applies to the form-valued linearization of a PDE of the form F (x, y, u, u x , u y , u xx , u xy , u yy ) = 0 (1.2) is developed. Moving from a single equation (1.2) to a system of equations (1.1) introduces a set of nonlinear integrability conditions which must be satisfied in order for the system to be involutive. These conditions will be analogous to the integrability conditions described in [22], with the partial derivatives having been replaced by total derivatives. Once linearized, the system may be analyzed using a generalization of the well-known classical Laplace method, used for integrating a single, linear hyperbolic PDE. This method was adapted to the study of systems of linear hyperbolic PDEs in n independent and one dependent variable in [22]. It was subsequently extended to the vector-valued case in [35]. In this paper, a form-valued version of the generalized Laplace transform is used to investigate the systems (1.1). Following the detailed computation of the essential structure equations and Lie bracket congruences for the Laplace adapted coframe which are provided in Appendix B, we are in a position to state and prove the first main result of Section 3, which is a structure theorem for the (2, s) conservation laws for systems (1.1).
In order to motivate the statement of this theorem, we will first describe an analogous result concerning classical conservation laws. The reader may refer to [24] for a detailed exposition on this topic, as all computations and proofs will be omitted here. Consider a system of differential equations . . , u α I with 1 ≤ i ≤ n, 1 ≤ α ≤ q and |I| ≤ k are local jet coordinates. A conservation law for this system is an n-tuple P = P 1 x, u (k) , . . . , P n x, u (k) whose total divergence vanishes identically for all solutions of (1.3). That is, where R ij = −R ji such that, taking into account R and its prolongations, Then two conservation laws, P andP , are considered equivalent if their difference P −P is a trivial conservation law. It can be shown using integration by parts that if R is totally nondegenerate, i.e., the system of equations and all its prolongations are of maximal rank and locally solvable, then any conservation law P is equivalent to a conservation lawP whose total divergence can be written in the form where the multipliers Q ν = Q ν x, u (k) are C ∞ functions on the infinite jet bundle. The m-tuple Q = (Q 1 , . . . , Q m ) is referred to as the characteristic of the conservation law P . A characteristic is trivial if it vanishes for all solutions of the system of equations, and two characteristics are said to be equivalent if they differ by a trivial characteristic. As stated in [24], two conservation laws P andP are equivalent if and only if their characteristics are equivalent. So it is natural to suspect that the study of characteristics and conservation laws go hand in hand. Applying the Euler-Lagrange operator, to both sides of (1.5), one obtains the identity where L * Q is the adjoint of the formal Fréchet derivative of the differential operator Q, and likewise for L * ∆ . Since the system of equations specifies ∆ = 0, the identity (1.6) leads to the following conclusion: Theorem 1.2. The characteristic Q ν x, u (k) , 1 ≤ ν ≤ m, of any conservation law (1.4) for a totally non-degenerate system of differential equations (1.3) lies in the kernel of L * ∆ , that is L * ∆ (Q) = 0.
Then the first main result of this paper, Theorem 3.15, which we will state in full presently, can be seen as an analogous result for conservation laws of systems of the type (1.1) for higher vertical degrees. Essentially, it tells us that any (2, s) conservation law can be constructed from certain (0, s − 1) contact forms which satisfy an equation involving the adjoints of the linearized equations of the original system. Precisely, it says the following: Theorem 1.3. Let R be a second order hyperbolic system of type (1.1). Then for s ≥ 1 and ω ∈ Ω 2,s R ∞ a d H -closed form, there exist contact forms and the ρ ij satisfy the equation where the Ψ ij are maps from the space of (0, s − 1) forms to the space of (2, s) forms defined explicitly in Section 3.5, and the operators L * ij are the adjoints of the operators appearing in the linearized system.
The second noteworthy result of Section 3 concerns the cohomology of the constrained variational bi-complex which is defined in Section 2.3. Before stating the theorem, we will take a moment to introduce some terminology which will be helpful in the discussion of this result, and those to follow. Just as in the case of the classical Laplace method, there are generalized Laplace invariants which arise during the application of the generalized Laplace transform. The generalized Laplace invariants are relative invariants with respect to contact transformations, so their vanishing is a contact-invariant condition. When the generalized Laplace transformation is applied repeatedly to a particular system of equations, a sequence of generalized Laplace invariants is generated and this sequence may or may not terminate at some point. If, for example, the (i, j) Laplace invariants are zero after p ij applications of the X ij -Laplace transform, then p ij is referred to as a Laplace index of the system of equations, written ind(X ij ) = p ij . If the sequence of Laplace invariants never terminates, then we write ind(X ij ) = ∞.
Then the crux of Theorem 3.17 is that if each of these sequences of Laplace invariants fails to terminate, there will be no non-trivial horizontal cohomology of bi-degree (2, s) for each s ≥ 3, as stated below. Theorem 1.4. Let R be a second order hyperbolic system of type (1.1) and suppose that ind(X ij ) = ∞ for 1 ≤ i < j ≤ 3. Then, for all s ≥ 3, all type (2, s) conservation laws are trivial. That is, The proof of this result utilizes the fact that relative invariant contact forms can be constructed from nonzero solutions to the adjoint equation seen in Theorem 3.15. A contact form, ω, is said to be invariant relative to the characteristic vector field X if X(ω) = λω for some function λ on R ∞ , where X(ω) denotes the projected Lie derivative of ω with respect to X. The existence of such contact forms then contradicts the hypothesis that none of the Laplace indices is finite.
The second set of results from this paper is contained in Section 4 and addresses the topic of the Darboux integrability of systems of the form (1.1). Classically, a second order scalar hyperbolic partial differential equation (1.2) is Darboux integrable if there exist smooth, realvalued functions I,Ĩ, J, andJ such that dI ∧ dĨ = 0, dJ ∧ dJ = 0 and where X and Y are the characteristic vector fields for the equation. It is well known, see for example [21], that for any pair of monotone functions f 1 , f 2 ∈ C ∞ (R; R), the system is completely integrable in the sense of the Frobenius theorem. It is therefore evident that Darboux integrable equations may be solved via ordinary differential equation techniques.
Example 1.5. To illustrate the concepts above, we return to the Liouville equation In this example, the characteristic vector fields are X = D x and Y = D y (the total derivatives with respect to x and y, respectively). Then I = y andĨ = u yy − 1 2 u 2 y are invariant functions with respect to X, and J = x andJ = u xx − 1 2 u 2 x are invariant functions with respect to Y as required by the definition above. Letting f 1 = g 1 and f 2 = g 2 , it is possible to integrate the system (1.7) to obtain u = log 2g 1 (y)g 2 (x) (g 1 (y) + g 2 (x)) 2 .
The concept of Darboux integrability has been studied extensively, and extended to new settings. Two noteworthy examples are [4], where the Darboux integrability of equations (1.2) is shown to imply the existence of infinitely many conservation laws of type (1, s) for all s ≥ 0, and [3] where the definition of Darboux integrability is recast in a group-theoretic approach that applies to the general framework of exterior differential systems. The definition introduced in [3] equates Darboux integrability with the existence of what the authors refer to as a Darboux pair. This terminology is explained in Section 4.2, where we also prove the lemma quoted below, demonstrating that if certain characteristic invariant functions exist, systems of the form (1.1) will satisfy the definition of Darboux integrability given in [3]. Lemma 1.6. Let u ij = f ij x 1 , x 2 , x 3 , u, u i , u j , 1 ≤ i < j < 3 be a system of three hyperbolic equations in one dependent and three independent variables with characteristic vector fields X 1 , X 2 , X 3 . If there exist smooth, real-valued functions, I,Ĩ, J,J, and K,K, such that the following two conditions hold, then the system defines a Darboux pair and thus satisfies the notion of Darboux integrability defined in [3].
1. I andĨ are invariant with respect to two of the characteristic vector fields, say X i and X j , and I andĨ are functionally independent: dI ∧ dĨ = 0.
2. J,J, K, andK are all invariant with respect to X l , l = i, j, and are all functionally independent from each other: dJ ∧ dJ ∧ dK ∧ dK = 0.
The next lemma describes a method of constructing a contact form invariant with respect to a pair of characteristic vector fields by using the characteristic invariant functions described in Lemma 1.6. Lemma 1.7. Let I and J be functions on R and X 1 , X 2 and X 3 characteristic vector fields. If I and J are invariant with respect to both X 1 and X 2 , such that X 3 (I) = I and X 3 (J) = 1, then is an X 1 and X 2 invariant contact form.
It is then possible to utilize the invariant contact form given in the previous lemma in order to draw a connection between Darboux integrability and the Laplace indices of the system, as the following corollary indicates. Corollary 1.8. Let R be a system of equations of the form (1.1). If R satisfies the conditions described in Lemma 1.6, then at least one of the Laplace indices ind(X ij ) must be finite.
The connection between Darboux integrability and the termination of sequences of Laplace invariants has been established in many different contexts. In particular, [4] and [20] taken together show that a scalar second-order hyperbolic PDE in the plane is Darboux integrable if and only if the Laplace indices of each of the two generalized Laplace transforms associated to the equation are finite. In Section 4, we also describe an algorithm for generating infinitely many (1, s) and (2, s) conservation laws for a system that the conditions of Lemma 1.6. Let us summarize here the procedure for generating (1, s) conservation laws in particular.
First, a d-closed (s + 1)-form is constructed by taking the wedge product of the exterior derivatives of a sequence of characteristic invariant functions whose existence is guaranteed in the hypotheses of Lemma 1.6. Let this (s + 1) form be written as α = dI 1 ∧ dI 2 ∧ · · · ∧ dI s+1 . It is then shown that the (1, s) component of α will be d H closed. Additional characteristic invariant functions are found by applying an appropriate choice of characteristic vector field repeatedly to the invariant functions described in Lemma 1.6, and in this way infinitely many conservation laws may be constructed. The method presented for the construction of type (2, s) conservation laws is completely analogous. It is furthermore shown in Section 4 that these conservation laws, of both types, can be written in such a way that they may be shown to be nontrivial, allowing us to conclude with the following theorem. Theorem 1.9. If R is a system of equations satisfying the hypotheses of Lemma 1.6, then there exist infinitely many nontrivial type (1, s) and type (2, s) conservation laws for all s ≥ 0.
The paper concludes with Section 5 which provides a summary of our findings and suggestions for future research. The results presented in this paper are a part of the author's Ph.D. Thesis at McGill University [18].

The variational bi-complex
Presently we will establish the necessary definitions and notation pertaining to several key concepts, including infinite jet bundles, split exterior differentiation, and the variational bicomplex, as presented in [1].

Infinite jet bundles
Begin with a fibered manifold π : E → M, with adapted coordinates x i , u α , for 1 ≤ i ≤ n and 1 ≤ α ≤ q, over a connected base manifold M of dimension n. In this paper we will be concerned with local properties, although the infinite jet bundle has important global properties as well (see [1]). Given our focus, we will take M to be an open connected subset of R n and π to be the trivial bundle. Let J k (E) = ∪ x∈M J k x (E) denote the bundle of k-jets of local sections of E, with local coordinates on J k (E) consisting of x i , u α , u α i 1 , u α i 1 i 2 , . . . , u α i 1 i 2 ...i k where 1 ≤ i 1 ≤ i 2 ≤ · · · ≤ i k ≤ n and with the natural projection maps In local coordinates a p-form on J ∞ (E) will be a sum of the form where the coefficients a IJ α are C ∞ real-valued functions defined on J ∞ (E) and I and J are multiindices, I = i 1 · · · i r and J = j 1 · · · j t . The contact ideal on J ∞ (E) generated by the contact 1-forms forms an ideal in Ω * (J ∞ (E)), which we will denote by C(J ∞ (E)) and whose exterior derivatives are given by the structure equations Two particular types of vector fields on J ∞ (E), total vector fields and vertical vector fields, will play an important role in the variational bi-complex: for which X ω = 0 for any contact form ω is referred to as a total vector field.
Note that in general total vector fields take the form

The bi-graded exterior derivative
We may now introduce a bi-grading of the p-forms ω on J ∞ (E) which will allow us to distinguish between independent and dependent variables when working with differential equations. Definition 2.3. A p-form ω is said to be of type (r, s) if r + s = p and ω(X 1 , . . . , X p ) = 0 whenever there are more than r total vector fields or more than s π ∞ M vertical vector fields among the vector fields X i . Using Definition 2.3, the de Rham complex on J ∞ (E) can be bi-graded as follows where in local coordinates a differential form ω of the type (r, s) will have the form ω = a iI β dx i 1 ∧ · · · ∧ dx ir ∧ θ β 1 I 1 ∧ · · · ∧ θ βs Is for real-valued functions a iI β on J ∞ (E). The bi-grading of forms in the de Rham complex induces a corresponding decomposition of the exterior derivative d : Ω n −→ Ω n+1 , given by Let π : E → M and ρ : E → N be two fibered manifolds and let Φ : J ∞ (E) → J ∞ (E ) be a smooth map. We define the projected pullback map, which will maintain a form's bi-graded type, to be the map Φ : Ω r,s (J ∞ (E )) → Ω r,s (J ∞ (E)) given by where π r,s is the projection map from Ω p (J ∞ (E)) to Ω r,s (J ∞ (E)). Likewise, for a total vector field X on J ∞ (E), and a type (r, s) form ω, we define X(ω) ∈ Ω r,s (J ∞ (E)) to be the projected Lie derivative, X(ω) = π r,s (L X ω).
The identity X(ω) = X d H (ω) + d H (X ω) follows from Cartan's formula and this, along with (2.1) and the antisymmetry of d H and d V , implies the following relations Finally, the proposition below gives additional properties of the projected Lie derivative that will be used in performing calculations, and are direct consequences of the definition of the projected Lie derivative and the properties of the Lie derivative.
Proposition 2.4. Let ω ∈ Ω r,s (J ∞ (E)). If X and Y are total vector fields on J ∞ (E) and Z is a π ∞ M vertical vector field on J ∞ (E) then the following two equations hold Using the bi-grading of the exterior derivative, we can now define the variational bi-complex for type (r, s) forms, which will provide a natural framework for our investigation of conservation laws for PDEs systems of the type (1.1).

2.3
The variational bi-complex for PDEs systems of the type (1.1) We will consider a system of three semi-linear PDEs defined on an open connected subset U of R 3 of the following form The equations (2.4) define a locus in J 2 (E) where E is the trivial bundle Let R denote an open contractible subset of this locus. Assume that the f ij are C ∞ functions in a neighborhood of R, and that π| R : R → U is a subbundle of the fiber bundle π : J 2 (E) → M . Note that for each ij pair, ∂F ij /∂u ij = 0 on R. Furthermore, in order for the system (2.4) to be locally solvable, we make the assumption that the following set of integrability conditions is satisfied, If the module of contact forms on J 2 (E) is pulled back to R, then a Pfaffian system, I, on R is obtained. The differential system I is generated by the one forms Solutions of (2.4) are then local sections σ : U → R such that σ * (ω) = 0 for all ω ∈ I. That is, solutions are integral manifolds of I with the independence condition The k th prolongation, written R (k) , of R is the locus in J k+2 (E) defined by the equations for 1 ≤ i < j ≤ 3 and l + m + n ≤ k. For example, the first prolongation of R is where (j k s)(x) denotes equivalence classes of k-jets of local sections s of E.
Each prolongation R (k) is a C ∞ submanifold of J k+2 (E). As a consequence of involutivity R (k+1) fibers over R (k) , π k+1 k | R (k+1) : R (k+1) → R (k) . Thus we can define the inverse limit of the system of k th prolongations to be R ∞ , called the infinite prolongation of R with the projection maps π ∞ k : R ∞ → R k and π ∞ U : R ∞ → U . Denote by C(R ∞ ) the pullback of the contact ideal on J ∞ (E) to R ∞ , C(R ∞ ) = ι * [C(J ∞ (E))] where ι is the inclusion map ι : R ∞ → J ∞ (E). Then we can define the constrained variational bi-complex to be the pullback of the free variational bi-complex.
Notice that there are no additional columns to the right of the constrained bi-complex due to the fact that our system (2.4) involves exactly 3 independent variables.
We have the following coordinates on R ∞ and a basis for the contact ideal on R ∞ is given by We will call the basis dx 1 , . . the coordinate coframe on R ∞ . It will be the first of several coframes introduced in our study of the systems (2.4).
In this coordinate system the total derivatives, D i , are then expressed as and likewise for D 2 and D 3 . A smooth function g x 1 , x 2 , x 3 , u, u 1 , u 2 , u 3 , . . . , u 1 k , u 2 k , u 3 k , . . . on R ∞ has the structure equations If in addition ω is exact, meaning that there exists a function N such that dN = ω, then ω is said to be a trivial conservation law. If the M i are themselves (0, s) contact forms, for s ≥ 1, then a (1, s) form such that d H ω = 0 is referred to as a contact form valued conservation law of R.
The following lemma follows from the fact that d H is an anti-commuting differential.
Lemma 2.7. If ω is a type (r, s) form valued conservation law and there exists a type (r − 1, s) form γ for which ω = d H γ, then d H ω = 0 and ω is called a trivial conservation law.
Remark 2.8. We will not need to investigate type (3, s) conservation laws owing to the fact that the systems (2.4) which are the focus of our study have three independent variables and hence any type (3, s) form will be trivially d H -closed. Furthermore, due to Vinogradov's "two line theorem", conservation laws of type (0, s) are also trivially d H -closed [33]. Thus the focus of what follows will be type (1, s) and (2, s) form valued conservation laws, henceforth simply referred to as conversation laws.
The conservation laws described above can also be viewed in terms of the horizontal cohomology of the constrained variational bi-complex. The cohomology space is made up of cohomology classes whose representatives are type (1, s) conservation laws, and the cohomology space consists of cohomology classes whose representatives are type (2, s) conservation laws.

Conservation laws for a particular class of semi-linear PDEs
From this point on we will be considering involutive semi-linear second order systems of PDEs specifically of the form where 1 ≤ i, j ≤ 3 and i = j. 1 A basis for the space of total vector fields defined on the infinite prolongation R ∞ of the equation manifold R defined by (3.1) is given by any set of three linearly independent vector fields where the D i are the total vector fields on J ∞ (E) restricted to R ∞ , defined by (2.7). The system (3.1) is hyperbolic in the sense that its associated exterior differential system has three distinct characteristics. In particular, the characteristic equation for each F ij = 0 is just λµ = 0 which has the roots (λ, µ) = (1, 0) and (λ, µ) = (0, 1) leading us to associate the total vector fields D i and D j to it . Thus we may take as our basis for the space of total vector fields on R ∞ the characteristic vector fields D 1 , D 2 and D 3 . These vector fields have the particularly convenient property of commuting with each other, in other words: Many of the results in this paper can be expected to remain true for more general involutive systems of the form with the property that the universal linearization of (3.2), defined in Section 3.1, is of the form where θ is a contact form on R ∞ and the total vector fields X i are given by and are not assumed to commute pairwise. Also note that we are not using the Einstein summation convention anywhere in this paper.
Since the vector fields (3.3) form a basis for the space of total vector fields on R ∞ , the commutator of X i and X j can be written as a linear combination of these, Choosing to consider systems of the form (3.1), for which we may let X i = D i , and utilizing the fact that these characteristics commute so that the coefficients B k ij = 0 above, will significantly simplify the expressions for the universal linearization of (3.1), as well as the structure equations and Lie bracket congruences for the Laplace adapted coframe which are given in the following sections. Therefore we will assume henceforth that we are considering systems of the form (3.1) and have chosen the characteristic vector fields for (3.1) to be X i = D i , unless explicitly stated otherwise.

Universal linearization
The first step in analyzing systems of the form (3.1) will be to linearize each equation in the system.
On the equation manifold R ∞ , the contact forms θ, θ i , θ j and θ ij are not independent, but rather related to each other according to the equations obtained by taking d V of each of the three equations F ij = 0, According to (2.2), D j (θ I ) = θ Ij , so it is straightforward to write these equations in terms of the characteristic vector fields, X i : We will refer to equation (3.4) as the universal linearization of F ij , and denote it by L ij (F ij ), or simply L ij when the system whose linearization we are considering is clear from the context.
Eventually we will also want to have the freedom to rescale the contact form θ in order to manipulate the form of equation (3.4). To this end, define Θ = µθ for a non-vanishing function µ defined on R ∞ . Then the universal linearization (3.4) can be written equivalently as An example of a non-linear involutive system of the form (3.1), can be found in Vassiliou [31]. The linearization of this system is described in the example below.
Example 3.1. Consider the following involutive system of the form (3.1) The corresponding linearized system is

Characteristic coframe
Our next goal is to describe a set of one-forms which will form a coframe on R ∞ . Begin by defining the (1, 0) forms σ i to be dual to the characteristic vector fields X i . That is, σ i (X i ) = 1 and σ i (X j ) = 0 for i = j. Then, since we will take where the total vector fields X i act on ω by projected Lie differentiation, as defined in Section 2.3. In particular, there will be several occasions where we wish to compute the horizontal exterior derivative of a type (1, s) form. That is, if ω is given by To complete the coframe, define the remaining contact forms inductively: for i = 1, 2, 3. Before giving their structure equations and stating formally that the forms (3.12), along with the previously defined forms {σ 1 , σ 2 , σ 3 , Θ}, do indeed make up a coframe on R ∞ , we need to pause to give the following definition and lemma.
has adapted order k if it lies in the exterior algebra generated, over the smooth functions on R ∞ , by the one-forms where k is minimal.
The adapted order of a form ω is invariant under contact transformations on R ∞ , and may differ from its order as a form on R ∞ .
The following lemma will be used to establish the structure equations for the one forms comprising the coframe given below. Lemma 3.3. For k ≥ 1 and i = j, X i ξ k j restricted to R ∞ has adapted order less than or equal to k.
Proof . The proof is completed by induction on the adapted order, k. First we show the statement holds for k = 1, that is which is evidently of adapted order ≤ 1. Next assume that X i ξ l j has adapted order ≤ l for all l ≤ k. Then X i ξ l j can be written as a linear combination with a, b l , c l , d l are all functions on R ∞ . As mentioned previously in this section, we are free to choose the characteristics for systems of the form (3.1) to commute. Then which has adapted order ≤ k + 1 as needed.
We are now prepared to state and prove the following proposition.
Proposition 3.4. The set of 1-forms forms a coframe on the equation manifold R ∞ , called the characteristic coframe. The d H structure equations for this coframe are given by where µ k i and ν k i are contact forms of adapted order ≤ k.
Proof . Repeatedly applying the vector fields X i to the contact form Θ = µθ, and using equation (2.2), we see that where δ consists of contact forms of order < k. Then the set of one forms where at least two of the indices m, n, l are ≥ 1, clearly spans the contact ideal C(R ∞ ). By Lemma 3.3, each X m 1 X n 2 X l 3 (Θ) can be expressed as a linear combination of the forms in (3.13) of order ≤ m + n + l. The structure equations (3.14)-(3.16) are a result of the definition of the forms ξ j i given in (3.12), formula (3.10), and Lemma 3.3.

Generalized Laplace method for systems
We will now introduce the generalized Laplace method as a means of solving systems of the form (3.5). The geometric origins of this method can be found in Chern's work [12] on the Laplace transformation of submanifolds admitting conjugate nets of curves. It was subsequently used in [22] to generalize the classical Laplace method (also described in [22]) to involutive overdetermined systems of linear equations in n independent and 1 dependent variable of the form where a k lk , a l lk , and c lk are differentiable functions of the independent variables x 1 , . . . , x n and each set of functions is symmetric in its lower indices. The similarity between the form of the system (3.17) and that of the linearized system (3.5) is apparent, and provides motivation for the description of the generalized Laplace method suited to studying the systems (3.5) given below.
A linearized system of the form (3.5) must satisfy certain integrability conditions stemming from the fact that D k (θ ij ) = D j (θ ik ) for i, j, k distinct. Applying D k to each equation L ij (Θ) = 0 and equating the coefficients on like-ordered contact forms leads to the following relations which the coefficients of any compatible system must satisfy for 1 ≤ i = j = k ≤ 3, These correspond to the compatibility conditions defined in [22], the only difference being that the partial derivatives found in the expressions written in [22] have been replaced with total derivatives in (3.18).

Laplace invariants and the generalized Laplace transform
The classical Laplace method associates to a given linear hyperbolic PDE in the plane, F (u, u x , u y , u xx , u xy , u yy ) = 0, two Laplace invariants h(F ) and k(F ). The vanishing of either invariant allows for the integration of the equation by quadratures, and the failure of these invariants to vanish allows one to apply the classical Laplace method in order to obtain a transformed equation whose Laplace invariants may or may not vanish. We will now develop an analogous procedure to investigate the systems (3.5). Begin by defining the following n(n − 1) 2 expressions which are invariant under rescaling of Θ by a nonvanishing function on R ∞ , and will be referred to as the higher-dimensional Laplace invariants of the system Note that when there are only two independent variables, and hence the system (3.5) is replaced by a single equation, the invariants (3.19) reduce to the classical Laplace invariants. In this case, the invariants (3.20) would not be defined.
Example 3.5. Recall the nonlinear example of an involutive system of the form (3.1), given in equation (3.9) and restated below, The Laplace invariants for this system are as follows , We are now in a position to introduce the concept of the generalized Laplace transform. Suppose that Θ solves the system (3.5). Then for any ordered pair (i, j), define the (i, j) Laplace transform of Θ to be (3.21) and denote this by ξ ij = X ij,L lk (Θ), or simply X ij (Θ) if the system of differential operators, L lk , 1 ≤ l = k ≤ 3, being considered is clear from context. One should be aware that the order of the pair (i, j) is significant, as we can see by comparing the (i, j) Laplace transform of Θ, given by (3.21), with the (j, i) Laplace transform of Θ, which would be The following proposition shows that the transformed contact form ξ ij = X ij (Θ) will solve a system of the same form as (3.5).
Proposition 3.6. Let ξ ij be the (i, j)-Laplace transform of Θ, where Θ satisfies the system of equations (3.5). In the case that the (i, j) Laplace invariants H ij and H ijk are nonzero for all k, k = i, j, ξ ij will satisfy a system of equations of the same form as (3.5).
Proof . We can compute expressions for the total derivatives of ξ ij directly from the definition of the (i, j) Laplace transform and the system (3.5): and then applying X j to both sides of (3.25), Finally, expressions for Θ k and Θ jj are gotten by solving for Θ k and Θ jj in (3.24) and (3.23) respectively, and making the appropriate substitutions using (3.25) and (3.26): Now we may take the total derivatives of (3.22)-(3.24) to see that indeed ξ ij satisfies a system of the form The coefficients in the equation (3.27) can be given explicitly in terms of the coefficients of the original system as follows, And for k = i, j we havê This proves the proposition.
According to the preceding proposition then, when a given system has non-vanishing (i, j) Laplace invariants, the system of three total differential operators L lk , 1 ≤ l = k ≤ 3, is transformed under the (i, j) Laplace transform into another system of three total differential operators of the same form. Denote the operator that L lk is transformed into under the (i, j) Laplace transform by X ij (L lk ). Then the (i, j) Laplace transform of Θ, X ij (Θ) = ξ ij , will satisfy the system [X ij (L lk )](ξ ij ) = 0 for all 1 ≤ l = k ≤ 3. When we wish to emphasize that the (i, j) Laplace transform is being performed, we will write (3.27) as Just as in the case of the classical Laplace method, the generalized Laplace transform has an inverse when the (i, j) invariant is nonzero, as the next result shows.
Proposition 3.7. Let ξ ij = X ij (Θ) be the (i, j) Laplace transform of Θ, which satisfies the system (3.5). If H ij = 0 then an inverse Laplace transform of the (i, j) Laplace transform exists and is given by Proof . By the definition of the Laplace transform, X ji (ξ ij ) = X i (ξ ij ) +Â j ij ξ ij . Substituting (3.21) and (3.22) into this expression and using the fact thatÂ j ij = A j ij , which follows from (3.18), we see that X ji (ξ ij ) = H ij Θ. Since H ij = 0, this can be solved for Θ = 1 H ij X ji (ξ ij ), as we wished to show.

Laplace adapted co-frame for systems of nonlinear PDEs
We can now introduce another coframe on the equation manifold R ∞ which is constructed by utilizing the generalized Laplace transforms, X ij , defined in Section 3.3. This coframe will be denoted by Before describing the elements of this coframe explicitly, we will pause to establish some necessary terminology and notation. As described in Proposition 3.6, if the (i, j) Laplace invariants of a given system (3.5) are nonzero, then the system may be transformed into another system, written X ij (L lk )(ξ ij ) = 0, of the same form. The (i, j) Laplace invariants of that transformed system may likewise be computed. Denote these by H 1 ij = H ij (X ij (L lk )) and H 1 ijq = H ijq (X ij (L lk )). The process may be repeated so long as the (i, j) Laplace invariants of the last system do not vanish. Assuming all the necessary Laplace invariants are nonzero, let the system of differential operators obtained by applying the (i, j) Laplace transform n times to the original system (3.5) be denoted by X n ij (L lk ). Accordingly, we will write to refer to the (i, j) Laplace invariants of the system of differential operators obtained by applying the (i, j) Laplace transform to the original system (3.5) n times. With this notation then, we will write H 0 ij and H 0 ijq to represent the Laplace invariants of (3.5) itself. Using equations (3.28)-(3.33) and the notation established in (3.34), we can write out explicitly the coefficients of the operator X n ij (L lk ) defining the n th application of the (i, j) Laplace transform to the system of equations L lk = 0. Write and the coefficients X n ij A l lk , X n ij A k lk , and X n ij (C lk ) can be computed in terms of the coefficients of the original system (3.5) as follows. Take the coefficient X n ij A l lk and proceed by induction, first expressing it in terms of the coefficient on X i in the operator X n−1 ij (L lk ): Likewise, it is straightforward to see that The following definition establishes some further notation concerning the generalized Laplace invariants.
Definition 3.8. Let ind(X ij ) = p ij denote the number of times the (i, j) Laplace transform must be applied to the system (3.5) in order to obtain vanishing (i, j) Laplace invariants. That is, H(X p ij ij (L lk )) = 0. If the (i, j) Laplace invariants never vanish despite repeated applications of the (i, j) Laplace transform, then we will write ind(X ij ) = p ij = ∞.
With this notation in mind, we may now introduce the Laplace adapted coframe as the set of one-forms where, for j = 1, 2, 3, And for i = j, The d H structure equations for the Laplace adapted coframe (3.35) can be computed directly using equation (3.10) and, given their complexity, are relegated to Appendix B.
Define the vertical vector fields dual to the contact one forms Θ andξ n j by where δ j k and δ n l are both Kronecker delta functions. So, for example,ξ n 1 V n 1 =1 andξ n 2 V n 3 =0. Then we can express the d H structure equations above in terms of the Lie brackets of the characteristic vector fields X i and the vertical vector fields, U and V l k . These congruences will be needed to prove some important results in the coming sections and have been provided in Appendix B for the reader's reference.

Adjoint of the linearized operator L ij
We will now introduce the adjoint to the linearized operator L ij , which will play a pivotal role in subsequent sections. Let us begin with the following definition. Definition 3.9. For a total differential operator F, its formal adjoint operator, denoted by F * , is the total differential operator such that for every ρ ∈ Ω 0,s and ω ∈ Ω 0,r there exists γ ∈ Ω 2,s+r such that Using this definition, we will state the following proposition which describes the adjoint L * ij of the linearized operator L ij . The proof of this result can be found in [32].
Proposition 3.10. The differential operator L ij which defines the universal linearization (3.5) has the following adjoint operator

Characteristic invariant contact forms
Characteristic invariant contact forms play an important role in the construction of conservation laws, as we will see in Sections 3.5 and 4.3. First we will define invariant functions and contact forms, and then analogously, relative invariant contact forms. Then an important result, Proposition 3.14, regarding the existence of relative invariant contact forms will be stated and proved.
Definition 3.11. For a total vector field X on the equation manifold R ∞ , a function f defined on R ∞ is said to be an X invariant function if X(f ) = 0.
Likewise for contact forms, we have Definition 3.12. A type (0, s) contact form ω is said to be invariant with respect to the total vector field X, or equivalently called an X invariant contact form, if X(ω) = X d H ω = 0. Furthermore, if ω is invariant with respect to two distinct total vector fields, X and Y , then we will say that ω is an X and Y invariant contact form.
If a slightly weaker condition is satisfied, then we have a contact form which is a relative X invariant, as described in the following definition.
Definition 3.13. For X a total vector field, the type (0, s) contact form ω is said to be a relative X invariant contact form if X(ω) = λω for some function λ on R ∞ . If ω is invariant relative to two distinct total vector fields, X and Y , then we say that ω is a relative X and Y invariant contact form.
Given this notation, we may state and prove the following proposition.
Finally, suppose ind(X 13 ) = ind(X 23 ) = ∞. Then clearly if ω is a relative X 1 invariant form then the preceding argument shows ω must equal 0. The proofs of (3.42) and (3.43) follow by analogous arguments, permuting the roles of X 1 , X 2 , and X 3 as needed.
Then for s ≥ 1, we can define a map as follows where the ψ ij k are defined by We may now state the following theorem whose proof will employ an integration by parts type of argument to express any d H -closed (2, s) form in terms of solutions to the adjoint equation  The proof of Theorem 3.15 will make use of the following lemma, which concerns the forms on J ∞ (E) that lie in the kernel of the inclusion map ι : R ∞ → J ∞ (E), where R ∞ is the infinitely prolonged equation manifold defined by the system of equations F ij = 0. Lemma 3.16. If ω ∈ Ω p (J ∞ (E)) satisfies ι * ω = 0 then ω can be expressed as where α ij lk , α lkm ∈ Ω p (J ∞ (E)) and β ij lk , β lkm ∈ Ω p−1 (J ∞ (E)).
Proof . Utilizing the system of equations (3.1), there is a set of coordinates on J ∞ (E) and a corresponding basis of one forms on J ∞ consisting of (3.52) Using this coframe, any p-form ω on J ∞ (E) can be written as where ω 0 is a p-form generated by the forms in (3.51) and the β are (p − 1)-forms on J ∞ (E). Since all the terms involving d V expressions will pull back to zero on the equation manifold R ∞ , ι * ω = ι * ω 0 . The forms in (3.51) are still independent when pulled back to R ∞ , so ι * ω = 0 if and only if all the coefficients vanish when ι * ω 0 is written in terms of the basis elements (3.51). We can conclude that each of these coefficients is a C ∞ (J ∞ (E)) linear combination of the functions Thus w 0 can be written as The above expression for ω 0 , along with equation (3.53), yields the desired form in Lemma 3.16.
In what follows, the notation D ij is used as a shorthand to mean D i D j . Defining ρ ij = 1 s ι * (ρ ij ), we will next show that the forms ρ ij satisfy the following relationship involving the adjoint equations L * ij in the coordinate coframe on R ∞ and that ω can be written as . Next the interior Euler-Lagrange operator J : Ω 3,s (J ∞ (E)) → Ω 3,s−1 (J ∞ (E)), which is defined by will be applied to both sides of (3.54). Note that in equation (3.58) D ij simply refers to the operator D i D j . By [2, Theorem 2.6], we know that for any (2, s) formω on J ∞ (E), J(d Hω ) = 0. When the operator J is applied to the right hand side of (3.54), we use the fact that where Q consists of terms depending linearly on F ij , d V F ij , etc., to obtain Since, for example, the expression restricted to R ∞ is the adjoint of L 12 in the coordinate frame on R ∞ , we see that J(d Hω ) = 0 implies that the ρ ij satisfy the equation (3.55).
To write the expression for ω in the coordinate coframe as well, we make use of the homotopy operator h r,s H : Ω r,s (J ∞ (E)) → Ω r−1,s (J ∞ (E)) as defined in [2]: We are concerned with the case where n = r = 3, so we have the equation Ifω is a (3, s) form of the type dx 1 ∧ dx 2 ∧ dx 3 ∧ M with M ∈ Ω 0,s (R ∞ ), then Therefore the pullback of (3.60) to R ∞ gives us the expression (3.56) for ω in the coordinate coframe on R ∞ . Now to write (3.57) and (3.56) in terms of the Laplace adapted coframe, we let L 0 ij be the operator L ij defined previously but with µ = 1. Then L 0 ij = L ij and L 0 * ij (ρ ij ) = L * ij (ρ ij ) so we have L 0 * ij (ρ ij ) = 0. Given that A i ij = F u i , A j ij = F u j and σ i = dx i , we can look at (3.57) and see directly that the expression for Ψ ij (ρ ij ) in the coordinate coframe corresponds to the expression for Ψ ij (ρ ij ) given by equations (3.47)-(3.49) with µ = 1.
We can now make use of the preceding result to prove the following important result concerning the (2, s) cohomology of the variational bi-complex for s ≥ 3.
Theorem 3.17. Let R be a second order hyperbolic system of type (3.1) and suppose that ind(X ij ) = ∞ for 1 ≤ i < j ≤ 3. Then, for s ≥ 3, all type (2, s) conservation laws are trivial. That is, Proof . According to Theorem 3.15, we only need to show that there do not exist nonzero type (0, s − 1) solutions ρ ij to the adjoint equations L * ij (ρ ij ) = 0 as this would preclude the existence of any nontrivial conservation law. Begin by rewriting the adjoint equation as a system of first order equations We will proceed by showing that if there is a nonzero solution to the system (3.61)-(3.62), a contradiction to Proposition 3.14 results. To that end, let ρ ij be a nonzero solution to (3.61)-(3.62) of adapted order k. Because we have taken s ≥ 3, ρ ij is a contact form of degree greater than or equal to 2 and the adapted order of ρ ij is k ≥ 1. Thus V k l ρ ij = 0 for some l = 1, 2, 3, where we recall that V k l is the vertical vector, defined by (3.40), dual to the Laplace adapted coframe. For the sake of clarity, and without loss of generality, we will take ρ ij to be ρ 12 and V k l to be V k 3 . We will begin by demonstrating that V k 3 ρ 12 is an X 1 invariant (0, s−1) contact form. Apply formula (2.3) to see that where the last equality holds since ρ 12 has adapted order k and hence V k+1 3 ρ 12 = 0. Now, using the Lie bracket congruences in Proposition B.2, we see that Since interior product of the right hand side of equation (3.62) with V k+1 3 is zero, this implies that Finally, take the interior product of (3.61) with V k 3 to obtain and conclude that V k 3 ρ 12 is a relative X 1 invariant contact form. Since for systems of the form (3.1) we may choose the characteristic vector fields to commute, the preceding argument can be replicated with the roles of X 1 and X 2 reversed to show that V k 3 ρ 12 is a relative X 2 invariant contact form as well. However, this contradicts Proposition 3.14, which states that if all Laplace indices are infinite, no nonzero relative X 1 and X 2 invariant contact forms exist. Hence there cannot in fact exist any nonzero (0, s − 1) solutions ρ 12 to the adjoint equation L * 12 = 0, and likewise regarding the equations L * 13 = 0 and L * 23 = 0, and the proof of Theorem 3.17 is complete.
We will conclude this section by applying Theorem 3.15 to an example found in [22].
Example 3.18. Consider the following involutive system of the form (3.1), The universal linearization for this system is given by L ij (θ) = θ ij + θ i + θ j + θ = 0 (3.64) and the system of adjoint operators is then Taking the triple (ρ 12 , ρ 23 , ρ 13 ) = e x 1 +x 2 , e x 2 +x 3 , e x 1 +x 3 , which solves the adjoint equation 1≤i<j≤3 L * ij (ρ ij ) = 0, we can construct the following conservation law according to the structure theorem presented in Theorem 3.15: 64). So ω is in fact a (2, 1) conservation law as we wished to show.
Example 3.19. We refer the reader to [23] for the following additional example of an involutive linear system of the form (3.1). This system is integrable by repeated applications of the Laplace transformation

Systems of Darboux integrable equations
We will begin by stating the definition of Darboux integrability for a single PDE in two independent variables as described, for example, in [21]. Given a second order scalar hyperbolic equation, R, in two independent variables x and y, F (x, y, u, u x , u y , u xx , u xy , u yy ) = 0. (4.1) Let X and Y be the characteristic total vector fields and I the Pfaffian system associated to R.
It can be shown that the integral manifolds of this system correspond to the general solution of R.
The relationship between the property of Darboux integrability for equations (4.1) and the vanishing of generalized Laplace invariants was investigated in [4] and [20]. We will outline one of the main findings of these works here.
Let X = m x D x + m y D y and Y = n x D x + n y D y , with δ = m x n y − m y n x = 0, be distinct characteristic vector fields associated to equation (4.1) such that we have the factorization (m x λ − m y µ)(n x λ − n y µ) = κ F uxx λ 2 − F uxy λµ + F uyy µ 2 .
In investigating equations of the form (4.1), no assumption is made in [4] or [20] regarding the existence of commuting characteristic vector fields, so we will need to consider the commutator of X and Y , which we write as [X, Y ] = P X + QY . Let L(θ) denote the universal linearization of (4.1) obtained from the identity d V F = 0: where the coefficients in (4.2) are given by Note that the operator in (4.2) is essentially the form-valued version of what is considered in [22] where the classical Laplace method is discussed. Likewise, in this setting two Laplace invariants are defined. The first is Also similar to the classical case, if for example H 0 is nonzero, (4.2) may be transformed into another equation of the same form. On the equation manifold R ∞ , L(θ) = 0 holds identically and defining η 1 = Y (θ) + Aθ, this identity can be expressed as When H 0 = 0, we see that η 1 satisfies an equation of the same form as (4.2) by rewriting (4.4) as where Equation (4.5) is called the Y-Laplace transform of (4.2). The Laplace invariant H 1 = X(A 1 ) + A 1 B 1 − C 1 of (4.5) can then be computed and if H 1 is also nonzero, the equation can be transformed again. Clearly this process can continue with In this way a second sequence of Laplace invariants, K i = Y (E i ) + E i D i − G i is generated, again noting that K i+1 will be defined so long as K i = 0. Let p = ind(Y) be the Laplace index of the Y-Laplace transform and q = ind(X ) be the Laplace index of the X -Laplace transform. That is, p is the number of times that the Y Laplace transform must be applied before the generalized Laplace invariant of the transformed equation equals zero and likewise for q. If one of the sequences of Laplace invariants never vanishes, we write p = ∞ or q = ∞. which was shown to be Darboux integrable in Section 1. Here the characteristic vector fields X = D x and Y = D y commute, and the linearization of (4.6) is Computing the Laplace invariants for this equation we see that H 0 = e u , so we may then write the Y-Laplace transform of (4.7) as XY (η 1 ) − u y X(η 1 ) − e u η 1 = 0. From this equation we can compute H 1 = 0, so for the Liouville equation p = 1. Likewise, we can compute K 0 = e u and K 1 = 0 to conclude that q = 1. Both Laplace indices are finite as expected from Theorem 3.17.

Characteristic systems and invariant functions
In order to extend the ideas developed in the preceeding section to the case of three equations in three independent variables of the form we begin by defining characteristic Pfaffian systems which correspond to the characteristics defined by Cartan's structural classification as is detailed in Appendix A.
Definition 4.4. Define the characteristic Pfaffian system of order k, for a given characteristic vector field, X i , as where θ,ξ 1 1 ,ξ 1 2 ,ξ 1 3 , . . . ,ξ k 1 ,ξ k 2 ,ξ k 3 is either the Laplace-adapted or characteristic coframe on R ∞ , the infinite prolongation of the equation manifold defined by (4.8). We can likewise define a characteristic Pfaffian system of order k with respect to a pair of characteristic vector fields, X i and X j , as Since dI = X 1 (I)σ 1 + X 2 (I)σ 2 + X 3 (I)σ 3 + d V I, it is immediately apparent that I is an X i invariant function of order k if and only if dI ∈ C k (X i ), and that I is invariant with respect to both X i and X j if and only if dI ∈ C k (X i , X j ).
Associated to any Pfaffian system I, there is a flag of Pfaffian subsystems called the derived flag, and so forth. The contact formξ 1 i can not be eliminated however since dσ 2 and/or dσ 3 may contain the term σ i ∧ξ 2 i . Likewise, for (4.10), we note that d Hξ k (X i , X j ) and the rest of the argument continues in the same fashion.
The next result follows immediately from Lemma 4.5.
Corollary 4.6. If J k is an X i -invariant function of order k, then for l, j = i and where a k b k = 0. If I k is invariant with respect to X i and X j , then for l = i, j and a k = 0.

Darboux integrability for a system of three equations
We will first discuss the definition of Darboux integrability given by Anderson et al. in [3] where the authors investigate a class of differential systems for which superposition formulas can be constructed. For a system of the form we have the associated Pfaffian system I generated by ω 0 , ω 1 , ω 2 , ω 3 with structure equations dω 0 ≡ 0, (4.13) dω 1 ≡ σ 1 ∧ π 1 mod I, (4.14) and characteristic vector fields X i = D i . To apply the definition of Darboux integrability established in [3], there must exist a coframe, {θ 0 , θ 1 , θ 2 , θ 3 , σ 1 , π 1 , σ 2 , π 2 , σ 3 , π 3 } such that I is generated algebraically by the 1-forms and 2-forms whereΩ i ∈ Ω 2 (σ 1 , π 1 , σ 2 , π 2 ), for i = 1, 2, andΩ 1 ∈ Ω 2 (σ 3 , π 3 ). One then defines two Pfaffian systemsV = {θ 0 , θ 1 , θ 2 , θ 3 , σ 1 , π 1 , σ 2 , π 2 } andV = {θ 0 , θ 1 , θ 2 , θ 3 , σ 3 , π 3 }, which will be referred to as the singular differential systems for I with respect to the decomposition (4.17). For simplicity of notation, we will continue using this particular choice of decomposition in all that follows. Note that we could just as well have grouped dω 1 and dω 3 , or dω 2 and dω 3 , together to make an analogous decomposition. Let V (∞) denote the largest integrable sub-bundle of a given Pfaffian system, V, and recall that the rank of V (∞) will equal the number of functionally independent first integrals of V. Then according to [3], a differential system I is Darboux integrable if the associated singular systemsV andV are Pfaffian and define a Darboux pair, meaning that the following three properties are satisfied: where M , thought of as an open subset of R 10 , is the manifold on which the exterior differential system I is defined.
We can now show that systems of equations of the form (4.12) which possess certain characteristic invariant functions, described in the lemma below, will define a Darboux pair and thus be Darboux integrable in the sense of [3].
Lemma 4.7. Let u ij = f ij x 1 , x 2 , x 3 , u, u i , u j , 1 ≤ i < j < 3 be a system of three hyperbolic equations in one dependent and three independent variables with characteristic vector fields X 1 , X 2 , X 3 . If there exist smooth, real-valued functions, I,Ĩ, J,J, and K,K, such that the following two conditions hold, then the system defines a Darboux pair and thus satisfies the notion of Darboux integrability defined in [3].
1. I andĨ are invariant with respect to two of the characteristic vector fields, say X i and X j , and I andĨ are functionally independent: dI ∧ dĨ = 0.
2. J,J, K, andK are all invariant with respect to X l , l = i, j, and are all functionally independent from each other: dJ ∧ dJ ∧ dK ∧ dK = 0.
From structure equations (4.13)-(4.16) it is clear that the third property defining a Darboux pair is satisfied. Property 2 asserts thatV andV have no integrals in common, and if necessary this condition can be satisfied by restrictingV andV to a level set of their common integrals. Thus the essential property to check when constructing a Darboux pair is that We will now show that this condition is indeed satisfied by the systems described in Lemma 4.7.
As in the preceding sections, for systems of the form (4.12), we take the characteristic vector fields to be X i = D i which commute. However one would expect to be able to reproduce the analysis carried out in this paper for more general systems of the form which are still contained in the same class of equations described in Cartan's structural classification (see Theorem A.3) as the systems (4.12), but for which the characteristics X i do not commute. We will now describe how for systems with non-commuting characteristic vector fields, the vector fields X i may be rescaled to make them commute pairwise, given the existence of sufficiently many characteristic invariant functions.
Proposition 4.9. Let I and K be functions on R ∞ such that I is invariant with respect to both X i and X j , and K is invariant with respect to X l . Assume that X i (K), X l (I) and X j (K) are all nonzero and define the following rescaled characteristic vector fields .
Then the vector fieldX l commutes with bothX i andX j : Proof . Write the commutator of X i and X j as follows Taking into consideration the conditions of invariance imposed on I and K, we may write Now use (4.20) and (4.21) in the following computation to confirm (4.19) An analogous calculation using B j jl X j (K) + B i jl X i (K) = −X l X j (K) and B l jl X l (I) = X j X l (I) shows that X j ,X l = 0.
The formulation of the Darboux-adapted coframe, as well as the construction of conservation laws for Darboux integrable systems, will rely on our ability to obtain characteristic invariant contact forms from the invariant functions whose existence is guaranteed in Lemma 4.7. For example, we can write down a contact form which is invariant with respect to the characteristic vector fields X 1 and X 2 as follows, Lemma 4.10. Let I, J, and K be functions on R ∞ that are invariant with respect to the characteristic vector field X 3 such that X 1 (I) = 1, X 2 (I) = 0, X 1 (J) = 0, and X 2 (J) = 1. Also define K = X 1 (K) and K = X 2 (K). Then Proof . Note that K and K are also X 3 invariant since the characteristic vector fields X i all commute. Computing d H ω, we obtain From this calculation we see that X 3 (ω) = X 3 d H (ω) = 0 so that ω is X 3 -invariant.
Likewise,we can use X 1 and X 2 invariant functions to obtain X 1 and X 2 invariant contact forms in the following way.
Lemma 4.11. Let I and J be functions on R ∞ which are invariant with respect to both characteristic vector fields X 1 and X 2 , such that X 3 (I) = I and X 3 (J) = 1. Then is an X 1 and X 2 invariant contact form.
Proof . In the same spirit as the last proof, we compute d H ω and by grouping terms arrive at Given the preceding lemma, the following corollary follows immediately from Proposition 3.14.
Corollary 4.12. Let R be a system of equations of the form (4.12). If R satisfies the conditions of Lemma 4.7, then at least one of the Laplace indices p ij = ind(X ij ) must be finite.
The construction of the Darboux adapted coframe will proceed according to the orders of the various characteristic invariant functions in Lemma 4.7. Here we do not have a complete analysis of the possible orders of invariants, as was given by Goursat [19] in the case of a single hyperbolic second order equation in the plane. Thus we will carry out the construction under certain assumptions on the orders of the invariants to serve as a demonstration of how a coframe may be constructed in a particular case.
Given two first order X 1 and X 2 invariant functions, I 1 andĨ 1 , we may obtain a sequence of functions invariant with respect to X 1 and X 2 by repeatedly applying X 3 : where X 3 Ĩ 1 = 1.
The set of contact forms {α 1 , α 2 , . . .} can thus replace the branch of the Laplace coframe given by ξ 1 3 ,ξ 2 3 , . . . previously. The structure equations for the α i are given by where (4.23) is calculated in the following way To complete the coframe, we need to construct contact forms which are equivalent to ξ 1 1 ,ξ 2 1 , . . . and ξ 1 2 ,ξ 2 2 , . . . . We accomplish this by first generating X 3 invariant functions, repeatedly applying the total vector fields X 1 and X 2 to the X 3 invariant functions J and K from Lemma 4.7 to obtain the following two sequences of X 3 invariant functions: Then using Lemma 4.10, we can write two corresponding sequences of X 3 invariant contact forms: · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · Then a similar argument to that given above for the sequence of α i shows that and likewise, where at least one of b k and c k is nonzero, and at least one of e k and f k is nonzero. While Lemma 4.7 requires that dJ ∧ dJ ∧ dK ∧ dK = 0, this condition is not enough to guarantee that d V J i and d V K i are functionally independent. As a result, there is some ambiguity about the construction of the remaining two branches of the coframe, as it will depend on whether the coefficients onξ i 1 andξ i 2 in β i and γ i satisfy b k f k − c k e k = 0. In other words on whether β i and γ i are independent mod θ,ξ 1 3 ,ξ 2 3 , . . . ,ξ i−1

3
. If any pair β i and γ i does fail to be independent, then we will complete the coframe by taking the necessary contact forms from the Laplace adapted coframe.
This coframe would likely be useful in the pursuit of an eventual classification of Darbouxintegrable systems, similar to that presented in [20] for the case of a scalar second order hyperbolic PDE in the plane.

Generating infinitely many (1, s) and (2, s) conservation laws
We may assume that for a Darboux integrable system, [X i , X l ] = P i il X i + P l il X l = 0 and [X j , X l ] = P j jl X j + P l jl X l = 0. Recalling the structure equation we see that if X l commutes with X i and X j , then d H σ l = 0. With this in mind, we may now demonstrate a way to construct infinitely many (1, s) conservation laws for Darboux integrable systems.
Theorem 4.13. Let R be a Darboux integrable system of equations. Then there exist infinitely many non-trivial type (1, s) conservation laws for all s ≥ 0.
Proof . Fix distinct i, j, and l as in Lemma 4.7. For s = 0, we may construct a conservation law by taking ω = Iσ l where I is an X i and X j invariant function. Suppose that ω is exact, that is ω = Iσ l = d HĨ for some functionĨ. Then soĨ is X i and X j invariant and X l (Ĩ) = I. Thus ω will be a nontrivial conservation law so long as I = X l (Ĩ) forĨ another X i and X j invariant function. Additional (1, 0) conservation laws may be constructed by generating additional X i and X j invariant functions by applying X l to I repeatedly. Note that X l (I) is still X i and X j invariant since X l commutes with both X i and X j . For s ≥ 1, we first note that if ω is a d closed (s+1) form with vanishing (2, s−1) component, then the (1, s) component of ω will be d H closed. To see this, write dω as a sum of forms of different bidegrees, according to the direct sum decomposition of the space of (s + 2) forms, (4.24) Let dω = ν 0 + ν 1 + ν 2 + ν 3 where ν i ∈ Ω (i,s+2−i) . Since each ν i belongs to a different component of (4.24), dω = 0 implies that ν i = 0 for all i. In particular, ν 2 = 0. The form ω can itself be decomposed according to bidegree as well. Write ω = µ 0 + µ 1 + µ 2 + µ 3 where µ i ∈ Ω (i,s+1−i) . Now decompose dω into its horizontal and vertical components, dω = (d H + d V )(ω), to see that d H µ 1 + d V µ 2 = ν 2 = 0. If the (2, s − 1) component, µ 2 , of ω vanishes, then this last equation becomes d H µ 1 = 0 and we can conclude that the (1, s) component of ω is d H closed. Now, to construct a d H closed (1, s) form, we begin by defining α = dI 1 ∧ dI 2 ∧ · · · ∧ dI s+1 where each of the functions I n is an X i and X j invariant function. These functions can be obtained as before by taking I = I 1 and repeatedly applying the total vector field X l to I to generate the remaining functions. The (s + 1) form α will have vanishing (2, s − 1) component: since the functions I n are all X i and X j invariant, the only horizontal form that will appear in any of the dI n is σ l . Thus, by the preceding discussion, π 1,s (dI 1 ∧ dI 2 ∧ · · · ∧ dI s+1 ) will be d H closed as required.
Furthermore, we can construct a d H closed (1, s) form that is not exact. In other words, we can write down a nontrivial (1, s) conservation law. Let p l = min{ind(X il ), ind(X jl )}, and take k ≥ p l + s. Then for a nonzero form η ∈ Ω s−2 I (α p l +1 , α p l +2 , . . . , α k−1 ) where I is the ring of functions which are both X i and X j invariant, the (1, s) form is d H closed, but not d H exact. It is straightforward to see that ω is indeed d H closed since σ l is itself d H closed and α k+1 ∧ α k ∧ η is both X i and X j invariant. Now suppose that ω were d H exact, that is, that there exists a (0, s) form, γ such that d H γ = ω. Since d H γ = σ 1 ∧ X 1 (γ) + σ 2 ∧ X 2 (γ) + σ 3 ∧ X 3 (γ), comparing this with the above expression, (4.25), for ω we conclude that γ must also be X i and X j invariant. Moreover, X l (γ) = α k+1 ∧ α k ∧ η which implies that γ ∈ Ω s I (α p l +1 , α p l +2 , . . .). The adapted order of ω being k + 1 implies that the adapted order of γ would be k. Then γ can be written as γ = α k ∧ β + δ where β and δ each have adapted order ≤ k−1, and are both X i and X j invariant forms. But then X l (γ) = X l (α k ∧β +δ) = α k+1 ∧α k ∧η which indicates that β = α k ∧η, showing that β has adapted order k. This contradicts the original premise that β be of adapted order ≤ k − 1, from which we conclude that ω is not exact.
It is possible to construct infinitely many nontrivial (2, s) conservation laws in a similar way. By an argument analogous to that given above in the case of (1, s) conservation laws, if ω is a d closed (s + 2) form with vanishing (3, s − 1) component, then the (2, s) component of ω is d H closed. Consider the form dJ 1 ∧ dJ 2 ∧ · · · ∧ dJ s+2 where the functions J n = X i (J n−1 ) are all X l invariant. This (s + 2) form has vanishing (3, s − 1) component since each term dJ n can only contain the horizontal forms σ i and σ j . So is a d H closed (2, s) form. Construct a second d H closed (2, s) form by taking the (2, s) component of the (s + 2) form dK 1 ∧ dK 2 ∧ · · · ∧ dK s+2 where K n = X j (K n−1 ) and each K n is X l invariant. Let µ = π 2,s (dK 1 ∧ dK 2 ∧ · · · ∧ dK s+2 ).
Take I to be the ring of X l invariant functions. Then the following (2, s) form will be a nontrivial conservation law where η 1 ∈ Ω s−2 I (ν p i +1 , ν p i +2 , . . ., ν k−1 ) and η 2 ∈ Ω s−2 I (µ p j +1 , µ p j +2 , . . ., µ k−1 ), and k ≥ max{p j + s, p i + s}. Clearly d H ω = 0 since σ i ∧ σ j is d H closed and ν k+1 ∧ ν k ∧ η 1 and µ k+1 ∧ µ k ∧ η 2 are X l invariant. To see that ω is not exact, we will assume otherwise and produce a contradiction. Suppose that there exists a (1, s) form, γ, such that d H γ = ω. First note that by writing d H γ = σ 1 ∧ X 1 (γ) + σ 2 ∧ X 2 (γ) + σ 3 ∧ X 3 (γ) and looking at equation (4.26), it may be concluded that γ is X l invariant, l = i, j. Furthermore, any (1, s) form can be written as where each M i is a (0, s) form. Comparing equation (3.11) for d H γ with the expression (4. 26) for ω, we see that d H γ = ω implies that and so M i , M j ∈ Ω s I (ν p l +1 , ν p l +2 , . . .). Since the adapted order of ω is (k + 1), the adapted order of γ must be k and hence M i and M j may be written as where α 1 , α 2 , δ 1 , and δ 2 all have adapted order ≤ k − 1. Plugging these expressions for M i and M j into (4.27) yields (4.28) From equation (4.28), we deduce that α 1 = µ k ∧ η 2 and α 2 = ν k ∧ η 1 . Since this contradicts the fact that α 1 and α 2 were taken to have order ≤ k − 1, we can conclude that ω is indeed not exact and thus represents a nontrivial conservation law.

Conclusion and further research
In this paper we have proven a number of results concerning the existence and structure of conservation laws for involutive systems of partial differential equations of the form After defining the constrained variational bicomplex, which provides the framework in which we investigate form-valued conservation laws, and developing other technical tools such as the generalized Laplace transform, and the structure equations and Lie bracket congruences for the Laplace adapted coframe, the first main result is proved in Section 3.5. This result, Theorem 3.15, describes the structure of (2, s) conservation laws for s ≥ 1 in terms of solutions to the adjoint equations of the associated linearized system. As a consequence of this structure theorem, we obtain another key result which states that if all of the Laplace indices of the linearized system are infinite, then all (2, s) conservation laws for s ≥ 3 are trivial. That is, in this case the horizontal cohomology of the constrained variational bi-complex associated to the system is trivial: Conditions are given in Lemma 4.7 which will guarantee that a system of the form (1.1) will satisfy the definition of Darboux integrability introduced in [3]. It is then show that these same conditions will imply that at least one of the Laplace indices for the linearization of such a system must be finite. A procedure for using characteristic invariant contact forms to construct non-trivial (1, s) and (2, s) conservation laws for s ≥ 0 is then described in Section 4.3.
There are several natural directions in which one may try to extend the results presented here. First, we should emphasize that in this paper we have considered the very specific class of semi-linear systems of the form (5.1). In particular, the characteristic vector fields, X i , for such systems commute. It should be fairly easy to investigate which of the results presented here could be replicated for more general classes of nonlinear systems belonging to the different cases of Cartan's structural classification listed in Theorem A.3, specifically systems with non-commuting characteristics. Other directions in which to extend this work would be to consider overdetermined systems in one dependent and n independent variables, instead of only 3 independent variables, or to allow u to be a vector-valued function (see, for example, [35]). Furthermore, it remains to find further examples and applications which fall under the form of systems (5.1) which we consider here. Constructing involutive overdetermined systems is a project in its own right. Another interesting avenue of research would consist in classifying Darboux-integrable systems of the form (5.1) or a generalization of these. A potentially fruitful approach to this problem would utilize the definitions and classification scheme presented in [3], which is discussed in Section 4.
Beyond these points, there is the overarching question of how the horizontal cohomology of the constrained variational bi-complex H r,s can be interpreted for s ≥ 1. In other words, what does the horizontal cohomology mean in terms of the analytic properties of the solution space of the system R? Here we refer the reader to [2] for some examples of the ways in which the cohomology of the variational bi-complex has been used to study, for example, the inverse problem of the calculus of variations for autonomous systems of ODE, Riemannian structures, and Gelfand-Fuks cohomology. In [28], Tsujishita also describes several applications of the higher-degree cohomology of the variational bi-complex, including its role in determining deformation classes of the solution space of R. An interpretation of H n−1,2 (R ∞ ) cohomology, where n is the number of independent variables in the associated PDEs system, is also given in [36]. It is shown there that the existence of nontrivial cohomology classes [ω] ∈ H n−1,2 indicates the possible existence of variational principles for R. Specifically, if R is the Euler-Lagrange equations for some Lagrangian λ, then it can be shown that d V λ = d H η on R ∞ for η a form of bi-degree (n − 1, 1). Therefore ω = d V (η), which Zuckerman calls the universal conserved current for λ, is a d H closed (n − 1, 2) form. Assuming that λ is nontrivial, Vinogradov's two-line theorem [33] implies that ω is not d H exact, so that if H n−1,2 (R ∞ ) = 0, then R does not admit a variational principle.

A Structural classification
In [10], Cartan studies involutive Pfaffian systems defined on a suitable jet space, which arise when considering involutive systems of PDEs from the perspective of exterior differential systems. In the course of his analysis, the Pfaffian systems' structure equations are computed and simplified. Systems with the same reduced structure equations are then regarded as being structurally equivalent. Some of the main results of this work are summarized in [25], and it is this exposition which we will follow as we proceed to place systems of the form (1.1) in the context of Cartan's structural classification.
Let ∆ a (x i , u, u i , u ij ) = 0 with i, j = 1, 2, 3 be an involutive second order system of PDEs considered as a subset of the bundle of 2-jets, J 2 R 3 , R . As defined in Section 2, consider the contact system on J 2 R 3 , R generated by the 1-forms u ij dx j , i = 1, 2, 3.
Let ω i = dx i , π i j = du ij = π j i so that a local coframe for J 2 R 3 , R is given by ω i , θ 0 , θ i , π i j for i, j = 1, 2, 3. Then we have the following structure equations The set of forms θ 0 , θ i comprises a basis for the contact system on J 2 R 3 , R . Keeping in mind that our goal is to obtain a simplified set of structure equations, another basis may be found by transforming the forms θ 0 , θ i in the following way where the entries a i j are functions on the jet space J 2 R 3 , R and det A = 0. Observe that this transformation preserves the contact form θ 0 , as well as the Pfaffian system θ 1 , θ 2 , θ 3 .
Denote by A the group consisting of such matrices 2 , A. Then we will show presently that there is a prolongation A of A that acts on the coframe ω i , θ 0 , θ i , π i j , i = 1, 2, 3, such that the form of the structure equations for the corresponding contact system is preserved.
Let A = (a i j ) and A −1 = (α i j ) so thatθ i = Aθ i and θ i = A −1θi . A straightforward calculation using the definitions provided shows that dθ 0 = Moreover, the π i j are transformed as the coefficients of a quadratic form.
The group A will later be utilized to simplify the possible structure equations which must be considered in the structural classification.
We will now consider more specifically systems of three second order equations, ∆ a (x i , u, u i , u ij ) = 0, 1 ≤ a ≤ 3, which satisfy the maximal rank condition that the 3 × 13 Jacobian matrix J ∆ (x, u, u i , u ij ) = ∂∆ a ∂x i , ∂∆ a ∂u I of ∆ a with respect to x i , u, u i , and u ij is of rank 3 whenever ∆ a = 0. To further ensure that the PDEs system is non-degenerate, we require the (3 × 6) submatrix ∂∆ a ∂u ij to have rank 3. Let R be the codimension 3 submanifold of J 2 R 3 , R thus defined by ∆ a = 0, and C 2 (R) the corresponding contact system restricted to R. When the one forms ω i , θ 0 , θ i , π i j are restricted to R, we obtain 3 linear relations between them. By replacing π i j with π i j plus an appropriate linear combination of the θ i , including θ 0 , these linear relations may be simplified to the following We will henceforth assume that the Pfaffian system C 2 (R) is an involution on which in other words, there exists at least one Kähler-regular integral element of the form where the α i jl are symmetric with respect to all indices. This allows the relations (A.2) to be further reduced to We would like to know how many of the coefficients a i jl are left arbitrary by these 9 equations since Cartan's involutivity test tells us that C 2 (R) is involutive with ω 1 ∧ ω 2 ∧ ω 3 = 0 if and only if the number of arbitrary coefficients a i jl above equals 3 · 3 − 2σ 1 − σ 2 where σ 1 and σ 2 are the first two reduced Cartan characters of C 2 (R). Given that we are concerned with the case of systems of 3 equations in 3 independent variables, we can see that the Cartan characters σ 1 and σ 2 satisfy σ 1 ≤ 3 and σ 1 + σ 2 ≤ 3. Thus the number of arbitrary coefficients a i jl will be 9 − 2σ 1 − σ 2 ≥ 9 − 6 = 3. Since the a i jl are symmetric with respect to all indices, there are 10 terms, which we now see are related by at most 10 − 3 = 7 independent linear equations. Hence, at most 7 of the 9 equations (A.4) are linearly independent. Making the identifications π i j ←→ξ iξj and a i jl ←→ξ iξjξl , we can state the following result: Lemma A.2. Defining the quadratic forms A jk iξ iξj , k = 1, 2, 3, the 9 cubic formsξ i M k are related by at least 9 − 7 = 2 linear equations.
The lemma states that the 9 cubic formsξ i M k are related by at least two linear relations, thus there are 6 linear forms, l i and m i , i = 1, 2, 3, in the variablesξ i , such that these two linear relations can be expressed by the two equations, Given that these relations among the M i exist, one may take a linear combination, and find z 1 and z 2 such that one of the coefficients z 1 l i + z 2 m i vanishes. In other words, we can choose z 1 and z 2 such that the three linear forms z 1 l i + z 2 m i are linearly dependent. Equivalently, consider z 1 and z 2 as homogeneous coordinates on the Riemann sphere, P 1 (C), and define z := z 2 /z 1 . Then we have three linear forms l i + zm i , each of which is a linear combination of the three variablesξ i . Taking the coefficient of each of the three variablesξ i in each of the three expressions l i + zm i yields a 3 × 3 matrix. Setting the determinant of this matrix equal to zero will determine z in such a way that the three linear forms l i + zm i are linearly dependent, and produce a cubic equation in z with three roots in P 1 (C).
To make this more explicit, and simplified, the group A may be used to induce linear transformations of the variablesξ i which then induce corresponding linear transformations of z 1 and z 2 . This allows us to assume, for example, that one root of the equation mentioned above is z = 0. Then l 1 , l 2 , and l 3 are linearly dependent and can be taken to be l 1 =ξ 1 , l 2 =ξ 2 , l 3 = 0.
Furthermore, if we specify the m i by m 1 = a 1ξ1 + b 1ξ2 + c 1ξ3 , m 2 = a 2ξ1 + b 2ξ2 + c 2ξ3 , m 3 = a 3ξ1 + b 3ξ2 + c 3ξ3 , then the equation which z must satisfy is the cubic polynomial given by, The structural classification then breaks into 5 possible cases, according to whether the equation (A.5) has (i) 3 simple roots, (ii) one double and one simple root, (iii) one triple root, or (iv) vanishes identically. Case (iv) splits further into two sub-cases according to which of the following two forms the matrix in equation (A.5) can be reduced, Hence there are 5 possible contact invariant structures in total, as described in the following theorem.
Finally, we will demonstrate to which of the five classes given in Theorem A.3 systems of the type (1.1) belong. Recall that the system of three equations (1.1) is given by Then by Lemma A.2 there exist at least 2 linear relations between the 3 quadratic forms Write these equations as l 1ξ1ξ2 + l 2ξ2ξ3 + l 3ξ1ξ3 = 0, m 1ξ1ξ2 + m 2ξ2ξ3 + m 3ξ1ξ3 = 0.

B Structure equations and Lie bracket congruences
Proposition B.1. The d H structure equations for the Laplace adapted coframe (3.35) for a system of the form (1.1) of three hyperbolic equations with commuting characteristics are given by the following equations. Recall that p 1 = min{p 21 , p 31 }, p 2 = min{p 12 , p 32 }, and p 3 = min{p 13 , p 23 }. Furthermore, p k1 , p k2 , and p k3 will be used to denote which Laplace index achieves the minimum in the definitions of p 1 , p 2 , and p 3 respectively; whereas p l1 , p l2 , and p l3 will denote the Laplace index which does not achieve the minimum in the definitions of p 1 , p 2 , and p 3 . So, for example, if p 1 = p 31 then p k1 = p 31 and p l1 = p 21 , n