Extended two dimensional characteristic framework to study nonrotating black holes

# Extended two dimensional characteristic framework to study nonrotating black holes

## Abstract

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

## I Introduction

Relevant astrophysical applications of the characteristic formulation of numerical relativity w12 () 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 puncture (), 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 p57 (). The LEO code, a large scale computational framework based on the characteristic formulation gbf07 (), was inspired by global forecasting techniques s72 (), rip96 () 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 christo (); hs (); koikeetal (); LRR (). One dimensional numerical simulations of the EKG led to the discovery of critical phenomena c93 () 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 gw92 (). The Bondi mass and news function reflect the discretely self-similar behavior pha05 ().

Extending the work of Gómez et al. gpw94 () and Papadopoulos p02 (), 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 pha05 (), b14 ().

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 Setup

### ii.1 The EKG problem

In general, the field equations for a massless scalar field minimally coupled to gravity are

 Gab=−8πTab, (1)

with

 Tab=∇aΦ∇bΦ−12gab∇cΦ∇cΦ, (2)

and can be reduced to

 Rab=−8π∇aΦ∇bΦ, (3)

which have to be considered together with the wave equation

 □Φ=0, (4)

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 wt65 (); tw66 (), 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 bvm62 ()

 ds2 = (Vre2β−U2r2e2γ)dv2−2e2βdvdr − 2Ur2e2γdvdθ−r2(e2γdθ2+e−2γsin2θdϕ2),

where , , , This metric is twist-free, and therefore rotation is not permitted. Thus we get the following field hypersurface equations

 β,r=12r(γ2,r+ψ2,r), (6)
 (r2Q),r = 2r2[2(γ,rγ,θ+ψ,rψ,θ)+r2(βr2),rθ (7) − 1sin2θ(sin2θγ),rθ],
 U,r=e2(β−γ)Qr2, (8)
 V,r = e2(β−γ)[1+(3γ,θ−β,θ)cotθ+γ,θθ (9) − ψ2,θ−β,θθ−β2,θ−2γ,θ(γ,θ−β,θ)] + r2sinθ[sinθ(rU,r+4U)],θ − 14e2(γ−β)r4U2,r,

and the evolution equations

 − e2β□(2)+(rγ)=−(Vr),rγ+14r3e2(γ−β)U2,r (10) + 1re2(β−γ)(ψ2,θ+β,θθ+β2,θ−β,θcotθ) − 12r[r2(2γ,θU+U,θ−Ucotθ)],r − rsinθ(γ,rUsinθ),θ,
 − e2β□(2)+(rψ)=−(Vr),rψ−1r(r2ψ,θU),r (11) + 1rsinθ[sinθ(e2(β−γ)ψ,θ−r2Uψ,r)],θ,

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:

 β(v,r=rW,θ) = 0, (12a) U(v,r=rW,θ) = 0, (12b) V(v,r=rW,θ) = r−2M, (12c)

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

 γ(v=v0,r,θ)=0 (13)

and

 ψ(v=v0,r,θ)=λre−(r−r0)2/σ2Pℓ(θ), (14)

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

 ds2=(1−2M/r)dv2−2dvdr−r2(dθ2+sin2θdϕ2). (15)

On this fixed background the linear approximation of Eq. (11) is

 2(rψ),vr+[(1−2M/r)(rψ),r],r= 2Mψr2−1rsinθ[sinθψ,θ],θ. (16)

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

 ψ(v,r,θ)=1r∞∑ℓ=0χℓ(v,r)Pℓ(θ) (17)

where each of the satisfies the one–dimensional wave equation in the plane ,

 2χ,vr+[(1−2M/r)χ,r],r=[2Mr3+ℓ(ℓ+1)r2]χ. (18)

This last equation is the usual which govern the scalar perturbations of a Schwarzschild black hole n99 (), written here in characteristic coordinates . Equation (18) has been studied extensively n99 (); ns92 (); k04 (), its most salient feature being the existence of QNM, whose frequencies have been tabulated; see for example Ref. k04 (). 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 gbf07 (). 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 wt65 (). 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,

 E(v)=∫TvvdV, (20)

where is the volume element of the surface at constant . For a sphere at constant , represents the energy flux across the surface,

 P(v)=∫Trvr2dΩ, (21)

with the solid angle element. The relevant components the energy-momentum tensor a massless scalar field are

 Tvv=e−2βψ,r[Vrψ,r−2Uψ,θ]+2r2ψ2,θe−2γ (22)
 Trv=−2e−2βψ,v[Vrψ,r+ψ,v−Uψ,θ] (23)

In the case of a linear scalar perturbation on a Schwarzschild background, the energy content of a hypersurface at constant is given by

 E(v)=2π∫[(1−2Mr)(rψ,r)2+2ψ2,θ]drsinθdθ, (24)

and the power radiated at time across a surface of constant is

 P=−4πr2∫ψ,v[ψ,v+(1−2Mr)ψ,r]sinθdθ. (25)

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,

 Σ(v)=E(v)+∫vv0[Pin(v′)−Pout(v′)]dv′=const. (26)

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 gmw97 (). If is tangent to the generators of , we get

 nα=grvv,α. (27)

Thus, for the diverging slices of , given by , the outgoing normal to is

 lα = r,α−θ,αR,θ (28) − 12grv[grr+gθθR2,θ−2grθR,θ]v,α.

For the projection tensor into the tangent space of

 hαβ=gαβ−nαlβ−lαnβ, (29)

the expansion associated to the null vector can be written as

 Θl = 2hαβ∇αlβ, (30)

which explicitly is

 12e2γr2Θl = R,θ[cotθ+2(β,θ−γ,θ)+U,re2(γ−β)]+R2,θ[2(β,θ−γ,θ)−1r]+R,θθ (31) + r2[U,θ+Ucotθ−Vr2]e2(γ−β).

This is an elliptic equation and has to be solved numerically with the system evolution. The convergence to leads to , which locates the MTS.

As an indicator of the trapped horizon location we can estimate the MTS using the q-boundary, following the method detailed in Ref. gmw97 (). If const., Eq. (31) reduces to

 q≡12e2βr2Θl=r2(U,θ+Ucotθ)−V. (32)

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. gpw94 (), p02 () 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.

### iii.1 Regularization

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

 β,r=r2(^γ2,r(1−y2)2+ψ2,r), (33)
 (r2^Q),r = 2r2{2(1−y2)^γ,r[(1−y2)^γ,y−2y^γ] (34) + 2ψ,rψ,y+r2(β/r2),ry+4y^γ,r − (1−y2)^γ,ry},
 ^U,r=e2[β−^γ(1−y2)]^Qr2 (35)
 V,r = e2[β−^γ(1−y2)]{1+(1−y2)2^γ,yy−2(1−5y2)^γ (36) − 8y(1−y2)^γ,y+2yβ,y−(1−y2)(β,yy+β2,y) − (1−y2)ψ2,y−2(1−y2)[^γ,y(1−y2) − 2y^γ]2+2(1−y2)(^γ,y(1−y2)−2y^γ)β,y} + r2[(1−y2)(r^U,ry+4^U,y)−2yr^U,r−8y^U] − 14e2(^γ(1−y2)−β)r4^U2,r(1−y2)

and the evolution equations

 −e2β□(2)+(r^γ)= − (Vr),r^γ+14r3e2(^γ(1−y2)−β)^U2,r (37) + 1re2(β−^γ(1−y2))(ψ2,y+β,yy+β2,y) − r[(^γ,ry^U+^γ,r^U,y)(1−y2) − 4y^γ,r^U]−(^γ,y(1−y2)−2y^γ)2^U − ^U,y−r2[(^γ,yr(1−y2)−2y^γ,r)2^U + (^γ,y(1−y2)−2y^γ)2^U,r + ^U,yr)],
 −e2β□(2)+(rψ)= − (Vr),rψ (38) − 2(1−y2)ψ,y^U − r(1−y2)[ψ,ry^U+ψ,y^U,r] + 1r{(1−y2){2[β,y−^γ,y(1−y2) + 2y^γ]ψ,ye2[β−^γ(1−y2)] + e2[β−^γ(1−y2)]ψ,yy − r2(^U,yψ,r+^Uψ,ry)} − 2y(ψ,ye2[β−^γ(1−y2)] − r2^Uψ,r)}.

We now review the essentials of the numerical integration procedure.

### iii.2 Discretization

#### Hypersurface equations

Equation (33) is easily discretized to get at once

 βnj,i=βnj,i−1−Δrβ,r|nj,i−1/2, (39)

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

 β,r|nj,i−1/2=12ri−1/2(ψ2,r+^γ2,r)|j,i−1/2, (40)

where derivatives of the RHS are evaluated numerically, as indicated above.

Now, combining Eqs. (34) and (35) we get

 C^U,r+12r2^U,rr=H^Ue2(β−γ), (41)

where

 C=r2[2r+^γ,r(1−y2)−β,r],
 H^U=2(1−y2)^γ,r[(1−y2)^γ,y−2y^γ]+¯ψ,rψ,y+¯ψ,yψ,r +r2(βr2),ry+4y^γ,r−(1−y2)^γ,ry.

Using centered finite differences at , , for stability reasons previously established in Ref. gpw94 (), 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

 −Cnj+1/2,i−1/22Δr(^Unj,i−^Unj,i−2) +r2i−1/22Δr2(^Unj,i−2^Unj,i−1+^Unj,i−2) =[H^Ue2(β−γ)]nj+1/2,i−1/2, (42)

from which we obtain

 ^Unj,i = 1fa{Δr2[H^Ue2[β−^γ(1−y2)]]nj+1/2,i−1/2 (43) − (fb^Unj,i−1+fc^Unj,i−2)},

where

 fa=12(r2i−1/2−ΔrCnj+1/2,i−1/2),
 fb=−r2i−1/2,
 fc=12(r2i−1/2+ΔrCnj+1/2,i−1/2).

Next, to discretize Eq. (36) we define the mass aspect

 M(v,r,θ)=12(r−V). (44)

 M,r=12(1−HM^U−HMβ^γψ), (45)

where

 HM^U = r2[(1−y2)(r^U,ry+4^U,y)−2yr^U,r−8y^U] (46) − 14e2(^γ(1−y2)−β)r4^U2,r(1−y2),
 HMβ^γψ = e2[β−^γ(1−y2)]{1+(1−y2)2^γ,yy−2(1−5y2)^γ (47) − 8y(1−y2)^γ,y+2yβ,y−(1−y2)(β,yy+β2,y) − (1−y2)ψ2,y−2(1−y2)[^γ,y(1−y2) − 2y^γ]2+2(1−y2)(^γ,y(1−y2) − 2y^γ)β,y}.

In this way we obtain

 Mnj,i = Mnj,i−1 (48) − Δr(1−HM^U|nj−1/2,i−1−HMβ^γψ|nj,i−1/2).

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 .

#### Evolution equations

The discretization of the evolution equation (37) proceeds in detail as follows (see Refs. gw92 (), gpw94 (), p02 ()). The core of the evolution integration is the ingoing marching algorithm from to involving the two time levels and .

We integrate Eq. (37) over the null parallelogram formed by ingoing and outgoing radial null rays in the plane that intersect at vertices , , and , as depicted in Fig. 2. Thus, we have

 (49)

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 gw92 () to get

 ∫Ae2β□(2)+~γdvdr=2(~γQ−~γP+~γR−~γS). (52)

 ~γQ=~γP+~γS−~γR+14Δv(rQ−rP+rS−rR)HC. (53)

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

 drdv=12(1−2Mr). (54)

Thus, the displacements and are calculate as

 δi=12⎡⎣1−⎛⎝Mn−1j,i+1ri+1+Mnj,i−1ri−1⎞⎠⎤⎦Δv (55)

and

 δi−1=12⎡⎣1−⎛⎝Mn−1j,iri+Mnj,i−2ri−2⎞⎠⎤⎦Δv. (56)

The vertices coordinates are positioned by

 rP = ri−1+12δi−1, (57a) rR = ri−1−12δi−1, (57b) rQ = ri+12δi, (57c) rS = ri−12δi, (57d)

and the center coordinate by

 rC=12(rP+rS). (58)

The interpolate field at each corner of the null parallelogram is

 ~γj,P = ~γnj,i−1−12δi−1Δr(~γnj,i−1−~γnj,i−2), (59a) ~γj,R = ~γn−1j,i−1+12δi−1Δr(~γn−1j,i−~γn−1j,i−1), (59b) ~γj,S = ~γn−1j,i+12δiΔr(~γn−1j,i+1−~γn−1j,i), (59c) ~γj,Q = ~γnj,i−12δiΔr(~γnj,i−~γnj,i−1). (59d)

To get these formulas we have used linear Lagrange interpolations. Now, from Eq. (53) we obtain the following extrapolation formula:

 ~γnj,i=11−δi/2Δr{~γj,P+~γj,S−~γj,R+14Δv(rQ−rP+rS−rR)HC−12δiΔr~γnj,i−1}. (60)

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

 δ0=12⎛⎝1−2Mn−1j,0r0⎞⎠Δv. (61)

The only point in the world tube is at

 rS=r0−δ0. (62)

Thus the interpolation leads us to:

 ~γj,S=~γn−1j,0+δ0Δr(~γn−1j,1−~γn−1j,0). (63)

Then we approximate

 ~γnj,0=~γj,Q←~γj,S. (64)

Observe that is placed at . The second displacement is

 δ1=12⎛⎝1−2Mn−1j,1r1⎞⎠Δv. (65)

to calculate the vertices placed at

 rP = r0, (66a) rQ = r1, (66b) rR = r0−δ0