\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} .