The finite volume method on a Schwarzschild background

We introduce a class of nonlinear hyperbolic conservation laws on a Schwarzschild black hole background and derive several properties satisfied by (possibly weak) solutions. Next, we formulate a numerical approximation scheme which is based on the finite volume methodology and takes the curved geometry into account. An interesting feature of our model is that no boundary conditions is required at the black hole horizon boundary. We establish that this scheme converges to an entropy weak solution to the initial value problem and, in turn, our analysis also provides us with a theory of existence and stability for a new class of conservation laws.


Introduction
We design and study a finite volume scheme for a class of nonlinear hyperbolic equations posed on a Schwarzschild black hole background. This paper is the follow-up of earlier investigations by LeFloch and co-authors [2,10,13,14]. As is common in the mathematical theory of hyperbolic balance laws, we consider a (drastically) simplified version of the compressible Euler equations and we describe the fluid evolution by a single scalar unknown function, typically representing the velocity of the fluid. For relativistic problems the velocity is naturally bounded and, after normalization, we seek for solutions v : M → [−1 , 1] (1.1) defined on a "spacetime" M (explicitly described below in a global coordinate chart) and satisfying the following hyperbolic balance law Here, X α = X α (w, ·) is the so-called flux vector field parametrized by the real variable w ∈ [−1, 1] and defined on M, while q = q(w, ·) is a prescribed real-valued function. Structural conditions (even for smooth solutions, as specified later in this text) must be imposed on the vector field in order for the balance law to admit a well-posed initial value formulation.
Keywords and phrases. Hyperbolic conservation law, Schwarzschild black hole, weak solution, finite volume scheme, convergence analysis.
Our objectives in this paper are as follows: • Choosing M to be (the outer domain of communication of) a Schwarzschild black hole, we introduce a class of hyperbolic balance laws (1.2) and formulate the associated initial value problem. We then seek for weak solutions v : M → [−1, 1] possibly containing shock waves which must satisfy a suitable entropy condition (discussed below). • Next, we design a finite volume scheme that allows us to numerically approximate these weak solutions and we derive several fundamental properties of interest: maximum principle, entropy inequalities, etc. We establish the strong convergence of this scheme toward a weak solution of the initial value problem.
Our arguments are based on a generalization of DiPerna's theory of measure-valued solutions [6] and require us to cope with the effects of the curved black hole geometry.
An outline of this paper is as follows. In Section 2 we introduce the class of hyperbolic equations of interest and provide a motivation from pressureless fluid dynamics. In Section 3, we analyze the geometry of the curved characteristics in the black hole geometry and the class of steady state solutions which represent a fluid at rest. In Section 4, we discuss an alternative choice of slicing and which illustrate how the balance gets transformed under change of coordinates. In Section 5, we introduce our finite volume scheme and state the convergence theory. The entropy inequalities satisfied by the weak solutions and their discrete version are also derived, and the proof of convergence is completed.

The choice of coordinates
The domain of outer communication of a Schwarzschild black hole, denoted by M, can be described in the so-called Schwarzschild coordinates x = (t, x j ) = (t, x 1 , . . . , x n ) in which the spacetime metric reads Here, the time variable t and the radius r defined by The light speed is normalized to unit while the parameter M ∈ [0, +∞) represents the mass of the black hole. Moreover, g S 2 denotes the canonical metric on the unit (n − 1)-sphere S n−1 ⊂ R n . The spacetime hypersurface is the boundary of our spacetime and represents the horizon of the black hole, from which nothing can propagate in the (outer communication) domain r > 2M of interest. Recall that the apparent singularity at r = 2M in the expression of the metric (2.1a) is not a physical singularity but is solely due to our choice of coordinates.

The model of interest
Choosing the vector field in the left-hand side of the balance law (1.2) to be , 0, . . . , 0 and the source term to be we arrive at the following hyperbolic balance law: Here, the functions f = f (w) and h = h(w) are prescribed functions, while the unknown scalar field is v : R + × Ω → [−1, 1], defined for all t ≥ 0 and r ≥ 2M , and we work in the exterior of the ball with radius 2M , that is In our model the unknown v need not be spatially symmetric, so it convenient to rewrite (2.2) in Cartesian coordinates, i.e.
Finally, in order to eliminate the singularity 1 1−2M/r , we propose an equivalent form, as follows.
is referred to as a hyperbolic balance law on a Schwarzschild black hole.
At this juncture, it should be emphasized that further conditions (presented in Sect. 3) will be required on the flux function f in order for the interval [−1, 1] to be an invariant domain.
We always tacitly assume that an entropy U is normalized to satisfy U (0) = 0. Then, by definition an entropy solution to the equation (2.5) must satisfy, for all convex entropy pairs (U, F ), We prescribe an initial data v 0 at the time t = 0, that is, and we formalize our notion of solution as follows.
for all convex entropy pairs (U, F ) and all compactly supported test-functions φ ≥ 0.

Derivation from the relativistic Euler system
Our motivation for introducing the above class of balance laws comes from a formal derivation made from the Euler equations for a relativistic compressible fluid, which read ∇ α T αβ (ρ, u) = 0, (2.10) in which ∇ denotes the Levi-Civita connection associated with the Schwarzschild metric (2.1a). We are interested here in the energy-momentum tensor of a pressureless fluid, given by where ρ : M → (0, +∞) denotes the density of the fluid and the velocity field u = (u α ) is normalized to be future-oriented, unit and timelike u α u α = g ββ u α u β = −1 with u 0 > 0 and, therefore, By assuming spherical symmetry, we can derive from the above system a single equation satisfied by a suitably normalized component of the velocity field, denoted below by v ∈ (−1, 1).

Maximum principle
The method of characteristics allows us to obtain a first insight about the properties of (sufficiently regular solutions) to our balance law (2.4). It leads to ordinary differential equations along characteristic curves parametrized with respect to some time parameter (denoted by s below). We would like to deduce some properties of solution v by proposing the following assumption on the flux f and the source h.
We motivate our condition by the following analysis along characteristic curves. So, we consider the coupled system A straightforward computation shows that This equation tells us how the values of a solution evolve along characteristics, and we use it in order to establish a maximum principle. It is convenient to assume a strict inequality in the data. sup as long as they remain sufficiently regular.
Proof. Observe first that if v 0 = ±1 initially then it remains so for all times. It is sufficient to show that v ≤ 1, since exactly the same arguments apply to showing v ≥ −1.

Geometry of the characteristic curves
Along a characteristic we see that and, more explicitly, where r 0 = r(s 0 ) and u 0 = u(s 0 ) are given data at some time s 0 .
Concerning the global behavior of the characteristics in the special case f (w) = w 2 /2 − 1/2 and h(w) = 0, which is the Burgers equation posed on the Schwarzschild background, the weak solutions in the (u, r) plane can be expressed in terms of the initial data via a minimisation formulation based on characteristics; see [2].
Here, to proceed with the study of the characteristic curves and for the sake of definitness, we assume some specific signs about the functions f and h.
and, by solving for u, the ordinary differential equation (3.4a) for the radius function r(s) can be written as (3.10) Here, F are the inverse functions of F + and F − , respectively, and Note that F ± are single-valued functions within the domain of interest. We follow [2] and introduce the escape velocity (whenever it exists) which satisfies the property Replacing the radius r 0 by the escape velocity parameter u E 0 in (3.9), we obtain 14) The late-time behavior of u = u(s) can be checked to be described as follows: • Negative initial data. The function u(s) decreases (as follows from (3.3) and the assumption (3.6)). If u 0 ∈ (−1, 0] with initial data (s 0 , r 0 ), then u(s) remains negative and decreasing and r(s) decreases towards 2M . More precisely, we have The positivity of u 0 initially ensures dr/ds > 0, that is, the characteristic curve initially moves away from the black hole. However, u(s) keeps decreasing and eventually reaches 0 at some time s 0 . The dynamics then coincides with that for negative initial data. We conclude that lim • Positive initial data with u 0 ≥ u E 0 . In this case, the characteristic curve moves away from the black hole for all times and the asymptotic behavior is

Steady state solutions
Finally, let us consider solutions that are steady states representing a fluid at rest in the curved black hole geometry. This is a special class of solutions of interest, for instance, in designing (well-balanced) numerical schemes and in finding test cases. In view of the radially symmetric form of our equation (2.16) (but possibly for non-radially symmetric solutions), for a time-independent solution we obtain the ordinary differential equation and, once again, we get the same ordinary differential equation as (3.4b) Under Assumptions (3.1) and (3.3), for any given data (r 0 , u 0 ) we can distinguish between two cases: • Negative u 0 . Then u is increasing and • Positive u 0 . Then u is decreasing and

An alternative choice of time slicing
In this section, we illustrate the fact that coordinates can be chosen in many different manners. While, for Schwarschild spacetime, this leads to significantly more involved algebraic expressions, such alternative coordinates may allow one to cover a larger region of the spacetime. For definiteness, in this section we take n = 3. Hence, we now introduce a nonlinear hyperbolic equation posed in a larger domain of the Schwarzschild geometry, obtained by "crossing" the horizon and we study the interior of the black hole. We follow [5] and introduce the following metric: in which t denotes the time variable and R the radial variable with Here R 0 ∈ (0, M ] is a parameter that is fixed, and we observe that the above expression is identical to the metric (2.1a) in the limit R 0 → 0, for which the radial variables R and r would then coincide. This new slicing (4.1a) allows us to go inside of the black hole (when R 0 > 0), and we cover the region {r : r + R 0 − 2M > 0}, within which the metric remains of a definite Lorenztian signature.
In fact, we can transform (4.1a) (for a restricted domain of the variables, only) into the metric (2.1a), by setting In the following, it will be convenient to rely on the vector fields We rewrite (4.1a) in the matrix form After a tedious computation, the (non-vanishing) Christoffel symbols Γ µ αβ = 1 2 g µν ( ∂ α g βν + ∂ β g αν − ∂ ν g αβ ) are found to be

Characteristics and maximum principle
As we mentioned in the beginning of this section, the new metric g coincides with the Schwarzchild metric g when r is replaced by R = r + R 0 in (2.1a), hence it is not surprising to have the following result. Namely, if we replace r by R throughout Section 2, then Burgers equation (2.15) is equivalent to (4.6). This is easy to check with From the equation (2.15) we have  Similarly to what we did in Proposition 3.2, we can check the following result.
which follows from (4.13) and (4.14). From this, we obtain where ( u 0 , R 0 ) is the initial location. Thus, our previous conclusions concerning the characteristic curves are recovered here. This is of course not surprising, since we are treating the same differential equation expressed in different coordinates. This second formulation however may have some numerical advantage when the horizon, instead of being fixed as it is in the present model, is dynamical.

Formulation of the finite volume scheme
Having considered the formulation (2.5) (in Sect. 2) and the formulation (4.6) (in Sect. 4), we now study the numerical approximation of the general balance law (1.2) in a setting that, in principle, may encompass both formulations. For definiteness, we treat the outer domain of communication so that the space variable varies in a half-line and no boundary condition is required at the boundary. The hyperbolic model of interest reads in which M denotes the outer domain of communication of a Schwarzschild black hole with radius 2M , as we described earlier. Here X α = X α (w, ·) is a smooth vector field on M, depending upon the real variable w ∈ [−1, 1]. An hyperbolicity condition and a condition at the boundary will be made explicit below. We are going to formulate a finite volume scheme for the equation (5.1) and establish its convergence by generalizing the technique of proof in [1]. In contrast with this later work, the spacelike slices in M are noncompact and the flux vector X α (v, x) is no longer assumed to be geometry compatible (i.e. X α (v, x) does not satisfy divergence free condition), and at the (horizon) boundary, no boundary data is needed.
In the class of interest in the present paper, the flux vector field satisfies the following property, which implies that no boundary condition is needed: the spatial components of the vector field X(·, x) vanishes on the boundary, i.e. X a (·, x) = 0 on the boundary ∂M.
Following [1] we design a finite volume scheme for (5.1) as follows. We introduce a spacetime triangulation T h = K∈T h K of M such that the boundary ∂K of each element K is the union of three possible types of faces: • A face e + K is spacelike and we denote its future-oriented outward unit normal by n K,e + K .
• A face e − K (in the past of e + K ) is also spacelike, and we denote its past-oriented outward unit normal by n K,e − K .
• A vertical face denoted by e 0 is timelike, and whose inward unit normal is denoted by n K,e 0 , and the union of all of such faces is denoted by ∂ 0 K := ∂K \ {e + K , e − K }. By definition, for every pair of distinct elements K, K ∈ T h , the intersection K ∩ K is either a common face of K, K or a submanifold of co-dimension at least 2. We use K e to denote the unique neighbor of K sharing the same edge e. We denote by K ± the neighbors of K which share the same edge e ± . We also write K e 0 for the element that shares the same edge e 0 . Furthermore, the notation |e ± K |, |e 0 |, |K|, etc. stands for the volume of the corresponding set.
By integrating the equation (5.1) over an arbitrary element K ∈ T h and applying the divergence theorem, we obtain Here, n denotes the exterior and unit, normal vector field along the boundary face under consideration, while dV is the induced measure element on the boundary. Our finite volume scheme is based on the following approximation formulas, in which e 0 , etc. denotes an edge of K: • Discretization of the main variable: • Discretization of the flux: • Discretization of the source term: Here, the numerical flux f K,e : R 2 → R is chosen to satisfy the properties of consistency, conservation and monotonicity: • Consistency property: g(X(v, p), n K,e 0 (p)) dV (p), v ∈ R. (5.5a) • Conservation property: • Monotonicity property: and µ X K,e (v) := 1 |e| e g(X(v, p), n K,e (p)) dV e . (5.7) Finally, the finite volume approximations are defined by

Convergence and existence theory
Based on the geometric formulation of a finite volume method above, we can now proceed with the analysis of our model problem (2.5). We integrate (2.5) over an element K and by applying the divergence theorem where n denotes the outward unit normal vector. We apply the following approximations: • Discretization of the main variable: • Discretization of the flux: tn+1 tn e • Discretization of the source term: where, r K being the radial variable evaluated at the center of K, In the above, we denoted by e some edge of K, and the numerical flux f K,e : R 2 → R is chosen to satisfy the properties of consistency, conservation and monotonicity, that is, in our case • Consistency property: • Conservation property: • Monotonicity property: The finite volume approximations are then given by the explicit scheme v n+1 For the sake of stability, we impose the CFL stability condition p K := e∈∂K |e| being the perimeter of K, as well as the source stability condition (5.14) Now we are ready to state our convergence result.
Let f K,e be a family of numerical flux satisfying the consistency, conservation and monotonicity conditions in (5.11a)-(5.11c) and satisfies the CFL condition (5.13) and the stability condition (5.14). Then the discrete scheme (5.12) uniquely defines the family of approximate solution v n K . By defining a piecewise constant function v h :

Discrete entropy inequalities
Entropy inequalities play a key role in the proof of Theorem 5.1. Proof. We first assume max T h |v n K | ≤ 1 for all elements K ∈ T h . We rewrite (5.12) as v n+1 which with an obvious notation we express in the form In view of (5.12) and the consistency property of f K,e , we have v n+1 The following lemma provides a standard result concerning the existence of discrete entropy flux terms and an entropy inequality relating v n+1 K,e and v n+1 K . We omit the proof and refer to [1] and the references therein.
Lemma 5.4 (Discrete entropy inequalities). Let (U, F ) be a convex entropy pair. Then there exists a family of discrete entropy flux functions F K,e : R 2 → R satisfying the following conditions: • Consistency with the entropy flux F : • Conservation property: F K,e (u, w) = F Ke,e (w, u), u, w ∈ R. (5.20b) • Discrete entropy inequality: Equivalently, (5.20c) can be written in terms of v n+1 K,e and v n K as K,e ). The entropy dissipation estimate below will serve to establish the convergence result.
Proof. By Lemma 3.5 in [4] and the convex decomposition identity (5.19c), we have Next, we multiply by |e||K|/p K in (5.21), and sum up over all K and e, This leads us to (5.22).
For the proofs of the following lemmas, see [1] and the references therein for details. The local entropy inequalities read as follows.
The global entropy inequalities read as follows.

Measure-valued solutions and strong convergence
We are now in a position to complete our proof of Theorem 5.1. Based on the entropy inequalities we have established, we are able to pass the limit in the inequality (5.24) as h → 0. Then we associate with a subsequence of v h (which is uniformly bounded in [0, T ) × Ω → R for fixed T ) a Young measure ν : [0, T ) × Ω → Prob(R), which is a family of probability measures in R parametrized by (t, x) ∈ [0, T ) × Ω. We then show that the Young measure, describing all the weak-star limits of v h , is an entropy measure-valued solution in the sense of DiPerna. The strong convergence result follows from the DiPerna's uniqueness theorem, see [6].
The Young measure allows us to write, for every continuous function a : R → R, a(v h ) ν, a as h → 0, (5.25) in the L ∞ weak-star topology. As presented in [1], it suffices to show that ν is an entropy measure-valued solution to our balance law, in order to imply that ν t,x reduce to a Dirac mass δ v(t,x) if this is true at the initial time t = 0. The convergence in (5.25) then holds in a strong sense and v h converges to the entropy solution v to the Cauchy problem. ν t,x , U ∂ t φ(t, x) + ν t,x , F ∂ j 1 − 2M r x j r φ(t, x) for all non-negative test functions φ : [0, T ) × Ω → R + .
For all convex entropy pairs, we thus have (5.27) and the proof of of Theorem 5.1 is completed.