# Multiple Riemann wave solutions of the general form of quasilinear
hyperbolic systems

A.M. Grundland  Centre de Recherches Mathématiques, Université de Montréal,
Succ. Centre-Ville, CP 6128, Montréal (QC) H3C 3J7, Canada
Department of Mathematics and Computer Science, Université du Québec, CP 500,
Trois-Rivières (QC) G9A 5H7, Canada<EMAIL_ADDRESS>J. de Lucas
Department of Mathematical Methods in Physics, University of Warsaw,
ul. Pasteura 5, 02-093, Warszawa, Poland<EMAIL_ADDRESS>
###### Abstract.

The objective of this paper is to construct geometrically Riemann $k$-wave
solutions of the general form of first-order quasilinear hyperbolic systems of
partial differential equations. To this end, we adapt and combine elements of
two approaches to the construction of Riemann $k$-waves, namely the symmetry
reduction method and the generalized method of characteristics. We formulate a
geometrical setting for the general form of the $k$-wave problem and discuss
in detail the conditions for the existence of $k$-wave solutions. An auxiliary
result concerning the Frobenius theorem is established. We use it to obtain
formulae describing the $k$-wave solutions in closed form. Our theoretical
considerations are illustrated by examples of hydrodynamic type systems
including the Brownian motion equation.

###### Key words and phrases:

$k$-wave solution, Riemann invariant, conditional symmetry, Lie point
symmetry.

###### 1991 Mathematics Subject Classification:

35Q53 (primary); 35A30, 35Q58, 53A05 (secondary)

##  1\. Introduction

Riemann waves represent a very important class of solutions of first-order
nonlinear systems of partial differential equations (PDEs). They are
ubiquitous in mathematical physics since they are admitted in all
multidimensional hyperbolic systems and constitute their elementary solutions.
Their characteristic feature is that, in most cases, they are expressible in
implicit form.

In this section we consider the possibility of generating the idea of Riemann
wave solutions of a hyperbolic system of partial differential equations to
what may be considered a nonlinear superposition of Riemann waves. This fact
will be explained in detail in Sections 5 and 6. All functions are assumed to
be smooth. We now discuss a homogeneous hyperbolic quasilinear first-order
system of $q$ equations in the dependent (properly determined) variables
$u=(u^{1},\ldots,u^{q})$ and independent variables $x=(x^{1},\ldots,x^{p})$ of
the form

$\sum_{i=1}^{p}A^{i}(u)u_{i}=0,\qquad u_{j}=\frac{\partial u}{\partial
x^{j}},\qquad j=1,\ldots,p,$ (1.1)

where $A^{1},\ldots,A^{p}$ are $q\times q$ matrices whose coefficients are
functions depending on $u$. All considerations here are local, so it suffices
to search for solutions defined on a neighbourhood of $x=0$. A Riemann wave
solution for (1.1) is defined locally on a neighbourhood of $x=0$ as a
solution of the equation

$u=f(r(x,u)),$ (1.2)

for certain functions $f:\mathbb{R}\rightarrow\mathbb{R}^{q}$ and
$r:\mathbb{R}^{p+q}\rightarrow\mathbb{R}$ given by

$r(x,u)=\sum_{i=1}^{p}\lambda_{i}(u)x^{i},$

where the function $r$ is called a Riemann invariant (RI) associated with the
wave vector $\lambda=(\lambda_{1},\ldots,\lambda_{p})$. Namely, it is assumed
that

$\ker\left(\sum_{i=1}^{p}\lambda_{i}(u)A^{i}(u)\right)\neq 0.$

Equation (1.2) defines a unique function $u(x)$ on a neighbourhood of $x=0$
for any function $f:\mathbb{R}\rightarrow\mathbb{R}^{q}$ and

$\frac{\partial u^{\alpha}}{\partial
x^{i}}(x)=\phi(x)^{-1}\lambda_{i}(u(x))\frac{df^{\alpha}}{dr}(r(x,u(x))),\quad
i=1,\ldots,p,\quad\alpha=1,\ldots,q,$

where

$\phi(x)=1-\sum_{\alpha=1}^{q}\frac{\partial r}{\partial
u^{\alpha}}(x,u(x))\frac{df^{\alpha}}{dr}(r(x,u(x))).$

The solutions (1.2) have rank at most equal to one. They are the building
blocks for constructing more general types of solutions that describe
superpositions of many simple Riemann waves, which are much more interesting
from the physical point of view. Until now, the only way to approach this task
was through the method of characteristics or via group theory. The method of
characteristics relies on treating RIs as new dependent variables (which
remain constant along appropriate characteristic curves of the basic system),
which leads to the reduction of the dimension of the problem. The most
important theoretical results obtained with the use of the method of
characteristics (see e.g. [5, 16, 17, 18, 21, 22] and [25] for a review)
include the finding of some necessary and sufficient conditions for the
existence of so-called Riemann $k$-waves in multidimensional systems. It was
shown in [3, 20, 21], through the use of involutive differential equations in
the Cartan sense, that these solutions of (1.1) depend on $k$ arbitrary
functions of one variable.

Certain criteria for determining the elastic or nonelastic character of the
superposition between Riemann waves, particularly useful in physical
applications, were also found [17, 21, 22]. In applications to fluid dynamics
and field theory, many new and interesting results were obtained [2, 6, 8, 10,
22, 24]. The method of characteristics, like all other techniques for solving
PDEs, has its limitations. These limitations motivated the search for the
means of constructing larger classes of multiple wave solutions expressible in
terms of Riemann invariants, i.e.

$u=f\left(r^{1}(x,u),\ldots,r^{k}(x,u)\right),$ (1.3)

where $f:\mathbb{R}^{k}\rightarrow\mathbb{R}$ and

$r^{a}(x,u)=\sum_{i=1}^{p}\lambda^{(a)}_{i}(u)x^{i},\qquad 1\leq a\leq k.$

The natural way to do this was to look at these solutions from the point of
view of the group invariance properties of the system (1.1). The feasibility
and clear advantages of such an approach were demonstrated for certain fluid
dynamics equations [6, 13].

The objective of this paper is to develop a systematic way of constructing
multiple Riemann waves by making use of both approaches, namely the method of
characteristics and the group analysis of differential equations. The basic
feature of Riemann $k$-waves is that they remain constant on $(p-k)$\-
dimensional hyperplanes, where the set of linearly independent wave covectors
$\lambda^{(1)},\ldots,\lambda^{(k)}$ entering the superposition vanish. In the
context of group theory, this means that the graph of a solution,
$\Gamma=\\{(x,u(x)):x\in\mathbb{R}^{p}\\}$, is invariant under all vector
fields on $\mathbb{R}^{p+q}$ of the form

$X_{a}=\sum_{i=1}^{p}\eta^{i}_{a}(u)\partial_{x^{i}},\qquad a=1,\ldots,k,$
(1.4)

that satisfy the orthogonality conditions

$\sum_{i=1}^{p}\eta^{i}_{a}(u)\lambda^{(a)}_{i}(u)=0,\qquad\forall
a=1,\ldots,k.$

Then, $u(x)$ is a solution of (1.1) for some function $f$ because the
functions $\\{r^{1},\ldots,r^{k},u^{1},\ldots,u^{q}\\}$ constitute a complete
set of invariants of the Abelian algebra $L$ of the vector fields (1.4). The
implicit form of these solutions (1.3) leads to some limitation in the
application of the classical symmetry reduction method to this case. To
overcome these limitations, we rectify the set of vector fields
$X_{1},\ldots,X_{k}$ by a change of independent and dependent variables in
$\mathbb{R}^{p+q}$, choosing RIs as the new independent variables. The initial
system (1.1) expressed in these new coordinates, complemented by the
invariance conditions for the rectified vector fields $X_{1},\ldots,X_{k}$,
form an overdetermined quasilinear system. Thus, the solutions of this system
are invariant under the Abelian group corresponding to $L$. These vector
fields $X_{1},\ldots,X_{k}$ are conditional symmetries of the initial system
(1.1).

In this paper we explore several different aspects of the $k$-wave problem. We
investigate certain general properties of the general form of first-order
quasilinear hyperbolic systems and propose several new approaches to the
construction of its solutions. We do not obtain all solutions of this system
of PDES, although the analysis of the obtained results can provide some
interesting physical information concerning nonlinear superpositions of
Riemann waves.

The paper is organised as follows. Section 2 describes the basic theory of
Riemann invariants. In Section 3, the relation between inhomogeneous and the
corresponding homogeneous first-order systems of PDEs in their general form is
studied. Section 4 contains simple examples to illustrate the construction
introduced in Section 3. Section 5 contains a geometric formulation of the
Riemann $k$-wave problem, which represents a generalisation of the approach
introduced in Section 2. In Section 6 a new technique for generating $k$-wave
solutions is proposed and conditions for their existence are found. Obtaining
these results required certain modifications of the Frobenius theorem, which
are presented in the Appendix. In Section 7 we analyse the existence
conditions for $k$-wave solutions and provide their geometric interpretation.
Section 8 contains two examples of the application of our approach to
hydrodynamic-type systems including the construction of Riemann double waves
for the Brownian motion equation.

##  2\. Riemann invariants

This section reviews the Riemann invariants method for hyperbolic quasi-linear
systems of first-order differential equations, also called hydrodynamic
equations. Our approach provides some new elements of the geometric
description of the method, which illuminate the proposed procedure for
constructing multi-wave solutions.

Let us discuss a properly determined homogeneous quasi-linear hyperbolic
system of first-order partial differential equations (PDEs) of the form

$\sum_{i=1}^{p}\sum_{\beta=1}^{q}A^{i\alpha}_{\beta}(u)u^{\beta}_{i}=0,\qquad\alpha=1,\ldots,q,$
(2.1)

in $p$ independent variables, $x^{1},\ldots,x^{p}$, and $q$ dependent
variables $u^{1},\ldots,u^{q}$. Each variable $u^{\beta}_{i}$ stands for the
derivative of $u^{\beta}$ with respect to $x^{i}$. All following
considerations are local. Due to the invariance of (2.1) under translations of
the independent variables, the study of (2.1) can be reduced to the analysis
of what happens close to the origin $x=0$. If we define
$u_{i}=(u^{1}_{i},\ldots,u_{i}^{q})$, then (2.1) can be rewritten as (1.1),
where each $A^{i}$, for $i=1,\ldots,p$, is a $q\times q$ matrix whose entries
are functions that depend only on the dependent variables. Frequently, one
writes $x^{1},\ldots,x^{n},t=x^{n+1}$ with $n+1=p$. Then, (2.1) can be
rewritten as

$\sum_{i=1}^{n}A^{i}(u)u_{i}+A^{n+1}(u)u_{t}=0.$ (2.2)

If $A^{n+1}={\rm Id}_{q}$, where ${\rm Id}_{q}$ is a $q\times q$ identity
matrix, then (2.2) takes the so-called evolutionary form

$u_{t}=-\sum_{i=1}^{n}A^{i}(u)u_{i}.$ (2.3)

A wave covector is a non-zero vector function
$\lambda:\mathbb{R}^{q}\rightarrow\mathbb{R}^{p}\,\\!{}^{*}$ such that

$\lambda(u)=(\lambda_{1}(u),\ldots,\lambda_{p}(u)),\qquad\ker\left(\sum_{i=1}^{p}\lambda_{i}(u)A^{i}(u)\right)\neq
0.$

The function

$r:(x,u)\in\mathbb{R}^{p}\times\mathbb{R}^{q}\mapsto\sum_{i=1}^{p}\lambda_{i}(u)x^{i}\in\mathbb{R},$
(2.4)

where $x=(x^{1},\ldots,x^{p})\in\mathbb{R}^{p}$ and
$u=(u^{1},\ldots,u^{q})\in\mathbb{R}^{q}$, is called the Riemann invariant
associated with the covector $\lambda$.

Consider a function $f:\mathbb{R}\rightarrow\mathbb{R}^{q}.$ Let us first
analyse when the condition $u=f(r(x,u))$ determines a solution $u=u(x)$ of
(2.1). The latter holds if and only if
$g:(x,u)\in\mathbb{R}^{p}\times\mathbb{R}^{q}\mapsto
f(r(x,u))-u\in\mathbb{R}^{q}$ satisfies $g(x,u(x))=0$. The mapping $g$ is
regular at every point with $x=0$ since $g(0,u)=f(0)-u$ and the tangent
mapping at 0, namely $T_{0}g$, in the basis
$\\{\partial_{x_{1}},\ldots,\partial_{x_{p}},\partial_{u_{1}},\ldots,\partial_{u_{q}}\\}$,
is given by

$[T_{(x,u)}g]=\left[\frac{df^{\alpha}}{dr}\lambda_{i}(u)\bigg{|}D^{\alpha}_{\beta}\right],\qquad
D_{\beta}^{\alpha}=\sum_{i=1}^{p}\frac{df^{\alpha}}{dr}\frac{\partial\lambda_{i}}{\partial
u^{\beta}}(u)x^{i}-\delta^{\alpha}_{\beta}$

and

$[T_{(0,u)}g]=\left[\frac{df^{\alpha}}{dr}\lambda_{i}(u)\bigg{|}-\delta^{\alpha}_{\beta}\right].$

The differentiability of $g$ ensures that $D^{\alpha}_{\beta}$ is regular in
an open subset containing the points $(0,u)$ where $u\in\mathbb{R}^{q}$. By
the implicit function theorem, $u=f(r(x,u))$ determines $u=u(x)$ in an open
subset of $\mathbb{R}^{p}$ containing $x=0$. In particular, $u(0)=f(0)$.
Moreover,

$\frac{\partial u^{\alpha}}{\partial
x^{i}}(x)=\frac{df^{\alpha}}{dr}(r(x,u(x)))\left(\sum_{\beta=1}^{q}\sum_{j=1}^{p}x^{j}\frac{\partial
u^{\beta}}{\partial x^{i}}(x)\frac{\partial\lambda_{j}}{\partial
u^{\beta}}(u(x))+\lambda_{i}(u(x))\right)$ (2.5)

amounts to

$\sum_{\beta=1}^{q}\\!\left[\delta^{\alpha}_{\beta}\\!-\\!\sum_{j=1}^{p}\frac{df^{\alpha}}{dr}(r(x,u(x))\frac{\partial\lambda_{j}}{\partial
u_{\beta}}(u(x))x^{j}\right]\\!\frac{\partial u^{\beta}}{\partial
x^{i}}(x)\\!=\\!\frac{df^{\alpha}}{dr}(r(x,u(x)))\lambda_{i}(u(x)).$

The $q\times q$ matrix $\Psi$, whose entries depend on dependent and
independent variables, is defined by

$\Psi^{\alpha}_{\beta}(x,u)=\delta^{\alpha}_{\beta}-\sum_{j=1}^{p}\frac{df^{\alpha}}{dr}(r(x,u))\frac{\partial\lambda_{j}}{\partial
u^{\beta}}(u)x^{j}\Rightarrow\Psi^{\alpha}_{\beta}(0,u)=\delta^{\alpha}_{\beta},\,\,\alpha,\beta=1,\ldots,q.$

Thus, $\Psi(x,u)$ is an invertible matrix in an open subset of each point
$(0,u)$. Then,

$\frac{\partial u^{\alpha}}{\partial
x^{i}}(x)=\sum_{\beta=1}^{q}(\Psi^{-1})^{\alpha}_{\beta}(x,u(x))\frac{df^{\beta}}{dr}(r(x,u(x)))\lambda_{i}(u(x))$
(2.6)

and
$\sum_{\beta=1}^{q}\Psi^{\alpha}_{\beta}(x,u)\frac{df^{\beta}}{dr}(r(x,u))=I$
can be rewritten as

$I=\frac{df^{\alpha}}{dr}(r(x,u))-\sum_{\beta=1}^{q}\sum_{j=1}^{p}\frac{df^{\alpha}}{dr}(r(x,u))\frac{\partial\lambda_{j}}{\partial
u^{\beta}}(u)x^{j}\frac{df^{\beta}}{dr}(r(x,u))\\\
=\frac{df^{\alpha}}{dr}(r(x,u))-\sum_{\beta=1}^{q}\frac{df^{\alpha}}{dr}(r(x,u))\frac{\partial
r}{\partial u^{\beta}}(u)\frac{df^{\beta}}{dr}(r(x,u))\\\
=\left(1-\sum_{\beta=1}^{q}\frac{\partial r}{\partial
u^{\beta}}(x,u)\frac{df^{\beta}}{dr}(r(x,u))\right)\frac{df^{\alpha}}{dr}(r(x,u)).$

The latter yields

$\sum_{\beta=1}^{q}\Psi^{\alpha}_{\beta}(x,u)\frac{df^{\beta}}{dr}(r(x,u))=\left(1-\sum_{\beta=1}^{q}\frac{\partial
r}{\partial
u^{\beta}}(x,u)\frac{df^{\beta}}{dr}(r(x,u))\right)\frac{df^{\alpha}}{dr}(r(x,u)).$

Let us define

$\Phi(x,u)=1-\sum_{\beta=1}^{q}\frac{\partial r}{\partial
u^{\beta}}(x,u)\frac{df^{\beta}}{dr}(r(x,u)).$

Hence,

$\sum_{\alpha,\beta=1}^{q}(\Psi^{-1})^{\pi}_{\alpha}(x,u)\Psi^{\alpha}_{\beta}(x,u)\frac{df^{\beta}}{dr}(r(x,u))=\sum_{\alpha=1}^{q}(\Psi^{-1})^{\pi}_{\alpha}(x,u)\Phi(x,u)\frac{df^{\alpha}}{dr}(r(x,u))$

and

$\frac{df^{\pi}}{dr}(r(x,u))=\sum_{\alpha=1}^{q}(\Psi^{-1})^{\pi}_{\alpha}(x,u)\Phi(x,u)\frac{df^{\alpha}}{dr}(r(x,u)).$

Using the previous equality in (2.6), we obtain that

$\frac{\partial u^{\beta}}{\partial
x^{i}}(x)=\sum_{\alpha=1}^{q}(\Psi^{-1})^{\beta}_{\alpha}(x,u)\frac{df^{\alpha}}{dr}(x,u)\lambda_{i}(u(x))=\frac{1}{\Phi(x,u)}\frac{df^{\beta}}{dr}(x,u)\lambda_{i}(u(x)).$
(2.7)

Let us now consider a function
$\eta:u\in\mathbb{R}^{q}\mapsto(\eta^{1}(u),\ldots,\eta^{p}(u))\in\mathbb{R}^{p}$
such that

$\langle\lambda(u),\eta(u)\rangle=\sum_{i=1}^{p}\eta^{i}(u)\lambda_{i}(u)=0.$

A solution $u(x)$ of (2.7) and the function $\eta$ give rise to a vector field
on $\mathbb{R}^{p}$ of the form

$X_{\eta,u(x)}(x)=\sum_{i=1}^{p}\eta^{i}(u(x))\frac{\partial}{\partial
x^{i}}.$

In view of (2.7), we obtain

$(X_{\eta,u(x)}u^{\alpha})(x)=\sum_{i=1}^{p}\eta^{i}(u(x))\frac{\partial
u^{\alpha}}{\partial x^{i}}(x)\\\
=\sum_{i=1}^{p}\eta^{i}(u(x))\Phi^{-1}(x,u(x))\frac{df^{\beta}}{dr}(r(x,u(x)))\lambda_{i}(u(x))=0.$

Hence, the coordinates $u^{\alpha}(x)$ of any particular solution $u(x)$ of
the hydrodynamic equations (1.1) are constants of the motion of
$X_{\eta,u(x)}$ and therefore $u(x)$ is constant along its integral curves
(2.1).

##  3\. A general inhomogeneous first-order system

Let us study a more general case of (properly determined) hydrodynamic
equations (2.1) when an additional affine term appears explicitly and the
matrix functions $A^{i}$ depend on the dependent variables
$u^{1},\ldots,u^{q}$, and also on the independent variables
$x^{1},\ldots,x^{p}$, namely

$\sum_{i=1}^{p}\sum_{\beta=1}^{q}A^{i\alpha}_{\beta}(x,u)u^{\beta}_{i}=b^{\alpha}(x,u),\qquad\alpha=1,\ldots,q,$
(3.1)

for certain functions
$b^{1},\ldots,b^{q}:\mathbb{R}^{p+q}\rightarrow\mathbb{R}$, where
$x=(x^{1},\ldots,x^{p})\in\mathbb{R}^{p}$ and
$u=(u^{1},\ldots,u^{q})\in\mathbb{R}^{q}$.

Let us show that an appropriate change of variables turns the inhomogeneous
system (3.1) with $b\neq 0$ into a homogeneous one. Assume, without loss of
generality, that $b^{1}\neq 0$ on an open subset $U\subset\mathbb{R}^{p+q}$.
Then, there always exists an invertible $q\times q$ matrix $M$ whose
coefficients are functions on $U$ such that
$Mb=(1,0,\ldots,0)^{T}\in\mathbb{R}^{q}$. For instance, let us assume

$M(x,u)=\frac{1}{b^{1}(x,u)}\left[\begin{array}[]{ccccc}1&0&0&\ldots&0\\\
-b^{2}(x,u)&b^{1}(x,u)&0&\ldots&0\\\ -b^{3}(x,u)&0&b^{1}(x,u)&\ldots&0\\\
\ldots&\ldots&\ldots&\ldots&0\\\ -b^{n}(x,u)&0&0&\ldots&b^{1}(x,u)\\\
\end{array}\right].$ (3.2)

Then,

$\sum_{i=1}^{p}\sum_{\beta=1}^{q}\mathcal{A}^{i\alpha}_{\beta}(x,u)u^{\beta}_{i}=\delta_{1}^{\alpha},\qquad\alpha=1,\ldots,q,$
(3.3)

where we set $\mathcal{A}^{i}(x,u)=M(x,u)A^{i}(x,u)$ and $\delta_{1}^{\alpha}$
is the Kronecker delta $\delta^{\alpha}_{\beta}$ with $\beta=1$. Define also

$\tilde{u}=u-x^{p+1}(1,0,\ldots,0)^{T}\in\mathbb{R}^{q}.$ (3.4)

Then, a particular solution $u(x^{1},\ldots,x^{p})$ of (3.1) induces a new
function

$\tilde{u}(x^{1},\ldots,x^{p+1})=u(x^{1},\ldots,x^{p})-x^{p+1}(1,0,\ldots,0)^{T}\in\mathbb{R}^{q}$

and therefore

$\frac{\partial\tilde{u}}{\partial
x^{i}}(x^{1},\ldots,x^{p+1})\\!=\\!\frac{\partial u}{\partial
x^{i}}(x^{1},\ldots,x^{p}),\,\,\frac{\partial\tilde{u}}{\partial
x^{p+1}}(x^{1},\ldots,x^{p+1})\\!=\\!(-1,0,\ldots,0)^{T},$

for $i=1,\ldots,p$. Let us also define the $q\times q$ matrices, with
$i=1,\ldots,p$ and coefficients given by functions on $\mathbb{R}^{p+q}$, of
the form

$\tilde{A}^{i}(x,\tilde{u})=\mathcal{A}^{i}(x,\tilde{u}+x^{p+1}(1,0,\ldots,0)^{T}),\,\,i=1,\ldots,p,\,\,\forall\tilde{u}\in\mathbb{R}^{q},\,\,\forall
x^{p+1}\in\mathbb{R},$

and let $\tilde{A}^{p+1}$ be the $q\times q$ identity matrix. Then, under the
change of variables (3.4), the first-order inhomogeneous system (3.3) becomes
a quasilinear homogeneous first-order system of PDEs in $p+1$ independent
variables

$\sum_{i=1}^{p+1}\sum_{\beta=1}^{q}\tilde{A}^{i\alpha}_{\beta}(x,\tilde{u})\tilde{u}^{\beta}_{i}=0,\qquad\alpha=1,\ldots,q.$
(3.5)

Note that the above method allows us to recover solutions of the inhomogeneous
system (3.1) as a certain subset of the solutions of a homogeneous equation
(3.5). However, a solution of (3.5) does not need to give rise to a solution
of the inhomogeneous system (3.1) since the expression

$\widetilde{u}(x^{1},\ldots,x^{p+1})+x^{p+1}(1,0,\ldots,0)^{T}$

need not be a function depending only on $x^{1},\ldots,x^{p}$.

The advantage of our previous approach is that homogeneous systems of first-
order differential equations are often easier to solve than the corresponding
inhomogeneous ones. As a consequence of previous considerations, we will
discuss, in what follows, properly determined homogeneous quasi-linear first-
order systems of PDEs. It is frequently convenient to denote the independent
variables by $x^{1},x^{2},\ldots,x^{n},t=x^{p}$, where $n+1=p$, so that our
hydrodynamic system, after reduction to the homogeneous form (3.5), can be
written in the evolutionary form

$\tilde{u}_{t}+\sum_{i=1}^{n}\tilde{A}^{i}(x,\tilde{u})\tilde{u}_{i}=0.$

In what follows, we suppress the tilde above the $u$’s and consider
homogeneous systems in the dependent variables without tilde, namely (2.3),
and, instead of referring to (3.5), we will use

$\sum_{i=1}^{p}\sum_{\beta=1}^{q}A^{i\alpha}_{\beta}(x,{u}){u}^{\beta}_{i}=0,\qquad\alpha=1,\ldots,q,$
(3.6)

where we assume a generic number $p$ of independent variables, which will be
useful in what follows.

Consider an underdetermined quasilinear system (3.1). We can slightly
generalise and rewrite such a system as follows

$\sum_{\mu=1}^{p}\sum_{i=1}^{q_{h}}A_{i}^{s\mu}(x,u)\lambda_{\mu}\gamma_{1}^{i}+\sum_{\mu=1}^{p}\sum_{i=q_{h}+1}^{q}A^{s\mu}_{i}(x,u)\lambda_{\mu}\gamma^{i}_{2}=b^{s}(x,u),\quad
s=1,\ldots,m,$

where $1\leq q_{h}\leq q$ is a certain value used to divide the quasilinear
system (3.1) into two parts. The vector $\gamma=(\gamma_{1},\gamma_{2})$ has
$q$ components, while

$\gamma_{1}=(\gamma^{q_{1}})_{q_{1}=1,\ldots,q_{h}},\qquad\gamma_{2}=(\gamma^{q_{2}})_{q_{2}=q_{h}+1,\ldots,q}$

and we define

$(A_{1}\lambda,A_{2}\lambda)=\left(\left(\sum_{\mu=1}^{p}A^{s\mu}_{q_{1}}\lambda_{\mu}\right)^{s_{1}},\left(\sum_{\mu=1}^{p}A^{s\mu}_{q_{2}}\lambda_{\mu}\right)^{s_{2}}\right),$

where $q_{1}=1,\ldots,q_{h},q_{2}=q_{h}+1,\ldots,q$, to be an $m\times q$
matrix which is a proper decomposition of the matrix
$\sum_{\mu=1}^{p}A^{s\mu}_{j}\lambda_{\mu}$ for $s=1,\ldots,m$ and
$j=1,\ldots,q$ into two submatrices of dimensions $m\times q_{1}$ and $m\times
q_{2}$.

If the matrix $A_{1}\lambda$ is a nonsingular $m\times m$ matrix, then one can
determine a quantity

$\gamma_{1}=(A_{1}\lambda)^{-1}(b-(A_{2}\lambda)\gamma_{2})$

and there exists a correspondence
$\gamma\otimes\lambda\mapsto(\gamma_{2},\lambda)$, which determines a map on
the set of simple integral elements for which the domain is formed by these
elements that have the nonsingular matrix $A_{1}\lambda$. All possible
divisions of such a type determine an atlas composed from $(lm)$ maps covering
the set of simple integral elements.

Note that it is remarkable that one can divide the set of simple integral
elements $\gamma_{(s)}\otimes\lambda^{(s)}$ via strata numbered by rank of
matrices $\sum_{\mu=1}^{p}A^{s\mu}_{j}\lambda_{\mu}$. Regular elements form
strata of the highest dimension. One can determine atlases for other strata
analogously.

##  4\. Reduction to homogeneous first-order systems

As an application of the approach of the previous section, let us study the
partial differential equations

$u_{t}-(1+\beta^{2}x^{2})^{\kappa}u_{xx}=0,\qquad
u\in\mathbb{R},\qquad\beta\in\mathbb{R},\qquad\kappa\in\mathbb{N}\cup\\{0\\},$
(4.1)

which reproduce the heat equation and the equation of Brownian motion as
particular cases. We introduce

$A^{t}=\left[\begin{array}[]{cc}0&0\\\ 0&-1\end{array}\right],\quad
A^{x}=\left[\begin{array}[]{cc}0&(1+\beta^{2}x^{2})^{\kappa}\\\
1&0\end{array}\right],\quad b=\left[\begin{array}[]{c}a\\\
0\end{array}\right],\quad V=\left[\begin{array}[]{c}a\\\ c\end{array}\right],$

where $a=u_{t}$ and $c=u_{x}$. Then, equation (4.1) can be written as

$A^{t}V_{t}+A^{x}V_{x}=b.$

According to (3.2), in this case, the matrix $M$ takes the form

$M=\frac{1}{a}\left[\begin{array}[]{cc}1&0\\\ 0&a\end{array}\right].$

One obtains that the inhomogeneous system of partial differential equations

$MA^{t}V_{t}+MA^{x}V_{x}=Mb$

can be rewritten in the new variables

$\widetilde{V}=V-\left[\begin{array}[]{c}y\\\ 0\end{array}\right]$

as a homogeneous system

$\left[\begin{array}[]{cc}0&0\\\
0&-1\end{array}\right]\widetilde{V}_{t}+\left[\begin{array}[]{cc}0&\frac{(1+\beta^{2}x^{2})^{\kappa}}{\widetilde{a}+y}\\\
1&0\end{array}\right]\widetilde{V}_{x}+\left[\begin{array}[]{cc}1&0\\\
0&1\end{array}\right]\widetilde{V}_{y}=\left[\begin{array}[]{c}0\\\
0\end{array}\right].$ (4.2)

Let us consider another simple example: the wave equation

$(\partial^{2}_{t}-\partial^{2}_{x}+k^{2})u=0,\qquad k\in\mathbb{R},$ (4.3)

appearing in relativity problems concerning the existence of gravitational
waves. To study it as a linear homogeneous partial differential equation, we
introduce new coordinates $v_{0}=u_{t}$ and $v_{1}=u_{x}$. Hence, (4.3) is
equivalent to

$\partial_{t}v_{0}-\partial_{x}v_{1}=-k^{2}u,\qquad\partial_{t}u=v_{0},\qquad\partial_{x}u=v_{1}.$

Consequently, one can write

$A^{t}=\left[\begin{array}[]{ccc}1&0&0\\\ 0&0&1\\\ 0&0&0\\\
\end{array}\right],\,\,A^{x}=\left[\begin{array}[]{ccc}0&-1&0\\\ 0&0&0\\\
0&0&1\\\ \end{array}\right],\,\,b=\left[\begin{array}[]{c}-k^{2}u\\\ v_{0}\\\
v_{1}\end{array}\right],\,\,V=\left[\begin{array}[]{c}v_{0}\\\ v_{1}\\\
u\end{array}\right]$

and therefore

$A^{t}V_{t}+A^{x}V_{x}=b.$ (4.4)

We consider the matrix

$M=\frac{1}{-k^{2}u}\left[\begin{array}[]{ccc}1&0&0\\\ -v_{0}&-k^{2}u&0\\\
-v_{1}&0&-k^{2}u\end{array}\right]$

and the change of variables

$\widetilde{V}=V-\left[\begin{array}[]{c}\widehat{x}\\\ 0\\\
0\end{array}\right],$

which maps (4.4) onto

$\frac{1}{-k^{2}u}\left[\begin{array}[]{ccc}1&0&0\\\
-\widehat{v_{0}}-\widehat{x}&0&-k^{2}u\\\
-v_{1}&0&0\end{array}\right]\widetilde{V}_{t}-\frac{1}{k^{2}u}\left[\begin{array}[]{ccc}0&-1&0\\\
0&\widehat{v_{0}}+\widehat{x}&0\\\
0&v_{1}&-k^{2}u\end{array}\right]\widetilde{V}_{x}\\\
+\left[\begin{array}[]{ccc}1&0&0\\\ 0&1&0\\\
0&0&1\end{array}\right]\widetilde{V}_{\widehat{x}}=\left[\begin{array}[]{c}0\\\
0\\\ 0\end{array}\right].$

##  5\. On $k$-wave solutions

A particular solution $u:\mathbb{R}^{p}\rightarrow\mathbb{R}^{q}$ of the
quasi-linear homogeneous system (3.6) such that
$T_{x}u:T_{x}\mathbb{R}^{p}\rightarrow T_{u(x)}\mathbb{R}^{q}$ has rank $k$ at
every point of $\mathbb{R}^{p}$ is called a $k$-wave solution. To simplify the
notation, we will use the Einstein summation convention. In this notation,
indices will be summed between the proper limits for each type of variable,
e.g. the independent variables are summed over their indices ranging from $1$
to $p$, while dependent variables are summed over their indices from $1$ to
$q$.

Let us consider solutions of the hydrodynamic equations given by

$u(x)=f\left(\lambda^{(1)}_{j}(u(x))x^{j},\ldots,\lambda^{(k)}_{j}(u(x))x^{j}\right)$

for certain functions $f:\mathbb{R}^{l}\rightarrow\mathbb{R}^{q}$,
$\lambda^{(i)}:\mathbb{R}^{q}\rightarrow\mathbb{R}^{p}$, with $i=1,\ldots,k$
such that $\lambda^{(1)}\wedge\ldots\wedge\lambda^{(k)}\neq 0$ on
$\mathbb{R}^{q}$. The family $\lambda^{(1)},\ldots,\lambda^{(k)}$ is called a
$k$-wave covector and the functions
$r^{(s)}:\mathbb{R}^{p}\times\mathbb{R}^{q}\rightarrow\mathbb{R}$ of the form

$r^{(s)}(x,u):=\lambda^{(s)}_{j}(u)x^{j},\qquad s=1,\ldots,k,$

are referred to as $k$-order Riemann invariants. Thus, the condition

$u(x)=f(r(x,u(x))),\quad r(x,u(x))=(r^{(1)}(x,u(x)),\ldots,r^{(k)}(x,u(x)))$

gives rise to a matrix of partial derivatives with entries

$\frac{\partial u^{\alpha}}{\partial x^{i}}(x)=\frac{\partial
f^{\alpha}}{\partial
r^{(s)}}(r(x,u))\left[\lambda_{i}^{(s)}(u(x))+x^{j}\frac{\partial\lambda_{j}^{(s)}}{\partial
u^{\beta}}(u(x))\frac{\partial u^{\beta}}{\partial x^{i}}(x)\right]$ (5.1)

and then

$\left[\delta_{\alpha}^{\beta}-x^{j}\frac{\partial f^{\alpha}}{\partial
r^{(s)}}(r(x,u(x)))\frac{\partial\lambda_{j}^{(s)}}{\partial
u^{\beta}}(u(x))\right]\frac{\partial u^{\beta}}{\partial
x^{i}}(x)=\lambda_{i}^{(s)}(u(x))\frac{\partial f^{\alpha}}{\partial
r^{(s)}}(r(x,u(x))).$

If $\Psi(x,u)$ is an invertible $q\times q$ matrix whose entries are given by
the functions

$\Psi^{\alpha}_{\beta}(x,u):=\delta_{\alpha}^{\beta}-x^{j}\frac{\partial
f^{\alpha}}{\partial r^{(s)}}(r(x,u))\frac{\partial\lambda_{j}^{(s)}}{\partial
u^{\beta}}(u),$

we get that $\Psi(0,u)$ is the $q\times q$ identity matrix and $\Psi(x,u)$ is
invertible if $x$ is close enough to zero. On the other hand,

$\Psi^{\alpha}_{\beta}(x,u)\frac{\partial f^{\beta}}{\partial
r^{(s)}}(r(x,u))\\\ =\frac{\partial f^{\alpha}}{\partial
r^{(s)}}(r(x,u))-x^{j}\frac{\partial f^{\alpha}}{\partial
r^{(v)}}(r(x,u))\frac{\partial\lambda_{j}^{(v)}}{\partial
u^{\beta}}(u)\frac{\partial f^{\beta}}{\partial r^{(s)}}(r(x,u))\\\
=\frac{\partial f^{\alpha}}{\partial r^{(s)}}(r(x,u))-\frac{\partial
f^{\alpha}}{\partial r^{(v)}}(r(x,u))\frac{\partial r^{(v)}}{\partial
u^{\beta}}(r(x,u))\frac{\partial f^{\beta}}{\partial r^{(s)}}(r(x,u)).$

Therefore,

$\Psi^{\alpha}_{\beta}(x,u)\frac{\partial f^{\beta}}{\partial
r^{(s)}}(r(x,u))\\\ =\left[\delta_{s}^{v}-\frac{\partial f^{\beta}}{\partial
r^{(s)}}(r(x,u))\frac{\partial r^{(v)}}{\partial
u^{\beta}}(r(x,u))\right]\frac{\partial f^{\alpha}}{\partial
r^{(v)}}(r(x,u))=[\Phi(x,u)]^{v}_{s}\frac{\partial f^{\alpha}}{\partial
r^{(v)}}(r(x,u)),$

where

$[\Phi(x,u)]^{v}_{s}:=\delta_{s}^{v}-\frac{\partial f^{\beta}}{\partial
r^{(s)}}(r(x,u))\frac{\partial r^{(v)}}{\partial u^{\beta}}(r(x,u)).$

Therefore,

$\frac{\partial f^{\alpha}}{\partial
r^{(v)}}(r(x,u))=[\Phi^{-1}(x,u)]^{s}_{v}\Psi^{\alpha}_{\beta}(x,u)\frac{\partial
f^{\beta}}{\partial r^{(s)}}(r(x,u)).$

From this, one gets that

$\frac{\partial u^{\alpha}}{\partial
x^{i}}(x)=(\Psi^{-1})^{\alpha}_{\beta}(x,u)\lambda^{(s)}_{i}(u(x))\frac{\partial
f^{\beta}}{\partial r^{(s)}}(r(x,u))\\\
=(\Psi^{-1})^{\alpha}_{\beta}(r(x,u))\lambda_{i}^{(s)}[\Phi^{-1}(x,u)]^{s}_{t}\Psi^{\beta}_{\delta}(x,u)\frac{\partial
f^{\delta}}{\partial r^{(t)}}(r(x,u))\\\
=[\Phi^{-1}(x,u)]^{t}_{s}\lambda_{i}^{(s)}(u(x))\frac{\partial
f^{\alpha}}{\partial r^{(t)}}(r(x,u)).$

Consider a function of the form
$\zeta:u\in\mathbb{R}^{q}\mapsto(\zeta^{1}(u),\ldots,\zeta^{p}(u))\in\mathbb{R}^{p}$
so that $\langle\lambda(u),\zeta(u)\rangle=0$. Then,

$\zeta^{j}(u(x))\frac{\partial u^{\alpha}}{\partial
x^{j}}(x)=[\Phi^{-1}(x,u)]^{t}_{s}\zeta^{j}(u(x))\lambda_{j}^{(r)}(u(x))\frac{\partial
f^{\alpha}}{\partial r^{(t)}}(r(x,u))=0$

and the coordinates $u^{1}(x),\ldots,u^{q}(x)$ of a solution $u(x)$ are
invariant relative to the vector fields

$X_{\zeta,u(x)}=\zeta^{i}(u(x))\frac{\partial}{\partial x^{i}}.$

In fact,

$(X_{\zeta,u(x)}u^{\alpha})(x)=\zeta^{i}(u(x))\frac{\partial
u^{\alpha}}{\partial x^{i}}(x)\\\
=\zeta^{i}(u(x))[\Phi^{-1}(x,u(x))]^{t}_{s}\zeta^{j}(u(x))\lambda_{j}^{(r)}(u(x))\frac{\partial
f^{\alpha}}{\partial r^{(t)}}(r(x,u))=0.$

Hence, the coordinates $u^{\alpha}(x)$ of any particular solution $u(x)$ of
the hydrodynamic equations (3.6) are constants of the motion of
$X_{\zeta,u(x)}$ and $u(x)$ is constant along its integral characteristic
curves.

##  6\. Conditions for the existence of $k$-wave solutions

We hereafter assume that $b=0$ in (3.1). This section generalises the theory
developed in [11, 12, 14] to the case of $k$-waves for the hydrodynamic
equations (3.1) whose coefficients $A^{1},\ldots,A^{p}$ depend on dependent
and independent variables. This dependence was missing in the theory of
$k$-wave solutions for hydrodynamic equations devised in [14]. Meanwhile, in
[11] the same hydrodynamic equations as in our present work were studied, but
only two-waves, were dealt with.

Most results given in the previous sections can be extended to the most
general case of $k$-wave solutions to be analysed in this section. Proofs
given hereafter essentially follow the ideas described in the previous
sections, but they are technically much more complicated. Moreover, to keep
the calculations easy to read, we skip certain details like the dependence of
each type of function. Nevertheless, these details can be inferred by
analysing the particular cases given in the previous sections.

Recall that a $k$-wave solution of the hydrodynamic equations (3.6) is a
solution $u:\mathbb{R}^{p}\rightarrow\mathbb{R}^{q}$ such that the tangent
mapping $T_{x}u:T_{x}\mathbb{R}^{p}\rightarrow T_{u(x)}\mathbb{R}^{q}$ has
rank $k$ at every $x\in\mathbb{R}^{p}$. Then, the tangent map can be written
as

$T_{x}u=\sum_{\sigma=1}^{k}\xi^{\sigma}(x)\gamma_{(\sigma)}(x,u(x))\otimes\lambda^{(\sigma)}(x,u(x)),\quad
x\in\mathbb{R}^{p},\quad u(x)\in\mathbb{R}^{q},$ (6.1)

where we assume
$\xi:x\in\mathbb{R}^{p}\mapsto(\xi^{1}(x),\ldots,\xi^{k}(x))\in\mathbb{R}^{k}$,
$\lambda^{(\sigma)}:(x,u)\in\mathbb{R}^{p+q}\mapsto(\lambda^{(\sigma)}_{1}(x,u),\ldots,\lambda^{(\sigma)}_{p}(x,u))\in\mathbb{R}^{p*}$,
and $\gamma_{(\sigma)}:\mathbb{R}^{p+q}\rightarrow\mathbb{R}^{q}$ for
$\sigma=1,\ldots,k$. Since $T_{x}u$ has rank $k$ by assumption,
$\gamma_{(1)}\wedge\ldots\wedge\gamma_{(k)}\neq 0$ and
$\lambda^{(1)}\wedge\ldots\wedge\lambda^{(k)}\neq 0$ on $\mathbb{R}^{p+q}$. It
is worth stressing that $T_{x}u$ is denoted by $du$ in other works [11, 14].
Recall also that the functions $u^{1},\ldots,u^{q},\xi^{1},\ldots,\xi^{k}$ are
considered as unknown functions of $x^{1},\ldots,x^{p}$. A $k$-wave solution
of (6.1), understood as a system of partial differential equations relative to
$u(x)$, satisfies the hydrodynamic system (3.6).

From now on, the dependence of $\xi$, $\lambda^{(1)},\ldots,\lambda^{(k)}$ and
$\gamma_{(1)},\ldots,\gamma_{(k)}$, on their domain variables will be omitted
to simplify the notation. We also keep on using the Einstein summation
convention as in the previous sections. Moreover, Latin indices stand for the
coordinates in $\mathbb{R}^{p}$. The Greek index $\sigma$ is used to sum over
the $k$ components of a $k$-wave, while the remaining Greek indices are
employed to sum over the coordinates in $\mathbb{R}^{q}$. Finally, we stress
that we are interested in decompositions (6.1) such that, for fixed values of
the simple elements
$\gamma_{(1)}(x,u)\otimes\lambda^{(1)}(x,u),\ldots,\gamma_{(k)}(x,u)\otimes\lambda^{(k)}(x,u)$,
the values of $\xi^{1}(x_{0}),\ldots,\xi^{k}(x_{0})$ can be arbitrarily
changed at a point $x_{0}\in\mathbb{R}^{p}$ without changing the fact that the
associated function $u(x)$ is a solution of (3.6).

Recall that the form of the vector functions
$\gamma_{(1)},\ldots,\gamma_{(k)}$ depends on the covector functions
$\lambda^{(1)},\ldots,\lambda^{(k)}$ so that $u(x)$ is a solution of the
system (3.6). In particular, the wave relation

$A^{i\alpha}_{\beta}\xi^{\sigma}\gamma^{\beta}_{(\sigma)}\lambda^{(\sigma)}_{i}=0,\qquad\alpha=1,\ldots,q,$
(6.2)

holds. Since the above must be satisfied for every value of
$\xi^{1},\ldots,\xi^{k}$ at a certain fixed point in $\mathbb{R}^{p}$ and
fixed $\lambda^{(\sigma)},\gamma_{(\sigma)}$, it follows that

$A^{i\alpha}_{\beta}\gamma^{\beta}_{(\sigma)}\lambda^{(\sigma)}_{i}=0,\qquad\alpha=1,\ldots,q,\quad\sigma=1,\ldots,k.$
(6.3)

In other words, if we fix the values of $\lambda^{(\sigma)}$, then

$\gamma_{(\sigma)}\in\bigcap_{\alpha=1}^{q}\ker(\lambda^{(\sigma)}_{i}A_{\beta}^{i\alpha})\neq
0,\qquad\sigma=1,\ldots,k.$ (6.4)

Our main aim in this section is to extend to $k$-waves the result that the
single waves of (3.6) can be written as implicit rank-one solutions as given
in the following proposition. Its proof, with less detail than the one below,
can be found in [11].

######  Proposition 6.1.

(Simple Riemann wave solutions) Suppose that there exists a map
$\varphi:\mathbb{R}^{p+q}\rightarrow\mathbb{R}$ such that the set of
implicitly defined relations between the variables $u,x$, and $s$, namely

$u=f(s),\qquad s=\varphi(x,u),$ (6.5)

can be solved so that $s$ and $u$ can be given as a graph over an open subset
$D\subset E$, i.e. one has $s=s(x)$ and $u=f(x)$ for $x$ in an open subset
$A\subset\mathbb{R}^{p}$. Assume that $\gamma(x,u)=\gamma(x^{\prime},u)$
whenever $s(x)=s(x^{\prime})$. Suppose that $\Gamma:u=f(s)$ is a
characteristic curve obtained by solving the system of ODEs

$\frac{du}{ds}=\alpha(s)\widetilde{\gamma}(s,u),$

where $\alpha$ is an arbitrary function of $s$ and
$\widetilde{\gamma}(s(x),u)=\gamma(x,u)$ for $u\in\mathbb{R}^{q}$ and $x\in
A$. Then, (6.5) constitutes an exact (Riemann-type wave) solution of the
quasilinear system (3.6) with coefficients depending on $x$ and $u$.

###### Proof.

The proof is obtained directly by the differentiation of (6.5), i.e. there
exists $s(x)$, univocally determined, such that

$s(x)=\varphi(x,f(s(x)))\Rightarrow(ds)(x)=\frac{(d_{x}\varphi)(x,u(x))}{1-\sum_{\alpha=1}^{q}\frac{\partial\varphi}{\partial
u^{\alpha}}(x,u(x))\frac{df^{\alpha}}{ds}(s(x))}.$

This shows that

$du(x)=\frac{\alpha(s(x))\widetilde{\gamma}(s(x),u(x)))\otimes(d_{x}\varphi)(x,u(x))}{1-\sum_{\alpha=1}^{q}\frac{\partial\varphi}{\partial
u^{\alpha}}(x,u(x))\frac{df^{\alpha}}{ds}(s(x))},$ (6.6)

where $\widetilde{\gamma}(s(x),u(x))=\gamma(x,u(x))$. Hence,

$du(x)=f(x)\gamma(x,u(x))\otimes\lambda(x,u(x)),$ (6.7)

for a certain function $f(x)$ and $du$ becomes a simple integral element (6.7)
of the quasilinear system (2.3), where we exclude the gradient catastrophe
where

$\sum_{\alpha=1}^{q}\frac{\partial\varphi}{\partial
u^{\alpha}}(x,u(x))\frac{df^{\alpha}}{ds}(s(x))=1.$

∎

It is worth noting that the previous proposition implies that
$(d_{x}\gamma)(x,u)$ is proportional to $ds(x)$. Moreover, the solution makes
sense when the denominator in the value of $ds(x)$ is different from zero,
which is just the condition of solvability of the system (6.5). Note also that
$\gamma$ must take a constant value, for each fixed $u\in\mathbb{R}^{q}$, on
the submanifold where $s(x)$ takes a constant value. This is a difference
relative to the formalism of simple waves where $\lambda$ depends only on the
dependent variables.

Note that from the the relation (6.6), it follows that on the hypersurface
$\mathcal{M}$ given by the relations

$u=f(s),\qquad
s=\varphi(x,u),\qquad\sum_{\alpha=1}^{q}\frac{\partial\varphi}{\partial
u^{\alpha}}(x,u(x))\frac{df^{\alpha}}{ds}(s(x))=1,$

the gradient of the function $s(x)$ becomes infinite. The solution can lose
its sense (becomes infinite) on the hypersurface $\mathcal{M}$. Consequently,
some types of discontinuities, namely shock waves, can appear.

For $k$-waves, the tangent map $T_{x}u:T_{x}\mathbb{R}^{p}\rightarrow
T_{u(x)}\mathbb{R}^{q}$ can be understood as an element of
$(\mathbb{R}^{p})^{*}\otimes\mathbb{R}^{q}$. Moreover, each solution, $u(x)$,
can be understood as a zero-form on $\mathbb{R}^{p}$ with values in
$\mathbb{R}^{q}$. As such, the differential of $u$, assuming that it has rank
$k$, takes the form (6.1).

To simplify the notation in what follows, we shall employ two differentials.
Given $f\in C^{\infty}(\mathbb{R}^{p+q},\mathbb{R}^{q})$, namely a function on
$\mathbb{R}^{p+q}$ taking values in $\mathbb{R}^{q}$, we write

$d_{x}f=\frac{\partial f}{\partial x^{i}}dx^{i},\qquad
df=d_{x}f+\frac{\partial f}{\partial u^{\alpha}}d_{x}u^{\alpha}(x).$

The second differential, $d^{2}u(x)$, must vanish and

$d_{x}^{2}u(x)=\xi^{\sigma}d\gamma_{(\sigma)}\wedge\lambda^{(\sigma)}+\gamma_{(\sigma)}\otimes(d\xi^{\sigma}\wedge\lambda^{(\sigma)}+\xi^{\sigma}d\lambda^{(\sigma)})=0.$
(6.8)

Moreover, we define

$X^{i}_{u^{\alpha}}=\frac{\partial X^{i}}{\partial u^{\alpha}},\quad
X^{i}_{Y}=\frac{\partial X^{i}}{\partial u^{\alpha}}Y^{\alpha},\quad\forall
X^{i}\in C^{\infty}(\mathbb{R}^{p+q}),\quad\forall Y\in\mathbb{R}^{q}.$

By using (6.1), we obtain

$\begin{gathered}d\gamma_{(m)}=d_{x}\gamma_{(m)}+\gamma_{(m),\alpha}d_{x}u^{\alpha}=d_{x}\gamma_{(m)}+\xi^{\sigma}\gamma_{(m),\gamma_{(\sigma)}}\otimes\lambda^{(\sigma)},\\\
d\lambda^{(m)}=du^{\alpha}\wedge\lambda^{(m)}_{,u^{\alpha}}+d_{x}\lambda^{(m)}=d_{x}\lambda^{(m)}+\xi^{\sigma}\lambda^{(\sigma)}\wedge\lambda^{(m)}_{,\gamma_{(\sigma)}},\end{gathered}$

for $m=1,\ldots,k$. Substituting the above expressions in (6.8), we get

$\frac{1}{2}\xi^{\sigma}\xi^{\sigma^{\prime}}(\gamma_{(\sigma)},\gamma_{(\sigma^{\prime})})_{u}\otimes\lambda^{(\sigma^{\prime})}\wedge\lambda^{(\sigma)}+\xi^{\sigma}\xi^{\sigma^{\prime}}\gamma_{(\sigma)}\otimes\lambda^{(\sigma^{\prime})}\wedge\lambda_{,\gamma_{(\sigma^{\prime})}}^{(\sigma)}\\\
+\xi^{\sigma}d_{x}\gamma_{(\sigma)}\wedge\lambda^{(\sigma)}+\gamma_{(\sigma)}\otimes
d\xi^{\sigma}\wedge\lambda^{(\sigma)}+\gamma_{(\sigma)}\otimes\xi^{\sigma}d_{x}\lambda^{(\sigma)},$
(6.9)

which has to be satisfied whenever (6.1) holds, where
$(\gamma_{(\sigma)},\gamma_{(\sigma^{\prime})})_{u}$, which can be understood
as a commutator, is an element of
$C^{\infty}(\mathbb{R}^{p+q},\mathbb{R}^{q})$ with coordinates

$(\gamma_{(\sigma)},\gamma_{(\sigma^{\prime})})^{\beta}_{u}=\gamma_{(\sigma^{\prime})}^{\alpha}\frac{\partial\gamma_{(\sigma)}^{\beta}}{\partial
u^{\alpha}}-\gamma_{(\sigma)}^{\alpha}\frac{\partial\gamma_{(\sigma^{\prime})}^{\beta}}{\partial
u^{\alpha}}=\gamma^{\beta}_{(\sigma),\gamma_{(\sigma^{\prime})}}-\gamma^{\beta}_{(\sigma^{\prime}),\gamma_{(\sigma)}},\quad\forall\beta=1,\ldots,q.$

For each fixed value of $x\in\mathbb{R}^{p}$, the functions
$\gamma_{(1)},\ldots,\gamma_{(k)}$ can be understood as the coordinates of the
$x$-dependent vector fields on (see [4] for further details on this notion)
$\mathbb{R}^{q}$ of the form

$\chi_{\sigma}:=\gamma_{(\sigma)}^{\beta}(x,u)\frac{\partial}{\partial
u^{\beta}},\qquad\sigma=1,\ldots,k.$ (6.10)

Hence, $(\gamma_{(\sigma)},\gamma_{(\sigma^{\prime})})_{u}$ can be understood,
for each fixed value $x\in\mathbb{R}^{p}$, as the Lie bracket of the vector
fields on $\mathbb{R}^{q}$ related to $\gamma_{(\sigma)}$ and
$\gamma_{(\sigma^{\prime})}$.

Let us consider an $x$-parametrised differential one-form on $\mathbb{R}^{q}$
satisfying $\omega(x,u)\in{\rm
span}\\{\gamma_{1}(x,u),\ldots,\gamma_{r}(x,u)\\}^{\circ}$, i.e. a mapping
$\omega:\mathbb{R}^{p+q}\rightarrow T^{*}\mathbb{R}^{q}$ such that
$\omega(x,\cdot):\mathbb{R}^{q}\rightarrow T^{*}\mathbb{R}^{q}$ is a standard
differential one-form on $\mathbb{R}^{q}$ for each $x\in\mathbb{R}^{p}$ and
$\omega(x,\cdot)$ vanishes on $\gamma_{(1)}(x),\ldots,\gamma_{(r)}(x)$.
Composing (6.9) with $\omega$, one obtains

$\frac{1}{2}\xi^{\sigma}\xi^{\sigma^{\prime}}\langle\omega,(\gamma_{(\sigma)},\gamma_{(\sigma^{\prime})})_{u}\rangle\lambda^{(\sigma^{\prime})}\wedge\lambda^{(\sigma)}+\xi^{\sigma}\langle\omega,d_{x}\gamma_{(\sigma)}\rangle\wedge\lambda^{(\sigma)}=0.$
(6.11)

By our initial assumption, the above must hold for all possible values of the
functions $\xi^{1},\ldots,\xi^{k}$ at any fixed $x\in\mathbb{R}^{p}$. Then,

$(\gamma_{(\sigma)},\gamma_{(\sigma^{\prime})})_{u}=\sum_{\sigma^{\prime\prime}=1}^{k}\tau^{\sigma^{\prime\prime}}_{\sigma\sigma^{\prime}}\gamma_{(\sigma^{\prime\prime})},\qquad\forall\sigma,\sigma^{\prime}=1,\ldots,k,$
(6.12)

for certain functions $\tau^{\sigma^{\prime\prime}}_{\sigma\sigma^{\prime}}\in
C^{\infty}(\mathbb{R}^{p+q})$, where
$\sigma,\sigma^{\prime},\sigma^{\prime\prime}=1,\ldots,k$. Since
$\gamma_{(1)},\ldots,\gamma_{(k)}$ can be understood as $x$-dependent vector
fields on $\mathbb{R}^{q}$ of the form (6.10) and in view of (6.12), it turns
out that, for every fixed $x\in\mathbb{R}^{p}$, the vector fields
$\chi_{1}(x),\ldots,\chi_{k}(x)$ span an involutive distribution on
$\mathbb{R}^{q}$. Since $\\{\gamma_{(1)},\ldots,\gamma_{(k)}\\}$ are linearly
independent by construction at every $(x,u)\in\mathbb{R}^{p+q}$, the
$x$-dependent vector fields on $\mathbb{R}^{q}$ given by
$\\{\chi_{1},\ldots,\chi_{k}\\}$ span a distribution $\mathcal{D}_{x}(u)={\rm
span}\\{\chi_{1}(x,u),\ldots,\chi_{k}(x,u)\\}$, for $u\in\mathbb{R}^{q}$, of
rank $k$ on $\mathbb{R}^{q}$ for every $x\in\mathbb{R}^{p}$. Hence,
$\mathcal{D}_{x}$ can be integrated on $\mathbb{R}^{q}$ for every
$x\in\mathbb{R}^{p}$ and there exists a family of $x$-parametric foliations of
$\mathbb{R}^{q}$ of $k$-dimensional leaves.

Expression (6.1) demonstrates that the tangent space to each solution $u(x)$
of the hydrodynamic equations at $x\in\mathbb{R}^{p}$ belongs to
$\mathcal{D}_{x}$. It is worth stressing that the distribution
$\mathcal{D}_{x}$ may be different at different values of $x$.

Recalling (6.11) once again, one obtains

$\langle w,d_{x}\gamma_{(m)}\rangle\wedge\lambda^{(m)}=0,\qquad m=1,\ldots,k.$
(6.13)

Let us assume that

$d_{x}\gamma_{(m)}=\beta^{\sigma}_{m}\gamma_{(\sigma)}\otimes\lambda^{(\sigma)}+\pi_{(m)}\otimes\lambda^{(m)},\qquad
m=1,\ldots,k,$ (6.14)

where $\pi_{(m)}$ is any $x$-parametrised vector field on $\mathbb{R}^{q}$, we
sum over the $\sigma$ index and we do not sum over the $m$, and the
$\beta^{\sigma}_{m}$ are $x$-parametrised functions on $\mathbb{R}^{q}$. Note
that condition (6.14) is not necessary for (6.13) to hold, since it can be
proved that a term $\gamma_{(k)}(x,u)\otimes\vartheta(x,u)$, where
$\vartheta(x,u)$ is any $u$-parametrised differential one-form on
$\mathbb{R}^{p}$, is a solution of (6.14).

By substituting (6.14) in (6.9), one obtains

$0=\frac{1}{2}\xi^{\sigma}\xi^{\sigma^{\prime}}(\gamma_{(\sigma)},\gamma_{(\sigma^{\prime})})_{u}\otimes\lambda^{(\sigma^{\prime})}\wedge\lambda^{(\sigma)}+\xi^{\sigma}\xi^{\sigma^{\prime}}\gamma_{(\sigma^{\prime})}\otimes\lambda^{(\sigma)}\wedge\lambda_{,\gamma_{(\sigma)}}^{(\sigma^{\prime})}\\\
+\xi^{\sigma}\beta^{\sigma^{\prime}}_{\sigma}\gamma_{(\sigma^{\prime})}\otimes\lambda^{(\sigma^{\prime})}\wedge\lambda^{(\sigma)}+\gamma_{(\sigma)}\otimes
d\xi^{\sigma}\wedge\lambda^{(\sigma)}+\gamma_{(\sigma)}\otimes\xi^{\sigma}d_{x}\lambda^{(\sigma)}.$
(6.15)

Let us now consider a set of $x$-dependent differential forms on
$\mathbb{R}^{q}$ such that

$\langle\omega_{\sigma^{\prime}},\gamma_{(\sigma^{\prime})}\rangle=\delta_{\sigma^{\prime}}^{\sigma},\qquad\sigma^{\prime},\sigma=1,\ldots,k.$

Applying $\omega_{s}$ to the differential two-form with values in
$\mathbb{R}^{q}$ given by (6.15), one obtains

$0=\frac{1}{2}\xi^{\sigma}\xi^{\sigma^{\prime}}\tau_{\sigma\sigma^{\prime}}^{s}\lambda^{(\sigma^{\prime})}\wedge\lambda^{(\sigma)}+\xi^{s}\xi^{\sigma}\lambda^{(\sigma)}\wedge\lambda_{,\gamma_{(\sigma)}}^{(s)}\\\
+\xi^{\sigma}\beta^{s}_{\sigma}\lambda^{(s)}\wedge\lambda^{(\sigma)}+d\xi^{s}\wedge\lambda^{(s)}+\xi^{s}d_{x}\lambda^{(s)},\quad
s=1,\ldots,k.$ (6.16)

The above expression is not summed over $s$, but it represents $k$ different
equalities for $s=1,\ldots,k$. Multiplying the above expression exteriorly by
$\lambda^{(s)}$, we obtain

$0=\frac{1}{2}\xi^{\sigma}\xi^{\sigma^{\prime}}\tau_{\sigma\sigma^{\prime}}^{s}\lambda^{(\sigma^{\prime})}\wedge\lambda^{(\sigma)}\wedge\lambda^{(s)}\\\
+\xi^{s}\xi^{\sigma^{\prime}}\lambda^{(\sigma^{\prime})}\wedge\lambda_{,\gamma_{(\sigma^{\prime})}}^{(s)}\wedge\lambda^{(s)}+\xi^{s}d_{x}\lambda^{(s)}\wedge\lambda^{(s)},\qquad
s=1,\ldots,k.$ (6.17)

Since the above expression must hold for any set of $\xi^{1},\ldots,\xi^{k}$,
one obtains that

$\tau_{\sigma\sigma^{\prime}}^{s}=0,\qquad(\sigma-\sigma^{\prime})(\sigma-s)(\sigma^{\prime}-s)\neq
0.$ (6.18)

This is a very important result since, as shown in the Appendix, it allows us
to rescale $\gamma_{({1})},\ldots,\gamma_{(k)}$ so that they commute among
themselves and it leads to very important applications, as noted in Section 8.
Moreover, (6.17) gives that

$0=\lambda^{(\sigma)}\wedge\lambda^{(s)}_{,\gamma_{(\sigma)}}\wedge\lambda^{(s)},\qquad\forall\sigma,s=1,\ldots,k.$
(6.19)

Additionally,

$d_{x}\lambda^{(s)}\wedge\lambda^{(s)}=0,\qquad\forall s=1,\ldots,k.$ (6.20)

In view of the integrability condition for Pffafian systems, the forms
$\lambda^{(s)}$ can be reparametrised by multiplying them by a function so
that the new reparametrised differential one-forms, which we still denote by
$\lambda^{(1)},\ldots,\lambda^{(k)}$, will satisfy

$\lambda^{(s)}=d_{x}\varphi^{s}(x,u),\qquad s=1,\ldots,k,$ (6.21)

for certain functions $\varphi^{1}(x,u),\ldots,\varphi^{k}(x,u)$, which are
defined up to a function on $u$ for each index $s=1,\ldots,k$. As
$x\in\mathbb{R}^{p}$, the functions $\varphi^{1},\ldots,\varphi^{k}$ can be
defined on the domain of $\lambda^{(1)},\ldots,\lambda^{(k)}$ as long as it is
simply-connected as a consequence of the Poincaré Lemma. The functions
$\varphi^{1},\ldots,\varphi^{k}$ are called the potential wave functions. It
is relevant to stress that every parametrisation (6.1) satisfying the previous
conditions (6.18)–(6.20) is such that the redefinition of each $\lambda^{(s)}$
due to multiplication by a function so as to satisfy (6.21) leads to a
rescaling of the associated $\gamma_{(s)}$ by the same function so as to keep
a parametrisation (6.1) of the same $k$-wave. Nevertheless, the new functions
$\gamma_{(1)},\ldots,\gamma_{(k)}$ still satisfy (6.18)–(6.20).

Consider the non-empty submanifolds of $\mathbb{R}^{p}$ of the form

$N_{\vec{a},u}:=\\{x\in\mathbb{R}^{p}:\varphi^{\sigma}(x,u)=a^{\sigma},\sigma=1,\ldots,k\\},$

which depend on $u\in\mathbb{R}^{q}$ and
$\vec{a}=(a^{1},\ldots,a^{k})\in\mathbb{R}^{k}$. Using the expression (6.14),
we get

$d_{x}\gamma_{(\sigma)}|_{TN_{\vec{a},u}}=(\beta^{\sigma^{\prime}}_{\sigma}\gamma_{(\sigma^{\prime})}\otimes\lambda^{(\sigma^{\prime})})|_{TN_{\vec{a},u}}=0,\qquad\sigma=1,\ldots,k.$

Therefore, each $\gamma_{(\sigma)}$ is a constant on each $N_{\vec{a},u}$ and

$\gamma_{(\sigma)}=\gamma_{(\sigma)}(\varphi^{1}(x,u),\ldots,\varphi^{k}(x,u),u),\qquad\sigma=1,\ldots,k.$

Recall that the vector fields $\chi_{1}(x),\ldots,\chi_{k}(x)$ span a
foliation on $\mathbb{R}^{q}$ of $k$-dimensional leaves for each value of $x$.
Note that at points of a submanifold
$N_{\vec{a},u}\times\\{u\\}\subset\mathbb{R}^{p+q}$, the vector fields
$\chi_{1},\ldots,\chi_{k}$ always span the same distribution on
$\mathbb{R}^{q}$.

Note that the condition (6.19) implies that if $s,\sigma=1,\ldots,k$ and
$\sigma\neq s$, then $\lambda^{(s)}$ can be written as

$\lambda^{(s)}_{,\gamma_{(\sigma)}}=\zeta_{s1}\lambda^{(s)}+\zeta_{s2}\lambda^{(\sigma)}\Rightarrow\lambda^{(s)}_{,\gamma_{(\sigma)}}|_{TN_{\vec{a},u}}=0$
(6.22)

for some functions $\zeta_{s1},\zeta_{s2}\in C^{\infty}(\mathbb{R}^{p+q})$.

If the condition (6.12) holds, then the distribution $\mathcal{D}_{x}$ is
integrable and there exists a foliation of $\mathbb{R}^{q}$ into
$k$-dimensional surfaces for every $x$. If, in addition, (6.19) and (6.21) are
satisfied, then one can consider an integral submanifold
$\mathcal{S}\subset\mathbb{R}^{q}$ of $\mathcal{D}_{x}$ and consider a
parametrisation, $f=f(\tau^{1},\ldots,\tau^{k})$, of $\mathcal{S}$ so that

$u=f(\tau^{1},\ldots,\tau^{k}),\qquad\tau^{\sigma}=\varphi^{\sigma}(x,f(\tau^{1},\ldots,\tau^{k})),\qquad\sigma=1,\ldots,k.$
(6.23)

For each $x$, the tangent vectors to the leaves of $\mathcal{D}_{x}$ take
values in the distribution $\mathcal{D}_{x}$. Hence, using (6.18) the explicit
parametrisation of the manifold $S$ requires us to solve the following system
of PDEs

$\frac{\partial
f}{\partial\tau^{\alpha}}=\mu_{\alpha}^{\alpha^{\prime}}\gamma_{(\alpha^{\prime})},\qquad\alpha=1,\ldots,k,$
(6.24)

for certain functions
$\mu_{\alpha}^{\alpha^{\prime}}=\mu^{\alpha^{\prime}}_{\alpha}(\tau^{1},\ldots,\tau^{k})$
for $\alpha,\alpha^{\prime}=1,\ldots,k$. Moreover, one can write

$d_{x}u(x)=\frac{\partial f}{\partial\tau^{\alpha}}\otimes
d_{x}\tau^{\alpha}.$

If $u(x)$ is a particular solution parametrised by
$f(\tau^{1},\ldots,\tau^{k})$, then

$d_{x}u(x)=\mu_{\alpha}^{\alpha^{\prime}}\gamma_{(\alpha^{\prime})}\otimes
d_{x}\tau^{\alpha}=\xi^{\alpha^{\prime}}\gamma_{(\alpha^{\prime})}\otimes\lambda^{(\alpha^{\prime})}$

and

$d_{x}\tau^{\alpha}\in\langle\lambda^{(1)},\ldots,\lambda^{(k)}\rangle,\qquad\alpha=1,\ldots,k.$
(6.25)

Since $d_{x}\tau^{1}\wedge\ldots\wedge d_{x}\tau^{k}\neq 0$ and
$\lambda^{(1)}\wedge\ldots\wedge\lambda^{(k)}\neq 0$, one obtains that
$\lambda^{(1)},\ldots,\lambda^{(k)}$ are functions of the form

$\lambda^{(\alpha)}=\lambda^{(\alpha)}(x,f(\tau^{1},\ldots,\tau^{k})),\qquad\alpha=1,\ldots,k.$
(6.26)

The conditions (6.19) and (6.20), in parametric form, are

$\frac{\partial\lambda^{(q)}}{\partial\tau^{m}}=\alpha_{q}^{q}\lambda^{(q)}+\alpha_{m}^{q}\lambda^{(m)},\quad\lambda^{(q)}=d_{x}\varphi^{q}(x,f(\tau^{1},\ldots,\tau^{k})),$
(6.27)

where $q,m=1,\ldots,k,$ with $q\neq m$, and the coefficients $\alpha^{q}_{m}$
depend on the variables $\tau^{1},\ldots,\tau^{k}$.

######  Theorem 6.2.

(A Riemann $k$-wave solution) Let us assume the following conditions:

  1. (1)

$S$ is a $k$-dimensional surface in $\mathbb{R}^{q}$ parametrised by a
function $f$ satisfying (6.24),

  2. (2)

There exists a set of functions $\varphi^{1}(x,u),\ldots,\varphi^{k}(x,u)$
such that

$\lambda^{(\alpha)}=d_{x}\varphi^{\alpha}(x,f(\tau^{1},\ldots,\tau^{k})),\qquad\alpha=1,\ldots,k,$
(6.28)

where the characteristic covectors $\lambda^{(1)},\ldots,\lambda^{(k)}$
satisfy the conditions (6.27) and are related to the $k$-dimensional
distribution spanned by the vector fields associated with
$\gamma_{(1)},\ldots,\gamma_{(k)}$ through the wave relation (6.3),

  3. (3)

The differential $d_{x}f$ has the factorial form

$\frac{\partial f}{\partial\tau^{\alpha}}\otimes
d_{x}\tau^{\alpha}=\gamma_{(\alpha)}\otimes\lambda^{(\alpha)}$ (6.29)

for the $k$ functionally independent differentials
$d_{x}\tau^{1},\ldots,d_{x}\tau^{k}$.

If the set of implicitly defined relations between $u,x,\tau$ of the form

$\begin{cases}u=f(\tau^{1},\ldots,\tau^{k}),\\\
\tau^{\alpha}=\varphi^{\alpha}(x,f(\tau^{1},\ldots,\tau^{k})),\qquad&\alpha=1,\ldots,k,\end{cases}$
(6.30)

can be solved on a certain open set of $\mathbb{R}^{p+q+k}$ so that $u$ and
$\tau$ become functions of $x$, then $u(x)=f(\tau^{1}(x)\ldots,\tau^{k}(x))$,
given by (6.30), constitutes a $k$-wave solution for the hyperbolic
homogeneous system of PDEs (3.6).

###### Proof.

Using our assumptions, we obtain that

$d_{x}u(x)=\frac{\partial f}{\partial\tau^{\alpha}}\otimes
d_{x}\tau^{\alpha}=\mu_{\alpha}^{\alpha^{\prime}}\gamma_{(\alpha^{\prime})}\otimes
d_{x}\tau^{\alpha}=\gamma_{(\alpha)}\otimes\lambda^{(\alpha)}.$

Consider the function $\Psi:\mathbb{R}^{p+q+k}\rightarrow\mathbb{R}^{q+k}$
given by

$\Psi^{\sigma}(x,\tau,u)=\tau^{\sigma}-\varphi^{\sigma}(x,u),\qquad\sigma=1,\ldots,k,$
$\Psi^{\beta+k}(x,\tau,u)=u^{\beta}-f^{\beta}(\tau^{1},\ldots,\tau^{k}),\qquad\beta=1,\ldots,q.$

Then, the condition $\Psi(x,\tau,u)=0$ allows us to determine
$\tau=\tau(x),u=u(x)$ via the Implicit Function Theorem, when the determinant
of the $(q+k)\times(q+k)$ matrix

$A:=\left[\begin{array}[]{c|c}{\rm Id}_{q\times q}&-\frac{\partial
f^{\alpha}}{\partial\tau^{\sigma}}\\\
\hline\cr-\frac{\partial\varphi^{\sigma}}{\partial u^{\alpha}}&{\rm
Id}_{k\times k}\end{array}\right]$

is assumed to be different from zero. By multiplying $A$ by

$\left[\begin{array}[]{c|c}{\rm Id}_{q\times q}&\frac{\partial
f^{\alpha}}{\partial\tau^{\sigma^{\prime}}}\\\ \hline\cr 0&{\rm Id}_{k\times
k}\end{array}\right]$

on the left, we obtain that

$\det A=\det\left[\begin{array}[]{c|c}{\rm Id}_{q\times
q}-\frac{\partial\varphi^{\sigma}}{\partial u^{\alpha}}\frac{\partial
f^{\alpha}}{\partial\tau^{\sigma^{\prime}}}&0\\\
\hline\cr-\frac{\partial\varphi^{\sigma^{\prime}}}{\partial u^{\sigma}}&{\rm
Id}_{k\times k}\end{array}\right],$

where we assume

$\det\left({\rm Id}_{q\times q}-\frac{\partial\varphi^{\sigma}}{\partial
u^{\alpha}}\frac{\partial f^{\alpha}}{\partial\tau^{\sigma}}\right)\neq 0$

to exclude the gradient catastrophe. Recall that from the second relation in
(6.30), one has that

$d_{x}\tau^{\sigma^{\prime}}=\left[\delta^{\sigma^{\prime}}_{\sigma}-\frac{\partial\varphi^{\sigma}}{\partial
u^{\alpha}}\frac{\partial
f^{\alpha}}{\partial\tau^{\sigma^{\prime}}}\right]\lambda^{(\sigma)},$

where $d_{x}\tau^{1}\wedge\ldots\wedge d_{x}\tau^{q}\neq 0$ because of
assumption 3. Hence, $\det A\neq 0$ since
$\lambda^{(1)}\wedge\ldots\wedge\lambda^{(k)}\neq 0$. Using the definition of
characteristic vectors, we obtain that $f(\tau^{1}(x),\ldots,\tau^{k}(x))$
satisfies equation (3.6), which completes the proof. ∎

Surfaces of the form $u=f(\tau^{1},\ldots,\tau^{k})$ are called hodograph
surfaces immersed in $\mathbb{R}^{q}$ and the function
$f(\tau^{1},\ldots,\tau^{k})$ is called a hodograph function. When the
conditions (6.24) and (6.27) are also satisfied, then we say that the surface
is involutive. The functions $\tau^{1},\ldots,\tau^{k}$ are called generalised
Riemann invariants.

Let us prove a relevant consequence of the existence of solutions of (6.1),
understood as a first-order system of partial differential equations in normal
form for $u(x)$. As $u(x)=f(\varphi(x,u))$, where
$\varphi(x,u)=(\varphi^{1}(x,u),\ldots,\varphi^{k}(x,u))$, one has that

$d_{x}u^{\gamma}(x)=\frac{\partial
f^{\gamma}}{\partial\varphi^{a}}d_{x}\varphi^{a}+\frac{\partial\varphi^{a}}{\partial
u^{\beta}}d_{x}u^{\beta}(x)\frac{\partial f^{\gamma}}{\partial\varphi^{a}}.$

Then

$d_{x}u^{\gamma}(x)=(\Psi^{-1})^{\gamma}_{\beta}\frac{\partial
f^{\beta}}{\partial\varphi^{j}}d_{x}\varphi^{j},$

where

$\Psi^{\alpha}_{\beta}=\delta^{\alpha}_{\beta}-\frac{\partial\varphi^{a}}{\partial
u^{\beta}}\frac{\partial f^{\alpha}}{\partial\varphi^{a}}.$

If $\langle\lambda^{(a)}(x,u),\vartheta(x,u)\rangle=0$ for every
$a=1,\ldots,k$, then

$\langle
d_{x}u^{\gamma}(x),\vartheta(x,u)\rangle=(\Psi^{-1})^{\gamma}_{\beta}\frac{\partial
f^{\beta}}{\partial\varphi^{i}}\langle d_{x}\varphi^{i},\vartheta\rangle=0.$

Hence, according to (6.26), the directions in
$\bigcap_{a=1}^{k}\ker\lambda^{(a)}_{x}$ give curves along which $u(x)$ takes
a constant value.

From the proof of Theorem 6.2, we obtain

$\mu^{\alpha^{\prime}}_{\alpha}d\tau^{\alpha}=\lambda^{(\alpha^{\prime})},\qquad\forall\alpha^{\prime}=1,\ldots,k.$

Hence, $\lambda^{(\sigma)}=\lambda^{(\sigma)}(x,f(\tau^{1},\ldots,\tau^{k}))$
for $\alpha=1,\ldots,k$. From the integrability conditions (6.19), it follows
that

$\frac{\partial\lambda^{(\sigma)}}{\partial\tau^{\sigma^{\prime}}}=g^{\sigma}_{\sigma^{\prime}}\lambda^{(\sigma)}+g^{\sigma^{\prime}}_{\sigma}\lambda^{(\sigma^{\prime})},\qquad\sigma\neq\sigma^{\prime},$

are conditions together with (6.27) for the existence of $k$-waves solutions
of (3.6).

##  7\. A geometric approach to hydrodynamic equations

Let us provide a geometric interpretation of the hydrodynamic equations (3.6)
for which we assume that the solutions satisfy

$T_{x}u=\xi^{a}(x)\gamma_{(a)}(x,u)\otimes\lambda^{(a)}(x,u).$ (7.1)

It is worth recalling that $T_{x}u$ stands for the tangent map to
$u:x\in\mathbb{R}^{p}\mapsto u(x)\in\mathbb{R}^{q}$ at $x\in\mathbb{R}^{p}$.
Recall again that, in the literature about hydrodynamic equations, $T_{x}u$ is
more commonly denoted by $du$. Unfortunately, $du$ is misleading from a
geometric point of view, as $d$ can be understood as a generalisation of the
exterior differential on a manifold whose meaning need not be unique or even
properly defined on a map $u:\mathbb{R}^{p}\rightarrow\mathbb{R}^{q}$.

Substituting (7.1) into (3.6), skipping the dependence on independent and/or
dependent variables, and taking into account our previous remarks, we obtain

$A^{i\alpha}_{\beta}\xi^{a}\gamma_{(a)}^{\beta}\lambda^{(a)}_{i}=0,\qquad\alpha=1,\ldots,q.$
(7.2)

Since we want our solutions to remain valid for every value of the functions
$\xi^{1},\ldots,\xi^{k}$ at a point, it follows that

$A^{i\alpha}_{\beta}\gamma_{(a)}^{\beta}\lambda^{(a)}_{i}=0,\qquad\alpha=1,\ldots,q,\quad
a=1,\ldots,k.$ (7.3)

Recall that the $\lambda_{i}^{(a)}$, with $i=1,\ldots,p$, are the coordinates
of each $\lambda^{(a)}$ while $\gamma_{(a)}^{\beta}$, for $\beta=1,\ldots,q$,
are the coordinates of each $\gamma_{(a)}$. Let us rewrite (7.2) more
geometrically by defining some new associated geometric structures. First, let
us define an $\mathbb{R}^{q}$-parametrised family of differential one-forms on
$\mathbb{R}^{p}$ taking values in $\mathbb{R}^{k}$ of the form

$\lambda\\!\\!\\!\\!\lambda=\lambda^{(a)}_{i}dx^{i}\otimes e_{(a)},$ (7.4)

where $\\{e_{(1)},\ldots,e_{(k)}\\}$ is a basis of $\mathbb{R}^{k}$. The space
of such differential forms can be denoted by
$C^{\infty}(\mathbb{R}^{q})\otimes\Omega^{1}(\mathbb{R}^{p})\otimes\mathbb{R}^{k}$.
Additionally, we define a family
$\mathfrak{X}\\!\\!\\!\\!\mathfrak{X}_{a}^{1},\ldots,\mathfrak{X}\\!\\!\\!\\!\mathfrak{X}_{a}^{q}$,
with $a=1,\ldots,k$, of $\mathbb{R}^{q}$-parametrised vector fields on
$\mathbb{R}^{p}$ taking values in $(\mathbb{R}^{k})^{*}$ given by

$\mathfrak{X}\\!\\!\\!\\!\mathfrak{X}_{a}^{\alpha}=A^{i\alpha}_{\beta}\gamma_{(a)}^{\beta}\frac{\partial}{\partial
x^{i}}\otimes e^{(a)},\qquad\alpha=1,\ldots,q,\quad a=1,\ldots,k,$ (7.5)

where $\\{e^{(1)},\ldots,e^{(k)}\\}$ stands for the dual basis of
$\\{e_{(1)},\ldots,e_{(k)}\\}$. It is worth noting that the above expression
is not summed over the values of the subindex $a$. Then, it makes sense to
define the contraction of $\lambda\\!\\!\\!\lambda$ and
$\xi^{a}\mathfrak{X}\\!\\!\\!\\!\mathfrak{X}_{a}^{\alpha}$, where the sum is
over $a=1,\ldots,k$, which makes the hydrodynamic equations look like

$\left\langle\lambda\\!\\!\\!\\!\lambda,\xi^{a}\mathfrak{X}\\!\\!\\!\\!\mathfrak{X}_{a}^{\alpha}\right\rangle=0,\qquad\forall\alpha=1,\ldots,q.$
(7.6)

Since we expect our solution to be valid for all possible values of
$\xi^{1},\ldots,\xi^{k}$ at an arbitrary point, it follows that

$\langle\lambda\\!\\!\\!\\!\lambda,\mathfrak{X}\\!\\!\\!\\!\mathfrak{X}_{a}^{\alpha}\rangle=0,\qquad\forall\alpha=1,\ldots,q,\quad\forall
a=1,\ldots,k.$ (7.7)

We are especially interested in the case when $\lambda\\!\\!\\!\\!\lambda$ is
an exact differential form, i.e. there exists a $\mathbb{R}^{q}$-parametrised
zero-form $\Upsilon\\!\\!\\!\\!\\!\Upsilon$ on $\mathbb{R}^{p}$ taking values
in $\mathbb{R}^{k}$ such that

$\lambda\\!\\!\\!\\!\lambda=d\Upsilon\\!\\!\\!\\!\\!\Upsilon.$ (7.8)

Let us propose a method for obtaining $\lambda\\!\\!\\!\\!\lambda$ while we
ensure that a $k$-wave solution exists. Consider the distribution on
$\mathbb{R}^{p}$, parametrised by elements of $\mathbb{R}^{q}$, spanned by
$\mathfrak{X}\\!\\!\\!\\!\mathfrak{X}_{a}^{1},\ldots,\mathfrak{X}\\!\\!\\!\\!\mathfrak{X}_{a}^{q}$
for $a=1,\ldots,k$. Assume that $\gamma_{(1)},\ldots,\gamma_{(k)}$ are chosen
so that the $\mathfrak{X}\\!\\!\\!\\!\mathfrak{X}_{a}^{\alpha}$, for
$\alpha=1,\ldots,q$ and $a=1,\ldots,k$, span an involutive distribution on
$\mathbb{R}^{p}$ for every value $u\in\mathbb{R}^{q}$. In this case, we can
choose, for every $u\in\mathbb{R}^{q}$, certain $\mathbb{R}^{q}$-dependent
functions $\Upsilon_{1},\ldots,\Upsilon_{k}$ on $\mathbb{R}^{p}$ that are
constant on the leaves of the involutive distribution (at least locally at a
generic point). Hence, (7.7) follows. This has to be done in such a way that
the $\gamma_{(1)},\ldots\gamma_{(k)}$ span an integral distribution on
$\mathbb{R}^{q}$ for every point in $\mathbb{R}^{p}$. We recall that a
$k$-vector field is said to be integrable when its components span an
integrable distribution. Hence, one can say that we are interested in an
integrable $k$-vector field on $\mathbb{R}^{q}$ of the form

$\gamma\\!\\!\\!\gamma=\gamma_{(a)}\otimes e^{(a)},$ (7.9)

for every fixed value of $x\in\mathbb{R}^{p}$. Moreover, we expect that, in
this case, we can consider for every value of $u\in\mathbb{R}^{q}$ the leaves
of the induced distribution on $\mathbb{R}^{p}$ spanned by
$\mathfrak{X}\\!\\!\\!\\!\mathfrak{X}_{a}^{1},\ldots,\mathfrak{X}\\!\\!\\!\\!\mathfrak{X}_{a}^{q}$
for $a=1,\ldots,k$. Hence, we can define a family of functions vanishing on
them. This ensures that $d\Upsilon\\!\\!\\!\\!\\!\Upsilon$ vanishes on the
$\mathfrak{X}\\!\\!\\!\\!\mathfrak{X}^{1}_{a},\ldots,\mathfrak{X}\\!\\!\\!\\!\mathfrak{X}^{q}_{a}$
with $a=1,\ldots,k$. This is to be done in such a way that the vector fields
$\gamma\\!\\!\\!\\!\gamma$ span an involutive distribution on $\mathbb{R}^{q}$
for every fixed $x$.

######  Proposition 7.1.

Let $\mathbb{R}^{k}\ni\tau\mapsto f(\tau)\in\mathbb{R}^{q}$ be an immersed
submanifold $S$ in $\mathbb{R}^{q}$. Consider the
$\mathbb{R}^{q}$-parametrised vector fields on $\mathbb{R}^{p}$ taking values
in $(\mathbb{R}^{k})^{*}$ of the form

$\mathfrak{X}\\!\\!\\!\\!\mathfrak{X}_{a}^{\alpha}=A_{\beta}^{i\alpha}\gamma_{(a)}^{\beta}\frac{\partial}{\partial
x^{i}}\otimes e^{(a)},\qquad\alpha=1,\ldots,q,\qquad a=1,\ldots,k,$

where the $A^{i\alpha}_{\beta}$ appear in (3.6). Suppose that there exists a
$\mathbb{R}^{q}$-parametrised zero-form
$\Upsilon\\!\\!\\!\\!\\!\Upsilon=\sum_{a=1}^{k}T^{a}\otimes e_{(a)}$ on
$\mathbb{R}^{p}$ taking values in $\mathbb{R}^{k}$ such that
$\lambda\\!\\!\\!\\!\lambda=d\Upsilon\\!\\!\\!\\!\\!\Upsilon$ and assume that
the $\gamma_{(a)}=\gamma_{(a)}(x,u)$ take a constant value on the submanifolds
$N_{\vec{a},u}$. Assume also that each $\partial f/\partial T^{s}$ is
proportional to $\gamma_{(s)}$ and the conditions (6.19) are satisfied. Then,
$u(x)=f(T^{1}(x),\ldots,T^{k}(x))$ constitutes an exact Riemann $k$-wave
solution of the system (6.1).

Note that the $T^{s}$ are first-integrals of the vector fields
$\mathfrak{X}\\!\\!\\!\\!\mathfrak{X}^{\alpha}_{a}$. As shown in the examples
of the next section, the existence of the functions $T^{s}$ imply, in some
cases, that the $\mathfrak{X}\\!\\!\\!\\!\mathfrak{X}^{\alpha}_{a}$ span an
involutive distribution.

##  8\. Examples

Now we shall present a few examples which illustrate our theoretical
considerations.

$\bullet$ Let us return to the Brownian motion equation (4.1) with $\kappa=1$.
In this case, the matrix related to $T_{x}u$ is given by

$[T_{x}u]=\left[\begin{matrix}\frac{\partial\widetilde{a}}{\partial
t}&\frac{\partial\widetilde{a}}{\partial
x}&\frac{\partial\widetilde{a}}{\partial y}\\\
\frac{\partial\widetilde{c}}{\partial t}&\frac{\partial\widetilde{c}}{\partial
x}&\frac{\partial\widetilde{c}}{\partial
y}\end{matrix}\right]=\xi\left[\begin{matrix}\gamma^{1}\\\
\gamma^{2}\end{matrix}\right]\otimes\left[\begin{matrix}\lambda_{1}&\lambda_{2}&\lambda_{3}\end{matrix}\right]=\xi\left[\begin{matrix}\gamma^{1}\lambda_{1}&\gamma^{1}\lambda_{2}&\gamma^{1}\lambda_{3}\\\
\gamma^{2}\lambda_{1}&\gamma^{2}\lambda_{2}&\gamma^{2}\lambda_{3}\end{matrix}\right],$

where we recall that $a=u_{t}$ and $c=u_{x}$. Meanwhile, we are interested in
solving the equation

$\lambda_{1}\left[\begin{array}[]{cc}0&0\\\
0&-1\end{array}\right]\left[\begin{array}[]{c}\gamma^{1}\\\
\gamma^{2}\end{array}\right]_{t}+\lambda_{2}\left[\begin{array}[]{cc}0&\frac{1+\beta^{2}x^{2}}{\widetilde{a}+y}\\\
1&0\end{array}\right]\left[\begin{array}[]{c}\gamma^{1}\\\
\gamma^{2}\end{array}\right]_{x}+\lambda_{3}\left[\begin{matrix}1&0\\\
0&1\end{matrix}\right]\left[\begin{array}[]{c}\gamma^{1}\\\
\gamma^{2}\end{array}\right]_{y}\\!\\!=\\!\\!\left[\begin{matrix}0\\\
0\end{matrix}\right].$

and we therefore analyse the vector fields

$X_{1}=\frac{1+\beta^{2}x^{2}}{\widetilde{a}+y}\gamma^{2}\frac{\partial}{\partial
x}+\gamma^{1}\frac{\partial}{\partial y},\qquad
X_{2}=-\gamma^{2}\frac{\partial}{\partial
t}+\gamma^{1}\frac{\partial}{\partial x}+\gamma^{2}\frac{\partial}{\partial
y}.$

Note that, indeed, if we want $\varphi(t,x,y,u)$ to exist, the vector fields
$X_{1},X_{2}$ must be such that their Lie bracket belongs to the distribution
on $\mathbb{R}^{3}$ spanned, at every $u\in\mathbb{R}^{2}$ by $X_{1},X_{2}$.
Consider then

$[X_{1},X_{2}]=-X_{1}\gamma^{2}\frac{\partial}{\partial
t}+\left(X_{1}\gamma^{1}-X_{2}\frac{1+\beta^{2}x^{2}}{a}\right)\frac{\partial}{\partial
x}+(X_{1}\gamma^{2}-X_{2}\gamma^{1})\frac{\partial}{\partial y}$

If $\gamma^{2}=0$ and $\gamma^{1}\neq 0$, then the distribution generated by
$X_{1},X_{2}$ for every value of $\widetilde{a},\widetilde{c}$ is involutive
as well as the distribution spanned by $\gamma$. Indeed, the obtained
distribution reads

$\left\langle\frac{\partial}{\partial x},\frac{\partial}{\partial
y}\right\rangle\Rightarrow\Upsilon=t.$

Hence, $\lambda=d_{(t,x,y)}t=dt$. Recall that $\gamma^{1},\gamma^{2}$ must be
constant on the submanifolds where $\Upsilon$ takes a constant value. Then,
$\gamma^{1}=\gamma^{1}(t)$. Since $\gamma$ is a unique vector, we do not have
any additional integrability condition. Solutions can be obtained by setting

$u=\left(\int\gamma^{1}(t)dt,\widetilde{c}_{0}\right)^{T}.$

It is worth noting that this is not the only possible choice for $\Upsilon$.
Hence, we obtain a simple wave.

$\bullet$ Let us analyse a simple hydrodynamical type equation given by

$\frac{\partial u}{\partial t}+xu\frac{\partial u}{\partial
x}+yu^{2}\frac{\partial u}{\partial y}=0,$ (8.1)

whose solutions are of the form $u:\mathbb{R}^{3}\rightarrow\mathbb{R}$. As
previously, let us parametrise the tangent map to the solution $u(p)$ in the
form

$d_{p}u(p)=\gamma(t,x,y,u)\otimes\lambda(t,x,y,u)$

to look for simple waves of (8.1). Then,
$\lambda=(\lambda_{1},\lambda_{2},\lambda_{3})$ and its associated $\gamma$
satisfy the relation

$(\lambda_{1}+xu\lambda_{2}+yu^{2}\lambda_{3})\gamma=0.$ (8.2)

Let us determine the values of
$\lambda=\lambda_{1}dt+\lambda_{2}dx+\lambda_{3}dy$ that give a solution of
(8.2) for any arbitrary associated $\gamma$. In this case,

$\lambda_{1}=-xu\lambda_{2}-yu^{2}\lambda_{3}.$

Hence, $\lambda$ takes the form

$\lambda=(-xudt+dx)\lambda_{2}+(-yu^{2}dt+dy)\lambda_{3}$

for arbitrary values of $\lambda_{2},\lambda_{3}$. To apply our techniques, we
require $d_{p}\lambda=0$. In particular, let us assume $\lambda_{2}=m/x$ and
$\lambda_{3}=k/y$ for any pair of constants $m,k\in\mathbb{R}$ such that
$m^{2}+k^{2}\neq 0$. Hence,

$\lambda=\left(-udt+\frac{dx}{x}\right)m+\left(-u^{2}dt+\frac{dy}{y}\right)k$

satisfies $d_{p}\lambda=0$ at $p=(t,x,y)$. Hence, we can ensure
$\lambda=d_{p}R$ for

$R=-t(um+u^{2}k)+\ln|x^{m}y^{k}|.$

Assuming that $\gamma$ is a constant on the submanifolds of $\mathbb{R}^{3}$,
where $R$ takes a constant value, a solution of (8.1) can be written as
$u=f(R)$, where $f(R)$ is a solution of

$\frac{df}{dR}=\gamma(R,u).$ (8.3)

Since $\gamma$ is an arbitrary function of the form $\gamma(R,u)$, a solution
can be written as

$u=f(R)$

for an appropriate function $\gamma(R,u)$. With respect to $R$, let us recall
that equation (6.5) can be determined in an implicit manner by the equation

$R=-t(f(R)m+f(R)^{2}k)+\ln|x^{m}y^{k}|$

for a solution $f$ of (8.3). In other words, the above expression determines
in an implicit manner the value of $R$ as a function of $t,x,y$. Then, in an
implicit manner,

$u=f(-t(um+u^{2}k)+\ln|x^{m}y^{k}|)$

describes the solutions of (8.1) determined by generalised Riemann invariants.
Let us provide a more particular case. Consider $\gamma=c$, which takes a
constant value on the submanifolds of $\mathbb{R}^{3}$, where $R$ takes a
constant value. Hence, $f(R)=cR$ and

$R=-t(cRm+c^{2}R^{2}k)+\ln|x^{m}y^{k}|,$

which, assuming that $k^{2}+c^{2}\neq 0$, gives that

$R=-\frac{\sqrt{4c^{2}kt\log\left(x^{m}y^{k}\right)+(cmt+1)^{2}}+ctm+1}{2c^{2}kt},$

and

$u=f(R(t,x,y))=-\frac{\sqrt{4c^{2}kt\log\left(x^{m}y^{k}\right)+(ctm+1)^{2}}+ctm+1}{2ckt},$

which can be proved to be an explicit solution of (8.1).

$\bullet$ Let us analyse the hydrodynamic system of partial differential
equations with two dependent and three independent variables of the form

$\begin{cases}\frac{\partial u_{1}}{\partial t}+x\frac{\partial
u_{2}}{\partial x}+yu_{1}\frac{\partial u_{2}}{\partial y}=0,\\\
\frac{\partial u_{2}}{\partial t}+xu_{2}\frac{\partial u_{1}}{\partial
x}+y\frac{\partial u_{1}}{\partial y}=0.\\\ \end{cases}$ (8.4)

We have the parametrization

$[T_{p}u]=\left[\begin{matrix}\frac{\partial u_{1}}{\partial t}&\frac{\partial
u_{1}}{\partial x}&\frac{\partial u_{1}}{\partial y}\\\ \frac{\partial
u_{2}}{\partial t}&\frac{\partial u_{1}}{\partial x}&\frac{\partial
u_{1}}{\partial y}\end{matrix}\right]=\xi\left[\begin{matrix}\gamma^{1}\\\
\gamma^{2}\end{matrix}\right]\otimes\left[\begin{matrix}\lambda_{1}&\lambda_{2}&\lambda_{3}\end{matrix}\right]=\xi\left[\begin{matrix}\gamma^{1}\lambda_{1}&\gamma^{1}\lambda_{2}&\gamma^{1}\lambda_{3}\\\
\gamma^{2}\lambda_{1}&\gamma^{2}\lambda_{2}&\gamma^{2}\lambda_{3}\end{matrix}\right].$

The hydrodynamic equation (8.4) in matrix form can be written as

$\lambda_{1}\left[\begin{array}[]{cc}1&0\\\
0&1\end{array}\right]\left[\begin{array}[]{c}\gamma^{1}\\\
\gamma^{2}\end{array}\right]_{t}+\lambda_{2}\left[\begin{array}[]{cc}0&x\\\
u_{2}x&0\end{array}\right]\left[\begin{array}[]{c}\gamma^{1}\\\
\gamma^{2}\end{array}\right]_{x}+\lambda_{3}\left[\begin{matrix}0&u_{1}y\\\
y&0\end{matrix}\right]\left[\begin{array}[]{c}\gamma^{1}\\\
\gamma^{2}\end{array}\right]_{y}=\left[\begin{matrix}0\\\
0\end{matrix}\right],$

where

$\gamma=\gamma^{1}\frac{\partial}{\partial
u_{1}}+\gamma^{2}\frac{\partial}{\partial u_{2}}.$

Here, consider the problem of the construction of a Riemann double wave. Let
us consider the vector fields (7.5) for (8.4) of the form
$\mathfrak{X}\\!\\!\\!\\!\mathfrak{X}^{1}=X^{1}\otimes e^{(1)}$ and
$\mathfrak{X}\\!\\!\\!\\!\mathfrak{X}^{2}=X^{2}\otimes e^{(1)}$, where

$X_{1}=\gamma^{1}\frac{\partial}{\partial
t}+x\gamma^{2}\frac{\partial}{\partial
x}+yu_{1}\gamma^{2}\frac{\partial}{\partial y},\qquad
X_{2}=\gamma^{2}\frac{\partial}{\partial
t}+xu_{2}\gamma^{1}\frac{\partial}{\partial
x}+y\gamma^{1}\frac{\partial}{\partial y}.$

More conveniently, let us define $\Delta=\gamma^{1}/\gamma^{2}$ and we verify
that the vector fields

$Y_{1}=\Delta\frac{\partial}{\partial t}+x\frac{\partial}{\partial
x}+yu_{1}\frac{\partial}{\partial y},\qquad Y_{2}=\frac{\partial}{\partial
t}+xu_{2}\Delta\frac{\partial}{\partial x}+y\Delta\frac{\partial}{\partial
y},$

span an involutive distribution of rank two, which implies that they have a
common constant of the motion. Note that

$[Y_{1},Y_{2}]=-(Y_{2}\Delta)\frac{\partial}{\partial
t}+(Y_{1}(u_{2}x\Delta)-Y_{2}x)\frac{\partial}{\partial
x}+[Y_{1}(y\Delta)-Y_{2}(u_{1}y)]\frac{\partial}{\partial y}.$

Hence,

$[Y_{1},Y_{2}]=-(Y_{2}\Delta)\frac{\partial}{\partial
t}+(Y_{1}(xu_{2}\Delta)-Y_{2}x)\frac{\partial}{\partial
x}+y(Y_{1}\Delta)]\frac{\partial}{\partial y}.$

Consider $u_{1}>0$. This implies that if $\Delta=\pm\sqrt{u_{1}}$, we obtain
that

$(u_{1}u_{2}-1)x\frac{\partial}{\partial x}=\pm\sqrt{u_{1}}Y_{2}-Y_{1}$

is proportional to $[Y_{1},Y_{2}]$. Note also that such values of
$\Delta_{\pm}$ allow us to define $\gamma^{1},\gamma^{2}$ to be constant.
Hence, they will be constant on the regions where
$\tau_{\pm}=\tau_{\pm}(t,x,y)$ take a constant value, which allows us to
satisfy condition (6.14). Hence, a common constant of the motion for
$Y_{1},Y_{2}$ can be found. In particular, it is immediate that the functions
corresponding to the relation (7.8) become

$\Upsilon^{\pm}=t\mp\frac{\ln|y|}{\sqrt{u_{1}}}$

and are constants of the motion. Therefore, again by using (7.8), we define

$\lambda^{\pm}=d_{x}\Upsilon^{\pm}=dt\mp\frac{dy}{y\sqrt{u_{1}}}.$

To find a simple wave, we consider a characteristic of the
$\mathbb{R}^{3}$-parametrised vector field on $\mathbb{R}^{2}$ of the form

$\gamma_{\pm}=\gamma_{\pm}^{1}\frac{\partial}{\partial
u_{1}}+\gamma_{\pm}^{2}\frac{\partial}{\partial u_{2}}$

Since $\Delta_{\pm}=\pm\sqrt{u_{1}}=\gamma^{1}_{\pm}/\gamma_{\pm}^{2}$, we
choose $\gamma_{\pm}^{1}=\pm\sqrt{u_{1}}$ and $\gamma_{\pm}^{2}=1$. Then, we
look for a solution of the characteristic system

$\frac{du_{1}}{d\tau}=\pm\sqrt{u_{1}},\qquad\frac{du_{2}}{d\tau}=1,$

which reads

$u_{1\pm}(\tau)=\left(\frac{\pm\tau_{\pm}-\tau_{0}}{2}\right)^{2},\qquad
u_{2\pm}(\tau)=\tau_{\pm}-\tau_{0}^{\prime},$

for arbitrary $\tau_{0},\tau_{0}^{\prime}\in\mathbb{R}$. The functions
$\tau_{\pm}$ can be written as functions of $t,x,y$ by solving the implicit
system

$f_{\pm}(\tau_{\pm})=(u_{1\pm}(\tau),u_{2\pm}(\tau)),\qquad\tau_{\pm}=\Upsilon_{\pm}(t,x,y,u_{1},u_{2})$

which give, for instance, the following solutions

$\tau_{\pm}=\frac{t\pm\tau_{0}-\sqrt{(t\mp\tau_{0})^{2}-8\ln|y|}}{2}.$

Consequently, we obtain the explicit solutions of our original system of PDEs
(8.4) given by

$u_{1\pm}(t,x,y)=\frac{1}{4}\left(\pm\frac{t\pm\tau_{0}-\sqrt{(t\mp\tau_{0})^{2}-8\ln|y|}}{2}-\tau_{0}\right)^{2},$
$u_{2\pm}(t,x,y)=\frac{t\pm\tau_{0}-\sqrt{(t\mp\tau_{0})^{2}-8\ln|y|}}{2}-\tau_{0}^{\prime},$

for arbitrary $\tau_{0},\tau_{0}^{\prime}\in\mathbb{R}$. It is worth remarking
that, if $p$ stands for the variables in $\mathbb{R}^{3}$, then

$d_{p}u_{\pm}(t,x,y)\\!\\!=\\!\\!\left[\begin{array}[]{c}\\!\\!\pm\frac{du_{1}}{d\tau^{\pm}}(\tau^{\pm}(t,x,y))\\\
\\!\\!\frac{du_{2}}{d\tau^{\pm}}(\tau^{\pm}(t,x,y))\end{array}\right]\otimes
d_{p}\tau^{\pm}(t,x,y)\\\
=\\!\left[\begin{array}[]{c}\\!\\!\pm\sqrt{u_{1}(t,x,y)}\\\
1\end{array}\right]\otimes d_{p}\tau^{\pm}(t,x,y),$

where $p=(t,x,y)$. But then,

${\small
d_{p}\tau^{\pm}\\!=\\!\frac{\sqrt{(t\mp\tau_{0})^{2}-8\ln|y|}-t\pm\tau_{0}}{2\sqrt{(t\mp
t_{0})^{2}-8\ln|y|}}\left(dt\\!-\\!\frac{4dy}{y(t\mp\tau_{0}-\sqrt{(t\mp\tau_{0})^{2}-8\ln|y|})}\right)}$

and hence

$d_{p}\tau_{\pm}=\frac{\sqrt{(t\mp\tau_{0})^{2}-8\ln|y|}-t\pm\tau_{0}}{2\sqrt{(t\mp
t_{0})^{2}-8\ln|y|}}\lambda^{\pm}(t,x,y,u_{1\pm}(t,x,y),u_{2\pm}(t,x,y)),$

as established in the proof of Proposition 6.5. Therefore, $u_{\pm}(t,x,y)$
are simple wave solutions for (8.4).

Finally, observe that
$\gamma_{\pm}=\pm\sqrt{u_{1}}\partial_{u_{1}}+\partial_{u_{2}},\lambda^{\pm}$
are such that

$[\gamma^{+},\gamma^{-}]=0,\qquad\gamma^{+}\wedge\gamma^{-}\neq
0,\qquad\lambda^{+}\wedge\lambda^{-}\neq 0,$

which allows us to obtain a superposition of simple waves giving rise to a
double-wave solution with

$T_{p}u(t,x,y)=\xi^{+}\gamma_{+}\otimes\lambda^{+}+\xi^{-}\gamma_{-}\otimes\lambda^{-},$

for certain coefficients
$\xi^{+},\xi^{-}:\mathbb{R}^{3}\rightarrow\mathbb{R}$. To obtain a $2$-wave,
we are going to consider the system of partial differential equations

$\frac{\partial
u_{1}}{\partial\tau^{+}}=\frac{1}{2}\sqrt{u_{1}},\qquad\frac{\partial
u_{1}}{\partial\tau^{-}}=-\frac{1}{2}\sqrt{u_{1}},\qquad\frac{\partial
u_{2}}{\partial\tau^{+}}=\frac{1}{2},\qquad\frac{\partial
u_{2}}{\partial\tau^{-}}=\frac{1}{2}.$

This is indeed the coordinate representation of the system

$\frac{\partial
u}{\partial\tau^{+}}=\frac{1}{2}\gamma_{+},\qquad\frac{\partial
u}{\partial\tau_{-}}=\frac{1}{2}\gamma_{-}.$

It is simple to see that a possible solution is given by

$u_{1}=\left(\frac{\tau^{+}-\tau^{-}}{4}\right)^{2},\qquad
u_{2}=\frac{1}{2}(\tau^{+}+\tau^{-}).$

This allows us to define a mapping $f:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$
of the form

$f(\tau^{+},\tau^{-})=\left(\left[\frac{\tau^{+}-\tau^{-}}{4}\right]^{2},\frac{\tau^{+}+\tau^{-}}{2}\right).$

Using $\tau_{+}=t+2\sqrt{-\ln|y|}$ and $\tau_{-}=t-2\sqrt{-\ln|y|}$, we obtain
the particular double-wave solution

$u^{\prime}(t,x,y)=(u_{1}(t,x,y),u_{2}(t,x,y))=(-\ln|y|,t),$

which is indeed a double-wave solution of (8.4) obtained by a superposition of
two single waves. In fact,

$d_{p}\tau^{+}=\lambda^{+},\qquad d_{p}\tau^{-}=\lambda^{-}$

and

$[T_{p}u^{\prime}]=\left[\begin{matrix}0&0&-1/y\\\
1&0&0\end{matrix}\right]=\frac{1}{2}\gamma_{+}\otimes\lambda^{+}+\frac{1}{2}\gamma_{-}\otimes\lambda^{-}.$

Appendix (A modification of Frobenius theorem for integration).

Let us prove a result that seems to be missing in the literature despite its
relation to the Frobenius Theorem and interest in the theory of hydrodynamic
equations. Such a result was employed without proof in [11] to study the
correspondence principle for $k$-waves. More specifically, it was stated in
[11] that the vector fields $\gamma_{(1)},\ldots,\gamma_{(r)}$ can be
redefined to commute among themselves. Nevertheless, in the theory of
$k$-waves, the vector fields $\gamma_{(1)},\ldots,\gamma_{(r)}$ must satisfy
the relation (6.3). Hence, we are only interested in rescaling them by
functions, which ensures that the rescaled vector fields still satisfy (6.3).
Note that the Frobenius Theorem is not sufficient to prove the existence of an
appropriate rescaling of $\gamma_{(1)},\ldots,\gamma_{(r)}$ leading the latter
to commute among themselves. In fact, the Frobenius Theorem proves that a
family of vector fields $X_{1},\ldots,X_{r}$ on a manifold $N$ satisfying
$X_{1}\wedge\ldots\wedge X_{r}\neq 0$ and spanning an involutive distribution
of rank $k$ gives rise to an Abelian Lie algebra of vector fields $\langle
Y_{1},\ldots,Y_{r}\rangle$ spanning the same distribution as
$X_{1},\ldots,X_{r}$ and such that

$Y_{i}=\sum_{j=1}^{r}f_{ij}X_{j}$

for a family of functions $f_{ij}\in C^{\infty}(N)$ with $i,j=1,\ldots,r$.
Nevertheless, the generalised theory of $k$-waves demands that a simple
rescaling of the vector fields $X_{1},\ldots,X_{r}$ by a function must give
rise to the family $Y_{1},\ldots,Y_{r}$ of vector fields that commute among
themselves. Meanwhile, the generalised theory of $k$-waves also assumes a more
restrictive form of the commutator between the vector fields
$\gamma_{(1)},\ldots,\gamma_{(r)}$. This restricted form is employed in the
following theorem, which constitutes a modification of the Frobenius Theorem
that is useful for our purposes. Note that, in the following theorem, for
clarity, the Einstein convention is not used.

######  Theorem 8.1.

(Modified Frobenius Theorem by rescaling) Consider a family
$X_{1},\ldots,X_{r}$ of vector fields on a $q$-dimensional manifold $N$ such
that

  1. (1)

$X_{1}\wedge\ldots\wedge X_{r}\neq 0$

  2. (2)

$[X_{i},X_{j}]=h_{i,j}^{i}X_{i}+h_{i,j}^{j}X_{j}$

for certain functions $h_{i,j}^{k}\in C^{\infty}(N)$ with $i,j,k=1,\ldots,r$.
Then, there exist functions $f_{1},\ldots,f_{r}\in C^{\infty}(N)$ such that
$[f_{i}X_{i},f_{j}X_{j}]=0$ for all $i,j=1,\ldots,r$.

###### Proof.

We proceed by induction. The result is immediate for $r=2$. Assume that our
theorem is true for $r<q$ vector fields on $\mathbb{R}^{q}$ and let us prove
that our theorem remains true for $r+1$ vector fields. Let
$X_{1},\ldots,X_{r},X_{r+1}$ be vector fields on $\mathbb{R}^{q}$ satisfying
the conditions 1 and 2. The first $r$ vector fields can be rescaled, by the
induction hypothesis, in order to obtain $Y_{1},\ldots,Y_{r}$ so that
$[Y_{i},Y_{j}]=0$ with $i,j=1,\ldots,r$. Then, there exist coordinates
$\chi_{1},\ldots,\chi_{r},\eta_{1},\ldots,\eta_{s}$ on $\mathbb{R}^{q}$ such
that

$X_{i}=\frac{\partial}{\partial\chi_{i}},\qquad i=1,\ldots,r.$

Let us determine whether there exists a function $f_{r+1}\in C^{\infty}(N)$
such that

$[X_{i},f_{r+1}X_{r+1}]\wedge X_{i}=0,\qquad i=1,\ldots,r.$ (8.5)

Condition 2 and (8.5) imply that

$[X_{i},f_{r+1}X_{r+1}]=(X_{i}f_{r+1})X_{r+1}+f_{r+1}(h_{i,r+1}^{r+1}X_{r+1}+h_{i,r+1}^{i}X_{i}),\,\,i=1,\ldots,r,$

for certain functions $h_{i,r+1}^{r+1},h_{i,r+1}^{i}$ on $N$ with
$i=1,\ldots,r$. Hence, to satisfy the conditions (8.5), the function $f_{r+1}$
must be such that

$X_{i}\ln|f_{r+1}|=-h_{i,r+1}^{r+1},\qquad i=1,\ldots,r.$ (8.6)

Let us prove that the above system of partial differential equations admits a
solution. Its compatibility condition reads

$\frac{\partial h_{i,r+1}^{r+1}}{\partial\chi_{j}}=\frac{\partial
h_{j,r+1}^{r+1}}{\partial\chi_{i}},\qquad i,j=1,\ldots,r.$

To show that the above compatibility condition holds, let us consider the
Jacobi identities

$[X_{i},[X_{j},X_{r+1}]]+[X_{j},[X_{r+1},X_{i}]]+[X_{r+1},[X_{i},X_{j}]]=0,\qquad
i,j=1,\ldots,r.$

Since $[X_{i},X_{j}]=0$ by the induction hypothesis, the Jacobi identities
reduces to

$[X_{i},[X_{j},X_{r+1}]]-[X_{j},[X_{i},X_{r+1}]]=0,\qquad i,j=1,\ldots,r.$

Using assumption 2, we obtain that

$[X_{i},[X_{j},X_{r+1}]]=\frac{\partial
h_{j,r+1}^{r+1}}{\partial\chi_{i}}X_{r+1}+h_{j,r+1}^{r+1}[X_{i},X_{r+1}]+\frac{\partial
h^{j}_{j,r+1}}{\partial\chi_{i}}X_{j}$

for any $i,j=1,\ldots,r$, and, moreover,

$0=[X_{j},[X_{i},X_{r+1}]]-[X_{i},[X_{j},X_{r+1}]]=\left(\frac{\partial
h_{i,r+1}^{r+1}}{\partial\chi_{j}}-\frac{\partial
h_{j,r+1}^{r+1}}{\partial\chi_{i}}\right)X_{r+1}\\\
+h_{i,r+1}^{r+1}[X_{j},X_{r+1}]+\frac{\partial
h^{i}_{i,r+1}}{\partial\chi_{j}}X_{i}-h_{j,r+1}^{r+1}[X_{i},X_{r+1}]-\frac{\partial
h^{j}_{j,r+1}}{\partial\chi_{i}}X_{j}\\\ =\left(\frac{\partial
h_{i,r+1}^{r+1}}{\partial\chi_{j}}-\frac{\partial
h_{j,r+1}^{r+1}}{\partial\chi_{i}}\right)X_{r+1}-\left(h_{j,r+1}^{r+1}h_{i,r+1}^{i}-\frac{\partial
h^{i}_{i,r+1}}{\partial\chi_{j}}\right)X_{i}\\\
+\left(h^{r+1}_{i,r+1}h_{j,r+1}^{j}-\frac{\partial
h^{j}_{j,r+1}}{\partial\chi_{i}}\right)X_{j}.$

Since $X_{1}\wedge\ldots\wedge X_{r+1}\neq 0$, it follows that

$\frac{\partial h_{j,r+1}^{r+1}}{\partial\chi_{i}}=\frac{\partial
h_{i,r+1}^{r+1}}{\partial\chi_{j}},\qquad i,j=1,\ldots,r,\qquad i\neq j.$
(8.7)

These equalities imply that the system (8.6) has a solution $f_{r+1}$ and
therefore, $X_{r+1}$ can be rescaled so that
$\widetilde{X}_{r+1}=f_{r+1}X_{r+1}$ satisfies

$[X_{i},\widetilde{X}_{r+1}]\wedge X_{i}=0,\qquad i=1,\ldots,r.$ (8.8)

Let us show that it is possible to rescale $X_{1},\ldots,X_{r}$ so that they
still commute among themselves and, in addition, commute with $X_{r+1}$. In
fact, $[f_{i}X_{i},f_{j}X_{j}]=0$ for every $i,j=1,\ldots,r$ if and only if

$f_{i}=f_{i}(\chi_{i},\eta),\qquad i=1,\ldots,r.$

The conditions (8.8) for $\widetilde{X}_{r+1}$ imply that the latter takes the
form

$\widetilde{X}_{r+1}=\sum_{i=1}^{r}f^{i}_{r+1}(\chi_{i},\eta_{1},\ldots,\eta_{s})\frac{\partial}{\partial\chi_{i}}+\sum_{\mu=1}^{s}f^{\mu}_{r+1}(\eta_{1},\ldots,\eta_{s})\frac{\partial}{\partial\eta_{\mu}}$

for some functions of the form
$f_{r+1}^{i}(\chi_{i},\eta_{1},\ldots,\eta_{s})$ and
$f^{\mu}_{r+1}(\eta_{1},\ldots,\eta_{s})$ for $i=1,\ldots,r$ and
$\mu=1,\ldots,s$. If $[f_{i}X_{i},\widetilde{X}_{r+1}]=0$, the function
$f_{i}$ must satisfy the partial differential equation

$0=\widetilde{X}_{r+1}f_{i}-f_{i}\frac{\partial
f^{i}_{r+1}}{\partial\chi_{i}}=f^{i}_{r+1}(\chi_{i},\eta)\frac{\partial
f_{i}}{\partial\chi_{i}}+\sum_{\mu=1}^{s}f^{\mu}_{r+1}(\eta)\frac{\partial
f_{i}}{\partial\eta_{\mu}}-\frac{\partial
f^{i}_{r+1}}{\partial\chi_{i}}f_{i}=0.$ (8.9)

For each $i=1,\ldots,r$, the equation (8.9) is a linear PDE. Due to the
dependence of its coefficients, it can be considered to be defined on
$\mathbb{R}^{s+1}$ and it always admits a solution by the method of
characteristics. Hence, we obtain that

$[f_{i}X_{i},f_{j}X_{j}]=0,\qquad i,j=1,\ldots,r+1.$

By the induction hypothesis, our theorem follows. ∎

##  9\. Acknowledgements

A.M. Grundland was partially supported by an Operating Grant from NSERC
(Canada) and the research grant ANR-11LABX-0056-LMHLabEX LMH (Fondation
Mathématique Jacques Hadamard, France). J. de Lucas and A.M. Grundland
acknowledge partial support from a HARMONIA project (numer
2016/22/M/ST1/00542) funded by the National Science Centre (Poland). J. de
Lucas would like very much to thank the staff of the Centre de Recherches
Mathématiques (CRM) of the University of Montreal for their hospitality and
assistance during his stay at the CRM during the Covid pandemic.

We would like to dedicate this paper to the memory of our friend Pavel
Winternitz (Université de Montréal), who passed away in February 2021. He was
always a source of inspiration and encouragement.

## References

  * [1] R. Abraham, J.E. Marsden, Foundations of mechanics, Addison-Wesley Publishing Co., Reading, Mass, 1978.
  * [2] G. Boillat, Sur la propagation des ondes, Traité de Physique Théorique et de Physique Mathématique, XXIII Gauthier-Villars, Paris, 1–42 (1965).
  * [3] E. Cartan, Les système différentiels extériurs et leurs applications géometrique, Hermann, Paris, 1946.
  * [4] J.F. Cariñena and J. de Lucas, Lie systems: theory, generalisations, and applications, Dissertationes Math. 479, 1–162 (2011).
  * [5] R. Courant and D. Hilbert, Methods of mathematical physics. Vol. II: Partial differential equations, Interscience Publishers (John Wiley & Sons), New York, 1962.
  * [6] P.W. Doyle and A.M. Grundland, Simple waves and invariant solutions of quasilinear systems, J. Math. Phys. 37, 2969–2978 (1996).
  * [7] B.A. Dubrovin and S.P. Novikov, Hamiltonian formalism of one-dimensional systems of hydrodynamic type and the Bogoliubov-Whitham averaging method, Math. Dokl. 270, 781–785 (1983).
  * [8] E.V. Ferapontov and K.R. Khusnutdinova, On the Integrability of (2+1)-Dimensional Quasilinear Systems, Comm. Math. Phys. 248, 187–206 (2004).
  * [9] A.M. Grundland and B. Huard, Conditional symmetries and Riemann invariants for hyperbolic systems of PDEs, J. Phys. A 40, 4093–4123 (2007).
  * [10] A.M. Grundland and V. Lamothe, Multimode solutions to first-order elliptic quasilinear systems, Acta Appl. Math. 146, 81–113 (2014).
  * [11] A.M. Grundland and B. Lawruk, Lagrange-Riemann method for integration of quasilinear partial differential equations, Analysis II, 1–20 (1991).
  * [12] A.M. Grundland, M.B. Sheftel, and P. Winternitz, Invariant solutions of hydrodynamic type equations, J. Phys. A 33, 81–93 (2000).
  * [13] A.M. Grundland and J. Tafel, Nonclassical symmetry reduction and Riemann wave solutions J. Math. Anal. Appl. 198, 879–892 (1996).
  * [14] A.M. Grundland and P.J. Vassiliou, On the solvability of the Cauchy problem for Riemann double waves by the Monge-Darboux method, Internat. J. Analysis Appl. 11, 221–279 (1991).
  * [15] A.M. Grundland and R. Żelazny, Simple waves in quasilinear systems. Part I and II, J. Math. Phys. 24, 2305–2314 and 2315–2329 (1983).
  * [16] R. Hyun-Ku, First-order partial Differential Equations, Theory and Applications of Hyperbolic Systems of Quasilinear Equations, Dover, New York, 1989.
  * [17] A. Jeffrey, Quasilinear hyperbolic systems and waves, Research Notes Math. 5 Pitman Publishing, London, 1976.
  * [18] F. John and S. Klainerman, Almost global existence to nonlinear wave equations in three space dimensions, Comm. Pure Appl. Math. 37, 443–455 (1984).
  * [19] N.G. Van Kampen, Stochastic differential applications, Phys. Reports, North-Holland Publ. Com., Amsterdam, 1976.
  * [20] Z. Peredzynski, Nonlinear plane $k$-waves and Riemann invariants, Bull. Acad. Polon. Sci. Ser. Tech. 19, 625–632 (1971).
  * [21] Z. Peredzynski, Geometry of interaction of Riemann waves, Advances in Nonlinear Waves, Research notes in Mathematics, Pitman Publishing, London, 1985.
  * [22] B.L. Roždestvenskiǐ and N.N. Janenko, Systems of quasilinear equations and their applications to gas dynamics, Translations Math. Monogr. 55 American Mathematical Society, Providence, 1983.
  * [23] W. Ślebodziński, Formes extérieures et leur applications, PWN, Warsaw, 1954.
  * [24] A. Trautman, On the propagation of information by waves, Pergamon Press, PWN, Warsaw, 1962.
  * [25] G.T. Whithan, Linear and Nonlinear waves, Willey, New York, 1974.