Extended two dimensional characteristic framework to study nonrotating black holes
We develop a numerical solver, that extends the computational framework considered in [Phys. Rev. D 65, 084016 (2002)], to include scalar perturbations of nonrotating black holes. The nonlinear Einstein-Klein-Gordon equations for a massless scalar field minimally coupled to gravity are solved in two spatial dimensions (2D). The numerical procedure is based on the ingoing light cone formulation for an axially and reflection symmetric spacetime. The solver is second order accurate and was validated in different ways. We use for calibration an auxiliary 1D solver with the same initial and boundary conditions and the same evolution algorithm. We reproduce the quasinormal modes for the massless scalar field harmonics , and . For these same harmonics, in the linear approximation, we calculate the balance of energy between the black hole and the world tube. As an example of nonlinear harmonic generation, we show the distortion of a marginally trapped two-surface approximated as a q-boundary and based upon the harmonic . Additionally, we study the evolution of the harmonic in order to test the solver in a spacetime with a complex angular structure. Further applications and extensions are briefly discussed.
pacs:04.25.D-, 04.30.Db, 04.70.Bw, 04.30.-w
Relevant astrophysical applications of the characteristic formulation of numerical relativity (1) require its adaptability and extension to a variety of scenarios. Although the Cauchy approach in numerical relativity has proven relatively successful in the simulation of binary black holes (2), the accurate prediction of wave forms from black hole-black hole, black hole-neutron star, black hole-boson star binaries as sources of gravitational radiation stands as formidable pending problems. The characteristic formulation of general relativity offers an alternative for the accurate prediction of waveforms from such astrophysical scenarios, but further improvements are mandatory to make it attractive and competitive.
One of the prime factors affecting the accuracy of any characteristic code is the introduction of a smooth coordinate system covering the sphere, which labels the null directions on the outgoing (ingoing) light cones. Interestingly, this is also an underlying problem in meteorology and oceanography (3). The LEO code, a large scale computational framework based on the characteristic formulation (4), was inspired by global forecasting techniques (5), (6) and showed great potential in handling 3D problems. The solver was tested solving the Einstein-Klein-Gordon (EKG). Despite its simplicity, analytical studies of this toy model for a self-gravitating massless scalar field show that it exhibits highly nonlinear physics (7); (8); (9); (10). One dimensional numerical simulations of the EKG led to the discovery of critical phenomena (11) and to reveal some features of the asymptotic spacetime structure. For instance, the Bondi mass and the scalar monopole moment satisfy an asymptotic relation at high amplitudes (12). The Bondi mass and news function reflect the discretely self-similar behavior (13).
Extending the work of Gómez et al. (14) and Papadopoulos (15), here we incorporate a massless scalar field and solve the 2D EKG system. We perform numerical validations that include tests of convergence, the simulation of the exponentially damped oscillation modes, called quasinormal modes (QNM) and the energy conservation (in the linear approximation) for the massless scalar field. Gravitational radiation waveforms and the nonlinear regime deserve especial attention and are postponed for a future study. However, we include a calculation of the distorted marginally trapped two-surface. The solver developed can be considered as an intermediate step, both in computational cost and dimensions. The 2D EKG is an interesting problem in itself, well suited to explore global issues (13), (16).
The paper is organized as follows. In Sec. II we setup the EKG system for a massless scalar field minimally coupled to gravity. In this section we also consider issues about QNM, energy conservation and marginally trapped surfaces. Sec. III is dedicated to details about the numerical implementation, and to the tests of convergence to second order of the nonlinear solver. In Sec. IV we present our results. Finally we conclude in Sec. V with some remarks.
ii.1 The EKG problem
In general, the field equations for a massless scalar field minimally coupled to gravity are
and can be reduced to
which have to be considered together with the wave equation
in order to complete the EKG system.
ii.2 Two–dimensional ingoing formulation
The initial-boundary value problem is formulated following the Winicour-Tamburino framework (17); (18), with ingoing light cones emanating from a timelike world tube (see Fig. 1). The characteristic initial value problem of the EKG system can be explicitly written using the ingoing line element in the case of axial and reflection symmetry (19)
where , , , This metric is twist-free, and therefore rotation is not permitted. Thus we get the following field hypersurface equations
and the evolution equations
where and .
ii.3 Scattering off a Schwarszhild black hole
On we can specify the boundary conditions in many ways. In this work we do the following choice:
where is the black hole mass. The geometry of is kept fixed at all times and is given by the Schwarzschild values. To be consistent we set the initial conditions
where , and are the Legendre polynomials. In this way the initial ingoing light cone at is distorted with the specification of an arbitrary outgoing massless scalar field. In general, with the nonlinear evolution, the black hole is distorted and the gravitational radiation is generated.
ii.4 Scalar field on a fixed background
The above system of equations (6)–(11) describes a self–gravitating scalar field. In the limit of small amplitudes, , the scalar field can be treated as a perturbation propagating on a fixed background. This considerably simpler model is contained in the fully nonlinear case, and is implemented in our code by integrating only Eq. (11). For a Schwarzschild background, the metric reads
On this fixed background the linear approximation of Eq. (11) is
For the simulations in the present work, we will be interested in solutions of the 2D scalar field on a fixed background, as discussed in the next subsection.
ii.5 QNM in a Schwarzschild background
The linear equation for the scalar field on a fixed background, Eq. (11), is separable; i.e. its solutions can be written in the form
where each of the satisfies the one–dimensional wave equation in the plane ,
This last equation is the usual which govern the scalar perturbations of a Schwarzschild black hole (20), written here in characteristic coordinates . Equation (18) has been studied extensively (20); (21); (22), its most salient feature being the existence of QNM, whose frequencies have been tabulated; see for example Ref. (22). In the present work we will use both the QNM equation, Eq. (18), and the linear Eq. (16), as tests to validate our numerical implementation. We do this in an incremental fashion, solving Eq. (18) for fixed values of , and comparing the effectiveness of the numerical integration scheme and of our boundary conditions in reproducing the QNM. To this end, we implement a purely radial code for Eq. (18) that employs the same numerical integration scheme that is used in the “linear” code [which solves Eqs. (16) and (18)], and in the full nonlinear code. In ingoing null coordinates, the slices at const. penetrate the event horizon at , effectively providing for an excision scheme, where the evolution can be stopped at a finite number of points inside the boundary, because the behavior of the field inside the horizon does not affect the solution outside. Evolutions in ingoing coordinates are carried out on a radial grid, for which boundary data are required at a fixed value of . Because of the presence of this outer boundary, simulations in ingoing coordinates can only be run for a limited time, typically , before outer boundary effects influence the signal extracted.
ii.6 Energy carried out by the scalar field
We calculate the balance of the scalar field energy contained between the inner and the outer boundary (4). The expressions we give here are valid in the linear case, where the background metric is the Schwarzschild metric. For a more general approach to this issue, the linkage integrals have to be calculated; specifically the asymptotic Killing vector field must be parallelly propagated from null infinity (17). Restricted to the background case, then, given a Killing vector field of the metric , , we can define the conserved quantity
In particular, selecting the timelike Killing vector , and for a surface of constant , is the energy contained on the surface,
where is the volume element of the surface at constant . For a sphere at constant , represents the energy flux across the surface,
with the solid angle element. The relevant components the energy-momentum tensor a massless scalar field are
In the case of a linear scalar perturbation on a Schwarzschild background, the energy content of a hypersurface at constant is given by
and the power radiated at time across a surface of constant is
In our simulations we place close enough to the Schwarzschild black hole, and the outer boundary such as . For the flux across the inner (outer) boundary, the integral as well as the spatial and time derivatives are to be taken as evaluated at . With these definitions, the following energy conservation law holds,
The expressions given above hold only in the limit in which is a Killing vector of the metric, so we use them as a criterion for code testing.
ii.7 Marginally trapped surfaces
In general, the constructed spacetime by the present approach contains a distorted black hole. This is made geometrically precise by the introduction of the concept of a marginally trapped two-surface (MTS) on a given ingoing light cone . A MTS is defined, in this context, as the two-parameter radial function on which the expansion of an outgoing null ray pencil vanishes (23). If is tangent to the generators of , we get
Thus, for the diverging slices of , given by , the outgoing normal to is
For the projection tensor into the tangent space of
the expansion associated to the null vector can be written as
which explicitly is
This is an elliptic equation and has to be solved numerically with the system evolution. The convergence to leads to , which locates the MTS.
The q-boundary is the slice with the largest const. on which . Such slice has , therefore is trapped, and trapped surfaces are inside the MTS. In the nonvacuum spherical symmetric case the MTS is given by , which determines the location of the apparent horizon at (corresponding in position with the event horizon for the vacuum Schwarzschild metric). In the absence of spherical symmetry vanishes at points for which . Thus, the q-boundary is everywhere trapped or marginally trapped and is a simple algebraic procedure for locating an inner boundary inside an event horizon.
Iii Numerical implementation
The computational algorithm is related to those developed in Refs. (14), (15) and, as we shall see, was shown to be second order accurate in the nonlinear regime. In the linear regime we get the expected QNM and the energy conservation. We briefly review some issues about the regularization and the discretization of equations.
The coordinate system consists of a radial coordinate and an angular coordinate. The numerical grid is uniformly spaced in both coordinates. Also we use the regularized variables , , . Thus, the hypersurface equations are given by
and the evolution equations
We now review the essentials of the numerical integration procedure.
Equation (33) is easily discretized to get at once
where indexes , and indicate discretization in , and , respectively. In general, any first order radial derivative is calculated as , because we proceed from to (), where is the inner boundary, , and the number of grid points in . Thus, the term is evaluated as
where derivatives of the RHS are evaluated numerically, as indicated above.
Using centered finite differences at , , for stability reasons previously established in Ref. (14), and dictated by the second order derivative in , we discretized Eq. (41); any other discretization in this same equation is staggered at , . Thus, we get
from which we obtain
Next, to discretize Eq. (36) we define the mass aspect
which leads to
In this way we obtain
Observe that the discretization is staggered (backward) for terms with respect to the other terms.
All the obtained formulae for the hypersurface equations are recursive and can be applied from to .
The discretization of the evolution equation (37) proceeds in detail as follows (see Refs. (12), (14), (15)). The core of the evolution integration is the ingoing marching algorithm from to involving the two time levels and .
where and is the RHS of Eq. (37) with the changed sign. Using the mean value theorem, we approximate
where the subscript indicates that the quantity is evaluated at the center of the null parallelogram . Now, easily we can get exactly
On the other hand, we use the conformal invariance (12) to get
Thus, the marching algorithm reads
The numerical implementation of this formula proceeds as follow. Referring again to Fig. 2, interpolations are required, and can be linear to keep a globally second order approximation. The ingoing null geodesic equation is given by
Thus, the displacements and are calculate as
The vertices coordinates are positioned by
and the center coordinate by
The interpolate field at each corner of the null parallelogram is
To get these formulas we have used linear Lagrange interpolations. Now, from Eq. (53) we obtain the following extrapolation formula:
The pending issue to use this formula, that evolves to the most advanced point, is the evaluation of . We do not show details here, but we make two comments in this respect: i) and its derivatives are calculated at (that is, at the center of the null parallelogram), and ; ii) , , , and its derivatives are calculated at , and .
The discretization of the evolution equation (38) for the scalar field proceeds in the same way.
Treatment of the initial-boundary conditions
The boundary conditions (12) are given at the first and second radial grid points to integrate the hypersurface equations. For these two points we implement the algorithm depicted in Fig. 3 to integrate the evolution equations. The first displacement is
The only point in the world tube is at
Thus the interpolation leads us to:
Then we approximate
Observe that is placed at . The second displacement is
to calculate the vertices placed at