\tolerance = 10000
\documentstyle[preprint,revtex]{aps}
\tightenlines
% Macro for lambda bar
\def\hlambda{{\mathchar'26\mkern-9mu\lambda}}
%
\begin{document}
\draft
\begin{title}
Numerically generated black hole spacetimes: \\ 
Interaction with gravitational waves
\end{title}
\medskip
\centerline{Andrew Abrahams$^{(1,2)}$, David Bernstein$^{(1)}$, 
David Hobill$^{(1,3)}$, 
Edward Seidel$^{(1)}$, and Larry Smarr$^{(1)}$}
\medskip
\begin{instit}
$^{(1)}$ National Center for Supercomputing Applications \\
605 E. Springfield Ave., Champaign, IL. 61820
\end{instit}
\smallskip
\begin{instit}
$^{(2)}$ Center for Radiophysics and Space Research \\
Cornell University, Ithaca, NY 14853
\end{instit}
\smallskip
\begin{instit}
$^{(3)}$ Department of Physics and Astronomy \\
University of Calgary, Calgary, Alberta, Canada T2N 1N4
\end{instit}
\begin{abstract}
In this paper we present results from a new two-dimensional
numerical relativity code used to study the 
interaction of gravitational waves with a 
black hole.  The
initial data correspond to a single black hole superimposed with
time-symmetric 
gravitational waves (Brill waves). 
A gauge invariant method is presented for extracting the gravitational
waves from the numerically generated
spacetime.  We show that the interaction between the gravitational
wave and the black hole excites the quasinormal modes of
the black hole.  An extensive comparison of these results is made to 
black hole perturbation theory.  For low amplitude initial 
gravitational waves,
we find excellent agreement between the theoretically predicted $\ell=2$ 
and $\ell=4$ waveforms and the waveforms generated by the code.  
Additionally, a code test is performed wherein the propagation of
the wave on the black hole background
is compared to the evolution predicted by perturbation theory.
\end{abstract}

\pacs{PACS \#'s 04.30.+x,  95.30.Sf}

%\begin{narrowtext}

\section{INTRODUCTION}\label{intro}

The study of black holes and their interactions 
is of great importance 
to physics and astronomy.  To the physicist black holes represent a rich 
laboratory for exploring general relativity.  Dynamic spacetimes 
containing black
holes will generally have strong time dependent gravitational fields,
complicated horizon structures, and gravitational waves.
To the astronomer the formation [1] or collision [2] of black holes may
provide one of the strongest 
sources of gravitational waves.  Coalescing black hole binary 
systems will emit radiation as the orbit decays; eventually 
the black holes will plunge together emitting a burst of gravitational
waves.
Radiation from such systems should be observable by LIGO [3].
Unfortunately, due to the complex nonlinear structure of Einstein's
equations, 
no general analytic techniques are available for obtaining strong-field,
highly dynamical solutions in the absence of spatial symmetries.
In such situations, the only recourse is numerical solution of the
Einstein equations.  In this study we employ the so-called 3+1 formalism 
developed by Arnowitt, Deser, and Misner[4] as a practical method for 
computer simulations of realistic relativistic systems.

We have recently developed an accurate, two-dimensional numerical
relativity code which solves the fully non-linear vacuum Einstein
equations.  This code is designed to compute axisymmetric spacetimes
consisting of black holes and gravitational waves.
To verify results, two completely 
independent codes that can be cross checked have been written 
(here we simply refer to them as ``the code").
The code has been thoroughly tested against results from a one 
dimensional code designed to evolve spherically symmetric spacetimes, and
it has 
also been tested for accuracy and stability in nonspherical spacetimes.
The code has also been written in such a way that different spatial gauges,

time slicing conditions, and  numerical evolution schemes can easily be 
implemented.  These features and tests are 
described in detail in Ref. [5].

In this paper we present an application of this numerical code to a 
spacetime consisting of a single black hole with a low amplitude 
gravitational wave present.  
The criterion by which we define the amplitude of the
wave is that the mass of the system measured by the area of the apparent
horizon
on the initial time slice is typically within about 10\% of the mass 
measured far from the black hole. This definition will be made more precise
 in Sec. IV.
In this regime the gravitational wave is not necessarily considered to be a

negligible perturbation on the black hole, but it is small enough that one 
would expect the evolution of the spacetime to be dominated by the
black hole.  One of the aims of the present paper is 
to examine the range of validity of black hole perturbation theory 
and to begin to explore the regime where nonlinear effects are seen. 
For certain sufficiently large perturbations of a single black hole, we
expect the evolution to resemble the late stages of a two black hole
spacetime 
after coalescence has occurred. 
Problems involving large amplitude waves, which may contain significantly
more 
mass than the black hole alone, as well as simulations of two black holes, 
will be considered in future papers.

In what follows we discuss the code and its application to black 
hole-gravitational wave interactions.  In Sec.\ref{code} we briefly
summarize 
the relevant mathematical and computational aspects of the 
code.  A discussion of the gauge-invariant 
extraction method used for calculating waveforms and for
code diagnostics then follows in Sec.\ref{extraction}.  The results from
numerical simulations of black hole-gravitational wave interactions are
presented in Sec.\ref{interaction} which is divided into the following
subsections:
Sec.\ref{excitation} which considers low amplitude perturbations and 
quasi-normal mode stimulation, Sec.\ref{perturbed} which 
presents a stringent code test where the full two-dimensional 
nonlinear code is compared with perturbation theory, and Sec.\ref{stronger}
which discusses some results in the regime where the initial 
gravitational wave has a moderate amplitude, so that nonperturbative 
effects begin to appear. 
In Sec.\ref{conclusions}
we summarize our results and discuss future work to be done with this code.
 
Finally, we provide some details of the waveform extraction method 
in the Appendix.

\section{CODE SUMMARY}\label{code}

Here we summarize the equations solved and the methods used
in the calculations.  Our code describes an 
axisymmetric, nonrotating, equatorial plane symmetric
spacetime with the metric written as
\begin{equation}
\label{comp metric}
g_{\alpha \beta}=\left(
\begin{array}{cc}
-\alpha^{2}+\beta^{i}\beta_{i} & 
\begin{array}{ccc}
\hspace{0.0in} \beta_{\eta} & \hspace{-0.04in} \beta_{\theta} & 
\hspace{0.16in} 0
\end{array}
\\
\begin{array}{c} \beta_{\eta} \\ \beta_{\theta} \\ 0 \end{array} 
&
\Psi^{4} 
\left(
\begin{array}{ccc}
A & C & 0 \\
C & B & 0\\
0 & 0 & D \sin^{2}\theta
\end{array} \right)
\end{array} \right)
\end{equation}
where $A,B,C,D,\alpha,\beta^{\eta},$ and $\beta^{\theta}$ are 
functions of our coordinates $\eta,$ $\theta,$ and $t$ ($\Psi$ is a 
function of $\eta $ and $\theta $ only).  We
use Greek indices to label 4-dimensional objects and latin indices
to label 3-dimensional objects on a given spacelike slice.  We use
units such that $G=c=1$ with $m$ our unit of length (time). 
Here $\eta$ is 
a logarithmic radial coordinate, chosen so that $\eta = 0$ labels the 
position of the isometry surface (or throat) of the black 
hole, $\theta$ is the usual angular polar
coordinate, and the spacetime is independent of the azimuthal coordinate
$\phi$
(there is an axial Killing vector $\partial/ \partial \phi$).  We write 
the extrinsic curvature in conformal form
\begin{equation}
K_{ab}=\Psi^{4} \hat{K}_{ab};
\end{equation}
 The conformal extrinsic curvature takes the form
\begin{equation}
\label{comp excurv}
\hat{K}_{ab}= \left(
\begin{array}{ccc}
H_a & H_c & 0 \\
H_c & H_b & 0 \\
0 &   0   & H_d \sin^2\theta 
\end{array} \right),
\end{equation}
and the 3-metric $\gamma_{ab}$ is given by 
\begin{equation}
\gamma_{ab} = \Psi^4 \left(
\begin{array}{ccc}
A & C & 0 \\
C & B & 0\\
0 & 0 & D \sin^{2}\theta
\end{array} \right).
\end{equation}
The 3-metric components $A, B, C$, and $D$, and their 
corresponding 
extrinsic curvature components are evolved according to the 3+1 
Einstein 
equations:
\begin{equation}
\label{metric evolution}
\partial_{t}\gamma_{ab}=-2\alpha K_{ab}+ D_{a}\beta_{b}+
D_{b}\beta_{a},
\end{equation}
and
\begin{equation}
\label{excurv evolution}
\partial_{t}K_{ab}=-D_{a}D_{b}\alpha+\alpha \left[ 
R_{ab}+(tr K)K_{ab}-2K_{ac}K^{c}{}_{b} \right] 
+\beta^{c}D_{c}K_{ab}+
K_{ac}D_{b}\beta^{c}+K_{cb}D_{a}\beta^{c},
\end{equation}
where $D_{a}$ is the three dimensional covariant derivative on a 
spacelike 
slice and $R_{ab}$ is the three dimensional Ricci tensor on that slice.  
These equations are expanded out into coordinate-component form and solved
as 
described in Ref [5].  

Due to the coordinate invariance of the Einstein equations, 
the lapse $\alpha$ and the shift vector $\beta^a$ can be chosen
arbitrarily as these quantities are related to our choice of
the spatial and time coordinates.
For this work, the shift vector $\beta^{a}$ 
is chosen so 
that the off diagonal component $C$ of the three metric vanishes.  
This has been found to be crucial in maintaining 
stability on the symmetry axis 
defined by $\theta = 0$.  The lapse $\alpha$ is determined
by using the maximal slicing condition  
\begin{equation}
tr K = 0.
\end{equation}

As described in Ref. [5], our spacetime consists of an 
Einstein-Rosen bridge [6] connecting two asymptotically flat
sheets.  We impose an isometry between the two sheets of the black hole
spacetime which yields  boundary conditions at the throat ($\eta = 0$).
Boundary conditions on the symmetry axis ($\theta = 0$) and
on the equator ($\theta = \pi /2$) are determined by symmetry. 
For the results presented in this paper we have placed the outer boundary
far enough from the black hole that gravitational waves do not reach the
edge of the grid.  Since at large distances from the black hole our radial 
coordinate $\eta$ is logarithmic
in the Schwarzschild radius, this procedure is not expensive in terms of
computational resources.  However, as we will discuss later, it
causes some problems with wave propagation at large radii.
We have recently experimented with an outgoing 
wave condition which allows us to move the outer boundary in 
closer to the black hole.  This condition, which improves the late time
behavior of the code, will be discussed in a future paper.

Initial data are obtained by assuming a three metric of the form:
\begin{equation}
ds^{2}=\Psi^{4} \left[ e^{2q} ( d \eta^{2} + d \theta^{2} )
+ \sin^2 \theta d \phi^{2} \right],
\end{equation}
where the function $q(\eta,\theta)$ must satisfy certain boundary 
conditions, but is otherwise arbitrary.  This form of the three metric 
was 
originally considered by Brill [7] in his study of time symmetric 
gravitational waves.  In our case, for $q = 0$, the hamiltonian 
constraint 
equation for $\Psi$ has a Schwarzschild solution (in isotropic
coordinates),
and for $q$ nonzero the solution corresponds to the superposition of a 
gravitational wave and a single black hole.  To facilitate calculation
of initial data, we assume time-symmetry, $K_{ab} = 0$, and solve
only the hamiltonian constraint.
We choose to take $q(\eta,\theta)$ in the form:
\begin{equation}
\label{Init eq}
q=A f(\theta) \left( \exp\left[-\left(\frac{\eta + 
\eta_0}{\sigma}\right)^{2}\right]+
\exp\left[-\left(\frac{\eta - \eta_0}{\sigma}\right)^{2}\right]\right).
\end{equation}

The function $q$ has three independent parameters $A,\eta_0$, and $\sigma$
which 
specify its amplitude, range, and width.  Its angular dependence is given
by the 
function $f,$  which we usually choose to be sin$^n \theta$, where 
$n$ is some even integer.  For a full discussion of code coordinate
conditions, boundary conditions and evolution schemes see Ref [5].

Given initial data representable by the parameters $ f(\theta), A, \eta_0$,
and $\sigma$,
the evolution equations, coupled with our choice of lapse $\alpha$
and shift vector $\beta ^i$, are used to advance the system forward in
time. 
Our code has been written so
that the numerical evolution scheme can be easily switched between the
MacCormack, Brailovskaya, and extrapolated leapfrog schemes [5].  The
results in this paper are obtained from the
leapfrog scheme.  

\section{WAVEFORM EXTRACTION}\label{extraction} 

A general problem associated with the extraction of physically relevant 
information in numerical relativity is that the dynamical variables
employed
in computer simulations
are generally tensor components, which are coordinate dependent, 
whereas one 
is primarily interested in calculating physically 
meaningful quantities such as the gravitational wave amplitudes 
$h_+$ and $h_ \times$ or the total energy radiated by the system. 
A further complication arises because one wishes to consider the full 
nonlinear evolution of the Einstein equations for strongly relativistic 
problems, yet the theory of gravitational radiation is best understood
in the context of linear perturbations of Minkowski space.  
One needs to identify quantities computed from the evolved metric 
functions
at finite radii with solutions for linearized 
gravitational waves in the radiation zone. 

The problem of radiation ``extraction" has been treated 
in number of ways in past studies.  One approach 
to extract gravitational wave information, developed by Smarr [2],
involved the calculation of a frequency averaged Bel-Robinson 
vector.  This was motivated by its 
formal similarity with the Poynting vector of electromagnetism. 
We will attempt to
evaluate the efficacy of the Bel-Robinson vector and other curvature based
quantities as
radiation indicators in a future paper.  Other
possible techniques have been proposed but not explored in detail
in full general relativistic simulations. Schutz [8] 
derived a method for calculating the quadrupole moment from projections of
Riemann tensor components in the near-zone.  Anderson and Hobill [9]
found a possible technique for matching an analytic solution on
null cones to an interior numerical solution on spacelike slices
enabling calculation of the Bondi ``news'' function.
One of the most practically useful advances
[10] was the construction of
metric variables in certain spatial gauges (quasi-isotropic and
radial) which asymptotically approach the gravitational wave 
quantities $h_+$ and $h_\times$.  Unfortunately, it is not always
possible to devise useful radiative variables for all gauge choices.
Also, by themselves, the radiative variables are not adequate for
radiation extraction in many circumstances; due to the finite 
spatial extent of the numerical grid these variables always contain
near-zone 
information and gauge effects.  

In this work we will use a gauge-invariant waveform extraction 
technique that can, 
in principle, be applied to any spacetime with a isolated wave source.
The method used in this paper, a slight variant of one of the 
techniques detailed by Abrahams and Evans [11], makes use of the fact that
the 
spacetime  
surrounding a perturbed star or black hole is basically spherical.  If the
energy 
carried by the gravitational waves is small compared with the mass 
of the 
star or black hole, as will usually be the case, the radiating system 
can be 
considered as one where gravitational wave perturbations are propagating on
a 
spherical 
spacetime background.  Thus, while the numerical code evolves the 
fully 
nonlinear spacetime, when extracting information from the 
numerically 
generated metric we may examine the system from the 
point of view of perturbation theory of spherical spacetimes [12].  
This perturbation theory can be applied to any 
spacetime dominated by a static monopole moment, 
whether it is a vacuum spacetime or one containing 
matter.  It is, of course, necessary to perform
the waveform extraction in a vacuum region of the spacetime. 
In the present case the entire spacetime is in vacuum. 

We assume that the numerically calculated spacetime metric
can be split as 
\begin {equation}
g_{\alpha \beta}  =   
\stackrel{o}{g}_{\alpha \beta} +  h_{\alpha \beta},
\end {equation} 
where  $\stackrel{o}{g}_{\alpha \beta}$ is spherically symmetric and 
diagonal:
\begin{equation}
\label{Background g eq}
\stackrel{o}{g}_{\alpha \beta}   =\left( \begin{array}{cccc}
-N^2 & 0 & 0 & 0 \\
0 & A^2 & 0 & 0 \\
0 & 0 & R^2 & 0 \\
0 & 0 & 0 & R^2 sin^2 \theta
\end {array} \right)
\end {equation}
with $N = N(t,r)$, $A = A(t,r)$, and $R = R(t,r)$, and where $t$ and 
$r$ are
the time 
and radial coordinates.  (Here $r$ may be any radial coordinate.)  The most
general spherical background 
metric 
would have a nonzero $\stackrel{o}{g}_{tr}$
component.  Although such a shift vector is present in our 
gauge, this term vanishes for a spherically symmetric 
spacetime; any contribution to this term is considered a nonspherical 
perturbation from the point of view of the waveform extraction 
method.  
In principle spherical contributions to this this term can be included 
as 
well.  We assume that we are dealing with axisymmetric, non-rotating,
systems which leads to a general  even-parity metric perturbation 
tensor of the form:
\begin{equation}
\label{Perturb g eq}
h_{\alpha\beta}=\left( \begin{array}{cccc}
N^2H_0Y_{\ell 0} & H_1Y_{\ell 0} & h_0 Y_{\ell 0,_\theta} & 0  \\
H_1Y_{\ell 0} & A^2H_2Y_{\ell 0} & h_1 Y_{\ell 0,_\theta}& 0 \\
h_0 Y_{\ell 0,_\theta} & h_1 Y_{\ell 0,_\theta} & R^2 \left(K + G 
\frac{\partial^2}{\partial \theta^2}\right)Y_{\ell 0} & 0  \\
0 & 0 & 0 & R^2 \left( K sin^2\theta + G sin\theta cos\theta 
\frac{\partial}{\partial \theta}\right) Y_{\ell 0}
\end {array} \right)
\end {equation}
where the Regge-Wheeler perturbation functions $(H_0, H_1, H_2, 
h_0, 
h_1, K, 
G)$ are all functions of the radial and time coordinates $(r,t)$ only, 
and the $Y_{\ell 0} = Y_{\ell 0}(\theta)$ are the usual spherical 
harmonic 
functions of order $\ell, m$, with $m = 0$.  For this study
of even-parity radiation, we are restricted 
to even values of $l$ by equatorial plane symmetry.  

Our first step will be to extract from the numerically 
generated metric $g_{\alpha \beta}$ the perturbed metric 
quantities $(H_0, H_1, H_2, h_0, h_1, K, G)$, 
as these represent the nonspherical part of the metric, and hence
contain all the information on any  gravitational waves present in the 
system.  
This extraction is done by projecting the full numerical metric 
onto the perturbed metric by performing integrals over a coordinate
two-sphere surrounding the source. The precise prescription 
will be given in the Appendix.  As these are metric quantities, 
dependent on gauge choice, they do not directly yield the actual 
wave perturbation of the spherical spacetime.
However, following the work of Moncrief and Gerlach and Sengupta [12], 
we can construct quantities from linear combinations of these metric 
coefficients and their derivatives which are 
gauge-invariant in the sense that 
they are invariant under general infinitesimal coordinate transformations 
which leave the background coordinate system fixed.  In the case 
where the background coordinate system is Schwarzschild, this 
construction of gauge-invariant quantities 
immediately gives a function $\psi$ 
which obeys the Zerilli equation describing the propagation of 
even-parity gravitational waves on a Schwarzschild background.  
If the background coordinate system is time dependent (e.g., a 
spherical vacuum spacetime sliced in an arbitrary way), 
generalizations of the Zerilli 
function and Zerilli wave equation 
can be found.  However, in our numerical code, the coordinate 
conditions 
are such that the spherical part of the metric closely approximates 
the Schwarzschild metric in the region where the waveform extraction is 
done.   (See Sec. \ref{perturbed}).  
In the analysis that follows we assume a static
Schwarzschild background metric.  In this case the areal radius
$R(t,r)$ in Eq. (\ref{Background g eq}) and (\ref{Perturb g eq}) 
is the same as the radial coordinate $r$, easily calculated from
the code coordinates.
The full set of gauge-invariant perturbation
equations on a general time dependent spherical background 
have also been worked out but will not be presented 
here.  See Gerlach and Sengupta or Seidel [12]
for details.    

The Zerilli function $\psi$ can be defined in the 
following way:  Once the Regge-Wheeler perturbation functions 
$(H_0, H_1, H_2, h_0, h_1, K, G)$ have been extracted from the 
metric (see the Appendix) one follows Moncrief [12] and  
defines the following quantities:
\begin {equation}
k_1  \equiv  K  +  Sr G,_r  -  2\frac{S}{r} h_1 ,
\end {equation}
\begin {equation}
k_2  \equiv  \frac{H_2} {2S}  -  \frac{1}{2S^{\frac{1}{2}}}
 \frac{\partial}{\partial r}\left(rS^{-\frac{1}{2}}K\right),
\end {equation}
where $S = 1 - 2M/r$.  Then the quantity defined by
\begin {equation}
\label{psi eq}
\psi  \equiv \sqrt{\frac{2(\ell -1)(\ell +2)}{\ell (\ell
+1)}}\frac{\left(4rS^2k_2+
\ell (\ell + 1) rk_1\right)}{\Lambda},
\end {equation}
where
\begin{equation}
\Lambda \equiv \left(\ell (\ell +1) -2+\frac{6M}{r}\right)
\end {equation}
satisfies the Zerilli equation, valid for any $\ell$ value, 
which is
\begin {equation}
\label{Zerilli eq}
\frac{{\partial}^2}{\partial t^2}\psi  -  \frac{{\partial}^2}{\partial 
r_*^2}
\psi  + 
V_Z(r)  \psi = 0.
\end {equation}
The gravitational scattering potential is
\begin {equation}
V_Z(r)  =   \left(1-\frac{2M}{r}\right) 
\left(\frac{1}{\Lambda^{2}}
\left(\frac{72M^3}{r^5} - \frac{12M}{r^3} (\ell-1) (\ell+2) \left(1-
\frac{3M}{r}\right)\right)
+
\frac{\ell (\ell-1) (\ell+2) (\ell+1)}
{r^2 \Lambda}\right),
\end {equation}
and $r_*$ is the tortoise coordinate defined by 
$r_* = r + 2M \ln\left(r/2M - 1\right)$, $t$ is the 
Schwarzschild time coordinate, and $r$ is the areal radius.  With our
choice of
normalization given by Eq. (\ref{psi eq}), the asymptotic energy flux
carried
by the gravitational wave is given by 
\begin {equation}
\frac{dE}{dt} = \frac{1}{32 \pi} \left({\frac{\partial \psi}{\partial
t}}\right)^2
\end {equation}
for all $\ell$ modes.  With this scaling, one can compare plots 
of $\psi$ for different runs and 
different $\ell$ modes and see directly the relative strengths of the
signals.

The function $\psi$ is the fundamental 
radiative quantity that can be extracted from the numerically 
generated metric 
$g_{\alpha \beta}$.  It represents the generally
small wave-like part of the metric, separated from the  
numerically dominant roughly static mass monopole part, that
will be radiative at large radii.
In addition to providing a clear way of tracking the
the propagation of gravitational waves, it can provide important 
tests of the accuracy of the code.  Since its evolution equation 
(\ref{Zerilli eq}) is known, it can be used as a check
on the code's ability to correctly propagate the waves as a perturbation
on a curved background.  Also, since it is the same function
ordinarily used in the analytic calculation of black hole normal 
mode frequencies [13], one can extract the oscillation frequencies 
of the black hole determined from $\psi$ and compare directly 
with analytically known results.  These procedures will be 
discussed in the next section. 

It is also possible to use $\psi$ calculated as a time-series at
a finite radius to eliminate linear near-zone field contamination
of the waveform and calculate a waveform valid in the local 
wavezone.  This is done [11]
by solving (\ref{Zerilli eq}) approximately in a weak-field region
for a given $l-$mode as a functional of the appropriate $l-$pole
moment and its time derivatives.  Solution of an ODE yields 
the asymptotically radiative ($r^{-1}$ dependent) part of the field.
Although we could perform this separation in the waveforms shown in
the next section, it is not of major consequence for a radiating
black hole because there is no weak field near-zone.  See discussion
in Sec.\ref{excitation}.

\section{GRAVITATIONAL WAVE - BLACK HOLE 
INTERACTIONS}\label{interaction}

Most previous studies concerned with actually solving the Einstein
equations for black hole spacetimes involve 
perturbation theory.  Generally, an analytic 
black hole solution such as Schwarzschild or Kerr is studied, and 
perturbations are explicitly added, either directly into the spacetime 
metric or with source terms which are considered small so the stress energy

contribution to the spacetime is neglected.  
The original aim was to examine the stability of black 
holes under small perturbations. Although it was relatively
easy to establish stability for Schwarzschild (see, for example Moncrief
[12]), 
it was not until recently
that rotating black holes were shown to be stable [14].  Through
these studies it was discovered that black holes also have a characteristic
set
of normal mode frequencies under infinitesimal perturbations, where
gravitational
waves are emitted as the black hole returns to a spherical, unperturbed
state.  Subsequently a number of studies[1,15] based on this 
perturbation formalism detailed the process of exciting the
black hole normal modes by orbiting particles, plane waves, or 
stellar collapse.


In the present study we are able to treat essentially arbitrary
amplitude Brill waves interacting with a black hole.
For low amplitude waves, we should be able to reproduce
the results known from black hole perturbation theory by
analyzing the spacetime using the techniques discussed above.
For larger amplitude waves, we can go beyond perturbation
theory and look for nonlinear effects.
For example, at linear order in perturbation theory the 
$\ell$-modes 
are completely decoupled, so that they may be treated separately.
At higher order the different
modes will mix.  This phenomenon will be explored in a future paper. 

\subsection{Excitation of Normal Modes}\label{excitation}

As discussed above the physics of a slightly 
nonspherical spacetime around a black hole can be described by a 
single wave equation, e.g., the Zerilli equation (\ref{Zerilli eq}).
This rather simple looking equation describes a 
rich physical system, with a  wide  spectrum of normal modes and wave 
scattering properties.  This physical information is encoded into the 
function $V(r)$, known as the curvature potential, which has
a maximum near $r = 3M$.
This potential barrier partially transmits and partially reflects
incident waves.  For waves of 
certain frequencies, a resonance occurs where the transmission coefficient
becomes singular.  This condition gives rise to a normal mode of the 
black hole, and is equivalent to waves ``originating" at the potential peak

and moving away from the black hole on one side, and down the black 
hole on the other side.  If the Zerilli function is assumed to have a
definite 
frequency $\omega$, then the Zerilli equation can be regarded as a 
time-independent Schr\"odinger equation with energy $\omega^2$.  The 
calculation of these quasinormal mode frequencies can be viewed as a 
solution to the energy eigenvalues of the Schr\"odinger equation when
outgoing 
boundary conditions are imposed.  For a nonrotating black hole, at each
$\ell$ 
value there are presumed to be an infinite number of normal mode
frequencies 
$\omega_n$, and these frequencies, which depend on the mass of 
the black hole, are independent of the azimuthal mode $m$.  
(The degeneracy is split if the black hole is rotating.)

If a black hole is perturbed in an arbitrary way, one expects 
its normal modes to be excited.  This effect has been seen in a number 
of studies over many years, where normal mode excitation is seen from 
perturbations due to processes ranging from infalling or orbiting particles
 
to nonspherical collapse of a star to a black hole [1, 15]. In the present
study 
we perturb the hole with a gravitational wave. 
This process should then produce a waveform far 
from the black hole with two distinct regimes:  1)  There will be a 
precursor waveform which contains a mixture of information from a 
number of processes:  a)  The initial distortion of the hole, which is 
related to the shape and amplitude of the incident gravitational wave, b)  
the reflection properties of the potential, and c) the normal mode 
excitation process itself;  2)  After the precursor emerges, 
one will be left with radiation composed of the spectrum of normal mode 
frequencies of the black hole as it rings down.  One might expect this to
be a 
complicated shape since in principle it may contain a superposition of a
large 
(possibly infinite) number of the normal mode frequencies.  
However, the fundamental frequency $\omega_0$ is significantly less damped
(i.e., 
it has a smaller imaginary part) than the higher overtone frequencies 
$\omega_n$, 
so any contribution from higher modes 
will be quickly damped out [13].  Therefore, one expects to see a precursor
waveform 
containing details of the perturbation itself followed by a waveform
dominated 
by the fundamental black hole normal mode frequency.


With the above considerations in mind, we present numerical results for
gravitational waves interacting with a single black hole. 
We have chosen a number of representative runs which
explore three different gravitational wave-black hole interactions (see 
Table 1).  For the first set of runs, labelled Run 1a-1e, we held the shape
and
position of the wave fixed while changing only its amplitude.  Cases 1a-1d
represent fairly low amplitude waves, and will be discussed in this
section, while case 1e is a much higher amplitude wave, and will
be discussed in Sec. \ref{stronger}.  In Run 2 we chose a much broader
initial Brill wave profile, while for Run 3 we chose to place a narrow
wave packet much farther from the hole, in the spirit of
simulating a particle falling into a black hole.  For each of these runs we
list in Table 1 the initial Brill wave parameters: the amplitude $A$, the
width $\sigma$, 
the radial position $\eta_0$, and the angular dependence parameter $n$
(see Eq. (\ref{Init eq})).  The grid resolution is listed as well.  We take
our
high resolution runs to have 201 radial $\times$
53 angular zones ($\Delta \eta = 
\Delta \theta = 0.03$) and low resolution
runs to have 101 radial by 27 angular zones ($\Delta \eta = 
\Delta \theta = 0.06$).  

To perform high resolution runs
to t = 100 m ($m$ is the scale parameter defining our units) requires
requires about 7000 time steps, or 8 hours, and 4 Megawords
of memory on the NCSA Cray-2, 
while low resolution runs require about 3500 time steps, or 1 hour, and 1
Megaword.  All the runs 
presented in this paper use the angular dependence parameter $n = 4$.
We have also used $n = 2$, but in this case there is very little admixture
of the
$\ell = 4$ wave mode present in the extracted signal.  The fundamental
transverse
variable has the angular dependence
$\sin \theta \partial _\theta( \sin^{-1}\theta \partial _\theta P_{\ell})$
for a given $\ell$-mode.  For the $n = 2$ initial data there can be
coupling to higher $\ell$-modes only through nonlinear effects.
This phenomenon will be studied in future work.


We also give the total 
mass of the initial slice, computed from the ADM  formula [4]
\begin{equation}
E =-\frac{1}{2 \pi} \oint_{\infty} \nabla_{a} \Psi dS^{a}.
\end{equation}
Evaluating this in $(\eta, \theta, \phi)$ coordinates over a 2--sphere of
radius $\eta$ and assuming axisymmetry gives
\begin{equation}
M_{\rm ADM} = - \frac{1}{2} \sqrt{2m} \; e^{\eta/2} \int_{0}^{\pi}
\left(\frac{\partial \Psi}{\partial \eta} -\frac{\Psi}{2} \right)  
sin \theta d\theta,
\end{equation}
Note that in practice this quantity can only be computed at a finite
radius (at the edge of our grid).
We have chosen the scale parameter $m$ to be unity.  When no
gravitational wave is present, the mass of the system is 
$M_{ADM} = m = 1$.  

Finally, the ratio of the mass of the apparent horizon, (i.e. the outermost
trapped 
surface) [16], to the
total mass, $M_{ADM}$,  provides a measure of the size of the initial
black hole and the strength of the perturbation.
The mass of the
apparent horizon, defined as $M_{AH} = \sqrt{\frac{Area}{16 \pi}}$,
must be less than or equal to that of the event horizon.  The
apparent horizon has the computational advantage that it
can be calculated from data on any slice of spacetime, whereas 
knowledge of the event horizon requires knowledge of the
global spacetime structure of the system.  
We use the method of Cook and York [17] to find
the apparent horizon numerically on each spacetime slice.
If this ratio, $M_{AH}/M_{ADM}$,
is very close to unity, essentially all the mass energy is in the black
hole initially.  A systematic study of the initial data sets 
is presented in [18].

In Fig. 1a-c we show the extracted Zerilli function $\psi$ on the initial 
hypersurface for the three distinct families of runs listed above. This
function shows
the shape and location of both the $\ell$ = 2 and $\ell$ = 4 gravitational 
wave on the initial time slice.  In all cases the numerical grid extended
out to about $r = 200 m$.  The extraction is done only
for $r > 2m$.
The function shown in Fig. 1a gives the shape of the initial wave for the
runs 1a-e;  in these runs only the amplitude of the initial wave was
changed.
In Fig. 1b we show the initial wave for Run 2, which has a much broader
profile.
Finally in Fig. 1c we give the initial profile for a wave considered in Run
3 which has been chosen
to be a localized wave packet, cleanly separated from the hole, so that we
could
study the collision of a wavepacket with a black hole.

As mentioned in Sec. \ref{extraction}, although 
the Zerilli function $\psi$ is gauge-invariant,
it does contain near zone information.  As shown earlier[11], this
near zone contribution to
the gauge-invariant field $\psi$ can easily be removed
if desired.  However, for black hole spacetimes, 
for which the characteristic reduced
wavelength is $\hlambda = 16.8M/2 \pi = 2.7M$, 
the near-zone effect will be small, about 5\%,
even at r = 15M (see Eq. 25 in Ref. [11c]).  For this reason, we have not
performed this 
correction to the waveforms presented in this paper.  However 
comparisons of these extraction methods with those obtained using the near
zone
methods confirm that the near zone effects are indeed small.
Furthermore, as we discuss in the next section, we have used the
Zerilli equation to evolve the function $\psi$, which includes the
near zone effect, as a test of our code.  For this integration we
must use $\psi$ itself without removing the near zone effect.
The waveforms presented here, which have been extracted
at r = 15m, approximate those which would be observed in the wave zone to
within a few percent.

In Figs. 2-7 we show the results of some of the evolutions listed in Table
1.  These
figures show the time development of the Zerilli function $\psi$ at a fixed
location, as calculated by our numerical gravitational wave 
``detectors."  These
detectors are typically located at 
$r = 15m$, $r=30m$ and $r = 60m$.  In practice,
these points are actually located at fixed positions in our 
radial coordinate $\eta$.  In
principle these coordinate positions are time dependent and fall towards
the 
black hole.  However, in 
practice this drifting of radial grid points is very small (less than 2\%),
so the detectors may be regarded as being at a fixed radius r.

Now we will examine in detail the series of simulations labelled Runs 1a-d
where the shape and location of the wave were held fixed;  only the 
amplitude was varied.
In Fig. 2a we show the $\ell$ = 2 Zerilli function $\psi$ extracted at a
radius 
$r = 15m$ as a function of time for Run 1b, while in Fig. 2b we show the
$\ell$
= 4 Zerilli function.  We consider this as our fiducial run which others
in this series are compared with.  In each figure we also show the 
appropriate normal mode for
a black hole of mass $m$ as a dashed line.  Note that the
time $t$ is also measured in units of $m$.  The normal mode curves are
matched to
the extracted waveforms
by simply adjusting the amplitude and phase of the fundamental mode of the
black hole, whose wavelength and damping time are known from perturbation
theory.
For the $\ell$ = 2 waveform, the matching was done at the peak near t =
42m, while for 
$\ell$ = 4 the matching was done at the peak near t = 35m.  On the whole,
this
match is excellent.  Note that this matching procedure is performed at only
one point; the agreement is not a result of a ``best fit"
performed over an extended region of the waveform.

 
However, the waveforms show
some differences from the known fundamental normal mode of the black hole. 
As discussed above, at very early times
we do not expect the waveform to match the black hole normal mode waveform.
This part of the waveform is dominated by the initial shape and amplitude
of the
Brill wave and by the excitation process of the normal mode.  After this
early precursor
phase, the fundamental normal mode of the black hole is fully excited by
about t = 30m.  From this time until about t = 60m, the 
match is extremely good, but after
t = 60m, the wavelength of the extracted wave has begun to increase.  
This lengthening of the wave can be understood as an artifact of our
spatial resolution.
Although this run was computed with 201 radial zones, as we showed in Ref
[5], this
radial resolution reproduces the Schwarzschild background less accurately
beyond
about $t$ =  50m. 
Therefore there will be a 
phase shift in the extracted
wave due to the fact that the wave is essentially propagating on a
background which
differs somewhat from true Schwarzschild.
Secondly, the inaccuracy of the background spacetime will
cause the effective gravitational scattering potential $V(r)$ to be
different from
the actual Schwarzschild case. 
%
Since the  
normal mode generation depends sensitively on the shape of the $V(r)$, we
can expect 
small deviations from the expected frequency.  In effect, the code is
generating the normal
modes of a somewhat different background spacetime.  Note that the late
time effects
cannot be due to an admixture of higher overtones of the 
fundamental normal mode frequency.  Although they are probably excited, 
they would be damped
to a  negligible level at late times.  It would be interesting 
to study the excitation of the higher modes at earlier times, say before 
t = 30m in these figures.  Higher overtones may partially account for the 
shape of this early part of the waveform,
although at very early times the waveform will be dominated by 
the initial Brill wave.  (See, e.g., Leaver [19] for a discussion of 
higher mode excitation).

Next, in order to test the effect of grid resolution, we considered 
identical initial data parameters on a grid with half as many radial and
angular zones.
In Figures 3a and 3b we show a comparison of waveforms extracted 
at r = 15m for Runs 1b and 1c.  For these simulations, all
parameters were the same except for grid resolution.  At early times the
waveforms 
agree fairly well with each other, but at later times there are two effects
which can be seen.  First,
there is a phase shift which occurs due to the difference in resolution. 
Note that the
effect of the lower resolution is to increase the wavelength.  This
lengthening
was already discussed above for high resolution, but the problem has been
exaggerated 
at this lower resolution.  Second, at late times in the low resolution
run there is a higher frequency
signal present in both the $\ell$ = 2 and $\ell$ = 4 extracted waves.  This
is a result
of backscatter of radiation off of the grid at larger radii.  We 
have confirmed that this is indeed the cause of this high
frequency component by examining the times at which
this signal is seen in various detectors.  For the lower resolution
run, the signal is seen first at the outer detector, and propagates
inward to the inner detector at the speed of light.  
Although the grid is
chosen to have $\Delta \eta $= const, because of the logarithmic nature of
this coordinate
the interval of the Schwarzschild-like radial coordinate $\Delta r$
increases
with increasing values of $\eta $.  Eventually the separation of the 
areal radius grid points becomes large enough that the outgoing waves
suffer
some reflection. The lower the resolution, the closer the reflection occurs
to the detectors.  This
effect is also present to a much smaller degree in the high resolution run.
Note also that at low resolution the $\ell = 4$ wave has a slightly lower
amplitude
than in the high resolution run, even at early times.  As one
would expect,
high resolution is clearly required for accurate generation and propagation of
gravitational
waves.

Next, we verify that these waves are in the linear
regime by checking the
scaling of the extracted waveform with the initial wave amplitude.
In Figures 4 a,b we show the $\ell$ = 2 and $\ell$ = 4 waves extracted for
Run 1a, along
with the normal mode match.  
In this run, all parameters
are the same as in Run 1b except that the amplitude of the Brill wave is
reduced by
a factor of 5, from $A = .05$ to $A = .01$.  Notice that the shapes of the
waveforms 
shown in these figures are nearly identical to those shown in Fig. 2a,b,
and that the 
normal mode matches are extremely good in both cases.  The fact that the
shapes are
identical clearly indicates that these runs are in the linear perturbation
regime.  This is borne out by a comparison of the amplitudes for the two
runs;  for both
$\ell$ = 2 and $\ell$ = 4, the waveforms scale by a factor of 5.0, as one
expects in the linear regime.

We now examine the effect of increasing the initial wave amplitude
while maintaining the same shape.  Consider now Run 1d, 
shown in Fig. 5a,b. It is again identical 
to that shown in Fig. 2a,b, except that now we increase the amplitude by a
factor of
10 to $A$ = 0.5.  Note that at first glance the $\ell$ = 2 waveform is very
well
fit by the normal mode match, 
even at late times, while the $\ell$ = 4 normal mode
match is not as good as that shown in Fig. 2b.  These effects can be 
understood by considering both small, nonlinear physical effects, 
and also small errors in the
calculating the background spacetime.  First of all, we 
note that although the amplitude
$A$ of the initial Brill wave has been scaled up by a factor 
of 10 (when compared to Run 1b), 
the resulting waveforms scale up by the smaller factor of 
9.1.  This indicates that
we are entering the nonlinear regime.  This 
effect can also be seen in the
calculation of the mass of the initial slice.  In Run 1b, 
the ADM mass of the initial
slice is $M_{ADM}$ = 0.98, while for Run 1d the mass is 
$M_{ADM}$ = 0.92.  Therefore, the black
hole has a smaller mass than in the previous cases, and it 
should radiate with its
own characteristic set of normal modes, which due to the slightly lower
mass, will
have slightly $shorter$ wavelengths.  Hence, we see clearly in Fig. 5b that
the
$\ell$ = 4 extracted wave has a slightly shorter wavelength than the
fundamental
$\ell$ = 4 normal mode for a black hole of unit mass.  At late
times the extracted wave ``catches up" with the matched normal mode; due to
the
resolution effects discussed previously the wavelength is artificially
lengthened
slightly at late times.  The same effect is present in Fig. 5a, but because
of the 
longer wavelength of the $\ell$ = 2 mode, the match simply appears better
over a
longer period of time.  In this case we see a conspiracy of a small
physical, nonlinear
effect and a small grid resolution effect which together cancel to first
order.  Note that
these effects would probably have gone unnoticed had we not analyzed the
$\ell$ = 4 mode
as well as the $\ell$ = 2 mode.


In the discussion to follow, we show how the extracted waveforms depend
upon
the initial shape parameters.
Consider the evolution of the initial data shown in Fig. 1b, listed as
Run 2 in Table 1.  In this case, the amplitude is again very small ($A$ =
0.05),
but the shape of the initial Brill wave has been changed dramatically. Here
we have created a very broad initial wave, centered on the isometry
surface of the black hole.
The waveform extraction for this run is shown in Figures 6a,b.  Note that
the precursor region of the waveform ($0 < t < 20m$) has a very different
shape
from all of the previous results, but the normal mode excitation is the
same
as before.  Thus, as expected, the fundamental normal mode of the black
hole
is excited for each $\ell$ value independent of the kind of perturbation, 
and it dominates the waveform after the initial
precursor.  The shape of the precursor, on the other hand, is very strongly
dependent 
on the shape of the perturbing gravitational wave.  In this particular case the
wavelength of the precursor is significantly longer than the normal mode
excitation,
while for the each of the previous cases it was of the same order.

Finally, in Fig. 7 we show the results of the evolution of the initial data
shown in Figure 1c, labeled Run 3.  In this case we created a narrow
wavepacket
cleanly separated from the black hole as initial data, in the spirit of
dropping
a particle into a black hole.  When this system is evolved, the initial
wavepacket
splits into two packets traveling in opposite directions.  One packet moves
to
the right, away from the hole, while the other moves toward and collides
with the
hole, perturbing it and causing it to radiate at its normal mode frequency.
 The
initial wavepacket is sharply peaked near $r = 10m$.  Since we have
extracted
the waveform at $r = 15m$,
the first narrow pulse in Figure 8, peaked near $t = 5m$, is 
actually the outgoing piece of this wavepacket.  The ingoing wavepacket
sets up
the normal mode excitation, which appears in Fig. 7 after $t = 20m$.  This
first broad
peak is clearly not completely characterized by the fundamental normal mode
of the hole; it has a complicated shape which contains information about
the
normal mode excitation process, the scattering properties of the potential,
and probably higher overtone modes of the hole.  By about $t = 50m$ the
system
is apparently oscillating at the fundamental normal mode frequency, but by
this
time because of grid resolution effects the waveform has widened somewhat
and
the normal mode match is less accurate.  At still later times, say for $t >
80m$,
we again see the high frequency component of the wave, which is a result
of backscattering of the initial outgoing wavepacket off of the 
coarser grid at larger radii.  With somewhat higher grid resolution
we will be able to study this normal mode excitation process in more
detail.
This will be particularly interesting in the context of much largeramplitude gravitational 
waves as it will be possible to regulate
the amount of wave energy one is putting into the perturbation
by starting the (ingoing) perturbing wave in the linear regime. 
This process will be studied in more detail in a future
paper.

\subsection{Evolution of Perturbations}\label{perturbed}

In the previous section our interest was in $interpreting$ 
physical results in the context of what is expected from black hole 
perturbation theory.  This is a code test only to the extent that one
expects 
normal modes to be excited, and the code should produce gravitational waves
of the 
correct wavelength and decay constant if the initial gravitational waves
have 
sufficiently small amplitude.  As discussed above the precursor waveform 
and the amplitude of the normal mode excitation depend in a complicated 
way on the details of the initial 
Brill wave and its interaction with the black hole.  Furthermore, it is 
conceivable that certain kinds of perturbations would preferentially excite

higher overtones of the black hole, leading to a more complicated waveform.
 
We have not 
attempted to use perturbation theory to predict such details of waveforms
to be 
expected from a given initial Brill wave.

In this section we are interested in taking initial Brill wave data 
which may be considered a perturbation on the black hole, evolving 
it with the full numerical code, and comparing the results to an
$evolution$ calculated by 
perturbation theory.  As described above in Sec. \ref{extraction}
and Sec. \ref{excitation}, the 
Einstein equations reduce to a system of background equations 
representing the spherically symmetric part of the system, and a set 
of equations representing the evolution of the perturbed quantities 
on the spherical background.  Our gauge choices are such that the 
spherically symmetric part of the field 
is essentially Schwarzschild in the region 
where we extract wave information, so we may approximate the 
system as the Zerilli function $\psi$ propagating on the 
Schwarzschild background through its wave equation Eq. (\ref{Zerilli eq}).   
Therefore we should be able to take initial data and evolve it with 
either the full two-dimensional nonlinear code or the Zerilli equation 
and get the same results.  The comparison of these two calculations 
provides a very strict code test.  It simultaneously tests both the 
accuracy of the waveform extraction method (since this procedure
provides the boundary data to be evolved and the comparison data at the 
outer radius) and the accuracy of the 
propagation of gravitational waves across the grid (since the waves are
propagated
by both the Zerilli equation and the nonlinear code).

In our code test, the integration of the Zerilli equation is most 
conveniently carried out on ingoing and outgoing null cones 
with coordinates given by:
\begin{equation}
u = t - r_{*} \mbox {   and   } v = t + r_{*},
\end{equation}
where $r_{*}$ is defined above.  In these coordinates the Zerilli equation
becomes
\begin{equation}
\frac{{\partial}^2 \psi}{\partial u \partial v}   + 4 V_Z(r)  \psi = 0.
\end{equation}
This is shown schematically in Fig. 8.  As the code runs, we 
extract data along a timelike line at along the left boundary at 
$r_{left} = 15m$ which we will 
use as boundary data for the Zerilli equation.  Similarly we 
extract data at each grid point along the first outgoing null cone 
labeled $u_{0}$ in Fig. 8 out to the end of the grid.  This procedure 
provides initial data for the Zerilli equation, which can then 
propagate the solution to any point in the interior region of the 
figure.  We then compare data evolved with the Zerilli equation 
in the interior along another timelike line at, say, $r_{right} = 30m$, 
to data extracted from 
the two-dimensional numerical code at that same radius.  
If the simulation runs from $t = t_0$ to $t = t_f$, note that we 
must provide data along the null cone $u_0$ all the way out
to $r_{edge}$ since in principle there may be ingoing waves
that were backscattered off
the background curvature.  (This backscattering is a physical
effect, unrelated to the backscattering due to a coarse grid discussed in
the
above section.)  Our scheme makes no
assumption that the waves are outgoing.  Furthermore, we can also
predict through the Zerilli equation the value of the Zerilli function
$\psi$ in the time interval $t_{f} < t < t_{p}$, which is beyond
the run time of the code.  Note that this test
has nothing to do with black hole normal modes;  it is a strict test of 
the code's ability to accurately evolve initial data regardless of 
the waveform.

In Figs. 9a and 9b we show the results of one such test.  
In this case we reexamine
our Run 1a from the previous section, with Brill wave parameters $A$ =
0.01,
$\sigma^2$ = 1, $\eta_0$ = 1, and $n$ = 4.  In this case we used the
Zerilli
equation to propagate extracted data from $r_{left} = 15m$ to $r_{right} =
30m$.
The wave extracted from our fully nonlinear numerical code is shown as a 
solid line, while the results from our perturbation theory test are shown
as
a dashed line.  The data from the Zerilli 
test begin at $t = 16.7m$, and as we
discussed above, we are able to propagate the solution beyond $t_f = 100m$
to
$t_p = 117m.$  These are typical of all of our results at this resolution.
This test indicates that both the extraction and propagation of waves in 
the system are very accurate. 

In Figure 10a,b we show results for the identical initial conditions, but
evolved
on a low resolution (101 x 27) grid.  In Fig 10a we see that the $\ell$ = 2

waves are 
in fairly good agreement with the result of integrating the Zerilli
equation
until about about $t = 60m$, after which time the signal is degraded.  The
high frequency component appearing at late times is again a result of the
backscatter off of the coarse grid at larger $r$.  Note that this
high frequency wave also appears in the Zerilli propagated waveform, but
with a time delay of about $30 m$.  This shows that this signal
really does represent ingoing waves coming from large 
$r$ back into the region of interest.  In Fig. 10b
we show the results for the $\ell$ = 4 wave..  It is immediately apparent
that the $\ell$ = 4 wave cannot be accurately resolved at this grid
resolution.
The reason the waveform shown in Fig. 11b appears to be so much worse than
the 
low resolution run shown in
Fig. 3b is that here all waveforms are computed at $r = 30m$.  At this
radius, not only
are the spatial zones spaced farther apart in the physical radius $r$, but
also the signal
arrives at a later time when the code is known to be less accurate. 
However, as shown
by Fig. 10a, at high resolution such effects are minimal.  Furthermore, had
we not gone
to the $\ell$ = 4 extraction, we would not have seen the extent of these
effects at
low resolution.

\subsection{Stronger Gravitational Waves}\label{stronger}

We will briefly examine results for an initial Brill wave which has an 
amplitude high enough that it contributes 
significantly to the total mass of the system.  Although
we defer a complete treatment of this problem to a future paper, we will
provide an example here to demonstrate the robustness of the code and to 
hint at what sorts of nonlinear effects one is likely to see in the strong 
wave regime.  Run 1e, listed in Table 1, is a high amplitude version of the

runs 1a-d discussed above. With the wave amplitude $A = 1.5$, the mass of
the 
initial slice is $M_{ADM} = 1.371$.  Since we expect that some of this
radiation will propagate away from the black hole during the evolution,
the final state of the system should be a black hole with a mass less than
the mass of the initial slice.  

We have also computed the mass of the apparent horizon, both on
the initial slice and during the evolution.  As shown in Table 1, initially
the apparent horizon has a mass which is only about 59\% that
of the total mass of the slice.  Although the event horizon 
has a larger mass than this, this is an indication that
initially the black hole contains less mass than is present in the
system.  As the system evolves, the mass of the apparent horizon
grows, leveling off after about $t = 30m$ to a mass of about 1.3m,
indicating that the mass of the black hole has increased significantly
as it absorbs gravitational wave energy. 


In Fig. 11a,b we show waveforms extracted at $r = 15m$ for both the $\ell$
= 2
and $\ell$ = 4 Zerilli functions.  For comparison, we have overlayed the
waveform extracted in Run 1d.  At first glance the waveforms appear to
be very similar to those shown in Figs.2-5.  However, note first of all
that
although the amplitude A of the initial Brill wave is larger than that for
Run 1d (also shown in Fig. 5a,b) by a factor of 3, the amplitude of the
extracted waves
is greater by less than a factor of 2.  While
the low amplitude runs scaled linearly with the amplitude of the
initial Brill wave, here it is clear that nonlinear effects are
important.  Next, we note that the scaling factors
for the two runs are different for the $\ell$ = 2 and $\ell$ = 4 waveforms.
Finally, we note that the wavelengths are significantly longer for this
higher amplitude run than in any of the previous cases studied so far.
In fact, for the $\ell$ = 2 waveform, we measure the wavelength of the
signal to be $\lambda_{extracted} = 22.8$, whereas for a black hole of unit
mass we would expect
a wavelength $\lambda_{unit mass} = 16.8$.  By taking the ratio of these
numbers, we obtain an estimate of the mass of the black hole in this system
of $M = 1.36m$, which is  consistent with the calculation
of the mass of the initial slice.  (In fact, since our study of low
amplitude waves in Sec.\ref{perturbed} shows a slight tendency for the
wavelength to increase artificially due to numerical effects,
the mass of the black hole is probably slightly below 1.36$m$.)
This indicates that most of the mass of
the Brill wave was swallowed by the black hole, which then radiates at its
(lower) quasinormal frequency.  Furthermore, by comparing the damping time
of the extracted wave with the known damping time for a black hole of unit
mass, we obtain an independent estimate of the black hole mass $M = 1.39m$.
A similar analysis of the extracted $\ell$ = 4 waveform shows that by comparing
wavelengths, we estimate the mass of the black hole to be $M = 1.31m$,
while
by comparing the damping times we find a mass $M = 1.2m$.  

These calculations
are all reasonably consistent with each other, and with the apparent
horizon mass measurements, and indicate that most of the radiation
present on the initial slice was captured by the black hole.  But since the
mass contained in this radiation was a significant fraction of the ``bare"
mass
of the black hole, the newly formed larger black hole radiates energy away
at a lower frequency.  This is a nonlinear effect in the sense
that much of the initial wave, which is essentially transverse in nature,
is effectively
converted into longitudinal field 
as the black hole absorbs the wave, increasing the
mass of the black hole background.
It should be noted that the radiation extraction method, which is
valid only at linear order in sphericity of the system, is still a good
approximation in this regime where $nonlinear$ physical effects are seen 
since the majority of the wave is going down the hole. 

\section{Conclusions and Future Directions}\label{conclusions}

There are several major results which we have reported in this paper.   
For low amplitude gravitational waves impinging on the black hole, we find 
that after a brief precursor waveform which depends on the shape of the 
ingoing wave pulse,  the fundamental normal modes of the black hole are 
strongly excited.  This excitation is independent of the shape of the
incoming 
gravitational wavetrain.  The waveforms extracted with the gauge 
invariant technique described in Sec.\ref{extraction} above agree, with 
high precision, with those predicted by black hole perturbation theory  
for both fundamental $\ell=2$ and fundamental $\ell=4$ modes.
The adherence to analytic results that we observe is an 
important test of the code in the strong-field regime.
With a systematic study of the waveforms, we have been able to
identify and understand very small discrepancies which arise
due to both physical and code effects.  Although these effects are
small, we believe that it is crucial to identify them as they
may be magnified when higher amplitude waves are considered.

For larger amplitude perturbations we have found cases where
apparently the incoming gravitational wave is adding non-negligible mass to
the black hole, thus changing the frequency of emitted normal modes.
We have also taken data on a null-cone, propagated it with the Zerilli
equation 
and compared it in detail with 
propagation of the same wave computed using the full nonlinear,
two-dimensional code.
This compelling test showed that the code accurately propagates 
gravitational waves in the weak field region of the spacetime.
It also demonstrates that the gauge invariant wave-extraction
method is useful for code-diagnostics as it facilitates the
comparison of numerical results against perturbation theory.

This code is a robust tool for studying black hole spacetimes in a
variety of contexts.  We have a number of projects underway.  In future
work we plan to use this code to study physics of large
amplitude gravitational waves interacting with a black hole.  We
have already evolved spacetimes where the mass of the system is
many times that of the bare black hole.  We are also performing a
systematic study of various scalar, tensor, and gauge-invariant
quantities to be used for tracking radiation in regimes where the
linearized waveform extraction may be inadequate.  Furthermore,
this code has also been modified to study the evolution of two 
colliding black holes.  These studies will be presented in future papers. 

\acknowledgements

We would like to thank Robert Bartnik, Greg Cook, and Wai-Mo Suen for
helpful discussions.  A.A. gratefully acknowledges the support of an NSF
Division of Advanced Scientific Computing fellowship.  The
calculations were performed on the Cray Y-MP and Cray-2 supercomputers
under a grant from the National Center for Supercomputing Applications.

\appendix{Projection of Metric Functions}

In this appendix we give details of the method used to project the
metric functions onto the spherical and perturbed metric.
The total spacetime metric $g_{\alpha \beta}$, which is generated 
numerically,
is written as a sum of a spherical and a nonspherical part:
\begin {equation}
g_{\alpha \beta}  =   
\stackrel{o}{g}_{\alpha \beta} +  h_{\alpha \beta},
\end {equation} 
where  
\begin {equation}
\stackrel{o}{g}_{\alpha \beta}=\left( \begin{array}{cccc}
-N^2 & 0 & 0 & 0 \\
0 & A^2 & 0 & 0 \\
0 & 0 & R^2 & 0 \\
0 & 0 & 0 & R^2 sin^2 \theta
\end {array} \right)
\end {equation}
represents the spherically symmetric part and
\begin{equation}
\label{pert metric}
h_{\alpha\beta}=\left( \begin{array}{cccc}
N^2H_0Y_{\ell 0} & H_1Y_{\ell 0} & h_0 Y_{\ell 0,_\theta} & 0 \\
H_1Y_{\ell 0} & A^2H_2Y_{\ell 0} & h_1 Y_{\ell 0,_\theta}& 0 \\
h_0 Y_{\ell 0,_\theta} & h_1 Y_{\ell 0,_\theta} & R^2 \left(K + G 
\frac{\partial^2}{\partial \theta^2}\right)Y_{\ell 0} & 0 \\
0 & 0 & 0 & R^2 \left( K sin^2\theta + G sin\theta cos\theta 
\frac{\partial}{\partial \theta}\right) Y_{\ell 0}
\end {array} \right)
\end {equation}
represents deviations from sphericity.  On the $(t,r)$ submanifold the 
metric functions can be easily projected using the orthogonality of the
spherical
harmonics $Y_{\ell m}$.  For example one can project the $g_{tt}$ component
of the metric onto $\stackrel{o}{g}_{tt}$ and $h_{tt}$ by performing 
the following integral over a 2-sphere of radius r:
\begin{equation}
\int g_{tt} Y_{\ell ' 0} d \Omega \ = 2 \pi \int_{0}^{2 \pi} \left(-N^2 
\sqrt{4 \pi} Y_{00}
+ N^2 H_0 Y_{\ell 0}\right) Y_{\ell ' 0},
\end{equation}
so that
\begin {equation}
N^2 = -\frac{1}{2} \int_{0}^{\pi} g_{tt} sin \theta d \theta .
\end{equation}
Similarly, one can compute
\begin{equation}
A^2 = \frac{1}{2} \int_{0}^{\pi} g_{rr} sin \theta d \theta .
\end{equation}
Furthermore, one can also compute the true spherical areal radius, in the
presence of a perturbation, by computing
\begin{equation}
R^2 = \frac{1}{2} \int_{0}^{\pi} g_{\theta \theta} sin \theta d \theta.
\end{equation}

The Regge-Wheeler perturbation functions $H_0$, $H_1$, and $H_2$ 
can also be computed in the same way, but in the case of a Schwarzschild
background, the functions $H_0$ and $H_1$ are not needed to compute the
gauge-invariant function $\psi$ (on a time dependent background they would
be required.)  The function $H_2$ can be computed from:
\begin{equation}
{H_2}^{(\ell )} = \frac{2 \pi}{A^2} \int_{0}^{\pi} g_{rr} Y_{\ell 0} 
sin \theta d \theta .
\end{equation}
The functions $H_0$ and $H_1$ could be computed in a similar manner if
needed.

%Derivatives of the $Y_{\ell m}$ are also orthogonal and have the 
%following normalization condition:
%\begin{equation}
%\int Y_{\ell \theta} Y_{\ell '\theta} d \Omega  = \left( \ell (\ell + 1)
%\right)\delta_{\ell \ell '} 
%\end {equation}
On a Schwarzschild background only $h_1$ is needed from this "vector" part
of the metric (the $t-[\theta \phi]$ and $r-[\theta \phi]$ part):
\begin{equation}
{h_1}^{(\ell)} = \frac{\pi}{3} \int_{0}^{\pi} Y_{\ell 0, \theta} g_{r 
\theta} sin \theta d \theta .
\end{equation}
In the gauge used by our code for this paper, the 3-metric is diagonal so
this term explicitly vanishes.  But note that the wave extraction 
method can easily incorporate a term like this if it is present. 

Finally one must compute the Regge-Wheeler functions $G$ and $K$.  
Unfortunately, the tensor spherical harmonics used on the $(\theta , \phi
)$ 
submanifold in Eq. (\ref{pert metric})
are not linearly independent, so a slight change of basis must be used 
here.  This part
of the metric can be written in terms of orthonormal tensor spherical 
harmonics [20], and then the functions $G$ and $K$ can be 
computed.  Explicit expressions for ${G}^{(\ell)}$ 
and ${K}^{(\ell)}$ are:
\begin{equation}
{G}^{(\ell)} = \frac{2 \pi}{R^2} \int_{0}^{\pi} 
\frac{ \left(g_{\theta \theta} sin^2 \theta - g_{\phi \phi} \right)
\left(-cos \theta Y_{\ell 0, \theta} + sin \theta Y_{\ell 0,\theta \theta}
\right)d\theta}
{\left(\ell (\ell+1)(\ell+2)(\ell-1)\right) sin^2 \theta} .
\end{equation}
 and
\begin{equation}
{K}^{(\ell)} = \frac{\ell (\ell+1)}{2} {G}^{(\ell)} + 
\frac{\pi}{R^2} \int_{0}^{\pi} 
\left(g_{\theta \theta} + \frac{g_{\phi \phi}}{sin^2\theta} \right)
 sin \theta Y_{\ell 0} d \theta .
\end{equation}

These Regge-Wheeler functions can also be written in terms of the
pure orbit harmonics used in [11].

%\newpage
\begin{references}
\bibitem{} a)  R.F. Stark and T. Piran,  {\it Phys. Rev. Lett.}
{\bf 55}, 891 (1985);  b)  C.T. Cunningham, R.H. Price and V. Moncrief, 
{\it Astrophys. J.} {\bf 224}, 643 (1978); {\bf 230}, 870 (1979); 
{\bf 236}, 674 (1980);  c) Y. Sun and R. Price, {\it Phys. Rev.} {\bf D
41}, 
2492, (1990);  d)  E. Seidel, {\it Phys. Rev D} {\bf 44}, 950 (1991).

\bibitem{} L. Smarr, in {\it Sources of Gravitational Radiation}, 
edited by L. Smarr, (Cambridge University Press, Cambridge, England, 1979).

\bibitem{} R. Vogt, to appear in {\it Proceedings of the Sixth Marcel
Grossmann Meeting on General Relativity}, (1992).

\bibitem{} R. Arnowitt, S. Deser and C.W. Misner, in {\it Gravitation: An 
Introduction to Current Research}, edited by L. Witten, (Wiley, New York,
1962).

\bibitem{} D.H. Bernstein, D.W. Hobill, E. Seidel and L.L. Smarr, in 
preparation.

\bibitem{} A. Einstein and N. Rosen, {\it Phys. Rev.}, {\bf 48}, 73,
(1935).

\bibitem{} D.R. Brill, {\it Ann. Phys.} {\bf 7}, 466, (1959).

\bibitem{} B. Schutz, in {\it Dynamical Spacetimes and Numerical
Relativity}, edited by J. Centrella, (Cambridge University Press, 
Cambridge, England, 1986).

\bibitem{} J.L. Anderson and D.W. Hobill {\it J. Comp. Phys.} , 
{\bf 75}, 283 (1988).

\bibitem{}  a) J.M. Bardeen and T. Piran, {\it Phys. Rep.} {\bf 96},
205 (1983);  b)  C.R. Evans, in {\it Dynamical Spacetimes and Numerical
Relativity}, edited by J. Centrella, (Cambridge University Press, 
Cambridge, England, 1986).

\bibitem{}  a)  A.M. Abrahams, in {\it Frontiers in Numerical Relativity}
edited by C.R. Evans, L.S. Finn, and D.W. Hobill
(Cambridge University Press, Cambridge, England, 1989);  
b) A. M. Abrahams and C. R. Evans, {\it Phys. Rev.} {\bf D37}, 318 (1988); 

c)  A. M. Abrahams and C. R. Evans, {\it Phys. Rev. } {\bf D 42}, 2585,
(1990).

\bibitem{} a)  T. Regge and J.A. Wheeler, {\it Phys. Rev.} {\bf 108}, 1063,

(1957);  b)  F.J. Zerilli, {\it Phys. Rev. Lett.}, {\bf 24}, 737, (1970);
c)  V. Moncrief, {\it Ann. Phys.}, {\bf 88}, 323, (1974);  d) U.H. Gerlach 

and U.K. Sengupta, {\it Phys. Rev. D} {\bf 19}, 2268 (1979);  e)  E.
Seidel,
{\it Phys. Rev.}, {\bf D42}, 1884, (1990).

\bibitem{}  a)  S. Chandrasekhar and S. Detweiler,  
{\it Proc. R. Soc. London A}, {\bf 344}, 441 (1975);  b) E. Leaver, {\it 
Proc. R. Soc. London A}, {\bf 402}, 285 (1985);  c)  S. Iyer and C.M. Will,
{\it Phys. Rev. D} {\bf 35}, 3621 (1987).

\bibitem{}  B. Whiting, {\it J. Math. Phys.} {\bf 30} {1301}, (1989).

\bibitem{}  a) M. Davis, R. Ruffini, W. Press, and R. Price, 
{\it Phys. Rev. Lett.} {\bf 27}, 1466, (1971); b) Y. Kojima and T.
Nakamura,
{\it Prog. Theor. Phys.} {\bf 71}, 79, (1984).

\bibitem{}  {\it The Large Scale Structure of Spacetime}, 
S. Hawking and G. Ellis, (Cambridge University Press, Cambridge, England,
1973).
 
\bibitem{} G. Cook and J. York, {\it Phys. Rev. D} {\bf 41}, 1077, (1990).

\bibitem{} D.H. Bernstein, D.W. Hobill, E. Seidel and L.L. Smarr, in 
preparation.

\bibitem{}  E. Leaver, {\it Phys. Rev. D} {\bf 34}, 384 (1986).

\bibitem{}  F. Zerilli, {\it J. Math. Phys.} {\bf 11}, 2203, (1970).

\end{references}
\newpage

%\centerline {Figure Captions}

\figure{The $\ell$ = 2 and $\ell$ = 4 Zerilli functions $\psi$ are shown 
on the initial time slice labeled by $t$ = 0.  (a)  The gravitational wave
is shown for Run 1b, described in Table 1.  Runs 1a-e have the same shape 
parameters, but the amplitude is varied.  (b)  The initial wave is shown
for 
Run 2.  In this case a much wider initial packet was chosen.  (c)  The
initial
wave is shown for Run 3.  In this case the parameters were chosen so that
a localized wavepacket would propagate in towards the hole.}

\figure{The gauge-invariant function $\psi$, extracted at r = 15m, is shown for
Run 1b.  (a)  The $\ell$ = 2 waveform is shown as a solid line, with the
$\ell$ = 2
fundamental normal mode frequency, shown as a dashed line, matched to 
the peak near t = 42m.  (b)  The $\ell$ = 4 waveform is shown as a solid
line, 
with the $\ell$ = 4
fundamental normal mode frequency, shown as a dashed line, matched to 
the peak near t = 35m.  See text for details.
}

\figure{A comparison is made of the extracted waveforms shown in Fig. 3
(Run 1b)
with a low resolution run with the same run parameters (Run 1c).  (a)  The
$\ell$ = 2
extraction is shown for high and low resolution.  (b)  The $\ell$ = 4
extraction is
shown for both resolutions.  In all cases, the extraction was done at r =
15m.
}

\figure{Extracted waveforms are shown for Run 1a, which is identical
to Run 1b, except that the amplitude has been reduced by a factor of 5 to
$A$ = 0.01.
(a)  The $\ell$ = 2 waveform, extracted at $r = 15m$ is shown as a solid
line, while
the normal mode match, performed on the peak near $t = 42m$, is shown as adashed line.
(b)  Similarly, the $\ell$ = 4 waveform, also extracted at $r = 15m$ is 
shown along with the $\ell$ = 4 normal mode match.
}

\figure {The extracted waveforms are shown for Run 1d, which is
identical to Run 1b except that the amplitude has been increased by a
factor of 10
to $A$ = 0.5.  
(a)  The $\ell$ = 2 waveform, extracted at $r = 15m$ is shown as a solid
line, while
the normal mode match, assuming unit black hole mass, is shown as a dashed
line.
(b)  Similarly, the $\ell$ = 4 waveform, also extracted at $r = 15m$ is 
shown along with the $\ell$ = 4 normal mode match.  In this run the ADM
mass
$M_{ADM}$ is somewhat lower than in the previous cases, causing a somewhat
shorter
wavelength than expected for a black hole of unit mass.  See text for
discussion.
}

\figure{The waveforms for Run 2, extracted at r = 15m, are
shown as solid lines, while the normal mode matches are shown as dashed
lines.  The shape of the initial Brill wave considered here, 
shown in Fig. 1b, is very different from the other cases considered
so far.
(a)  The fundamental $\ell$ = 2 mode is shown.  (b)  The fundamental $\ell$
= 4
mode is shown.
}

\figure{The $\ell$ = 2 waveform is shown for Run 3, extracted
at $r = 15m$.  This waveform results from the evolution of the initial data
shown in Fig. 1c, which represents a narrow wavepacket which collides with
the hole.
The first pulse in this figure, peaked near t = 5m, results from the
outgoing
part of the initial wavepacket, while the later development results from
the
response of the black hole to the ingoing wavepacket.  The dashed line
shows
the fundamental normal mode match to the peak near $t = 50m$.  See text for
details.
}

\figure{The method used to test the numerical code against perturbation
theory is shown schematically.  Data from the numerical code are extracted
on the left boundary labeled by $r_{left}$ and also along the first
outgoing null ray labeled by $u_0$.  These lines are shown in dark.  The 
code is run from $t = t_0$ until $t = t_f$.  The boundary data on the
dark lines are provided to a numerical code written to evolve the Zerilli
equation.  The data evolved using the Zerilli equation are compared with
data extracted from the full 2-D numerical code at the radius $r =
r_{right}$.
Note that the boundary data supplied to the Zerilli equation allow us to 
extend the solution beyond the evolution of the full numerical code, to a
time $t = t_p$.
}

\figure{(a) A waveform extracted at $r = 30m$ from the results of a 
fully nonlinear
evolution (solid line) is compared to that obtained through perturbation
theory (dashed line) as discussed in the text.  These results are
for an $\ell = 2$ extraction.  Note the excellent agreement.
(b)  The same test shown in Fig. 9a, but for an $\ell = 4$ extraction.
In both cases high grid resolution (201 x 53) was used.
}

\figure{(a) The same test shown in Figure 9a is performed here, but for
the low grid resolution (101 x 27).  The accuracy is significantly less than
in the high resolution run shown in Figure 9a, especially after $t = 60m$.
(b)  The $\ell$ = 4 code test is performed for the low resolution grid.  In
this
case the grid is clearly inadequate to propagate the higher frequency mode,
so that although the $\ell$ = 2 mode is reasonably propagated for part of
the run, the higher frequency components are grossly inaccurate.  
}

\figure{Extracted waveforms are shown for Run 1e as solid lines.  This run
had the same initial data shape as Runs 1a-d, but a significantly
larger amplitude. For comparison, we have overlayed the waveforms for
Run 1d as dashed lines. (a)  The $\ell = 2$ extracted waveform is shown.  
As discussed in the text, the addition of the gravitational wave
has created a larger mass black hole, which oscillates with a longer
wavelength.  The wavelength of this signal is measured to be $\lambda =
22.8m$,
from which the mass of the hole is estimated to be about 1.36m.  
(b)  The $\ell = 4$ extracted waveform is shown.  The wavelength of 
this signal is measured to be $\lambda = 10.2m$,  which is consistent with
the
above measurement for the mass of the black hole.
}
\widetext
\begin{table}
\begin{tabular}{|l||llllllll|} \hline
     & $A$  &  $\sigma^2$ & $\eta_0$  &  n  &  $M_{ADM}$ & $M_{AH}/M_{ADM}$
& Resolution\\ \hline \hline
Run 1a    &   0.01    &   1   &  1  &  4  &  .996 &  .999 & high \\
Run 1b    &   0.05    &   1   &  1  &  4  &  .979 &  .999 & high \\
Run 1c    &   0.05    &   1   &  1  &  4  &  .979 &  .999 & low \\
Run 1d    &   0.50    &   1   &  1  &  4  &  .921  &  .916 & high \\
Run 1e    &   1.50    &   1   &  1  &  4  &  1.371 &  .591 & high \\
Run 2     &   0.05    &   2   &  0  &  4  &  .960 &  1.000 & high \\
Run 3     &   0.01    &  .2   &  3  &  4  &  .9995 &  .999 & high \\
\end{tabular}
\caption{Run parameters for various low amplitude gravitational waves
incident on a black hole.  The radial parameters are the initial 
amplitude A, the width $\sigma^2$ of the 
gaussian, and the initial location $\eta_0$ of the gaussian. The angular
dependence of the gravitational wave is given by sin$^n\theta$, the total
ADM mass $M_{ADM}$ of the initial hypersurface, and the ratio
of the mass of the apparent horizon to the total ADM mass, are also listed,
as is the
grid resolution.  High resolution is 201 x 53 zones ($\Delta \eta = 
\Delta \theta = 0.03$), while low resolution
is 101 x 27 zones ($\Delta \eta = 
\Delta \theta = 0.06$).}
\end{table}
%\end{narrowtext}
\end{document}


.
