An application of Green-function methods to gravitational radiation theory
Giampiero Esposito$^{1,2}$
Previous work in the literature has studied gravitational radiation in black-hole collisions at the speed of light. In particular, it had been proved that the perturbative field equations may all be reduced to equations in only two independent variables, by virtue of a conformal symmetry at each order in perturbation theory. The Green function for the perturbative field equations is here analyzed by studying the corresponding second-order hyperbolic operator with variable coefficients, instead of using the reduction method from the retarded flat-space Green function in four dimensions. After reduction to canonical form of this hyperbolic operator, the integral representation of the solution in terms of the Riemann function is obtained. The Riemann function solves a characteristic initial-value problem for which analytic formulae leading to the numerical solution are derived.
LNSIM.tex
\documentclass[preprint,eqsecnum,prd,aps,epsfig,nofootinbib]{revtex4}
\usepackage{epsfig}
\usepackage{amsmath}
\begin{document}

\def\cstok#1{\leavevmode\thinspace\hbox{\vrule\vtop{\vbox{\hrule\kern1pt
\hbox{\vphantom{\tt/}\thinspace{\tt#1}\thinspace}}
\kern1pt\hrule}\vrule}\thinspace}

\title{An application of Green-function methods 
to gravitational radiation theory}
\author{Giampiero Esposito$^{1,2}$}
\thanks{Submitted to {\it Lecture Notes of S.I.M.}, volume
edited by D. Cocolicchio and S. Dragomir, with kind permission
by IOP to use material in \cite{espo01}.}
\address{${ }^{1}$Istituto Nazionale di Fisica Nucleare, Sezione
di Napoli, Complesso Universitario di Monte S. Angelo, Via Cintia,
Edificio N', 80126 Napoli, Italy\\
${ }^{2}$Universit\`a di Napoli Federico II, Dipartimento
di Scienze Fisiche, Complesso Universitario di Monte S. Angelo,
Via Cintia, Edificio N', 80126 Napoli, Italy}

\vspace{2cm}
\begin{abstract}
Previous work in the literature has studied gravitational radiation
in black-hole collisions at the speed of light. In particular, it had
been proved that the perturbative field equations may all be reduced to
equations in only two independent variables, by virtue of a conformal
symmetry at each order in perturbation theory. The Green function for the
perturbative field equations is here analyzed by studying the corresponding
second-order hyperbolic operator with variable coefficients, instead of
using the reduction method from the retarded flat-space Green function
in four dimensions. After reduction to canonical form of this hyperbolic
operator, the integral representation of the solution in terms of the
Riemann function is obtained. The Riemann function solves a characteristic
initial-value problem for which analytic formulae leading to the
numerical solution are derived.
\end{abstract}

\maketitle
\section{Introduction}

In the theory of gravitational radiation, it remains true what 
R. Sachs stressed in his 1963 Les Houches lectures, i.e. that we
understand the following three main features \cite{sach64}:
\vskip 0.3cm
\noindent
(i) We can give a description of radiation at large distances from
its sources in an asymptotically flat universe; this description is
geometrically elegant and sufficiently detailed to analyze all
conceptual experiments concerning the behaviour of test particles or
test absorbers in the far field.
\vskip 0.3cm
\noindent
(ii) How the exact theory relates the far field to the near field and
the sources or non-gravitational fields. 
\vskip 0.3cm
\noindent
(iii) We have approximation methods that make it possible to obtain
numerical results for the amount of radiation emitted in a particular
situation, or for scattering cross-sections, etc.; these approximation
methods do not have a profound geometrical nature, but are very
important in comparing theory with the experiment, in case the latter
succeeds in finding observational evidence in favour of gravitational
waves.
\vskip 0.3cm
Within this framework, it is often desirable to use Green-function 
methods, since the
construction of suitable inverses of differential operators lies still
at the very heart of many profound properties in classical and quantum
field theory. For example, the theory of small disturbances in local field
theory can only be built if suitable invertible operators 
are considered \cite{dewi65}. In a
functional-integral formulation, these correspond to the gauge-field and
ghost operators, respectively 
\cite{dewi84,espo01a}. Moreover, the Peierls bracket on the
space of physical observables, which is a Poisson bracket preserving the
invariance under the full infinite-dimensional symmetry group of the theory,
is obtained from the advanced and retarded Green functions of the theory via 
the supercommutator function \cite{dewi65,dewi84,espo01a,cart00} 
and leads possibly to a deeper approach
to quantization. Last, but not least, a perturbation approach to classical
general relativity relies heavily on a careful construction of Green functions
of operators of hyperbolic \cite{deat92a,deat92b,deat92c,deat96} 
and elliptic \cite{espo00a,espo00b} type. In particular,
following \cite{deat92a,deat92b,deat92c,deat96}, 
we shall be concerned with the axisymmetric collision of two
black holes travelling at the speed of light, each described in the
centre-of-mass frame before the collision by an impulsive plane-fronted shock
wave with energy $\mu$. One then passes to 
a new frame to which a large Lorentz
boost is applied. There the energy $\nu=\mu {\rm e}^{\alpha}$ 
of the incoming shock
$1$ obeys $\nu >> \lambda$, where $\lambda=\mu {\rm e}^{-\alpha}$ is the energy
of the incoming shock $2$ and ${\rm e}^{\alpha} \equiv \sqrt{{1+\beta}\over
{1-\beta}}$ ($\beta$ being the usual relativistic parameter). 
In the boosted frame, to the future of the strong
shock $1$, the metric can be expanded in the form \cite{deat92b,deat96}
\begin{equation}
g_{ab} \sim \nu^{2} \left[\eta_{ab}+\sum_{i=1}^{\infty}
\left({\lambda \over \nu}\right)^{i}h_{ab}^{(i)}\right],
\label{(1.1)}
\end{equation}
where $\eta_{ab}$ is the standard notation for the Minkowski metric. 
The task of solving the Einstein field equations becomes then a problem
in singular perturbation theory, having to 
find $h_{ab}^{(1)},h_{ab}^{(2)},...$
by solving the linearized field equations at first, second, ... order 
respectively in ${\lambda \over \nu}$, once that characteristic initial
data are given just to the future of the strong shock $1$. The perturbation
series (1.1) is physically relevant because, on boosting back to the
centre-of-mass frame, it is found to give an accurate description of
space-time geometry where gravitational radiation propagates at small 
angles away from the forward symmetry axis ${\hat \theta}=0$. The news
function $c_{0}$ (see appendix), 
which describes gravitational radiation arriving at
future null infinity in the centre-of-mass frame, is expected to have the
convergent series expansion \cite{deat92b,deat96}
\begin{equation}
c_{0}({\hat \tau},{\hat \theta})=\sum_{n=0}^{\infty}
a_{2n}({\hat \tau}/\mu)(\sin {\hat \theta})^{2n},
\label{(1.2)}
\end{equation}
with $\hat \tau$ a suitable retarded time coordinate, and $\mu$ the energy
of each incoming black hole in the centre-of-mass frame. 
In \cite{deat92b,deat96} a very
useful analytic expression of $a_{2}({\hat \tau}/\mu)$ was derived,
exploiting the property that perturbative field equations may all be reduced
to equations in only two independent variables, by virtue of a remarkable
conformal symmetry at each order in perturbation theory. The Green 
function for perturbative field equations was then found by reduction 
from the retarded flat-space Green function in four dimensions. 

However, a {\it direct} approach to the evaluation of Green functions
appears both desirable and helpful in general, and it has been our 
aim to pursue such a line of investigation. For this purpose, 
following hereafter our work in \cite{espo01}, reduction
to two dimensions with the associated hyperbolic operator is studied
again in section 2. Section 3 performs reduction to canonical form
with the associated Riemann function. Equations for the Goursat problem
obeyed by the Riemann function are derived in section 4, while the
corresponding numerical algorithm is discussed in section 5.
Some backgound material is described in the appendix.

\section{Reduction to two dimensions and the associated operator}

As is well known from the work in 
\cite{deat92b} and \cite{deat96}, the field equations for
the first-order correction $h_{ab}^{(1)}$ in the expansion (1.1) are
particular cases of the general system given by the flat-space wave
equation (here $u \equiv {1\over \sqrt{2}}(z+t),
v \equiv {1\over \sqrt{2}}(z-t)$)
\begin{equation}
\cstok{\ }\psi=2{\partial^{2}\psi \over \partial u \partial v}
+{1\over \rho}{\partial \over \partial \rho} \left(
\rho {\partial \psi \over \partial \rho}\right)
+{1\over \rho^{2}}{\partial^{2}\psi \over \partial \phi^{2}}=0,
\label{(2.1)}
\end{equation}
supplemented by the boundary condition
\begin{equation}
\psi(u=0)={\rm e}^{im \phi}\rho^{-n}f[8 \log(v \rho)-\sqrt{2}v],
\label{(2.2)}
\end{equation}
\begin{equation}
f(x)=0 \; \; \forall x < 0.
\label{(2.3)}
\end{equation}
Moreover, $\psi$ should be of the form 
${\rm e}^{i m \phi}\rho^{-n}\chi(q,r)$
for $u \geq 0$, where
\begin{equation}
q \equiv u \rho^{-2},
\label{(2.4)}
\end{equation}
\begin{equation}
r \equiv 8 \log(\nu \rho)-\sqrt{2}v.
\label{(2.5)}
\end{equation}
For the homogeneous wave equation (2.1) there is no advantage in eliminating
$\rho$ and $\phi$ from the differential equation. However, the higher-order
metric perturbations turn out to obey inhomogeneous flat-space wave equations
of the form
\begin{equation}
\cstok{\ }\psi=S,
\label{(2.6)}
\end{equation}
where $S$ is a source term equal to 
${\rm e}^{im \phi}\rho^{-(n+2)}H(q,r)$. This
leads to the following equation for 
$\chi \equiv {\rm e}^{-im \phi}\rho^{n}\psi$:
\begin{equation}
{\cal L}_{m,n}\chi(q,r)=H(q,r),
\label{(2.7)}
\end{equation}
where ${\cal L}_{m,n}$ is an hyperbolic operator in the independent 
variables $q$ and $r$, and takes the form \cite{deat92b,deat96}
\begin{eqnarray}
{\cal L}_{m,n}&=&-(2\sqrt{2}+32q){\partial^{2}\over \partial q \partial r}
+4q^{2}{\partial^{2}\over \partial q^{2}}
+64{\partial^{2}\over \partial r^{2}} \nonumber \\
&+& 4(n+1)q{\partial \over \partial q}-16n{\partial \over \partial r}
+n^{2}-m^{2}.
\label{(2.8)}
\end{eqnarray}
The proof of hyperbolicity of ${\cal L}_{m,n}$, with the associated
normal hyperbolic form, can be found in 
section 3 of \cite{deat92b}, and in \cite{deat96}. The advantage
of studying Eq. (2.7) is twofold: to evaluate the solution at some
space-time point one has simply to integrate the product of $H$ and the
Green function $G_{m,n}$ of ${\cal L}_{m,n}$:
\begin{equation}
\chi(q,r)=\int G_{m,n}(q,r;q_{0},r_{0})H(q_{0},r_{0})dq_{0}dr_{0},
\label{(2.9)}
\end{equation}
and the resulting numerical calculation of the solution is now 
feasible \cite{deat92c,deat96}.

If one defines the variables 
\begin{equation}
X \equiv \log(q)+{r\over 4}, \; 
Y \equiv \log(q)-{r\over 4},
\label{(2.10)}
\end{equation}
the operator ${\cal L}_{m,n}$ is turned into
\begin{equation}
T_{m,n} \equiv 16{\partial^{2}\over \partial Y^{2}}
+8n{\partial \over \partial Y}+n^{2}-m^{2} 
-{1\over \sqrt{2}}{\rm e}^{-(X+Y)/2}
\left({\partial^{2}\over \partial X^{2}}
-{\partial^{2}\over \partial Y^{2}}\right).
\label{(2.11)}
\end{equation}
The operator $T_{m,n}$ is the `sum' of an elliptic operator in the
$Y$ variable and a two-dimensional wave operator `weighted' with the
exponential ${\rm e}^{-(X+Y)/2}$, which is the main source of technical
complications in these variables.

\section{Reduction to canonical form and the Riemann function}

It is therefore more convenient, in our general analysis, to reduce first
Eq. (2.7) to canonical form, and then find an integral representation 
of the solution. Reduction to canonical form means that new coordinates
$x=x(q,r)$ and $y=y(q,r)$ are introduced such that the coefficients
of ${\partial^{2}\over \partial x^{2}}$ and 
${\partial^{2}\over \partial y^{2}}$ vanish. 
As is shown in \cite{deat92b,deat96}, this is achieved if 
\begin{equation}
{\partial x \over \partial r}={\partial y \over \partial r}=1,
\label{(3.1)}
\end{equation}
\begin{equation}
{\partial x \over \partial q}={1+8q \sqrt{2}+\sqrt{1+16q\sqrt{2}} \over
2\sqrt{2}q^{2}},
\label{(3.2)}
\end{equation}
\begin{equation}
{\partial y \over \partial q}={1+8q \sqrt{2}-\sqrt{1+16q\sqrt{2}} \over
2\sqrt{2}q^{2}}.
\label{(3.3)}
\end{equation}
The resulting formulae are considerably simplified if one defines
\begin{equation}
t \equiv \sqrt{1+16q \sqrt{2}}=t(x,y).
\label{(3.4)}
\end{equation}
The dependence of $t$ on $x$ and $y$ is obtained implicitly by 
solving the system \cite{deat92b,deat96}
\begin{equation}
x=r+\log \left({{t-1}\over 2}\right)-{8\over (t-1)}-4,
\label{(3.5)}
\end{equation}
\begin{equation}
y=r+\log \left({{t+1}\over 2}\right)+{8\over (t+1)}-4.
\label{(3.6)}
\end{equation}
This leads to the equation
\begin{equation}
\log{(t-1)\over (t+1)}-{2t \over (t^{2}-1)}
={(x-y)\over 8},
\label{(3.7)}
\end{equation}
which can be cast in the form
\begin{equation}
{(t-1)\over (t+1)}{\rm e}^{2t \over (1-t^{2})}
={\rm e}^{(x-y) \over 8}.
\label{(3.8)}
\end{equation}
This suggests defining
\begin{equation}
w \equiv {(t-1)\over (t+1)},
\label{(3.9)}
\end{equation}
so that one first has to solve the transcendental equation
\begin{equation}
w{\rm e}^{(w^{2}-1)\over 2w}={\rm e}^{(x-y)\over 8},
\label{(3.10)}
\end{equation}
to obtain $w=w(x-y)$, from which one gets
\begin{equation}
t={(1+w)\over (1-w)}=t(x-y).
\label{(3.11)}
\end{equation}
On denoting by $g(w)$ the left-hand side of Eq. (3.10), one finds that,
in the plane $(w,g(w))$, the right-hand side of Eq. (3.10) is a line
parallel to the $w$-axis, which intersects $g(w)$ at no more than
one point for each value of $x-y$. For
example, when $w=1$, $g(w)$ intersects the line taking the constant value
$1$, for which $x-y=0$. The function
$$
g:w \rightarrow g(w)=w{\rm e}^{(w^{2}-1)\over 2w}
$$
is asymmetric and has the limiting behaviour described by
\begin{equation}
\lim_{w \to 0^{-}}g(w)=-\infty, \; \;
\lim_{w \to 0^{+}}g(w)=0,
\label{(3.12)}
\end{equation}
\begin{equation}
\lim_{w \to -\infty}g(w)=0, \; \; 
\lim_{w \to +\infty}g(w)=\infty .
\label{(3.13)}
\end{equation}
Thus, in the lower half-plane, $g$ has an horizontal asymptote given
by the $w$-axis, and a vertical asymptote given by the line $w=0$,
while it has no asymptotes in the upper half-plane, since
$$
\lim_{w \to \infty}{g(w)\over w}=\infty
$$
in addition to (3.13). The first derivative of $g$ reads as
\begin{equation}
g'(w)={(w+1)^{2}\over 2w}{\rm e}^{(w^{2}-1)\over 2w}.
\label{(3.14)}
\end{equation}
One therefore has $g'(w)>0$ for all $w>0$, and $g'(w)<0$ for all
$w \in (-\infty,0) - \left \{ -1 \right \}$, and $g$ is monotonically
decreasing for negative $w$ and monotonically increasing for positive $w$.
The point $w=-1$, at which $g'(w)$ vanishes, is neither a maximum nor 
a minimum point, because 
\begin{equation}
g''(w)=\left({1\over 4w^{3}}+{1\over 2w}+1+{w\over 4}\right)
{\rm e}^{(w^{2}-1)\over 2w},
\label{(3.15)}
\end{equation}
\begin{equation}
g'''(w)=\left({1\over 8w^{5}}-{3\over 4w^{4}}
+{3\over 8 w^{3}}+{3\over 8w}+{3\over 4}+{w\over 8}\right)
{\rm e}^{(w^{2}-1)\over 2w}.
\label{(3.16)}
\end{equation}
These formulae imply that $g''(-1)=0$ but $g'''(-1)=-1 \not = 0$,
and hence $w=-1$ yields a flex of $g(w)$ (see Fig. 1).

\begin{figure}
\begin{center}
\epsfig{height=10truecm,file=figure1.eps}
\caption{Plot of the function $g: w \rightarrow 
g(w)=w{\rm e}^{(w^{2}-1)/2w}$.}
\end{center}
\end{figure}

In the $(x,y)$ variables, the operator ${\cal L}_{m,n}$ therefore reads
\begin{equation}
{\cal L}_{m,n}=f(x,y){\partial^{2}\over \partial x \partial y}
+g(x,y){\partial \over \partial x}+h(x,y){\partial \over \partial y}
+n^{2}-m^{2},
\label{(3.17)}
\end{equation}
where, exploiting the formulae
\begin{equation}
{\partial x \over \partial q}={64 \sqrt{2}\over (t-1)^{2}},
\label{(3.18)}
\end{equation}
\begin{equation}
{\partial y \over \partial q}={64 \sqrt{2}\over (t+1)^{2}},
\label{(3.19)}
\end{equation}
one finds
\begin{eqnarray}
f(x,y)&=& -(2\sqrt{2}+32q)\left({\partial x \over \partial q}
+{\partial y \over \partial q}\right)
+8q^{2}{\partial x \over \partial q}{\partial y \over \partial q}
+128 \nonumber \\
&=& 256 \left[1-{2t^{2}(t^{2}+1)\over (t-1)^{2}(t+1)^{2}}\right],
\label{(3.20)}
\end{eqnarray}
\begin{equation}
g(x,y)=4(n+1)q{\partial x \over \partial q}-16n 
=16 \left[1+{2(n+1)\over (t-1)}\right],
\label{(3.21)}
\end{equation}
\begin{equation}
h(x,y)=4(n+1)q{\partial y \over \partial q}-16n 
=16 \left[1-{2(n+1)\over (t+1)}\right].
\label{(3.22)}
\end{equation}
The resulting canonical form of Eq. (2.7) is
\begin{eqnarray}
L[\chi]&=& \left({\partial^{2}\over \partial x \partial y}
+a(x,y){\partial \over \partial x}
+b(x,y){\partial \over \partial y}+c(x,y)
\right)\chi(x,y) \nonumber \\
&=& {\widetilde H}(x,y)
\label{(3.23)}
\end{eqnarray}
where
\begin{equation}
a(x,y) \equiv {g(x,y)\over f(x,y)}
={1\over 16}{(1-t)(t+1)^{2}(2n+1+t)\over (t^{4}+4t^{2}-1)},
\label{(3.24)}
\end{equation}
\begin{equation}
b(x,y) \equiv {h(x,y)\over f(x,y)}
={1\over 16}{(t+1)(t-1)^{2}(2n+1-t)\over (t^{4}+4t^{2}-1)},
\label{(3.25)}
\end{equation}
\begin{equation}
c(x,y) \equiv {n^{2}-m^{2}\over f(x,y)}
={(m^{2}-n^{2})\over 256}
{(t-1)^{2}(t+1)^{2}\over (t^{4}+4t^{2}-1)},
\label{(3.26)}
\end{equation}
\begin{equation}
{\widetilde H}(x,y) \equiv {H(x,y)\over f(x,y)}
=-{H(x,y)\over 256}
{(t-1)^{2}(t+1)^{2}\over (t^{4}+4t^{2}-1)}.
\label{(3.27)}
\end{equation}
Note that $a(-t)=b(t), b(-t)=a(t), c(-t)=c(t),
{\widetilde H}(-t)={\widetilde H}(t)$.

For an hyperbolic equation in the form (3.23), we can use the Riemann
integral representation of the solution. For this purpose, recall 
from \cite{cour61} that, on denoting by 
$L^{\dagger}$ the adjoint of the operator 
$L$ in (3.23), which acts according to
\begin{equation}
L^{\dagger}[\chi]=\chi_{xy}-(a\chi)_{x}-(b \chi)_{y}+c \chi,
\label{(3.28)}
\end{equation}
one has to find a `function' $R(x,y;\xi,\eta)$ (actually a kernel)
subject to the following conditions ($(\xi,\eta)$ being the coordinates
of a point $P$ such that characteristics through it intersect a curve
$C$ at points $A$ and $B$, $AP$ being a segment with constant $y$,
and $BP$ being a segment with constant $x$, as is shown 
in Fig. 2):

\begin{figure}
\begin{center}
\epsfig{height=10truecm,file=figure2.eps}
\caption{Geometry of the characteristic initial-value problem
in two independent variables.}
\end{center}
\end{figure}

\vskip 0.3cm
\noindent
(i) As a function of $x$ and $y$, $R$ satisfies the adjoint equation
\begin{equation}
L^{\dagger}_{(x,y)}[R]=0,
\label{(3.29)}
\end{equation}
(ii) $R_{x}=bR$ on $AP$, i.e.
\begin{equation}
R_{x}(x,y;\xi,\eta)=b(x,\eta)R(x,y;\xi,\eta) \; {\rm on} \;
y=\eta,
\label{(3.30)}
\end{equation}
and $R_{y}=aR$ on $BP$, i.e.
\begin{equation}
R_{y}(x,y;\xi,\eta)=a(\xi,y)R(x,y;\xi,\eta) \; {\rm on} \;
x=\xi,
\label{(3.31)}
\end{equation}
(iii) $R$ equals $1$ at $P$, i.e.
\begin{equation}
R(\xi,\eta;\xi,\eta)=1.
\label{(3.32)}
\end{equation}
It is then possible to express the solution of Eq. (3.23) in the form
\begin{eqnarray}
\chi(P)&=& {1\over 2}[\chi(A)R(A)+\chi(B)R(B)] 
+\int_{AB}\biggr( \left[{R\over 2}\chi_{x}
+\left(bR-{1\over 2}R_{x}\right)\chi \right]dx \nonumber \\
& - & \left[{R\over 2}\chi_{y}
+\left(aR-{1\over 2}R_{y}\right)\chi \right]dy \biggr) 
+ \int \int_{\Omega}R(x,y;\xi,\eta){\widetilde H}(x,y)dx dy,
\label{(3.33)}
\end{eqnarray}
where $\Omega$ is a domain with boundary.

Note that Eqs. (3.30) and (3.31) are ordinary differential equations
for the Riemann function $R(x,y;\xi,\eta)$ along the characteristics
parallel to the coordinate axes. By virtue of (3.32), 
their integration yields
\begin{equation}
R(x,\eta;\xi,\eta)={\rm exp} \int_{\xi}^{x}b(\lambda,\eta)d\lambda,
\label{(3.34)}
\end{equation}
\begin{equation}
R(\xi,y;\xi,\eta)={\rm exp} \int_{\eta}^{y}a(\lambda,\xi)d\lambda,
\label{(3.35)}
\end{equation}
which are the values of $R$ along the characteristics through $P$.
Equation (3.33) yields instead the solution of Eq. (3.23) for arbitrary
initial values given along an arbitrary non-characteristic curve $C$,
by means of a solution $R$ of the adjoint equation (3.29) which depends 
on $x,y$ and two parameters $\xi,\eta$. Unlike $\chi$, the Riemann
function $R$ solves a characteristic initial-value problem.

\section{Goursat problem for the Riemann function}

By fully exploiting the reduction to canonical form of Eq. (2.7)
we have considered novel features with respect to the analysis in
\cite{deat92b,deat96}, because the Riemann formula 
(3.33) also contains the integral
along the piece of curve $C$ from $A$ to $B$, and the term
${1\over 2}[\chi(A)R(A)+\chi(B)R(B)]$. This representation of the 
solution might be more appropriate for the numerical purposes
considered in \cite{deat92c}, 
but the task of finding the Riemann function $R$
remains extremely difficult. One can however use approximate methods
for solving Eq. (3.29). For this purpose, we first point out that,
by virtue of Eq. (3.28), Eq. (3.29) is a canonical hyperbolic equation
of the form
\begin{equation}
\left({\partial^{2}\over \partial x \partial y}
+A{\partial \over \partial x}+B{\partial \over \partial y}
+C \right)R(x,y;\xi,\eta)=0,
\label{(4.1)}
\end{equation}
where
\begin{equation}
A \equiv -a,
\label{(4.2)}
\end{equation}
\begin{equation}
B \equiv -b,
\label{(4.3)}
\end{equation}
\begin{equation}
C \equiv c-a_{x}-b_{y}.
\label{(4.4)}
\end{equation}
Thus, on defining
\begin{equation}
U \equiv R,
\label{(4.5)}
\end{equation}
\begin{equation}
V \equiv R_{x}+BR,
\label{(4.6)}
\end{equation}
the equation (4.1) for the Riemann function is equivalent to the
hyperbolic canonical system \cite{cour61}
\begin{equation}
U_{x}=f_{1}(x,y)U+f_{2}(x,y)V,
\label{(4.7)}
\end{equation}
\begin{equation}
V_{y}=g_{1}(x,y)U+g_{2}(x,y)V,
\label{(4.8)}
\end{equation}
where
\begin{equation}
f_{1} \equiv -B=b,
\label{(4.9)}
\end{equation}
\begin{equation}
f_{2}=1,
\label{(4.10)}
\end{equation}
\begin{equation}
g_{1} \equiv AB-C+B_{y}=ab-c+a_{x},
\label{(4.11)}
\end{equation}
\begin{equation}
g_{2} \equiv -A=a.
\label{(4.12)}
\end{equation}
For the system described by Eqs. (4.7) and (4.8) with boundary data
(3.34) and (3.35) an existence and uniqueness 
theorem holds (see \cite{cour61}
for the Lipschitz conditions on boundary data), and we can therefore
exploit the finite differences method to find approximate solutions
for the Riemann function $R(x,y;\xi,\eta)$, and eventually $\chi(P)$
with the help of the integral representation (3.33).

\section{Concluding remarks}

The inverses of hyperbolic operators 
\cite{lera53} and the Cauchy problem for
hyperbolic equations with polynomial coefficients 
\cite{lera56} have always been
the object of intensive investigation in the mathematical literature.
We have here considered the application of such issues to axisymmetric
black hole collisions at the speed of light, relying on the work in
\cite{deat92a,deat92b,deat92c,deat96}. We have pointed out that, 
for the inhomogeneous equations (2.7)
occurring in the perturbative analysis, the task of inverting the 
operator (2.8) can be accomplished with the help of the Riemann integral
representation (3.33), after solving Eq. (4.1) for the Riemann function.
One has then to solve a characteristic initial-value problem for a
homogeneous hyperbolic equation in canonical form in two independent
variables, for which we have developed formulae to be used for the
numerical solution with the help of a finite differences scheme.
For this purpose one studies the canonical system (cf (4.7) and (4.8))
\begin{equation}
U_{x}=F(x,y,U,V),
\label{(5.1)}
\end{equation}
\begin{equation}
V_{y}=G(x,y,U,V),
\label{(5.2)}
\end{equation}
in the rectangle ${\cal R} \equiv \left \{ x,y:
x \in [x_{0},x_{0}+a], y \in [y_{0},y_{0}+b] \right \}$ with known
values of $U$ on the vertical side $AD$ where $x=x_{0}$, and known
values of $V$ on the horizontal side $AB$ where $y=y_{0}$. The 
segments $AB$ and $AD$ are then divided into $m$ and $n$ equal parts,
respectively. On setting ${a\over m} \equiv h$ and 
${b \over n} \equiv k$, the original differential equations become
equations relating values of $U$ and $V$ at three intersection points
of the resulting lattice, i.e.
\begin{equation}
{U(P_{r,s+1})-U(P_{rs})\over h}=F,
\label{(5.3)}
\end{equation}
\begin{equation}
{V(P_{r+1,s})-V(P_{rs})\over k}=G.
\label{(5.4)}
\end{equation}
It is now convenient to set $U_{rs} \equiv U(P_{rs}), 
V_{rs} \equiv V(P_{rs})$, so that these equations read as
\begin{equation}
U_{r,s+1}=U_{rs}+hF(P_{rs},U_{rs},V_{rs}),
\label{(5.5)}
\end{equation}
\begin{equation}
V_{r+1,s}=V_{rs}+kG(P_{rs},U_{rs},V_{rs}).
\label{(5.6)}
\end{equation}
Thus, if both $U$ and $V$ are known at $P_{rs}$, one can evaluate
$U$ at $P_{r,s+1}$ and $V$ at $P_{r+1,s}$. The evaluation at
subsequent intersection points of the lattice goes on along horizontal
or vertical segments. In the former case, the resulting algorithm is
\begin{equation}
U_{rs}=U_{r0}+h \sum_{i=1}^{s-1}F(P_{ri},U_{ri},V_{ri}),
\label{(5.7)}
\end{equation}
\begin{equation}
V_{rs}=V_{r-1,s}+kG(P_{r-1,s},U_{r-1,s},V_{r-1,s}),
\label{(5.8)}
\end{equation}
while in the latter case one obtains the algorithm expressed by the
equations
\begin{equation}
V_{rs}=V_{0s}+\sum_{i=1}^{r-1}G(P_{is},U_{is},V_{is}),
\label{(5.9)}
\end{equation}
\begin{equation}
U_{rs}=U_{r,s-1}+hF(P_{r,s-1},U_{r,s-1},V_{r,s-1}).
\label{(5.10)}
\end{equation}
Stability of such solutions is closely linked with the geometry of the
associated characteristics, and the criteria to be fulfilled are
studied in section 13.2 of \cite{gara64} 
(stability depends crucially on whether
or not ${h\over k} \leq 1$).

To sum up, one solves numerically Eq. (3.10) for $w=w(x,y)=w(x-y)$,
from which one gets $t(x-y)$ with the help of (3.11), which is a
fractional linear transformation.
This yields $a,b,c$ and $\widetilde H$ as functions of $(x,y)$
according to (3.24)--(3.27), and hence $A,B$ and $C$ in the equation
for the Riemann function are obtained according to (4.2)--(4.4), where
derivatives with respect to $x$ and $y$ are evaluated numerically.
Eventually, the system given by (4.7) and (4.8) is solved according
to the finite-differences scheme of the present section, with
\begin{equation}
F=f_{1}U+f_{2}V=f_{1}R+f_{2}(R_{x}+BR),
\label{(5.11)}
\end{equation}
\begin{equation}
G=g_{1}U+g_{2}V=g_{1}R+g_{2}(R_{x}+BR).
\label{(5.12)}
\end{equation}
Once the Riemann function $R=U$ is obtained with the desired accuracy,
numerical evaluation of the integral (3.33) yields $\chi(P)$, and
$\chi(q,r)$ is obtained upon using Eqs. (3.5) and (3.6) for the
characteristic coordinates.

Our steps are conceptually desirable since they rely on well
established techniques for the solution of hyperbolic equations
in two independent variables \cite{cour61,gara64}, 
and provide a viable alternative
to the numerical analysis performed in \cite{deat92c}, because
all functions should be evaluated numerically. Our method is not
obviously more powerful than the one used in 
\cite{deat92a,deat92b,deat92c,deat96}, but is well
suited for a systematic and lengthy numerical analysis, while its
analytic side provides an interesting alternative for the evaluation
of Green functions both in black hole physics and in other problems
where hyperbolic operators with variable coefficients might occur.
This task remains very important because a strong production of
gravitational radiation is mainly expected in the extreme events studied
in \cite{deat92a,deat92b,deat92c,deat96} and which motivated our 
paper. Any viable way of looking 
at mathematical and numerical aspects of the problem is therefore
of physical interest for research planned in the 
years to come \cite{alle00}.

\appendix
\section{The characteristic initial-value problem}

In our expository article, we find it appropriate to include
some background material, following, for example, the presentation
in section IV of \cite{sach64}. We therefore consider some 
four-dimensional region of space-time and choose in it a set of
null hypersurfaces $u={\rm constant}$; the corresponding ray
congruence with tangent vector $k_{a}=u_{,a}$ is assumed to have
expansion $\rho \equiv {1\over 2}k_{\; ;a}^{a} \not=0$, 
which can always be arranged in a space-time patch,
whereas outside of some patch the rays start to cross and hence our
construction breaks down globally. On completing the $k^{a}$ direction
to a quasi-normal tetrad $(k,m,t)$, one finds the following split
of the vacuum Einstein equations with Einstein tensor $G_{ab}$:

Main equations (6 equations)
\begin{equation}
k^{a}G_{ab}=0, \; G_{ab}t^{a}t^{b}=0,
\label{(A1)}
\end{equation}
trivial equation
\begin{equation}
G_{ab}t^{a}{\overline t}^{b}=0,
\label{(A2)}
\end{equation}
and 3 supplementary conditions
\begin{equation}
G_{ab}m^{a}t^{b}=0, \; G_{ab}m^{a}m^{b}=0,
\label{(A3)}
\end{equation}
where a single complex equation has been counted as two real
equations. Remarkably, if the main equations hold everywhere, then
the trivial equation holds everywhere and the supplementary conditions
hold everywhere if they hold at one point on each ray. The fulfillment
of the trivial equation is proved by writing, from the vacuum
Einstein equations, that
\begin{equation}
G_{\; b;a}^{a}=0,
\label{(A4)}
\end{equation}
and then exploiting the main equations (A1) jointly with the split
of $k_{a;b}$ as given in \cite{sach64}:
\begin{eqnarray}
k_{a;b}&=& zt_{a}{\overline t}_{b}+\sigma t_{a}t_{b}
+\Omega t_{a}k_{b} +\zeta k_{a}t_{b}+{\rm c.c.} \nonumber \\
&+& \xi k_{a}k_{b}.
\label{(A5)}
\end{eqnarray}
Hence one gets
\begin{equation}
0=k^{a}G_{\; a;b}^{b}=-k_{a;b}G^{ab}=-2\rho G_{ab}
t^{a}{\overline t}^{b}.
\label{(A6)}
\end{equation}
By hypothesis the expansion $\rho$ does not vanish, so that the 
trivial equation is, indeed, identically satisfied. The fulfillment of
(A3) everywhere is proved along similar lines. Thus, one can again
integrate the main equations (A1) first and worry about the
supplementary conditions (A3) later. 

Choose now the coordinate $x^{0}$ as the retarded time: $x^{0}=u$.
Let $r=x^{1}$ be a luminosity distance along the rays; let
$x^{\alpha}$ (with $\alpha=2,3$) be any other pair of 
coordinates constant along the rays. The line element in these 
coordinates takes therefore the form (no confusion should arise
with the $\beta$ of section 1)
\begin{equation}
ds^{2}=W\; du^{2}+2{\rm e}^{2\beta} \; dudr
-r^{2}h_{\alpha \gamma}\Bigr(dx^{\alpha}-U^{\alpha}du \Bigr)
\Bigr(dx^{\gamma}-U^{\gamma} \; du \Bigr),
\label{(A7)}
\end{equation}
where $W,\beta,h_{\alpha \gamma},U^{\alpha}$ depend on the
$x^{\alpha}$ coordinates.
Since $r$ is a luminosity distance, the determinant of 
$h_{\nu \mu}$ is independent of $r$. Bearing in mind that the 
luminosity distance is defined only up to a factor constant along
each ray, one can demand without loss of generality that
(here $\theta \equiv x^{2}, \; \phi \equiv x^{3}$)
\begin{eqnarray}
2h_{\mu \nu}dx^{\mu}dx^{\nu}&=& \Bigr({\rm e}^{2\gamma}
+{\rm e}^{2\delta}\Bigr)d\theta^{2}+4 \sinh(\gamma-\delta)
d\theta d\phi \; \sin \theta \nonumber \\
&+& \sin^{2}\theta \Bigr({\rm e}^{-2\gamma}
+{\rm e}^{-2\delta}\Bigr)d\phi^{2}.
\label{(A8)}
\end{eqnarray}
The metric corresponding to the line element (A7) contains only
six unknown functions of four variables, and our coordinate system 
is `rigid' enough for our purposes \cite{sach64}.

One can either analyze the field in the neighbourhood of some point, or
the field near infinity in an asymptotically flat space-time. Indeed,
if in Minkowski space-time one uses a retarded time $u=t-r$ and
spherical coordinates $r,\theta,\phi$ one finds for the line element
\begin{equation}
ds^{2}=du^{2}+2dudr-r^{2}\Bigr(d\theta^{2}+\sin^{2}\theta \;
d\phi^{2}\Bigr),
\label{(A9)}
\end{equation}
hence one is led to require that, if asymptotic flatness holds,
\begin{equation}
\lim_{r \to \infty}W=1, \;
\lim_{r \to \infty}(rU^{\alpha})=\lim_{r \to \infty}\beta
=\lim_{r \to \infty}\delta =\lim_{r \to \infty}\gamma=0,
\label{(A10)}
\end{equation}
where all limits are taken as $r$ approaches infinity with 
$u,\theta,\phi$ fixed. The second requirement in (A10), i.e.
that all quantities of interest admit a power-series expansion
in ${1\over r}$, e.g.
\begin{equation}
(1-i)(\delta+i \gamma)/2={c(u,\theta,\phi)\over r}
+{d(u,\theta,\phi)\over r^{2}}+{\ldots} ,
\label{(A11)}
\end{equation}
is indeed restrictive. Such a requirement can be drastically
weakened but not fully eliminated; moreover, it is closely related
to an outgoing radiation condition of the Sommerfeld type. It
should be stressed that all these requirements no longer hold when
$r$ becomes small to the extent that rays start to cross each other.

In the axially- and reflection-symmetric case considered by 
Bondi et al. \cite{bond62}, one has
\begin{equation}
\delta=\gamma, \; U^{3}=0, \;
{\partial g_{ab}\over \partial \phi}=0,
\label{(A12)}
\end{equation}
and the $\phi$-direction is a hypersurface-orthogonal Killing
direction. The line element acquires the simpler form
\begin{equation}
ds^{2}={V{\rm e}^{2\beta} \over r}du^{2}
+2{\rm e}^{2\beta}\; dudr -r^{2}\Bigr[{\rm e}^{2\gamma}
(d\theta-U \; du)^{2}+{\rm e}^{-2\gamma}
\sin^{2}\theta \; d\phi^{2}\Bigr],
\label{(A13)}
\end{equation}
where the peculiar form of the first coefficient is chosen to
simplify the resulting calculations.

Interestingly, two main equations are found to be identically
satisfied by virtue of axial symmetry, whereas the other four
turn out to be linear combinations of \cite{bond62}
\begin{equation}
R_{11}=-{4\over r}\Bigr(\beta_{1}-{1\over 2}r\gamma_{1}^{2}\Bigr)=0,
\label{(A14)}
\end{equation}
\begin{equation}
-2r^{2}R_{12}= {\Bigr[r^{4}{\rm e}^{2(\gamma-\beta)}U_{1}\Bigr]}_{1}
-2r^{2}\biggr(-\gamma_{12}+2\gamma_{1}\gamma_{2}-2\gamma_{1}
\cot \theta +\beta_{12}-2{\beta_{2}\over r}\biggr)=0,
\label{(A15)}
\end{equation}
\begin{eqnarray}
\; & \; & R_{22}{\rm e}^{2(\beta-\gamma)}-r^{2}R_{3}^{3}
{\rm e}^{2\beta}=2V_{1}+{1\over 2}r^{4}
{\rm e}^{2(\gamma-\beta)}U_{1}^{2} \nonumber \\
&-& r^{2}U_{12}-4rU_{2}-r^{2}U_{1} \cot \theta 
-4r U \cot \theta \nonumber \\
&+& 2 {\rm e}^{2(\beta-\gamma)}
\Bigr[\beta_{22}+\beta_{2}^{2}-1-(3\gamma_{2}-\beta_{2})
\cot \theta -\gamma_{22} 
+2\gamma_{2}(\gamma_{2}-\beta_{2})\Bigr]=0,
\label{(A16)}
\end{eqnarray}
\begin{eqnarray}
\; & \; & -r^{2}R_{3}^{3}{\rm e}^{2\beta}= 
2r(r\gamma)_{01}+(1-r \gamma_{1})V_{1}
-(r\gamma_{11}+\gamma_{1})V-r(1-r \gamma_{1})U_{2}
\nonumber \\
&-& r^{2}(\cot \theta-\gamma_{2})U_{1}
+r(2r \gamma_{12}+2\gamma_{2}+r \gamma_{1}\cot \theta
-3 \cot \theta)U \nonumber \\
&+& {\rm e}^{2(\beta-\gamma)}\Bigr[-1
-(3 \gamma_{2}-2 \beta_{2})\cot \theta -\gamma_{22}
+2 \gamma_{2} (\gamma_{2}-\beta_{2})\Bigr]=0.
\label{(A17)}
\end{eqnarray}
Equations (A14)--(A16) are called {\it hypersurface equations}
because they contain no $u$ derivatives, while Eq. (A17) is called
the {\it standard equation}.

Now if $\gamma$ is given for one value of $u$, Eq. (A14) and the
boundary conditions (A10) determine $\beta$ uniquely. Next Eq. (A15)
and the boundary conditions determine $U$ up to a function of
integration $-6N(u,\theta)$ that can be added to 
$r^{4}{\rm e}^{2(\gamma-\beta)}U_{1}$. Equation (A16) determines $V$
up to the additive function $-2M(u,\theta)$; last, Eq. (A17)
determines $\gamma_{0}$ up to an additive function 
${c_{0}(u,\theta)\over r}$. One can then differentiate Eqs. 
(A14)--(A17) with respect to $u$ and repeat the whole procedure.
To sum up, given $\gamma$ at one moment the main equations 
determine the future or past up to the three integration functions
just mentioned. In the general case, the results are completely
similar. One has to assign at one value of $u$ the two functions
$\gamma$ and $\delta$. The future is then determined up to five
integration functions: a term $-2{M(u,\theta,\phi)\over r}$
to be added to $W$; two functions $N^{\alpha}(u,\theta,\phi)$ which
occur in the $r^{-3}$ term for $U$; and two `news functions'
$c_{0}(u,\theta,\phi)$, where the complex function $c$ is given
by Eq. (A11).

As far as the supplementary conditions are concerned, the lemma just
given makes it clear that they should only involve the functions
$M,N$ and $c_{0}$, while a long calculation yields \cite{sach64}
\begin{equation}
M_{0}=-|c_{0}|^{2}+{1\over 2}(\sin \theta)^{-1}{\rm Re}
\left \{ {\overline \nabla}\Bigr[(1/\sin \theta)
{\overline \nabla}(c_{0}\sin^{2}\theta)\Bigr] \right \},
\label{(A18)}
\end{equation}
\begin{equation}
3(N^{2}+i\sin \theta N^{3})=-\nabla M
-[4c \; \cot \theta +(\nabla c)+3c \nabla]
{\overline c}_{0},
\label{(A19)}
\end{equation}
where $\nabla \equiv {\partial \over \partial \theta}
+i(\sin \theta)^{-1}{\partial \over \partial \phi}$. The desired
$M$ and $N^{\alpha}$ can be determined once that 
$c(u,\theta,\phi)$ and some initial values are given. In the
axially symmetric case, Eqs. (A18) and (A19) take the simpler form
\cite{bond62}
\begin{equation}
M_{0}=-c_{0}^{2}+{1\over 2}(c_{22}+3c_{2} \; \cot \theta -2c)_{0},
\label{(A20)}
\end{equation}
\begin{equation}
-3N_{0}=M_{2}+3c \; c_{02}+4c c_{0} \; \cot \theta
+c_{0}c_{2}.
\label{(A21)}
\end{equation}
The functions $\gamma$ and $\delta$ given on the initial hypersurface
$u={\rm constant}$, jointly with the two news functions $c_{0}$ given
at $r=\infty$ describe the two transverse degrees of freedom. Moreover,
one should specify $M$ and $N$ at the initial (or final) retarded time,
and these three functions of two variables must be related to the
longitudinal-timelike degrees of freedom of the gravitational field.
In the characteristic value problem for general relativity, the 
independent data appear therefore in a very explicit form \cite{sach64}.

\acknowledgments
The work of G. Esposito has been partially supported by PRIN
{\it SINTESI}. He is grateful to Decio Cocolicchio and Sorin
Dragomir for warm encouragement. IOP Publishing has kindly granted
permission to republish material from the author's paper in
\cite{espo01}.

\begin{references}
\bibitem{sach64}
R.K. Sachs, in {\it Relativity, Groups and Topology}, edited
by C. DeWitt and B.S. DeWitt (Gordon \& Breach, New York, 1964).
\bibitem{dewi65}
B.S. DeWitt, {\it Dynamical Theory of Groups and Fields}
(Gordon and Breach, New York, 1965).
\bibitem{dewi84}
B.S. DeWitt, in {\it Relativity, Groups and Topology II}, edited 
by B.S. DeWitt and R. Stora (North--Holland, Amsterdam, 1984).
\bibitem{espo01a}
G. Esposito, {\it Quantum Gravity in Four Dimensions}
(Nova Science, New York, 2001).
\bibitem{cart00}
P. Cartier and C. DeWitt--Morette, {\it J. Math. Phys.}
{\bf 41}, 4154 (2000).
\bibitem{deat92a}
P.D. D'Eath and P.N. Payne, {\it Phys. Rev.} D {\bf 46}, 658 (1992).
\bibitem{deat92b}
P.D. D'Eath and P.N. Payne, {\it Phys. Rev.} D {\bf 46}, 675 (1992).
\bibitem{deat92c}
P.D. D'Eath and P.N. Payne, {\it Phys. Rev.} D {\bf 46}, 694 (1992).
\bibitem{deat96}
P.D. D'Eath, {\it Black Holes: Gravitational Interactions}
(Clarendon Press, Oxford, 1996).
\bibitem{espo00a}
G. Esposito and C. Stornaiolo, {\it Class. Quantum Grav.}
{\bf 17}, 1989 (2000).
\bibitem{espo00b}
G. Esposito and C. Stornaiolo, {\it Found. Phys. Lett.}
{\bf 13}, 279 (2000).
\bibitem{espo01}
G. Esposito, {\it Class. Quantum Grav.} {\bf 18},
1997 (2001).
\bibitem{cour61}
R. Courant and D. Hilbert, {\it Methods of Mathematical
Physics. II. Partial Differential Equations} 
(Interscience, New York, 1961).
\bibitem{lera53}
J. Leray, {\it Hyperbolic Differential Equations}
(Princeton University Press, Princeton, 1953).
\bibitem{lera56}
J. Leray, {\it C. R. Acad. Sci. Paris} {\bf 242}, 953 (1956).
\bibitem{gara64}
P.R. Garabedian, {\it Partial Differential Equations}
(Chelsea, New York, 1964).
\bibitem{alle00}
B. Allen and A. Ottewill, {\it Gen. Rel. Grav.} 
{\bf 32}, 385 (2000).
\bibitem{bond62}
H. Bondi, M.G.J. van der Burg, A.W.K. Metzner, 
{\it Proc. Roy. Soc. Lond.} 
{\bf A 269}, 21 (1962).
\end{references}
\end{document}


astro-ph0504599.extracted_bib
\bibitem{sach64}
R.K. Sachs, in {\it Relativity, Groups and Topology}, edited
by C. DeWitt and B.S. DeWitt (Gordon \& Breach, New York, 1964).

\bibitem{dewi65}
B.S. DeWitt, {\it Dynamical Theory of Groups and Fields}
(Gordon and Breach, New York, 1965).

\bibitem{dewi84}
B.S. DeWitt, in {\it Relativity, Groups and Topology II}, edited 
by B.S. DeWitt and R. Stora (North--Holland, Amsterdam, 1984).

\bibitem{espo01a}
G. Esposito, {\it Quantum Gravity in Four Dimensions}
(Nova Science, New York, 2001).

\bibitem{cart00}
P. Cartier and C. DeWitt--Morette, {\it J. Math. Phys.}
{\bf 41}, 4154 (2000).

\bibitem{deat92a}
P.D. D'Eath and P.N. Payne, {\it Phys. Rev.} D {\bf 46}, 658 (1992).

\bibitem{deat92b}
P.D. D'Eath and P.N. Payne, {\it Phys. Rev.} D {\bf 46}, 675 (1992).

\bibitem{deat92c}
P.D. D'Eath and P.N. Payne, {\it Phys. Rev.} D {\bf 46}, 694 (1992).

\bibitem{deat96}
P.D. D'Eath, {\it Black Holes: Gravitational Interactions}
(Clarendon Press, Oxford, 1996).

\bibitem{espo00a}
G. Esposito and C. Stornaiolo, {\it Class. Quantum Grav.}
{\bf 17}, 1989 (2000).

\bibitem{espo00b}
G. Esposito and C. Stornaiolo, {\it Found. Phys. Lett.}
{\bf 13}, 279 (2000).

\bibitem{espo01}
G. Esposito, {\it Class. Quantum Grav.} {\bf 18},
1997 (2001).

\bibitem{cour61}
R. Courant and D. Hilbert, {\it Methods of Mathematical
Physics. II. Partial Differential Equations} 
(Interscience, New York, 1961).

\bibitem{lera53}
J. Leray, {\it Hyperbolic Differential Equations}
(Princeton University Press, Princeton, 1953).

\bibitem{lera56}
J. Leray, {\it C. R. Acad. Sci. Paris} {\bf 242}, 953 (1956).

\bibitem{gara64}
P.R. Garabedian, {\it Partial Differential Equations}
(Chelsea, New York, 1964).

\bibitem{alle00}
B. Allen and A. Ottewill, {\it Gen. Rel. Grav.} 
{\bf 32}, 385 (2000).

\bibitem{bond62}
H. Bondi, M.G.J. van der Burg, A.W.K. Metzner, 
{\it Proc. Roy. Soc. Lond.} 
{\bf A 269}, 21 (1962).
