\section{Examples}\label{examples}
\def\Fsigma{F}
\def\wavephi{\thetabar}
In this section we review the classifications of second-order Darboux-integrable \MA/ equations,
due to Goursat and Vessiot, and discuss the connection between our results and
the work of Zvyagin.  We will also give examples of a method for explicitly solving for \BT/s linking
these equations to the wave equation.

\subsection{The Goursat-Vessiot List}
Goursat \cite{G1899} studied non-linear PDE of the form
\begin{equation}\label{sform}
u_{xy} = \Fsigma(x,y,u,u_x,u_y)
\end{equation}
which are Darboux-integrable at the 2-jet level, classifying them up to complex contact transformations that preserve the
form \eqref{sform}.
%(Moreover, it is known from the work of Jur\'a\v{s} \cite{Juras} that any Darboux integrable
% Monge-Ampere equation is contact-equivalent to one of the form \eqref{sform}.)
Using Lie-theoretic techniques,
Vessiot \cite{V1942} reproduced Goursat's classification,
expanded to include linear equations, and showed that some of the equations on Goursat's list were equivalent under more
general contact transformations.  Recently, Biesecker \cite{Mattsthesis} re-proved Vessiot's classification using Cartan's method of equivalence, with respect to real contact transformations.
Retaining Goursat's numbering, the list is:
\begin{align}
(x+y) u_{xy} &= 2\sqrt{u_x u_y}; \tag{I}\\
u\, u_{xy} &= \sqrt{1+u_x^2}\sqrt{1+u_y^2}; \tag{II}\\
(\sin u) u_{xy} &= \sqrt{1+u_x^2}\sqrt{1+u_y^2}; \tag{III} \\
u\, u_{xy} &= \pm \phi(u_x) \psi(u_y), \tag{IV} \\
\intertext{where $\phi(t),\psi(t)$ satisfy the ODE $df/dt \pm t/f= K$ for
some nonzero constant $K$;}
(x+y) u_{xy} &= \gamma(u_x) \gamma(u_y), \tag{V} \\
\intertext{where $\gamma$ is implicitly defined by $\gamma(t)-1 =\exp(t - \gamma(t))$;}
u_x - u\dfrac{u_{xy}}{u_y} &= f\left(x,\dfrac{u_{xy}}{u_y}\right);\tag{VI}\\
u_{xy} &= e^u \sqrt{1+(u_x)^2};\tag{VII} \\
u_x - y\, u_{xy} &= f(x,u_{xy}); \tag{VIII}\\
u_{xy} &= e^u; \tag{IX} \\
u_{xy} &= u_x e^u,\tag{X};\\
u_{xy} &= \left(\dfrac1{u+x}+\dfrac1{u+y}\right)u_x u_y  \tag{XI}. \\
\intertext{(In (VI) and (VIII) the function $f$ is arbitrary.)
To Goursat's original list, Vessiot added representatives of the two equivalence classes of
Darboux-integrable linear equations:}
u_{xy} &= a(x,y)u_x + b(x,y)u_y -a(x,y)b(x,y)u, \tag{XII} \\
\intertext{where $h(x,y)=-a_x$ and $k(x,y)=-b_y$ must satisfy the system $(\ln h)_{xy} = 2h-k$, $(\ln k)_{xy} = 2k-h$
with $h\ne k$; and finally,}
u_{xy} &= \dfrac{2u}{(x+y)^2}.\tag{XIII}
\end{align}

In the above list, we have replaced Goursat's original versions of (VII) and (XI) by
simpler equations that Vessiot showed were equivalent to them by contact transformations;
see \cite{V1942}, part 2, pages 5 and 6, respectively.
Vessiot also observed that
(VI) is contact-equivalent to (X), and (VIII) may be reduced by a contact transformation
to the special case
\begin{equation}
u_{xy} = \dfrac{u_x}{x+y}. \qquad \tag{VIII*}
\end{equation}

\bigskip
Our Theorem \ref{ourtheorem}, together with Goursat's classification, implies the following
\begin{cor}\label{Cor-3} If a second-order \MA/ PDE for one function of two variables is
linked to the standard wave equation by a normal \BT/ with 1-dimensional fibers,
then the PDE is either equivalent to the wave equation by a contact transformation,
or equivalent to one of the equations (I)-(XIII) in the above list.
\end{cor}

\subsection{Zvyagin's List}\label{zsec}
Zvyagin \cite{Zvyagin2} investigated second-order \MA/ equations linked to the standard
wave equation by a \BT/, and asserted that all such transformations
that are non-holonomic are exhausted by a list of six examples in addition to Liouville's equation.
Zvyagin did not publish a proof of this classification, and did not give
explicit forms for the \MA/ equations for some of the transformations on
his list.   He did give an explicit transformation for Goursat-Vessiot equation (I):
\begin{equation}
\sqrt{p}-\sqrt{P} = \sqrt{\dfrac{Z-u}{x-y}},\quad
\sqrt{q}-\sqrt{Q}=-\sqrt{\dfrac{Z-u}{x-y}},\tag{Z.I}
\end{equation}
where, as in \S4, $Z$ is the solution to the wave equation, with
$x$- and $y$- derivatives $P$ and $Q$.  (The $x$- and $y$-coordinates
are preserved by the transformation.)

Corollary \ref{Cor-3} implies that every one of Zvyagin's
transformations must be identifiable with an equation on the Goursat-Vessiot
list.   We have calculated explicit forms for certain transformations
on Zvyagin's list, and we can identify the following
transformations as belonging to equations (II), (III), and (VII), respectively:
\begin{gather}
\label{G2bt}
Z P - u p = \sqrt{Z^2 - u^2} \sqrt{1+P^2},\quad
Z Q - u q = \sqrt{Z^2 - u^2}\sqrt{1+Q^2};
\tag{Z.II} \\[0.15in]
p = \dfrac{(\sinh(z) +\tfrac12 e^{-w})P + e^{(Z-w)/2}\sqrt{P^2-1}}{-\sin u}, \quad
q = \dfrac{(\sinh(z) -\tfrac12 e^{-w})Q + e^{-(w+Z)/2}\sqrt{Q^2-1}}{-\sin u},
\tag{Z.III}\\
\intertext{where $w$ is related to $u$ and $Z$ by
$\cos u = \cosh Z - \tfrac12 e^{-w}$;}
p = (1-2e^{u+Z})P - 2 e^{(u+Z)/2}\sqrt{e^{u+Z}-1}\sqrt{P^2+1},
\qquad
q = -Q - e^{(u-Z)/2}\sqrt{e^{u+Z}-1}.
\tag{Z.VII}
\end{gather}
Each of the above transformations, which preserve the $x$- and $y$-coordinates,
 may be verified as being non-holonomic.
To see how this is done, suppose that the transformation equations, when solved for $p$ and $q$, take the form
\begin{equation}\label{transeq}
p = f(x,y,u,Z,P), \qquad q = g(x,y,u,Z,Q).
\end{equation}
The Cartan system
%\footnote{The {\em Cartan system} for a given
%EDS $\I$ is the smallest Frobenius system that contains $\I$.}
for the \MA/ equation is spanned
by $dx,dy,du,dp$ and $dq$, while the Cartan system for the
wave equation is spanned by $dx,dy,dZ, dP,dQ$.
Recall from \S2 that a transformation is holonomic if the intersection
of these systems is Frobenius.  In light of the
transformation equations \eqref{transeq}, the intersection of
these two systems is spanned by $dx$, $dy$, $dp - f_u\,du$
and $dq - g_u\,du$.  The last two 1-forms are congruent modulo $dx$ and $dy$
to
$$\xi_1 =  f_P dP + f_Z dZ, \qquad \xi_2 = g_Q dQ + g_Z dZ,$$
respectively.  To check that the system $\{ dx, dy, \xi_1, \xi_2\}$ is not Frobenius,
compute
\begin{equation}
\label{holonomicxi}
\begin{aligned}
d\xi_1 \wedge \xi_1 \wedge \xi_2 &\equiv g_Q (f_{uZ} f_P - f_{Pu} f_Z)dP \& dQ \& du \& dZ, \\
d\xi_2 \wedge \xi_1 \wedge \xi_2 &\equiv f_P (g_{uZ} g_Q - g_{Qu} g_Z)dP \& dQ \& du \& dZ
\end{aligned}
\mod dx, dy.
\end{equation}
In each case, the coefficients on the right are nonzero, and we
conclude that the transformation is non-holonomic.

\subsection{Solving for \BT/s}\label{solving}
In this subsection, we will set up systems of PDE
whose solutions are \BT/s to the wave equation
for some examples on the Goursat-Vessiot list.
Although the existence of these transformations follows
from the arguments of \S4, here we will be
able to go further in writing down explicit formulas
for the transformations.   Because we will work with
specific \MA/ equations on the list, we can take advantage
of explicit formulas for the characteristic invariants.
(These invariants are computed, for example, in the dissertation
of M. Biesecker \cite{Mattsthesis}.)

The general approach is as follows.  We write a PDE on the list in the form
\begin{equation}\label{sform4}
s=\Fsigma(x,y,u,p,q).
\end{equation}
This form always has $x$ and $y$ as characteristic invariants, and we assume
these are the only functionally independent invariants up to first order
for the equation (i.e., we assume that the equation is not Monge-integrable).
The \BT/ must take these invariants to
corresponding characteristic invariants for the wave equation. By employing
a change of variables on the wave equation side, of the form $X \mapsto \phi(X)$, $Y \mapsto \psi(Y)$, and
interchanging $X$ and $Y$ if necessary, we may assume that the transformation has
$$x=X, \qquad y=Y.$$
Now suppose that the remaining equations defining the \BT/ take the form
\begin{equation}\label{pandqgeneral}
p = f(x,y,u,Z,P,Q), \qquad q = g(x,y,u,Z,P,Q).
\end{equation}

The \MA/ system on $\R^5$ encoding the PDE \eqref{sform4} is generated
algebraically by the contact form $\theta = du - p\, dx - q\, dy$ and the 2-forms
$$\Omega_1 = (dp - \Fsigma(x,y,u,p,q) dy) \& dx,
\qquad \Omega_2 = (dq - \Fsigma(x,y,u,p,q) dx) \& dy.$$
The defining property of the \BT/ is that substituting \eqref{pandqgeneral}
into $\Omega_1, \Omega_2$ must make them congruent to linear combinations
of $dP \& dx$ and $dQ \& dy$ (the 2-forms defining the \MA/ system
for the wave equation) modulo $\theta$ and the contact form on the
wave equation side,
$$\wavephi = dZ - P\,dx - Q\,dy.$$
In fact, $\Omega_1$ must become congruent to a multiple of $dP \& dx$ and
$\Omega_2$ congruent to a multiple of $dQ \& dy$.  Using \eqref{pandqgeneral},
we compute
\begin{align*}
\Omega_1 &\equiv \left((f_y + f_u g + f_Z Q-\Fsigma) dy + f_P dP + f_Q dQ\right) \& dx
\qquad \mod \theta, \wavephi,\\
\intertext{and}
\Omega_2 &\equiv \left((g_x + g_u f + g_Z P - \Fsigma) dx + g_P dP + g_Q dQ\right) \& dy
\qquad \mod \theta, \wavephi.
\end{align*}
We immediately conclude that $f_Q = g_P = 0$, so that
the transformation is of the form
\begin{equation}\label{pandq}
p = f(x,y,u,Z,P), \qquad q = g(x,y,u,Z,Q),
\end{equation}
and $f,g$ must satisfy
two additional first-order PDEs,
%\begin{equation}\label{genfygx}
%f_y =\Fsigma(x,y,u,f,g)-  f_u g - f_Z Q, \qquad
%g_x =\Fsigma(x,y,u,f,g) - g_u f - g_Z P.
%\end{equation}
\begin{align}
f_y &=\Fsigma(x,y,u,f,g)-  f_u g - f_Z Q, \label{genfy}\\
g_x &=\Fsigma(x,y,u,f,g) - g_u f - g_Z P. \label{gengx}
\end{align}

We derive additional first- and second-order PDEs that $f$ and $g$ must satisfy
by differentiating the conditions so far.
Taking derivatives with respect to $Q$ in \eqref{genfy} and $P$ in \eqref{gengx}
gives
\begin{equation}\label{genfzgz}
f_Z = \left(\Fsigma_q - f_u\right) g_Q,\qquad
g_Z = \left(\Fsigma_p - g_u\right) f_P,
\end{equation}
%\begin{align}
%f_Z = \left(\Fsigma_q - f_u\right) g_Q,\label{genfz} \\
%g_Z = \left(\Fsigma_p - g_u\right) f_P, \label{gengz}
%\end{align}
where the partials $\Fsigma_p = \di\Fsigma/\di p$ and $\Fsigma_q=\di\Fsigma/\di q$ are taken
and then evaluated with $p$ and $q$ given by \eqref{pandq}.   As we will
see in specific cases below, this will sometimes imply that $f$ and $g$
must be linear in $P$ and $Q$.
%Moreover, differentiating \eqref{genfz}
%with respect to $y$ and differentiating \eqref{genfy} with respect to
%$Z$, and using $f_{Zy} = f_{yZ}$, will generically let us determine $g_{Qy}$
%in terms of lower-order derivatives; we can similarly determine $f_{Px}$.

In what follows, let $J_1$ and $J_2$ denote the second-order characteristic
invariants for the given PDE (whose existence makes the equation
Darboux-integrable), expressed in
terms of $x,y,u,p,q$ and the second-order jet coordinates
$r$ and $t$.  (We make the convention that
$J_1$ is invariant along the characteristic curves where $x$ is constant,
and $J_2$ is invariant when $y$ is constant.)
Then the \BT/ must take $J_1$ and $J_2$
to second-order characteristic invariants for the wave equation.
In order to compute these additional constraints, we must take total
$x$- and $y$-derivatives in \eqref{pandq}
to deduce how the second-order jet coordinates $r$ and $t$
transform in terms of those of the wave equation:
\begin{equation}\label{randt}
r = f_x + f_u f + f_Z P + f_P R, \qquad
t = g_y + g_u g + g_Z Q + g_Q T.
\end{equation}
Requiring that, under these substitutions, $J_1$ transforms to be a function of only  $x,P,R$,
and $J_2$ transforms to be a function of only $y,Q,T$,
will lead to
additional second-order PDEs which $f$ and $g$ must satisfy.

We now turn to specific examples.
%
%For example, in the case that $J_1$ and $J_2$ are linear in
%$r$ and $t$, then requiring that the coefficient of $R$ in the transform
%of $J_1$ be a function of $x$ and $P$ only
%will yield additional partial differential equations that
%determine $f_{Pu},f_{Py},f_{Pz}$; we can similarly determine $g_{Pu}, g_{Px}$ and
%$g_{Qz}$.
%
%For each representative equation, we encode the overdetermined system of
%PDEs that $f$ and $g$ must satisfy as an EDS, and use Cartan's test for
%involutivity (see \cite{CFB}, Chapter 7) to investigate the space of
%solutions.  We will only consider solutions which are given by functions
%of $x,y,u,Z,P,Q$, and so we will require the integral submanifolds to
%satisfy the {\em independence condition} $dx \& dy \& du \& dZ \& dP \& dQ \ne 0$
%(i.e., this 6-form pulls back to be nonzero at every point of the
%submanifold).  In what follows, we will refer to integral elements\footnote{%
%For an EDS $\I$ defined on a manifold $\scM$, potential tangent spaces for integral submanifolds of the
%EDS are referred to as {\em integral elements}; these are subspaces in
%the tangent spaces of $\scM$ to which all forms in the ideal restrict to be zero.}
% as {\em admissible} if they are 6-dimensional and
%satisfy the independence condition.

\subsubsection*{Equation IX (Liouville's equation)}
In this case, $\Fsigma = e^u$, and the equations
\eqref{genfy} through \eqref{genfzgz} become
\begin{align}
f_y &= e^u +f_u(Q g_Q - g), \label{Lioufy}\\
g_x &= e^u + g_u( P f_P - f), \label{Liougx}\\
f_Z &= -f_u g_Q,\label{Lioufz} \\
g_Z &= -g_u f_P. \label{Liougz}
\end{align}
As mentioned in \S1, the characteristic invariants are
$$J_1 = r -\tfrac12 p^2, \qquad J_2 =t - \tfrac12 q^2.$$
Under \eqref{pandq}, the first invariant
transforms as
$$r-\tfrac12 p^2 = f_x + f_u f + f_Z P + f_P R - \tfrac12 f^2.$$
Requiring that this be a function of $x,P,R$ only immediately implies
that $f_P$ can depend on $x$ and $P$ only, and that the remaining
terms have no dependence on $u,Z$ or $y$.  This gives us 6
additional second-order PDEs for $f$:
\begin{align}
\label{Liousecond1}
f_{Pu}=f_{Py}=f_{PZ}&=0, \quad & \di_u, \di_y, \di_Z(f_x +f_u f + f_Z P  - \tfrac12 f^2) &=0. \\
\intertext{Similarly, we also get}
\label{Liousecond2}
g_{Qu}=g_{Qy}=g_{QZ}&=0, \quad & \di_u, \di_x, \di_Z (g_y +g_u g + g_Z Q  - \tfrac12 g^2) &=0.
\end{align}
Note that some of these second-order equations are redundant,
in light of the derivatives of \eqref{Lioufz} and \eqref{Liougz}.

Next, we derive additional equations by differentiation.
Note that \eqref{Lioufy} shows that $f_u$ cannot be identically zero;
then, taking a $Q$ derivative of \eqref{Lioufz} shows that  $g_{QQ}=0$.
Similarly, $f_{PP}=0$, so that $f$ and $g$ are linear in $P$ and $Q$.  Thus, we may set
$$f(x,y,u,P,Z) = f^0(x)P + f^1(x,y,u,Z), \qquad
g(x,y,u,Q,Z) = g^0(y)Q + g^1(x,y,u,Z).
$$
In particular, taking the terms in \eqref{Liousecond1}, \eqref{Liousecond2} that are linear
in $P$ and $Q$ respectively gives
\begin{equation}\label{Lioufugu}
\di_u\left((f^0-g^0) f^1_u - f^0 f^1\right) = 0,
\qquad
\di_u\left((g^0-f^0) g^1_u - g^0 g^1\right) = 0.
\end{equation}
Furthermore, equating the $Z$-derivative of \eqref{Lioufy}
with the $y$-derivative of \eqref{Lioufz}, and using the $u$-derivatives of these equations to
determine $f_{yu}$ and $f_{Zu}$, gives the compatibility condition
\begin{align*}
(g^0-f^0) f^1_u g^1_u &= g^0_y f^1_u + e^u g^0;\\
\intertext{we similarly derive}
(f^0-g^0) f^1_u g^1_u &= f^0_x g^1_u + e^u f^0.
\end{align*}
Adding and differentiating with respect to $u$, and using the values
for $f^1_{uu}$ and $g^1_{uu}$ given by \eqref{Lioufugu},
shows that $f^0=-g^0=k$ for some nonzero constant $k$, and $f^1_u g^1_u = \tfrac12 e^u$.
Integrating the remaining equations shows that the most general form for the
transformation is
\begin{equation}\label{Liouvillegeneral}
\begin{aligned}
p &= k P + 2\exp\left(\dfrac{u+k Z + v(x) + w(y)}{2}\right)+v'(x),\\
q &= -k Q + \exp\left(\dfrac{u-k Z - v(x)-w(y)}{2}\right) - w'(y),
\end{aligned}
\end{equation}
where $v(x),w(y)$ are arbitrary functions.

Using the calculation \eqref{holonomicxi}, it is easy to verify that none
of these transformations is holonomic.

%\begin{rem}
%The process of generating compatibility conditions
%for overdetermined systems of PDE, by differentiating and equating mixed partials,
%can be carried out more systematically by rewriting the system of PDE as a Pfaffian
%system, and computing the exterior derivatives of the generator 1-forms.
%It is this approach which we will use in the next two examples.
%Moreover, applying Cartan's Test for involutivity (see \cite{CFB}, Chapter 7)
%to the resulting EDS will tell us
%when we can stop differentiating:  if the system is involutive then
%no further compatibility conditions arise through differentiation,
%and solutions exist that may be constructed by applying
%the Cartan-K\"ahler Theorem.
%\end{rem}
\bigskip
In the next two examples, we will analyze the system of PDEs that $f$ and $g$ must satisfy
using the techniques of exterior differential systems.

\bigskip
\subsubsection*{Equation XIII}%Biesecker 5
This PDE,
$$s = \Fsigma(u,x,y) := \dfrac{2u}{(x+y)^2},$$
has second-order characteristic invariants
$$J_1 =r + \dfrac{2p}{x+y}, \qquad J_2 = t + \dfrac{2q}{x+y}$$
in the $x$- and $y$-directions respectively
(see \cite{Mattsthesis}, Appendix A).
Substituting for $p$ and $r$ from \eqref{pandq} and \eqref{randt}
yields
$$J_1 = \left(f_x+ f_u f + f_Z P+\dfrac{2f}{x+y}\right) + f_P R,$$
so that $f_P$ and the expression in parentheses must be functions of $x$ and $P$ only.
Similarly, we have
$$J_2 = \left(g_y + g_u g + g_Z Q+\dfrac{2g}{x+y}\right) + g_Q T,$$
hence $g_Q$ and the expression in parentheses must be functions of $y$ and $Q$ only.

In this case, \eqref{genfy} through \eqref{genfzgz} specialize to
\begin{align}
f_y &= \Fsigma - (g - Q g_Q)f_u, \label{b5fy} \\
g_x &= \Fsigma -(f - P f_P) g_u, \label{b5gx} \\
f_Z &= -g_Q f_u, \label{b5fz} \\
g_Z &= -f_P g_u. \label{b5gz}
\end{align}
If $f_u$ were identically zero, then $f_Z$ would also be identically zero, but
then $f_y = \Fsigma= 2u/(x+y)^2$ would give a contradiction.  So, we may assume that
$f_u$ and (similarly) $g_u$ are nonzero on an open dense set.
It then follows from \eqref{b5fz} that $g_{QQ}=0$ and from \eqref{b5gz} that $f_{PP}=0$,
i.e., $f$ and $g$ are again linear in $P$ and $Q$.

Differentiating \eqref{b5gx},\eqref{b5gz} with respect to $x$ and $Z$, and equating mixed partials, enables us
to solve for $f_{Px}$ as
\begin{equation}\label{b5fpx}
f_{Px} = f_u (f_P - g_Q) - \dfrac{2 f_P}{g_u(x+y)^2},
\end{equation}
while from \eqref{b5fy},\eqref{b5fz} we similarly obtain
\begin{equation}\label{b5gqy}
g_{Qy} = g_u(g_Q - f_P) - \dfrac{2 g_Q}{f_u(x+y)^2}.
\end{equation}

To encode the PDEs that $f$ and $g$ must satisfy as an exterior differential system,
we will use $x,y,u,Z,P,Q$ as independent variables, and use $f,f_x,f_u,f_P,r_1$
and $g,g_y,g_u,g_Q,t_1$ as dependent variables.
(The role of the coefficients $r_1$ and $t_1$ will
be made clear below.)  We will regard these variables as coordinates on $\R^{16}$.
As stated above, we restrict to the
open subset $\scrU\subset \R^{16}$ where $f_P,f_u,g_Q$ and $g_u$ are nonzero.

The generator 1-forms are $\psi_1$ through $\psi_6$,
where
\begin{align*}
\psi_1 &= -df + f_x dx + f_u du + f_P dP
%-g_Q f_u dZ + \left(\Fsigma- f_u (g - Q g_Q)\right)dy,
+ f_Z dZ + f_y dy,
\\
\psi_2 &= -dg + g_y dy + g_u du + g_Q dQ
%-f_P g_u dZ + \left(\Fsigma -g_u (f - P f_P)\right)dx,
+ g_Z dZ + g_x dx,
\\
\psi_3 &= -df_P + f_{Px} dx, \\
\psi_4 &= -dg_Q + g_{Qy} dy,
\end{align*}
with $f_y,f_Z,g_x,g_Z$ given by equations \eqref{b5fz} through \eqref{b5gx} and
$f_{Px},g_{Qy}$ given by \eqref{b5fpx} and \eqref{b5gqy}.
The remaining generators $\psi_5,\psi_6$
encode the rest of the condition that the second-order characteristic
invariants be preserved.  Differentiating the first term in $J_1$ gives
\begin{multline*}
d\left(f_x + f_z P + f_u f + \dfrac{2f}{x+y}\right) \equiv
d(f_x) + P d(f_Z) + \left(d(f_u)-\dfrac{2\,dy}{(x+y)^2}\right) f\\
+ \left(f_u+\dfrac{2}{x+y}\right)
\left(f_u \theta -g_Q f_u \wavephi + \Fsigma dy\right)
\mod \psi_1, dx, dP.
\end{multline*}
Let $\eta_1$ be the 1-form on the right; then for any \BT/ $\eta_1$ must
be a linear combination of $dx$ and $dP$.  In fact, since only the first
term in $\eta_1$ can contain $dP$, the coefficient of $dP$ in $\eta_1$ must be $f_{Px}$.  Thus,
our remaining generators are
$$
\psi_5 = \eta_1 - f_{Px} dP - r_1 dx,\qquad
\psi_6 = \eta_2 - g_{Qy} dQ - t_1 dy,
$$
where, based a similar calculation of $dJ_2$, we set
$$\eta_2 = d(g_y) + Q d(g_Z) + \left(d(g_u)-\dfrac{2\,dx}{(x+y)^2}\right)g
+ \left(g_u+\dfrac{2}{x+y}\right)
\left(g_u \theta -f_P g_u \wavephi + \Fsigma dx\right).
$$

We seek to construct integral manifolds of the given
differential ideal, i.e., submanifolds of $\scrU$ to which
all the forms in the ideal pull back to be zero.
An {\em integral element} for an EDS is an infinitesimal
version of an integral manifold, i.e., a subspace in
the tangent space to $\scrU$ at some point, to which all
the forms in the ideal restrict to be zero.  Because
we want integral manifolds which are graphs of functions
of $x,y,u,P,Q,Z$, we will only consider integral elements
which are 6-dimensional, and to which the differentials
$dx,dy,du,dP,dQ,dZ$ restrict to be linearly independent;
we will call these {\em admissible} integral elements.

Applying Cartan's Test to the Pfaffian system
generated by $\psi_1, \ldots, \psi_6$ shows that
it has last nonzero Cartan character $s_1 =4$, but is not involutive,
as the space of admissible integral elements
has 2-dimensional fiber at each point.  However, the system becomes
involutive after one prolongation, and this establishes the existence of
the required \BT/s.  The last nonzero Cartan character of the involutive
prolongation is $s_1=2$.  By the Cartan-K\"ahler Theorem (see \cite{CFB}, Chapter 7)
 we conclude that
6-dimensional integral submanifolds, satisfying the independence condition,
exist through every point of $\scrU$, and that the construction of such
submanifolds depends on a choice of $2$ functions of one variable.

The additional 1-forms that generate the prolongation include
\begin{align*}
\psi_7 &= d(f_u) +\left(f_u^2+2\dfrac{f_u}{x+y}\right)dx+\left(f_u g_u-\dfrac{2}{(x+y)^2}\right)dy,\\
\psi_8 &= d(g_u) + \left(g_u^2 + 2\dfrac{g_u}{x+y}\right) dy +\left(f_u g_u - \dfrac{2}{(x+y)^2}\right)dx,
\end{align*}
which are actually defined on the original manifold $\R^{16}$.
These forms vanish on all integral elements of the original system, and if they had
been included in the ideal, it would have been involutive with $s_1 = 2$.

The vanishing of $\psi_7,\psi_8$ implies that $f_u$ and $g_u$ are
functions of $x$ and $y$ only.
Moreover, forms $\psi_7,\psi_8$ define a smaller Pfaffian system, involving only $f_u,g_u$ as functions of $x$ and $y$, which satisfies the Frobenius condition.  This means that $f_u(x,y)$ and $g_u(x,y)$ can be determined by solving systems of ODE.  Once these are determined, substituting the solutions into \eqref{b5fpx} and \eqref{b5gqy} gives a Frobenius system which may be solved for the functions $f_P(x)$ and $g_Q(y)$.  Then $f$ and $g$ may be determined by integrating first-order PDE, with $f$ including an arbitrary function of $x$ and $g$ an arbitrary function of $y$.

For example, by observing that $f_u + g_u$ must satisfy a Riccati equation as a function of $x+y$, we are led to a solution $$f_u = \dfrac{y}{(x+y)x}, \qquad g_u = \dfrac{x}{(x+y)y}.$$ Substituting these into \eqref{b5fpx},\eqref{b5gqy} leads to $f_P + g_Q = k(x+y)/(xy)$ for a constant $k$.  It is simplest to choose $k=0$ with $f_P = 1$ and $g_Q = -1$.  Integrating then gives the solution $$f =  P+\dfrac{y(u+Z)}{x(x+y)},\qquad g = -Q +\dfrac{x(u-Z)}{y(x+y)}.$$

\begin{prop} All \BT/s between (XIII) and the wave equation are holonomic.\end{prop}
\begin{proof}  As noted in \S\ref{zsec}, the holonomic condition is equivalent to
the Pfaffian system on $\B^6$ spanned by $dx,dy$ and
$$dp - f_u du \equiv (-g_Q f_u) dZ + f_P dP, \quad dq -g_u du\equiv (-f_P g_u) dZ + g_Q dQ, \mod dx, dy$$
being Frobenius.
It is straightforward to check that $d(f_P dP - g_Q f_u dZ)$ and $d(g_Q dQ - f_P g_u dZ)$ are
zero modulo $dx,dy$ and the 1-forms of the above EDS.  (For example,
$d(f_P) \equiv 0$ modulo $dx, dy, \psi_1, \ldots, \psi_8$, and
the same is true for $d(g_Q)$ and $d(f_u)$.)
\end{proof}

\bigskip
\subsubsection*{Equation IV}  This PDE has the form
$$s =\Fsigma(u,p,q) := \pm \dfrac{\alpha(p)\beta(q)}{u},$$
where $\alpha$ and $\beta$ are arbitrary solutions of the ODE $df/dt \pm t/f = K$ for some fixed $K\ne 0$.
(We will take the plus sign in these equations, the computation
for the other sign being completely analogous.)  In this case,
\eqref{genfzgz} takes the form
\begin{align}
f_Z = \left( \dfrac{\alpha(f)}u \left(K -\dfrac{g}{\beta(g)}\right) - f_u\right) g_Q, \label{G4fz} \\
g_Z = \left( \dfrac{\beta(g)}u \left(K -\dfrac{f}{\alpha(f)}\right) - g_u\right)f_P. \label{G4gz}
\end{align}
(From now on, instead of writing $\alpha(f)$ and $\beta(g)$,
$\alpha$ and $\beta$ will be understood to be composed with
$f(x,y,u,Z,P)$ and $g(x,y,u,Z,Q)$ respectively.)

Unlike in previous examples,
here it is not valid to conclude that $f$ and $g$ are linear in $P$ and $Q$.
In fact, differentiating
\eqref{G4fz} with respect to $Q$ gives
\begin{equation}\label{G4gqq}
0=(u \beta f_u -(K\beta-g)\alpha)g_{QQ}
+\dfrac{(\beta^2 -(K\beta - g)g)}{\beta^2}\alpha g_Q^2,
% \left(\dfrac{\alpha}u \left(K -\dfrac{g}{\beta}\right) - f_u
% \right)g_{QQ}
%+\dfrac{ g_Q^2\alpha}{u \beta}\left(\dfrac{g}{\beta}\left(K - \dfrac{g}{\beta}\right)-1\right),
\end{equation}
enabling us to determine $g_{QQ}$.  (If the coefficient in front were identically
zero, then $\beta(q)$ would be identically equal to a constant times $q$,
which contradicts $K\ne 0$.)  Similarly, differentiating \eqref{G4gz} yields
\begin{equation}\label{G4fpp}
%0 = \left( \dfrac{\beta}u \left(K -\dfrac{f}{\alpha}\right) - g_u\right)f_{PP}
%+\dfrac{f_P^2 \beta}{u \alpha}\left(\dfrac{f}{\alpha}\left(K - \dfrac{f}{\alpha}\right) - 1\right).
0 =(u \alpha g_u - (K\alpha -f)\beta) f_{PP}
+ \dfrac{(\alpha^2 -(K\alpha - f)f)}\beta f_P^2.
\end{equation}

The characteristic invariants for (IV) are
$$J_1 = \dfrac{r}{\alpha} - \dfrac{\alpha}{u}, \qquad J_2 = \dfrac{t}{\beta} - \dfrac{\beta}{u}.$$
Substituting for $p$ and $r$ from \eqref{pandq} and \eqref{randt} gives
\begin{equation}\label{B4J1after}
J_1 = \dfrac{f_x + f_u f + f_Z P + f_P R}{\alpha} - \dfrac{\alpha}{u},\
\end{equation}
so that $f_P /\alpha$ must be a function of $x$ and $P$ only.
Setting the derivatives of this with respect to $u,y$, and $Z$ equal to zero yields
\begin{equation}\label{G4fpu}
f_{Pu} = \left(K -\dfrac{f}{\alpha}\right)  \dfrac{f_P f_u}{\alpha}, \qquad
f_{Py} = \left(K -\dfrac{f}{\alpha}\right)  \dfrac{f_P f_y}{\alpha},\qquad
f_{PZ} = \left(K -\dfrac{f}{\alpha}\right)  \dfrac{f_P f_Z}{\alpha},
\end{equation}
where $f_y$ is given by \eqref{genfy} and $f_Z$ is given by \eqref{G4fz}.
Similarly, from the $T$ coefficient in $J_2$ we get that $g_Q/\beta$ must be a function of
$y$ and $Q$ only, and hence
\begin{equation}\label{G4gqu}
g_{Qu} = \left(K -\dfrac{g}{\beta}\right)  \dfrac{g_Q g_u}{\beta}, \qquad
g_{Qx} = \left(K -\dfrac{g}{\beta}\right)  \dfrac{g_Q g_x}{\beta},\qquad
g_{QZ} = \left(K -\dfrac{g}{\beta}\right)  \dfrac{g_Q g_Z}{\beta}.
\end{equation}
We may also differentiate \eqref{G4fz} and  \eqref{G4gz} to obtain
equations for $f_{Px}$ and $g_{Qy}$.

We encode the various first- and second-order partial differential
equations for $f$ and $g$ derived so far into an exterior differential system generated by 1-forms $\psi_1, \ldots, \psi_6$, as we did for equation (XIII).
Unlike the previous example,
we do not need to prolong, but instead obtain integrability conditions which take the form
\begin{equation}\label{G4intcond2}
\alpha g_u = \beta f_u,
\end{equation}
and
\begin{equation}\label{G4intcond3}
 (u \beta f_u -(K\beta-g)\alpha) g_Q = (u \alpha g_u - (K\alpha -f)\beta) f_P.
\end{equation}
(Note that, by using \eqref{G4fz}, \eqref{G4gz}, this implies that
$\alpha g_Z = \beta f_Z.$)
With these conditions incorporated into the EDS, it becomes involutive
with last nonzero character $s_1=2$.

Solutions of this system may be obtained by observing that the quantities
\begin{equation}\label{defgammadelta}
\lambda = \dfrac{f_u}{\alpha}, \quad \mu = \dfrac{f_Z}{\alpha},
\quad \gamma = \alpha - (\lambda f + \mu P)u,
\quad \delta = \beta - (\lambda g + \mu Q)u
\end{equation}
must be functions of $x,y,u$ and $Z$ only, and satisfy the following compatible
system of first-order PDE:
\begin{align*}
\dfrac{\di\lambda}{\di u} &=-u\lambda^3 + K \lambda^2 - \dfrac{2}{u}\lambda,
&
\dfrac{\di \mu}{\di u} &=-\dfrac{(u^2\lambda^2 - K u\lambda+1)}{u}\mu,
&
\dfrac{1}{\gamma}\dfrac{\di \gamma}{\di u} &=
\dfrac{1}{\delta}\dfrac{\di \delta}{\di u} = (K-u\lambda)\lambda,
\\
\dfrac{\di\lambda}{\di Z} &= \dfrac{\di\mu}{\di u},
&
\dfrac{\di\mu}{\di Z} &= (K-u\lambda) \mu^2,
&
\dfrac{1}{\gamma}\dfrac{\di \gamma}{\di Z} &=
\dfrac{1}{\delta}\dfrac{\di \delta}{\di Z} = (K-u\lambda)\mu,
\\
\dfrac{\di\lambda}{\di x} &= -\dfrac{(u^2\lambda^2 - K u\lambda+1)}{u^2}\gamma,
&
\dfrac{\di\mu}{\di x} &= (K-u\lambda) \dfrac{\gamma \mu}{u},
&
\dfrac{\di\delta}{\di x} &= (K - u\lambda) \dfrac{\gamma\delta}{u},
\\
\dfrac{\di\lambda}{\di y} &= -\dfrac{(u^2\lambda^2 - K u\lambda+1)}{u^2}\delta,
&
\dfrac{\di\mu}{\di y} &= (K-u\lambda) \dfrac{\delta \mu}{u},
&
\dfrac{\di \gamma}{\di y} &= \dfrac{\di\delta}{\di x}.
\end{align*}
A solution $(\gamma, \delta, \lambda,\mu)$ to this PDE system
may be constructed by integrating successively in the $u$-direction,
the $Z$-direction, the $x$-direction and the $y$-direction.  (Note that
the $x$-dependence of $\gamma$ and the $y$-dependence of $\delta$ are
given by arbitrary functions.)  Once $\gamma$ and $\delta$ are
known, they implicitly determine $f$ and $g$.

Using \eqref{holonomicxi}, one can check that
the resulting \BT/s are holonomic if and only if, in the above system,
$u^2\lambda^2 -Ku\lambda+1=0.$  Thus, holonomic transformations
exist, and depend on fewer arbitrary constants but the same number
of arbitrary functions.  For example, if $K=2$, then a solution to the above
system is given by
$$\lambda =\dfrac1{u}, \qquad \mu = \dfrac{-1}{Z -v(x)-w(y)},\qquad
\gamma = \dfrac{v'(x)u}{Z -v(x)-w(y)}, \qquad \delta = \dfrac{w'(y) u}{Z-v(x)-w(y)}.
$$
Then, using \eqref{defgammadelta},  a holonomic \BT/  is implicitly
defined by
$$\alpha(p) - p = \dfrac{(v'(x)-P)u}{Z-v(x)-w(y)},
\qquad
\beta(q) - q = \dfrac{(w'(y)-Q)u}{Z-v(x)-w(y)}.
$$

\subsection{Summary}
Besides equations (IV) and (XIII) discussed above,
we have also investigated the exterior differential
system for \BT/s to the wave equation for
equations (V), (VII), (IX), (XI) and (XII).
Even if explicit formulas are not available,
in each case we use the Cartan-K\"ahler Theorem to determine
(in terms of the last nonzero Cartan character) the size of the
solution set, in both the holonomic and non-holonomic cases.
The results are summarized in the table below.

\bigskip
\begin{tabular}{|c|c|c|c|}\hline
Equation & Monge-Integrable & Holonomic BTs & Non-holonomic BTs \\ \hline
I 	& no & yes, $s_1=2$	& yes, $s_1=2$ \\
II 	& no & no          &	 yes, $s_1=2$ \\
III 	& no & no 			& yes, $s_1=2$ \\
IV	& no & yes, $s_1=2$ & yes, $s_1=2$ \\
V	& no & yes, $s_1=2$ & yes, $s_1=2$ \\
VI	& yes& yes, $s_1=3$ & yes, $s_3=1$ \\
VII	& no & no & yes, $s_1=2$ \\
VIII	& yes& yes, $s_1=3$ & yes, $s_3=1$ \\
IX	& no & no			& yes, $s_1=2$ \\
X 	& yes& yes, $s_1=3$ & yes, $s_3=1$ \\
XI	& no  & yes, $s_1=2$ & no \\
XII	& no  & yes, $s_1=2$ & no \\
XIII	& no 	 & yes, $s_1=2$ & no \\ \hline
\end{tabular}

\medskip
Note that the approach described in \S\ref{solving}
is not feasible for the Monge-integrable equations (VI, VIII and X),
but for completeness we include them in the table,
together with the results from
the analysis in \S\ref{DtoBmonge}.  The Cartan character for
the system for holonomic \BT/s for such equations is variable,
depending on whether one considers the cases (1), (2), or (3), as
 described in \S\ref{DtoBmonge}.

It is interesting to note that equations (I), (IV) and (V) have both holonomic
and non-holonomic transformations, in roughly the same degree of generality.
In fact, it is possible that these two kinds of \BT/s linking the same pair
of equations may be closely related.  In our previous paper \cite{CI}, we pointed out that
the transformation (Z.I) is a composition of two simpler transformations, a holonomic
\BT/ to the wave equation, and a contact transformation from the wave equation to itself.
It is possible that, more generally, the non-holonomic transformations for these
equations are obtainable from holonomic transformations in this way.
