Conclusions

This chapter has presented the work related to FO-CCZ4, a new formulation of the Einstein equations based on the conformal and covariant Z4 (CCZ4) system ofAlic et al.(2012), which is strongly hyperbolic and first order in both space and time. The system consists of 58 evolution equations for the vector of state variable given by

Q^T :=

˜

γij,lnα, βⁱ,lnφ,A˜ij, K,Θ,Γˆⁱ, bⁱ, Ak, Bⁱ_k, Dkij, Pk

.

To the best of our knowledge, this is the first time that a first-order, strongly hyperbolic formulation of the CCZ4 system has been proposed. The chapter also summarizes the results of numerical tests conducted with the FO-CCZ4 formulation when discretized with a path-conservative ADER-DG numerical method.

The guiding principle in order to obtain a strong hyperbolic formulation of the equation has been the approximate symmetrization of the sparsity pat-tern of the system matrix, via the appropriate use of various typologies of con-straints, with the goal of avoiding the appearance of Jordan blocks that cannot be diagonalized. A second technique employed to obtain the FO-CCZ4 formu-lation as presented in this chapter is the use of first-order ordering constraints in such a way to reduce the evolution equations for the lapseα, the shiftβⁱ, the conformal metric˜γij and the conformal factorφto a system of pure ordinary differential equations (as opposed to PDEs) (Alcubierre,2008). In other words, whenever differential terms with respect toα,βⁱ,φandγ˜ijappear, they are re-placed by the corresponding auxiliary variablesAk,Bⁱ_k,PkandDkij, thus be-coming algebraic source terms. This choices greatly simplified the hyperbolic-ity analysis of the FO-CCZ4 system, since the eigenvalues and eigenvectors as-sociated withα,βⁱ,φand˜γijbecome trivial. A further advantage is that for the rest of the analysis a reduced system of partial differential equations relative to only 47 dynamic variables (namelyU = ( ˜Aij, K,Θ,Γˆⁱ, bⁱ, Ak, Bⁱ_k, Dkij, Pk)) can be considered. What’s more the system matrix of the reduced system is a function ofα,βⁱ,φand˜γijonly, which not only substantially simplifies the hy-perbolicity analysis but also leads to the important result that all fields of our FO-CCZ4 system are linearly degenerate. This in turn implies that no shock waves can be generated from smooth initial data by evolving it with the FO-CCZ4 system.

In contrast with the first-order Z4 system proposed in Bona et al.(1997);

Alic et al.(2009), the FO-CCZ4 system is written in a fully non-conservative form, which is another key idea of the present approach. It’s worth pointing out that the above mentioned simplifications are not possible if a conserva-tive formulation of the system is sought,e.g. the ones proposed inBona et al.

(1997);Alic et al.(2009). This follows by the fact that the Jacobian∂F/∂Qof the fluxF(Q)would in general depend also on the dynamical variables and the quasi-linear form of the system would also contain differential terms inα, βⁱ,φand˜γij. This by construction not the case in the present non-conservative formulation.

Thanks to the choices outlined above we have been able to provide a proof of strong hyperbolicity for completely general lapse and spatial metric for the simple gauge choiceβⁱ = 0, by a direct explicit calculation of all eigenvalues

and all left and right eigenvectors of the system. While it would be desirable to prove hyperbolicity also for more realistic choices of the shift condition,e.g.

for the Gamma-driver gauge, at the moment we have only reached a proof for general lapse and shift conditions and general spatial metrics only as long as one of the three independent off-diagonal components of the three-metric is zero. We conjecture that the formulation is strongly hyperbolic in all cases, but a general analytical proof is at present missing and left to future work.

We have implemented numerically the FO-CCZ4 system via a discretiza-tion based on a family of high order fully-discrete one-step ADER discontin-uous Galerkin (DG) schemes, coupled with an ADER-WENO finite-volume limiter (the latter being necessary in order to deal with the physical singular-ities arising in the case of black hole spacetimes). The non-conservative na-ture of the formulation is naturally treated within the DG framework by the use of path-conservative schemes, first proposed by Castro and Parés in the finite-volume context (Castro et al.,2006;Pares,2006) and later extended also to ADER-DG schemes inDumbser et al.(2009,2010). In order to ensure the positivity of the numerical solution ofαandφ, we have evolve the logarithms of these quantities.

Following a well established practice, we have applied the FO-CCZ4 for-mulation to a series of standard test cases, i.e. a subset of the Apples-with-Apples tests suggested in Alcubierre et al. (2004b), namely the gauge-wave test, the robust stability test and the linear-wave test bed. Besides providing ev-idence that the new system is able to reproduce accurately and with moderate costs the known analytic solutions, we have carried out numerical convergence studies of the method on the gauge-wave test in the highly nonlinear regime, as well as on further tests involving Schwarzschild and Kerr black holes using 3D Cartesian Kerr-Schild coordinates. We have also provided numerical evi-dence that the FO-CCZ4 formulation coupled with ADER-DG schemes (with the ADER-WENO finite-volume subcell limiter) is able to perform a long time integration of a single puncture black hole with the usual Gamma driver and 1 + loggauge conditions. Finally we have also shown some first preliminary results for two moving puncture black holes. To the best of our knowledge, the numerical results shown here represent the first simulations of the 3+1 Einstein equations ever done with high order DG and WENO finite-volume schemes on three-dimensional grids.

Future research will concern the extension of the present algorithms to dy-namic AMR (adaptive mesh refinement), as well as the extraction of the gravi-tational waveforms generated by binary black-hole mergers (seeCentrella et al.

(2010);Bishop and Rezzolla(2016) for reviews) and binary neutron-star merg-ers (seeBaiotti and Rezzolla(2017) for a review), and also the characterization of black hole horizon properties.

High-order numerical schemes for relativistic

hydrodynamics

95

The Entropy Limited

Hydrodynamics Scheme

5.1 Introduction

In order to study astrophysical systems involving compact objects, such as black holes and neutron stars, large-scale general-relativistic hydrodynamical numerical simulations have been shown to be a very powerful tool (Font,2008;

Shibata and Taniguchi,2011; Rezzolla and Zanotti, 2013; Martí and Müller, 2015;Shibata,2016;Baiotti and Rezzolla,2017;Paschalidis,2017). Performing such simulations however is a very non-trivial task, which requires dealing with a plethora of different physical, mathematical and computational issues.

One of the most challenging of such issues, which can lead to significant dif-ferences on the outcome of said simulations, especially when the resolution employed is not very high, is the choice of the numerical method for which is employed in the solution of the relativistic hydrodynamics equations.

As already mentioned in chapter3, in this context the most commonly used methods are generally known as high-resolution shock-capturing (HRSC) tech-niques. HRSC methods have been shown in general to be very effective in take care of shocks waves and suppressing spurious oscillations in the numerical solution of PDEs, and have been employed with varying degree of success in astrophysical simulations. Recently much effort has gone into improving these schemes (e.g. by employing innovative mesh refinement techniques such as inDeBuhr et al.(2015)) or moving beyond them; one promising and popular alternative is that of discontinuous Galerkin (DG) methods, which have been introduced in section3.4and used in obtaining the results of chapter4. Both

“standard” HRSC schemes and their improvements however potentially suf-fer from a few shortcomings. First they are in general complex to derive and implement, or to extend and modify (e.g.in order to increase the formal order of accuracy. This does not apply to DG methods however); they often depend on a large number of a priori unknown coefficients, requiring some degree of optimisation (e.g.a typical example being WENO methods); they may lead to load imbalance in parallel implementations as a result of their complexity.

In this chapter an alternative approach, different from the HRSC mindset, is proposed, able to address some of these shortcomings. These alternative

97

scheme is named “entropy-limited hydrodynamics” (ELH) and we formulate it in a finite-differences framework. This is a variant of the “flux-limiting” FD scheme described in chapter3, where the underlying concept prescribing how the limiter should be activated and driven is relatively straightforward: to de-termine which gridpoints are in need of the low-order contribution, we employ a “shock detector”, which not only marks regions of the computational domain in need of limiting, but also determines the relative ratio of the high and low-order fluxes.

Such a shock detector is offered by the entropy viscosity function described by Guermond et al. (we refer primarily to Guermond et al.(2011), but see also Guermond and Pasquetti(208);Zingan et al.(2013)), in which the local production of entropy is used to identify shocks. Since entropy is produced only in the presence of shocks, this choice results in a stable method, which is nonetheless able to recover high-order in regions of smooth flow. We have extended the definition of the entropy viscosity function from the classical to the relativistic regime, and rather than a prefactor to additional viscous terms in the hydrodynamical equations, we choose to employ it to drive the lower-order flux in the flux limiting scheme. Therefore in contrast to the approach of Guermond et al.(2011) the underlying equations of relativistic hydrodynamics are not modified in this approach by introducing additional entropy-related terms.

In the following the method itself and the details of our implementation are described, followed by a summary of the results of tests we conducted in order to gauge its behaviour against a standard HRSC method, namely the MP5 scheme (Suresh and Huynh,1997). This chapter is structured as follows:

in sections 5.2 and 5.3the ELH method and its present implementation are described, while the results of the numerical tests are presented in section5.4.

Conclusions and an outlook are collected in section5.5.