# Entropy based flux limiting scheme for conservation laws

Georgios Doulis1,2 Sebastiano Bernuzzi2 Wolfgang Tichy3 1Institut für
Theoretische Physik, Goethe-Universität Frankfurt, 60438 Frankfurt am Main
2Theoretisch-Physikalisches Institut, Friedrich-Schiller-Universität Jena,
07743 Jena  3Department of Physics, Florida Atlantic University, Boca Raton,
FL 33431

###### Abstract

The construction of high-resolution shock-capturing schemes is vital in
producing highly accurate gravitational waveforms from neutron star binaries.
The entropy based flux limiting (EFL) scheme is able to perform fast
converging binary neutron star merger simulations reaching up to fourth-order
convergence in the gravitational waveform phase. Here, we extend the
applicability of the EFL method beyond special/general relativistic
hydrodynamics to scalar conservation laws and show how to treat systems
without a thermodynamic entropy. This is an indication that the method has
universal applicability to any system of partial differential equations that
can be written in conservation form. We also present some further very
challenging special/general relativistic hydrodynamics applications of the EFL
method and use it to construct eccentricity reduced initial data for a
specific neutron star binary and show up to optimal fifth-order convergence in
the gravitational waveform phase for this simulation.

###### pacs:

04.25.D-, 04.30.Db, 95.30.Sf, 95.30.Lz, 97.60.Jd

##  I Introduction

Binary neutron star (BNS) coalescences are one of the confirmed [1, 2] sources
of gravitational waves (GW) detected by the LIGO-VIRGO interferometers.
Numerical relativity (NR) is able to produce merger waveforms that replicate
these GWs accurately enough for BNS parameter estimation. The production of
highly accurate numerical waveforms plays a vital role in understanding the
properties of neutron star (NS) binaries during the inspiral, merger and
postmerger regime. In this spirit, several high-resolution shock-capturing
(HRSC) schemes [3] have been developed over the years aiming in increasing the
accuracy and the convergence properties of the numerically produced waveforms.
Going beyond second-order convergence is currently the major challenge for the
existing HRSC schemes, which has been pointed out long ago, but addressed only
in a few places [4, 5, 6, 7]. Here, we focus on a specific family of flux
altering HRSC methods: the flux-limiters [8], which have been proven very
successful in dealing with shocks and suppressing spurious oscillations.

In [9] a fast converging HRSC scheme that utilises the local entropy
production in the design of an entropy based flux-limiter (EFL) was proposed.
The EFL scheme developed from an entropy method that started as an artificial
viscosity method [10, 11] and grew with time to a flux-limiting scheme [12].
In [12], for the first time, entropy was used in the design of a flux-limiter.
The EFL scheme, apart from overcoming the shortcomings of the method developed
in [12], is considerably more automated and simpler in its implementation.
These features make the EFL method universally applicable to different special
and general relativistic scenarios. The EFL scheme was successfully used in
[9] to evolve several different BNS systems and managed to deliver up to
fourth-order convergent waveforms in the dominant $(2,2)$-mode and in the
subdominant $(3,2)$\- and $(4,4)$-mode.

In the present work we extend the applicability of the EFL method to scalar
conservation laws without a thermodynamically defined entropy. To do so, an
auxiliary entropy functional has to be defined [11], which can be used to flag
the presence of non-smooth features in the solution space. This entropy can be
used, like in [9], to define the local weights entering the usual definition
of the entropy flux-limiter. The weights combine an unfiltered high-order flux
used in regions of smooth flow with a stable low- or high-order flux used in
regions of non-smooth flow.

We also use the EFL method for the first time in the construction of BNS
initial data. (Notice that now we have a well defined thermodynamic entropy,
therefore the use of the well established EFL machinery of [9] is more than
sufficient.) In [9], the EFL scheme was exclusively used for the temporal
evolution of initial data computed with the Lorene library [13]. Here, we
involve the EFL method in the construction of eccentricity reduced initial
data with the SGRID code [14] and study the effect of this action on the
accuracy and the convergence properties of the produced merger waveforms.

The article is structured as follows. Theoretical and numerical aspects of the
EFL method are discussed in Sec. II and the details of its numerical
implementation in Sec. III. The first application of the EFL method to scalar
conservation laws is presented in Sec. IV. Sec. V includes our results for
some standard benchmark tests of special relativity and in Sec. VI the
performance of our method is tested against three-dimensional general
relativistic non-static/stationary single NS configurations. Sec. VII contains
our results of the first BNS simulation where the EFL method is not only used
for the temporal evolution, but also in the construction of the initial data.
Finally, we conclude in Sec. VIII.

Throughout this work we use geometric units. We set $c=G=1$ and the masses are
expressed in terms of solar masses $M_{\odot}$.

##  II EFL Method

In the present work, the entropy flux limiting scheme developed in [9] is
extended to systems without a thermodynamically defined entropy. Therefore,
here, we consider any system of PDEs that can be written in conservation form:

$\partial_{t}\textbf{Q}+\partial_{i}\textbf{F}^{i}(\textbf{Q})=\textbf{S},$
(1)

where the summation is performed over the spatial dimensions, Q is the vector
of the conserved variables, $\textbf{F}^{i}$ the vector of the physical
fluxes, and S is the vector of the sources. From the plethora of system that
can be written in the form (1), in the following, we are going to study the
scalar transport and Burgers equations and the equations of special and
general relativistic hydrodynamics [15].

The basic idea of the EFL method consists of expressing the numerical fluxes
resulting from the spatial discretisation of (1) as a superposition of an
unstable high-order and a stable low- or high-order flux, where the weight
dictating the transition between the two fluxes is computed from the local
entropy production of the system under investigation. Thus, the entropy of the
system is used as a “shock detector” indicating when to switch from the
unstable high-order scheme to the stable one.

For systems with thermodynamic entropy the implementation of the EFL method
follows the exposition in [9]. In the following, we will show how to extend
the EFL machinery of [9] to systems without a physically motivated definition
of entropy, like the ones of Sec. IV.

In accordance with [9], we start by approximating the spatial derivative of
the $x$ component, $\textbf{F}^{x}$, of the physical flux appearing in (1)
with the conservative finite-difference formula111For clarity and without loss
of generality, from now on the presentation is restricted to one dimension,
say $x$. A multidimensional scheme is obtained by considering fluxes in each
direction separately and adding them to the r.h.s.

$\partial_{x}F^{x}_{i}=\frac{\hat{f}_{i+1/2}-\hat{f}_{i-1/2}}{h},$ (2)

where $F^{x}$ is any one of the components of $\textbf{F}^{x}$ with
$F^{x}_{i}=F^{x}(x_{i})$, $\hat{f}_{i\pm 1/2}$ are the numerical fluxes at the
cell interfaces and $h$ is the spatial grid spacing.

In order to construct a flux-limiter, the numerical fluxes in (2) have to be
split into two contributions in the usual way [3]:

$\hat{f}_{i\pm 1/2}=\theta_{i\pm 1/2}\hat{f}^{\,\mathrm{HO}}_{i\pm
1/2}+(1-\theta_{i\pm 1/2})\hat{f}^{\,\mathrm{LO}}_{i\pm 1/2}.$ (3)

The quantities entering the r.h.s of (3) read:

  * •

$f^{\,\mathrm{HO}}$ is a high-order (HO) unstable numerical flux used in
regions where the numerical solution is smooth.

  * •

$f^{\,\mathrm{LO}}$ is a high- or low-order (LO) stable numerical flux used in
regions where the numerical solution develops non-smooth features. Notice that
in order to be consistent with previous work [12, 9], we keep the notation LO
although in the following we use exclusively a stable high-order scheme to
approximate $f^{\,\mathrm{LO}}$.

  * •

$\theta\in[0,1]$ is a continuous weight function indicating how much from the
numerical fluxes $f^{\,\mathrm{HO}}$ and $f^{\,\mathrm{LO}}$ to use at every
instance.

In particular, the flux $\hat{f}^{\,\mathrm{HO}}$ is built using the Rusanov
Lax-Friedrichs flux-splitting technique and performing the reconstruction on
the characteristic fields [16, 7]. A fifth-order central unfiltered stencil
(CS5) is always used for reconstruction. The flux $\hat{f}^{\,\mathrm{LO}}$ is
approximated by the LLF central scheme with reconstruction performed on the
primitive variables [17]. Primitive reconstruction here is performed with the
fifth-order weighted-essentially-non-oscillatory finite difference schemes
WENO5 [18] and WENOZ [19].222To produce the results of FIG. 5 the third-order
convex-essentially-non-oscillatory (CENO3) algorithm [20, 21] is used to
perform primitive reconstruction. With this choice, originally introduced in
[9], we generalise the traditional notion of a flux-limiter where
$\hat{f}^{\,\mathrm{LO}}$ is always a LO monotone flux [3, 22].

The weight function $\theta$ is computed from the so-called entropy production
function $\nu$, a quantity that depends on the local entropy production.
Explicitly, the relation between $\theta$ and $\nu$ is

$\theta_{i\pm 1/2}=1-\frac{1}{2}(\nu_{i}+\nu_{i\pm 1}).$ (4)

Below, we summarise how to compute $\nu$ for systems with and without
thermodynamically defined entropy.

Systems with thermodynamic entropy: In order to quantify the relation between
$\nu$ and the entropy produced by the systems of Secs. V, VI, and VII, we
define the specific entropy (entropy per unit mass) of any piecewise
polytropic EoS333For a more general EoS the specific entropy $s$ can be taken
from the EoS. as

$s=\ln\left(\frac{p}{\rho^{\Gamma}}\right),$ (5)

where the pressure is computed in accordance with the EoS in use.

Following [12], we employ the second law of thermodynamics to define the
entropy residual:

$\mathcal{R}=\nabla_{\mu}(s\,\rho\,u^{\mu})\geq 0,$ (6)

which provides a quantitative estimation of the rate of the entropy produced
by the system under study. Using the continuity equation and writing the
4-velocity $u^{\mu}$ in terms of the fluid 3-velocity $\upsilon^{i}$, the
above expression can be written [12] in terms of the time and spatial
derivatives of the specific entropy as

$\mathcal{R}=\frac{\rho
W}{\alpha}\left(\partial_{t}s+(\alpha\,\upsilon^{i}-\beta^{i})\partial_{i}s\right).$
(7)

Following the argumentation in [9], we suppress the multiplication factor
$\frac{\rho W}{\alpha}$ and replace $\mathcal{R}$ by

$R=\partial_{t}s+(\alpha\,\upsilon^{i}-\beta^{i})\partial_{i}s,$ (8)

which amounts to a rescaling of $\mathcal{R}$ so that the coefficient of
$\partial_{t}s$ is equal to one.

Systems without thermodynamic entropy: The systems of Sec. IV do not have an
EoS which can be used to define the entropy. For such systems (1) takes the
scalar form

$\partial_{t}Q+\partial_{i}F^{i}(Q)=S.$ (9)

Following [10, 11], we define an entropy for (9) by introducing the so-called
entropy pair: $E(Q)$ and $\mathcal{F}^{i}(Q)$. The convex function $E$ is
called entropy functional and $\mathcal{F}^{i}(Q)$ is the associated entropy
flux given by the expression

$\mathcal{F}^{i}(Q)=\int\partial_{Q}E(Q)\,\partial_{Q}F^{i}(Q)\,dQ.$ (10)

The entropy pair can be used to define an entropy residual:

$R=\partial_{t}E(Q)+\partial_{i}\mathcal{F}^{i}(Q).$ (11)

For convex physical fluxes $F^{i}$, like the ones of Sec. IV, the entropy pair
generated by the choice $E(Q)=\frac{1}{2}Q^{2}$ has been proven [10, 11]
enough for our purposes here. With the above choice of the entropy functional,
the entropy residual (11) can be straightforwardly computed for the systems of
Sec. IV.

Finally, we define the so-called entropy production function, for systems with
and without thermodynamic entropy, in terms of the rescaled entropy residual
$R$ defined by (8) or (11),

$\nu_{E}=c_{E}|R|,$ (12)

where $c_{E}$ is a tunable constant used to scale the absolute value of $R$.
In all our three-dimensional general relativistic simulations of Secs. VI and
VII we do not have to tune $c_{E}$, its value is set always to unity, i.e.
$c_{E}=1$. For the one-dimensional simulations of Secs. IV and V the tunable
constant $c_{E}$ is set either to 10 or 50 depending on the case studied.

Keeping in mind that the parameter $\theta$ cannot exceed unity, we have to
impose a maximum value of $\nu_{\rm max}=1$ for the entropy production
function in order to ensure that the rhs of (4) does not exceed the range
$[0,1]$. Accordingly, the entropy production function entering (4) is given by

$\nu=\min\left[\nu_{E},1\right].$ (13)

##  III Numerical implementation

The finite differencing code BAM [23, 17, 24, 7] is used to solve numerically
the system of equations of the form (1) discussed in the present work. The
aforepresented EFL method is implemented into BAM and is part of its
infrastructure. BAM uses the method-of-lines with Runge-Kutta (RK) time
integration and finite differences for the approximation of spatial
derivatives. The value of the Courant-Friedrich-Lewy (CFL) condition is set to
0.25 for all runs.

For the numerical implementation of the EFL method the BAM routines computing
the numerical fluxes had to be modified in order to accommodate the hybrid
flux (3). In order to compute the entropy production function we have to
approximate the time and spatial derivatives of the entropy quantities in (8)
and (11). We use finite differences to do so. Specifically, the spatial
derivatives are approximated with a standard centered finite-difference
stencil of order $p+1$ or higher, where $p$ is the order of the stencil used
to approximate the physical fluxes. (In the present work we use $p=5$.) With
this restriction it is ensured that the entropy production function $\nu$
converges to zero in the smooth regions faster than the overall convergence of
the scheme. The time derivative is also approximated with finite differences.
We employ a third-order one-sided stencil by using, at every Runge-Kutta
iteration, the current value of the specific entropy and the values at the
three previous timesteps. The fact that we manage to achieve higher than
third-order convergence in the majority of our simulations can be possibly
attributed to the dominance of the spatial error over the time discretization
error.

The numerical domain contains a mesh made of a hierarchy of cell-centered
nested Cartesian boxes and consists of $L$ refinement levels $l=0,\dots,L-1$
ordered with increasing resolution. Each refinement level is made out of one
or more equally spaced Cartesian grids with grid spacing $h_{l}$. There are
$n$ points per direction on each grid plus a certain number of buffer points
on each side. (For simplicity, we always quote grid sizes without buffer
points.) The resolution between two consecutive levels is doubled such that
the grid spacing at level $l$ is $h_{l}=h_{0}/2^{l}$, where $h_{0}$ is the
grid spacing of the coarsest level. The inner levels move in accordance with
the moving boxes technique, while the outer levels remain fixed. The number of
points in one direction of a moving level can be set to a different value than
the number of points of a fixed level. The coordinate extent of a grid at
level $l$ entirely contains grids at any level greater than $l$. The moving
refinement levels always stay within the coarsest level. For the time
evolution of the grid the Berger-Oliger algorithm is employed enforcing mass
conservation across refinement boundaries [24, 25]. Restriction and
prolongation is performed for the matter fields with a fourth-order WENO
scheme and for the metric fields with a sixth-order Lagrangian scheme.
Interpolation in Berger-Oliger time stepping is performed at second-order.

The derivatives of the metric components are approximated by fourth-order
accurate finite-differencing stencils. In addition, sixth-order artificial
dissipation operators are employed to stabilize noise from mesh refinement
boundaries. The general relativistic hydrodynamic equations (1) are solved by
means of a high-resolution-shock-capturing method [17] based on primitive
reconstruction and the aforedescribed high-order entropy limited scheme for
the numerical fluxes. In the present work spacetime is dynamically evolved
using either the BSSNOK [26, 27, 28] or the Z4c [29, 30] evolution scheme.

Vacuum regions are simulated with the introduction of a static, low-density,
cold atmosphere in the vacuum region surrounding the star [17]. The atmosphere
density is defined as

$\rho_{\rm atm}=f_{\rm atm}\,\mathrm{max}\,\rho(t=0).$ (14)

All grid points with rest-mass density below a threshold value $\rho_{\rm
thr}=f_{\rm thr}\rho_{\rm atm}$ are set automatically to $\rho_{\rm atm}$.
Transition to low-density regions is one of the main sources of error in NS
simulations. This is a common feature in all current numerical relativity
implementations of NS dynamics. To deal with this challenging feature they
also make use of similar assumptions and algorithms at low densities as those
employed here.

The equations of state (EoS) used in the present work are the following. For
the special relativistic tests of Sec. V we use a $\Gamma$-law EoS,

$p=(\Gamma-1)\rho\,\epsilon,$ (15)

with $\Gamma$ the adiabatic index. The neutron star matter of the single
neutron star evolutions of Sec. VI is also modelled by a $\Gamma$-law EoS of
the form (15). Finally, the matter of the neutron stars comprising the binary
of Sec. VII is described by a SLy EoS [31]. The latter is implemented by a
piecewise polytrope fit [32], and thermal effects are modeled by an additive
pressure contribution given by the $\Gamma$-law EoS with $\Gamma=1.75$ [33,
34, 17].

##  IV Scalar 1D tests

###  IV.1 Transport equation with non-smooth initial data

We start from the simplest one-dimensional scalar equation that describes the
evolution of shocks:

$\partial_{t}u(t,x)+\partial_{x}u(t,x)=0\quad\mathrm{with}\quad
u(0,x)=u_{0}(x).$ (16)

In order to test the convergence properties of the EFL method we start with
smooth initial data of the form $u_{0}(x)={\rm e}^{-40x^{2}}$. The exact
solution is infinitely smooth and is deduced from the initial data:
$u(t,x)=u_{0}(x-t)$. Convergence studies are performed on an equidistant grid
with range $[-1,1]$ that is progressively refined with $n=2^{k}\times 200$
grid points, where $k=0,\ldots,5$. Different high-order RKp (with $3\leq p\leq
5$) schemes are used for the temporal evolution. Our convergence analysis
indicates that optimal convergences $p$ is achieved in all cases independently
of the scheme used for the LO flux.

Figure 1: Transport equation with non-smooth square wave initial data. The
numerical solution (top panel) and the corresponding entropy production
profile (bottom panel) are depicted at different stages of the evolution. The
solution is computed with the WENOZ reconstruction scheme on a grid of 800
points with resolution $\Delta x=2.5\times 10^{-3}$. Dotted black lines depict
the exact solution.

Next, we consider non-smooth square wave initial data:

$u_{0}(x)=\begin{cases}1,&\text{if\,\, $-0.5\leq x\leq 0.5$},\\\
0,&\text{otherwise}.\end{cases}$ (17)

The grid set-up is the same as above. The initial data (17) is evolved with a
RK4 routine and the CFL condition is set to $0.125$. The WENOZ scheme is used
for the LO flux and the value of the tunable constant is set to $c_{E}=10$.
FIG. 1 depicts the numerical solution $u$ of (16) and the corresponding
entropy production profile $\nu$ at different stages of the evolution on a
grid composed by 800 grid-points ($\Delta x=2.5\times 10^{-3}$). It is
apparent that the EFL method locates successfully the position of the shocks
and tracks them efficiently during the evolution. The fact that around the
shocks of the square wave $\nu=1$ guarantees the use of the non-oscilatory LO
flux there. Thus, the EFL method is able to suppress the emergence of
oscillations around the vicinity of the shocks while achieving HO converenge
on the rest of the computational domain where the solution is smooth. In
addition, the behaviour of the entropy production profiles presented in FIG. 2
entails that with increasing resolution the entropy production gets better
localised around the shocks and thus reduces the impact of the LO flux only to
areas very close to the shocks.

Figure 2: Entropy production profiles for the transport equation with
increasing resolution. The profile of $\nu$ at $t=0.2$ for a low $n=100$
(solid blue line), an intermediate $n=400$ (dashed green line) and a high
$n=3200$ (dotted red line) resolution. The solution is computed with the WENOZ
reconstruction scheme.

Finally, TABLE 1 contains the results of the convergence analysis for the
numerical solution of FIG. 1 at $t=0.2$. The observed convergence rates are in
accordance with their theoretically expected values [3] of $1$ and $0.5$ for
the norms $L_{1}$ and $L_{2}$, respectively. During the evolution the
convergence rates are kept constant. Using lower that fifth order
reconstruction schemes for the LO flux results in slightly smaller convergence
rates, something expected as the lower accuracy of the LO flux drives down the
overall accuracy of the entropy flux limiter. The use of other RK time
integration routines, like RK3 and RK5, did not significantly affect the
results of TABLE 1.

Table 1: Convergence analysis of the numerical solution of FIG. 1 at $t=0.2$. $L_{1}$ and $L_{2}$ are normalised norms and the convergence rate is calculated as the $\log_{2}$ of the ratio of two successive normalised norms. n | $L_{1}$ | Conv. | $L_{2}$ | Conv.
---|---|---|---|---
200 |  $8.8\times 10^{-3}$ | – |  $4.8\times 10^{-2}$ | –
400 |  $5.1\times 10^{-3}$ | 0.802 |  $3.7\times 10^{-2}$ | 0.398
800 |  $2.9\times 10^{-3}$ | 0.802 |  $2.8\times 10^{-2}$ | 0.396
1600 |  $1.7\times 10^{-3}$ | 0.821 |  $2.1\times 10^{-2}$ | 0.397
3200 |  $9.3\times 10^{-4}$ | 0.834 |  $1.6\times 10^{-2}$ | 0.401
6400 |  $5.2\times 10^{-4}$ | 0.836 |  $1.2\times 10^{-2}$ | 0.412

###  IV.2 Burgers equation

Next, we consider the inviscid Burgers equation in one-dimension:

$\partial_{t}u(t,x)+u(t,x)\partial_{x}u(t,x)=0\quad\mathrm{with}\quad
u(0,x)=u_{0}(x).$

This is the simplest non-linear scalar equation that develops shocks from
originally smooth initial data:

$u_{0}(x)={\rm e}^{-40x^{2}}.$ (18)

During the evolution the smooth initial profile becomes steeper and steeper
and at around $t\simeq 0.1843$ a shock is formed.

Figure 3: Burgers equation. The numerical solution (top panel) and the
corresponding entropy production profile (bottom panel) are depicted at
different stages of the evolution. The solution is computed with the WENOZ
reconstruction scheme on a grid of 1600 points with resolution $\Delta
x=1.25\times 10^{-3}$. Dotted black lines depict the exact solution.

The exact solution is derived from the initial data (18) in the implicit form
$u(t,x)=u_{0}(x-ut)$. The grid set-up is the same as in the case of the
transport equation, see subsection IV.1. An RK4 routine is used for the
temporal evolution and the CFL condition is set to 0.125. The WENOZ scheme is
used for the LO flux and the value of the tunable constant is set to
$c_{E}=10$. FIG. 3 depicts the numerical solution of Burgers equation (top
panel) for the initial data (18) and the corresponding entropy production
profile (bottom panel) at different stages of the evolution on a grid of 1600
grid-points ($\Delta x=1.25\times 10^{-3}$). It is evident that the EFL method
recognizes when the solution is smooth–the entropy production is zero–and, in
addition, locates successfully the position of the emerging shock and tracks
it efficiently during the evolution. It can be also confirmed by inspection
that the EFL solutions are oscillation-free. The localisation of the entropy
production profiles with increasing resolution first observed in FIG. 2 is
also observed here.

Figure 4: Convergence to zero with increasing resolution of the error between
the numerical and the exact solution of Burgers equation at $t=1$. The
vertical dashed black line depicts the position at which the shock is first
formed. Dotted lines show results scaled to fifth-order.

FIG. 4 depicts the numerical error with increasing resolution over the whole
computational domain at $t=1$. Optimal fifth-order convergence is observed at
the parts of the domain where the solution is smooth—dotted lines show results
scaled to fifth-order. As expected close to the location of the shock the
convergence drops to first-order. The position of the right-moving shock can
be easily identified by the lower order convergence spike at $x\approx 0.6$.
It is clear that the movement of the shock through the computational grid does
not contaminate the rest of the numerical domain with its low order
convergence. It is evident that the points of the computational domain
affected by the shock recover the optimal fifth-order convergence almost
instantly. Therefore, lower converegence is a local characteristic of our
solution confined to the neighborhood of the shock.

Figure 5: Convergence rates with respect to the $L_{1}$-norm for the numerical
solution of Burgers equation. The behaviour of the convergence rates with time
for two different LO schemes is shown: EFL-WENOZ (solid blue line) and EFL-
CENO3 (dotted green line). The vertical dashed black line depicts the instant
at which the shock is first formed.

In FIG. 5 the convergence rates with respect to the $L_{1}$-norm for the
numerical solution of Burgers equation are computed at different stages of the
evolution for two different choices of the LO flux. At early stages of the
evolution the convergence rates maintain their optimal high-order values. When
the shock starts forming—see vertical dashed black line—the rates decrease
rapidly and, after the shock has been fully formed, reach a plateau at around
1. Notice that in the smooth region before the formation of the shock, the use
of a lower than fifth-order reconstruction scheme, i.e. CENO3, leads to
slightly smaller convergence rates. The use of other RK routines does not
alter the convergence rates of FIG. 5 significantly.

##  V Special relativistic 1D tests

###  V.1 Simple wave

Figure 6: Simple wave. The numerical solution of the velocity (top panel) and
the corresponding entropy production profile (bottom panel) are depicted at
different stages of the evolution. The solution is computed with the WENOZ
reconstruction scheme on a grid of 1600 points with resolution $\Delta
x=1.875\times 10^{-3}$. Dotted black lines depict the exact solution.
t=0.6t=0.84
---
t=0.96t=1.2
Figure 7: Behaviour of the simple wave’s numerical error with increasing
resolution at different instances of the evolution. The vertical dashed black
line depicts the position at which the shock is first formed. Dotted lines
show results scaled to fifth-order.

Relativistic simple waves are exact analytical solutions of special-
relativistic hydrodynamics (SRHD) that demonstrate some of the main features
of non-linear wave propagation, such as steepening, breaking, and shock
formation. Simple wave solutions can be written analytically in an implicit
form using the method of characteristics. The study of relativistic simple
waves sheds light on the problem of relativistic shock formation due to non-
linear steepening. Like Burgers equation the relativistic simple wave starts
off from smooth initial data and at some point of its evolution develops a
shock. These analytic solutions have been discussed in [35, 36]. Here, we use
the simple wave described in [7], therein the initial velocity profile is of
the form

$\upsilon=a\,\Theta(|x|-X)\sin^{6}\left(\frac{\pi}{2}\left(\frac{x}{X}-1\right)\right),$
(19)

where $\Theta(x)$ is the Heaviside function, $a=0.5$ and $X=0.3$. For the
construction of the initial data a polytropic EoS, $p=K\rho^{\gamma}$, with
$K=100$ and $\gamma=5/3$ is employed; this data is then evolved using a
$\Gamma$-law EoS (15) with $\Gamma=5/3$. During the evolution the smooth
initial profiles of all primitive variables become steeper and steeper and at
around $t_{sh}\simeq 0.63$ they form a shock [35]. Numerical solutions are
computed on the one-dimensional domain $x\in[-1.5,1.5]$ with the RK4 time-
integrator and a CFL factor of 0.125. The tunable constant is set to
$c_{E}=10$ and the LO flux is evaluated with the WENOZ reconstruction scheme.

The top panel of FIG. 6 illustrates the evolution of the simple wave from its
initial smooth configuration (19) up to the instance $t=1.2$ when the shock
has been fully formed. The solution is computed on a grid composed of 1600
points ($\Delta x=1.875\times 10^{-3}$). By inspection, the numerical
solutions reproduce the correct physics even after the formation of the shock
$t\geq t_{sh}$. No spurious oscillations are visible on the solutions, a sign
that shows the effectiveness of the EFL method in suppressing oscillatory
effects. The bottom panel of FIG. 6 depicts the corresponding profiles of the
entropy production function $\nu$. Clearly, the EFL method successfully
locates the initial formation of the shock at $t=t_{sh}$ and effectively
tracks its location at later times $t\geq t_{sh}$. During the “smooth period”
$t<t_{sh}$, the entropy production is, as expected, trivial. Notice that
shortly after the formation of the shock, i.e. at $t=0.75$, the solution is
oscillation-free although the contribution of the stable LO flux in the hybrid
numerical flux (3) is only $\sim 33\%$. This fact points out that the EFL
method can guarantee the smoothness of the solutions for some time after the
shock formation without switching entirely to the LO flux. This property of
the EFL method has a direct effect on the behaviour of the convergence rates
of FIG. 8. As expected, optimal fifth-order convergence is maintained up to
the formation of the shock—denoted here by a vertical dashed line. With the
shock formed the convergence of our numerical solution decreases slowly up to
roughly $t\sim 0.85$—notice that the convergence rates are still higher than
fourth-order. Beyond this point the convergence decreases rapidly and reaches
a plateau at around 1. Clearly, there is a direct connection between the
increasing use of the LO flux and the decreasing order of convergence.

Figure 8: Convergence rates of the simple wave with respect to the
$L_{1}$-norm. The behaviour of the convergence rates with time for the LO
scheme EFL-WENOZ is shown. The vertical dashed black line depicts the instant
at which the shock is first formed.

FIG. 7 depicts the absolute error between the numerical and the exact solution
with increasing resolution at different stages of the evolution. For times
$t<t_{sh}$ before the formation of the shock the numerical error decreases to
zero with the expected fifth-order convergence, see top-left panel of FIG. 7.
With the formation of the shock at $t=t_{sh}$ the expected lower-order
convergence spike appears. Shortly after, the numerical error convergences to
zero with fifth-order over the whole computational domain except at the
location of the shock where the convergence drops to first-order, see top-
right panel of FIG. 7. Notice that, at this stage, the convergence of the
parts of the computational domain impacted by the shock is quickly restored to
fifth-order. At late times, see bottom panels of FIG. 7, this is not anymore
possible and the low-order convergence spreads through the computational grid
at the opposite direction of the movement of the shock. This is in stark
contrast with the behaviour of the numerical error of Burgers equation, see
FIG. 4. The coupled non-linear nature of the SRHD equations is responsible for
this difference as the EFL method operates in both cases in a similar fashion:
the LO flux dominates the hybrid numerical flux (3) and its use has been
restricted around the location of the shock.

###  V.2 Sod shock-tube

Figure 9: Sod shock-tube. Profiles of the rest-mass density (top panel) and of
the corresponding entropy production (bottom panel) at different stages of the
evolution. The solution is computed with the WENO5 reconstruction scheme on a
grid of 1600 points with resolution $\Delta x=6.25\times 10^{-4}$. Dotted
black lines depict the exact solution.

We move on now to a standard benchmark Riemann problem used in SRHD: the
relativistic version of Sod’s shock-tube [37]. For a simple ideal fluid EoS of
the form (15) with adiabatic index $\Gamma=1.4$, the discontinuous initial
data for the pressure $p$, the rest-mass density $\rho$, the velocity
$\upsilon$, and the specific energy $\epsilon$ are the following

$\displaystyle(p_{L},\rho_{L},\upsilon_{L},\epsilon_{L})$
$\displaystyle=(1,1,0,2.5),$ (20)
$\displaystyle(p_{R},\rho_{R},\upsilon_{R},\epsilon_{R})$
$\displaystyle=(0.1,0.125,0,2).$

The initial discontinuity at $x=0$ splits during the evolution into a shock
wave followed by a contact discontinuity, both travelling to the right, and a
rarefaction wave travelling to the left. The numerical solutions are computed
on the one-dimensional grid $x\in[-0.5,0.5]$ with the RK4 scheme and a CFL
factor of 0.25. The tunable constant is set to $c_{E}=50$ and the WENO5
reconstruction scheme is used to evaluate the LO flux.

The top panel of FIG. 9 shows how the initial discontinuity at $x=0$ evolves
with time on a grid of 1600 points ($\Delta x=6.25\times 10^{-4}$). Even at
early stages of the evolution the three characteristic features of Sod’s
shock-tube are clearly visible: the left-moving rarefaction wave and the
right-moving shock and contact discontinuities. It is evident that all the
features of the Sod shock-tube are reproduced quite accurately and oscillatory
effects are absent from the solutions. At the bottom panel of FIG. 9 the
corresponding entropy production profiles are shown. Again the EFL method
locates and tracks the right-moving discontinuities extremely accurately.
Notice that the entropy production peaks are very well localised around the
position of the discontinuities, restricting in this way the use of the stable
LO flux only to these problematic regions and enabling the use of the HO flux
in the rest of the computational domain where the solution is smooth. At early
stages of the evolutions a third entropy production peak related to the
rarefaction wave is clearly visible. The initial very steep slope of the
rarefaction triggers the production of this peak. With time, as the
rarefaction develops, the amplitude of the entropy production peak reduces and
tends to zero when the rarefaction wave has been fully developed. This
comprises another manifestation of the sensitivity and accuracy of the EFL
method.

Figure 10: Percentage of the stable LO flux used at $t=0.5$ over the whole
computational domain.

A more quantitative picture of the use of the LO flux at each point of the
computational domain is given in FIG. 10. As mentioned in Sec. II the entropy
production takes values in the interval $\nu\in[0,1]$. In the extreme cases of
$\nu=0$ or $\nu=1$ the HO or LO flux is exclusively used by the EFL method.
For any other value of $\nu$ the resulting numerical flux (3) is a mix of
these two fluxes. FIG. 10 shows exactly that: the percentage of each flux used
at each point of the computational domain at $t=0.5$. As expected the stable
LO flux is exclusively used at the locations of the contact and shock
discontinuities, in between them the use of the LO flux drops drastically to
$\sim 0.5\%$ and at any other point of the computational domain does not
exceed $0.01\%$. Notice that the positions of the head and tail of the
rarefaction wave are clearly indicated by the two $0.01\%$ entropy production
peaks—the EFL method is so sensitive that recognises the slight non-smoothness
of these two features compared to their neighbouring points.

Similar results we get also for the other two standard SRHD shock-tube tests:
the relativistic blast wave 1 and 2 described in [38].

##  VI Single star evolutions

We proceed further by applying the EFL method to three-dimensional general
relativistic single NS spacetimes. In addition to the static TOV and rotating
NS tests presented in [9], we study here a couple of more involved single NS
spacetimes that unveil the ability of the EFL method to cope extremely well
with non static/stationary spacetimes. The treatment of the rapidly declining
gradient of the hydrodynamical variables at the surface of the NS remains also
here the most challenging feature of our simulations. In addition, the non
static/stationary nature of these spacetimes increases further the level of
complexity we encounter. Our results are compared with those of [39, 40, 41,
17]. The EFL method is also compared to the best performing schemes currently
implemented in BAM: i) a second-order scheme (LLF-WENOZ) that uses the LLF
scheme for the fluxes and WENOZ for primitive reconstruction [17] and ii) a
“hybrid” algorithm (HO-LLF-WENOZ) that employs the high-order HO-WENOZ scheme
above a certain density threshold $\rho_{\mathrm{hyb}}$ and switches to the
standard second-order LLF-WENOZ method below $\rho_{\mathrm{hyb}}$ [7].

In the following, we evolve a stable TOV star that moves along the
$x$-direction [39, 17] and an unstable migrating NS [40, 41, 17], both in a
dynamically evolved spacetime. The NS matter is here described by a
$\Gamma$-law EoS with $\Gamma=2$. The grid is composed of five fixed
refinement levels. Simulations are performed at resolutions $n=(64,96,128)$
points leading to a grid spacing $h$ that depends on the specific setting of
the NS under investigation. For each NS configuration the resolution is
explicitly given in TABLE 2. It is ensured that the NS is entirely covered by
the finest box at any given resolution. Radiative (absorbing) boundary
conditions are used for all single star simulations.

Table 2: Grid configurations of single star simulations. Columns (left to right): name of simulation, $L$: number of fixed refinement levels, $n$: number of points per direction, $h_{L-1}$: resolution per direction in the finest level $l=L-1$, $h_{0}$: resolution per direction in the coarsest level $l=0$. Name | $L$ | $n$ | $h_{L-1}$ | $h_{0}$
---|---|---|---|---
TOVboost | 5 | 64 | 0.500 | 8.000
5 | 96 | 0.333 | 5.333
5 | 128 | 0.250 | 4.000
TOVmig | 5 | 64 | 0.125 | 2.000
5 | 96 | 0.083 | 1.333
5 | 128 | 0.063 | 1.000

###  VI.1 Boosted TOV star

The Tolmann-Oppenheimer-Volkoff (TOV) initial data used here is constructed
using a $\Gamma=2$ polytrope model with $K=100$, gravitational mass $M=1.4$,
baryonic mass $M_{b}=1.506$ and central rest-mass density $\rho_{c}=1.28\times
10^{-3}$. The resulting TOV star is boosted in the $x$-direction at a speed of
$v=0.5$ corresponding to a Lorentz factor of $W=1.15$. The spacetime is
dynamically evolved and the BSSNOK scheme is used for the evolution of the
metric.

Figure 11: Boosted TOV star. Profiles of the rest-mass density (top panel) and
of the corresponding entropy production (bottom panel) at different instances
of the movement of the TOV star along the $x$-direction. The solution is
computed with the WENOZ reconstruction scheme and $n=128$. Dotted black lines
depict the exact solution.

The one-dimensional profiles of the rest-mass density (top panel) and entropy
production function $\nu$ (bottom panel) are depicted in FIG. 11. The top
panel shows the evolution of the TOV star along the $x$-direction from its
initial position (dashed black lines) up to $t=3000M_{\odot}$. The numerical
solution reproduces quite accurately the exact one (dotted black lines) even
at very late times. At the bottom panel the corresponding entropy production
profiles are depicted. As expected, a local peak of the entropy production
function $\nu$ is observed around the surface of the TOV star. There, the
gradient of the hydrodynamical variables experiences a violent decline which
leads to the production of large values of $\nu$. In the interior of the NS
the entropy production function $\nu$ is vanishing. It is evident from FIG. 11
that the EFL scheme locates and tracks extremely accurately the surface of the
star during its movement along the $x$-axis. This, in turn, triggers the use
of the stable numerical flux in the neighbourhood of the surface. The use of
the stable scheme around the surface of the NS guarantees its long term
stability.

Figure 12: Central rest-mass density evolution of a boosted TOV star in a
dynamical spacetime with $n=128$. Beside the EFL-WENOZ scheme, the LLF-WENOZ
and HO-LLF-WENOZ schemes are also shown.

After securing the proper flagging of the problematic regions and the correct
implementation of the EFL method, we proceed further by checking its
performance during the evolution of the TOV star by monitoring the dynamical
behaviour of the central rest-mass density. The oscillation of the central
rest-mass density $\rho_{\rm max}$ for the EFL method is presented, together
with the results for the LLF-WENOZ and HO-LLF-WENOZ methods [7], in FIG. 12.
Clearly, the EFL method is performing better than the LLF-WENOZ and HO-LLF-
WENOZ schemes as the profile of the central rest-mass density is more stable
over time and the amplitude of the oscillations is smaller.

Figure 13: Evolution of the $L_{1}$ distance $||\rho(t)-\rho(0)||_{1}$ of a
boosted TOV star in a dynamically evolved spacetime. The EFL-WENOZ scheme is
compared to the LLF-WENOZ and HO-LLF-WENOZ schemes. Dashed lines show results
scaled to second order.

The non-dispersive character of the boosted TOV star enables us to check the
convergence properties of the EFL method as the exact solution can be read off
from the initial data. Here we consider the $L_{1}$-norm of the difference
between the three-dimensional evolution profile of the rest-mass density and
the corresponding exact solution (initial data) and study its behaviour with
time. The $L_{1}$-distance from the exact solution for the three schemes used
here are depicted in FIG. 13. The convergence rate of all the schemes
considered is approximately second-order in agreement with the result of [17]
and the fact that the error at the stellar surface dominates the evolution.

###  VI.2 Migrating neutron star

We turn now to the study of migrating neutron stars and focus on the migration
of unstable neutron stars to stable ones. This is a quite challenging test as
during the process of migration oscillations of large amplitude are produced
by the pulsating NS. The rapid periodic expansion and contraction of the
surface of the NS is the ideal arena to test the EFL method’s performance in
an extreme dynamical setting where the location of the surface is constantly
changing in a non-trivial way.

The unstable initial configuration we consider here is constructed using a
$\Gamma=2$ polytrope model with $K=100$, gravitational mass $M=1.447$,
baryonic mass $M_{b}=1.535$ and central rest-mass density
$\rho_{c}=7.993\times 10^{-3}$ that is larger than the central rest-mass
density of the stable model with the maximum possible mass. Therefore the NS
is initially in an unstable equilibrium. Numerical truncation errors cause the
NS to move away from this initial equilibrium and to expand, evolving
gradually to smaller central rest-mass densities, until it reaches a new
stable NS configuration.

The star is evolved with the $\Gamma$-law EoS and the metric components with
the BSSNOK scheme. The spacetime is dynamically evolved.

Figure 14: Migration test. Profiles of the rest-mass density (top panel) and
of the corresponding entropy production (bottom panel) at different stages of
the initial rapid expansion of the migrating TOV star, see FIG. 15. The
solution is computed with the WENOZ reconstruction scheme and $n=128$.

We start by inspecting the one-dimensional profiles of FIG. 14 which depict
the rest-mass density (top panel) and the entropy production function (bottom
panel). The top panel shows the behaviour of the rest-mass density during the
initial rapid expansion of the NS. The NS evolves from its initial unstable
configuration at $t=0$ to $t=175M_{\odot}$ where the central rest-mass density
acquires its minimum value $\rho_{c}=3.139\times 10^{-4}$, see FIG. 15. In
just 0.8 ms the NS has expanded in such a degree that its central rest-mass
density is about 4% of its initial value. As can be seen on the bottom panel,
the EFL method reacts extremely well to this rapid change of the location of
the NS surface. The peaks on the profiles of the entropy production $\nu$ are
localised correctly around the moving surface of the NS ensuring that the
stable numerical flux is used there. This, in turn, guarantees that the
numerical solution is free of spurious oscillations and long term stable.

Figure 15: Evolution of the central rest-mass density for the migration test
with $n=128$. The results of the EFL method are compared to the second-order
scheme LLF-WENOZ used in [17, 24].

We continue by checking the behaviour of the central rest-mass density with
time. FIG. 15 presents the evolution of the central rest-mass density for both
the EFL and LLF-WENOZ method. The resulting oscillating behaviour is due to
the strong non-linear pulsation of the originally unstable NS around a stable
configuration. At every pulsation shock waves are formed that dissipate
kinetic into thermal energy and gradually damp the oscillations. The results
of the EFL method are comparable to the corresponding ones of the LLF-WENOZ
scheme and the ones in [40, 41, 17].

---
Figure 16: Two-dimensional hybrid plots depicting the entropy production and
rest-mass density profiles across different stages of the pure EFL BAM:95
simulation. Left: inspiral. Middle: merger. Right: post-merger.

##  VII Binary neutron star evolutions

###  VII.1 Initial data and numerical setup

We turn now to general relativistic simulations of neutron star binaries with
the EFL method. In the following, we study the dynamics of a specific BNS
configuration: BAM:95 [42]. This configuration is similar to BAM:97 [43],
which was studied in the EFL framework in [9]. The main difference between
these two BNS configurations lies in the way their initial data have been
constructed: BAM:97 uses the Lorene library [13] while BAM:95 uses initial
data constructed with the SGRID code [14].

We chose this BNS simulation because it enables us to study the effect of
using the EFL method in the construction of the initial data. In all the BNS
simulations presented in [9] the EFL method was used only in the dynamical
evolution of Lorene initial data. BAM:95 gives us the opportunity to use the
EFL method in the construction of SGRID initial data through the eccentricity
reduction process.

The initial data for BAM:95 is constructed with the pseudospectral SGRID code
[44, 45, 46, 47], which makes use of the constant rotational velocity approach
[14, 48] to construct spinning BNS configurations in hydrodynamical
equilibrium. In accordance with [44, 49], eccentricity reduced initial data
are constructed through an iterative procedure of monitoring and varying the
binary’s initial radial velocity and eccentricity [44]. As an initial guess, a
quasi-equilibrium configuration in the usual quasi-circular orbit is employed
with residual eccentricity of the order of $e\sim 10^{-2}$. The iterative
procedure consist of the following steps:

  * i)

Evolve the data for 2-3 orbits.

  * ii)

Measure the eccentricity $e$ from the proper distance as described in [44].

  * iii)

Re-compute the initial data with adjusted radial velocity and eccentricity.

The actual initial data we evolve can be found in TABLE 3. As mentioned above,
it has been constructed with the SGRID code [14] and is characterised by the
Arnowitt-Deser-Misner (ADM) mass-energy $M_{\rm ADM}$, the angular momentum
$J_{0}$, the baryonic mass $M_{b}$ and the dimensionless GW circular frequency
$M\Omega_{0}$.

Table 3: BAM:95 initial data. Columns: name; EoS; number of orbits; binary mass; rest-mass; ADM mass; angular momentum; GW frequency.  Name | ID | EoS | orbits | $M$ | $M_{b}$ | $M_{\mathrm{ADM}}$ | $J_{0}$ | $2M\Omega_{0}$
---|---|---|---|---|---|---|---|---
BAM:95 | SGRID | SLy | 10 | 2.700 | 2.989 | 2.678 | 7.686 | 0.038

In the present work, we chose to experiment with the use of different
numerical schemes, for computing the numerical fluxes, both in the
construction of the SGRID initial data and in the evolution of these data with
BAM. The main reason for doing this is to study if the effect of using the EFL
method will lead to higher accuracy and better convergence properties of our
numerical solutions. In [42] the HO-LLF method was used both in the
construction of the initial data for BAM:95 and in its subsequent evolution
with BAM. Here, we use the EFL method to conduct two additional studies of the
BAM:95 configuration. In the first, the HO-LLF method is used for the initial
data construction and the EFL method for its BAM evolution. In the second,
both the construction of the initial data and its BAM evolution is conducted
with the EFL method. The aforementioned different combinations of numerical
schemes used to study BAM:95 in the present work are summarised in TABLE 4.

Table 4: BAM:95 case studies. Columns: name of BAM:95 case study; method used in the construction of the initial data; method used to evolve the initial data with BAM; reconstruction scheme used in the simulations; eccentricity. Name | SGRID | BAM | Reconstruction | $e[10^{-3}]$
---|---|---|---|---
pure HO-LLF | HO-LLF | HO-LLF | WENOZ | 0.4
hybrid EFL | HO-LLF | EFL | WENOZ | 0.4
pure EFL | EFL | EFL | WENOZ | 0.7

---

pure HO-LLFhybrid EFLpure EFL

Figure 17: BAM:95 waveforms. The amplitude (blue line), the real part (green
line) and the instantaneous frequency $M\omega_{22}$ (red line) of the GW
signal obtained from the ten-orbit BAM:95 simulation using WENOZ
reconstruction for the three case studies of TABLE 4: pure HO-LLF (top),
hybrid EFL (middle), and pure EFL (bottom). The time of merger is defined as
the first peak of the amplitude $A_{22}$ and is indicated by a black vertical
line.

The initial data for BAM:95 is evolved in the three different combinations of
TABLE 4 with the WENOZ reconstruction scheme. WENOZ reconstruction is used
because the results in [9] strongly indicate that WENOZ is the best performing
reconstruction scheme from those available. For each case study three
different grid resolutions are considered. The grid specifications for all
BAM:95 runs are the same and are reported in TABLE 5. The atmosphere setting
for all simulations are $f_{\rm atm}=10^{-11}$ and $f_{\rm thr}=10^{2}$. The
metric is evolved with the Z4c scheme. Standard radiative boundary conditions
are used for all BNS simulations. The behaviour and the convergence properties
of the conserved quantities of the present BAM:95 runs are similar to those of
BAM:97 presented in [9]. Therefore, in the following, for the sake of
presentational brevity and clarity, we solely focus on the qualitative
behaviour of the entropy production and the gravitational wave analysis of our
BAM:95 runs.

Table 5: Grid configurations for all the BAM:95 case studies. Columns (left to right): name of BNS simulation; refinement levels; minimum moving level index; number of points per direction in fixed levels; number of points per direction in moving levels; resolution per direction in the finest level $l=L-1$; resolution per direction in the coarsest level $l=0$. Name | $L$ | $l^{\mathrm{mv}}$ | $n^{\mathrm{fix}}$ | $n$ | $h_{L-1}$ | $h_{0}$
---|---|---|---|---|---|---
BAM:95 | 7 | 2 | 160 | 64 | 0.235 | 15.040
7 | 2 | 240 | 96 | 0.157 | 10.027
7 | 2 | 320 | 128 | 0.118 | 7.520

###  VII.2 Qualitative behaviour of the entropy production

---

pure HO-LLFhybrid EFLpure EFL

Figure 18: GW phase difference convergence rate study for the ten-orbit BAM:95
simulation. Left: pure HO-LLF case. Middle: hybrid EFL case. Right: pure EFL
case.

The entropy production function $\nu$ plays central role in our method. Hence,
it is of great interest to study its behaviour during the evolution of BNS
merger simulations. In the following, we discuss the two-dimensional entropy
production profiles of the ten-orbit simulation BAM:95.

In FIG. 16 we present two-dimensional hybrid plots depicting the entropy
production function $\nu$ and rest-mass density $\rho$ profiles at different
stages of the pure EFL simulation. (The hybrid EFL case shows similar
behaviour.) The different panels show (from left to right) selected snapshots
of the inspiral, merger and post-merger stages, respectively. Taking a closer
look at the $\nu$ profile, we observe that during the inspiral and merger the
EFL method is, as expected, only activated at the surface of each NS. The same
can be also observed in the post-merger regime, however, here also regions in
the exterior of the remnant NS are flagged by the EFL method. Judging from the
corresponding rest-mass density plots, the observed flagging in the exterior
is actually capturing the outward dynamics of the spiral density waves emitted
by the rotating remnant of the two inspiraling NSs. This fact shows the
effectiveness of the EFL method in locating any non-smooth features developing
in the numerical domain.

###  VII.3 Gravitational wave analysis

We study now the effect of using the EFL method on the gravitational waveforms
(GWs). Following [7, 9], we compute the GWs from the curvature scalar
$\Psi_{4}$ on coordinate spheres that are at distance $r$ from the origin of
the computational domain. First, spin weighted spherical harmonics are used to
expand $\Psi_{4}$ into its modes $\psi_{\ell m}$ and then the multipolar modes
$h_{\ell m}$ are reconstructed by solving $\ddot{h}_{\ell m}=\psi_{\ell m}$ in
the frequency domain [50]. In terms of the amplitude $A_{\ell m}$ and phase
$\phi_{\ell m}$ the different modes of $h_{\ell m}$ are given by the
expression

$\displaystyle R\,h_{\ell m}=A_{\ell m}e^{-i\phi_{\ell m}}.$ (21)

In the present section, we plot our results against the retarded time

$\displaystyle u=t-R_{\ast}=:t-R(r)-2M\log\left(\frac{R(r)}{2M}-1\right),$
(22)

where $M$ is the total gravitational mass of the BNS system.
$R(r)=r(1+M/2r)^{2}$ is the radius in Schwarzschild coordinates and $r$
corresponds to the radius in isotropic coordinates which is the extraction
radius in our simulations. The time of merger $u_{\rm mrg}$ is estimated from
the first peak of $A_{22}$ of the dominant $(\ell,m)=(2,2)$ mode

The waveforms obtained from the ten-orbit inspiral BAM:95 simulation for the
three case studies of TABLE 4 are presented in FIG. 17. The wave train of all
three waveforms after a first peak, indicating the time of merger, shows a
more complicated structure that involves multiple denser distributed peaks and
a slow amplitude decay. The instantaneous GW frequency
$\omega_{22}=-\mathcal{I}(\dot{h}_{22}/h_{22})$ is also plotted (red line),
which undergoes a sharp frequency increase near merger (a feature
characteristic of a chirp-like signal). Notice that from top to bottom the
time of merger is increasing, with the difference between pure HO-LLF and pure
EFL being around $u_{\rm mrg}\approx 150M$. Clearly, the use of the EFL method
elongates the inspiral and increases the merger time.

Next, self-convergence studies are performed based on simulations using
$(n_{i})=(64,96,128)\equiv(\textrm{LOW, MID, HIG})$ points per direction on
the highest AMR level. The phase difference between pairs of resolutions is
defines by the expression

$\displaystyle\Delta\phi^{(n_{i},n_{j})}_{\ell m}=\phi^{(n_{i})}_{\ell
m}-\phi^{(n_{j})}_{\ell m}\,.$ (23)

The experimental convergence rate $p$ is determined through the rescaling of
the resulting differences by a factor $s$ that expresses the theoretically
expected rate of decrease of these differences with increasing resolutions,
provided of course that our simulations are converging. The definition of the
factor $s$ follows [51]:

$\displaystyle s(p,n_{i},n_{j},n_{k})$
$\displaystyle=\frac{1-(n_{i}/n_{j})^{p}}{(n_{i}/n_{j})^{p}-(n_{i}/n_{k})^{p}},$
(24)

where $n_{i}<n_{j}<n_{k}$.

FIG. 18 shows a self-convergence study of the waveform phase differences for
the BAM:95 simulations. The three panels (from left to right) correspond to
results obtained with the pure HO-LLF, hybrid EFL, and pure EFL schemes,
respectively. The pure HO-LLF result [42] (left panel) serves as a reference
for comparison with the case studies involving the EFL method. Solid lines in
FIG. 18 represent phase differences between runs with consecutively increasing
resolutions; dashed lines correspond to rescaled differences of MID-HIG
differences, where the scaling factor is computed from (24) using an integer
convergence rate $p$; the gray shaded regions mark the differences in merger
times between runs with increasing resolution.

As expected the phase differences between runs with increasing resolutions
decrease for all simulations. A feature that indicates that all schemes are
capable of providing results that converge to some order. This observation is
further supported by the fact that the difference between merger times is
reducing, indicated by the narrowing of the gray shaded regions, with
increasing resolution. Notice that, in accordance with FIG. 17, the use of the
EFL method clearly increases the time of merger.

---

pure EFLpure EFL

Figure 19: GW phase difference convergence rate study for higher modes of the
pure EFL ten-orbit BAM:95 simulation with WENOZ reconstruction. Left:
$(3,2)$-mode. Right: $(4,4)$-mode.

All convergence series demonstrate a clear convergence behaviour: the pure HO-
LLF second-order, the hybrid EFL third-order, and the pure EFL fifth-order
convergence. Comparing the results obtained with the EFL method to the ones
obtained with the HO-LLF method, we find that the use of the EFL method leads
to higher convergence rates and to generically smaller absolute differences
MID-HIG at merger. The pure HO-LLF case (left panel) yields a clear second-
order convergence, consistent with the results in [42], and phase difference
MID-HIG at merger a factor $\sim 2$ larger than the pure EFL case (right
panel). We also note that in the early post-merger phase the convergence rate
for the EFL cases is higher (third-order) and the difference MID-HIG smaller.
An observation that unveils the robustness of the EFL method and its ability
to capture satisfactorily the violent dynamics of the remnant NS.

Given this clear higher-order convergence pattern for the dominant
$(2,2)$-mode of the pure EFL case, it is natural to also investigate the
convergence of other (subdominant) modes of $h_{\ell m}$ with $\ell>2$. FIG.
19 shows a convergence study of the $(\ell,m)=(3,2)$ and $(4,4)$ modes. Also
in this case, the phase differences show a consistent decrease with increasing
resolution and a clear fifth-order convergence behaviour. Furthermore, a
third-order convergence pattern clearly holds through merger and in the early
post-merger regime.

##  VIII Conclusions

The main purpose of the present work is to extend the use of the EFL method:
i) to scalar PDE systems and ii) to the construction of BNS initial data sets.
We address also the question whether the use of the EFL scheme in the
construction of BNS initial data has a positive impact on the convergence
properties of the produced gravitational waveforms. Based on our results, we
answer this in the affirmative. In addition, the EFL scheme is used to
successfully simulate a couple of very demanding single NS spacetimes: the
migrating and boosted TOV NS.

We extend the EFL method presented in [9] and show how to generalised it to
scalar PDE systems in Sec. II through the introduction of an auxiliary entropy
functional. The resulting EFL scheme successfully passes in Sec. IV a quite
demanding set of scalar benchmark problems, with results comparable to the
ones in the literature [3].

In Sec. V we revisit a couple of benchmark problems in special relativistic
hydrodynamics that were also studied in [9]. Here, we focus on the study of
some quantitative and qualitative properties of the entropy production
function that provide a better understanding of the workings of the EFL
method. Our results i) show the direct connection between the increasing use
of the LO flux and the decreasing order of convergence and ii) demonstrate the
extreme sensitive of the EFL method in recognising the slightest non-
smoothness features present in our solutions.

Next, the EFL method is tested in Sec. VI against three-dimensional general-
relativistic non static/stationary single NS configurations. As expected, the
EFL scheme accurately locates the surface of the NS, see FIG. 11 and FIG. 14,
and enables the use of the LO flux in this region while the interior remains
mainly resolved by the HO flux. EFL simulations give results comparable to
those obtained with standard WENO schemes [17].

Finally, in Sec. VII, the entropy method is applied to a specific BNS merger
simulation. The new feature that we add with the present work to the BNS
simulations carried out in [9] is that the EFL method is used in the
construction of the SGRID initial data through the process of eccentricity
reduction. As in [9], the EFL scheme can be used to successfully evolve this
initial data and the properties observed in the single star tests carry over
to this simulation where the spacetime is non-stationary and the neutron stars
are moving freely in the computational grid. As shown in FIG. 16, the entropy
limiter locates the surface of the inspiraling NSs extremely accurately and it
converges to zero in regions of smooth flows. Further, it captures the
collisional shocks at merger and the outward dynamics of the spiral density
waves, thus being robust also for postmerger evolutions.

A convergence study of the gravitational waveforms obtained from this
simulation shows that when the EFL method is used both in the construction of
the initial data and in the temporal evolution, then it can deliver fifth-
order convergent waveforms at current production resolutions (subsection
VII.3). Such a convergence is measured in the $(2,2)$ dominant mode of the
strain but also in the next subdominant modes $(3,2)$ and $(4,4)$. To our
knowledge, these are the first results in which fifth-order convergence is
demonstrated. The estimated phase error in the waveform is about a factor
${\sim}2$ smaller than the error in the state-of-the-art high-order WENOZ
scheme used in the same BAM code, at the same resolution. Thus, we conclude
that the use of the EFL scheme in the construction of the initial data can
have positive impact on the accuracy and the convergence properties of BNS
simulations and can be used for high-quality waveform production.

Naturally the question arises of whether the entropy pair introduced in Sec.
II and used to define the entropy of systems without EoS can be also used for
systems with a well defined thermodynamic entropy like the ones of Sec. V,
Sec. VI and Sec. VII. Such a possibility would highly simplify the use of the
EFL scheme as the same entropy pair could be used for a wide variate of
different conservation laws. The investigation of this possibility is
currently in progress and results will be published as soon as available.

###### Acknowledgements.

We would like to thank members of the Jena group for fruitful discussions and
invaluable input. GD acknowledges funding from the European High Performance
Computing Joint Undertaking (JU) and Belgium, Czech Republic, France, Germany,
Greece, Italy, Norway, and Spain under grant agreement No 101093441 (SPACE)
and from the EU Horizon under ERC Starting Grant, grant agreement no.
BinGraSp-714626. SB acknowledges funding from the EU Horizon under ERC
Starting Grant, grant agreement no. BinGraSp-714626, and ERC Consolidator
Grant, grant agreement no. InspiReM-101043372. WT acknowledges support by the
National Science Foundation under grant PHY-2136036. Computations where
performed on the ARA and DRACO clusters at Friedrich Schiller University Jena.
The ARA cluster is partially funded by the DFG grants INST 275/334-1 FUGG and
INST 275/363-1 FUGG, and ERC Starting Grant, grant agreement no.
BinGraSp-714626. The authors also gratefully acknowledge the Gauss Centre for
Supercomputing e.V. (www.gauss-centre.eu) for funding this project by
providing computing time on the GCS Supercomputer SuperMUC-NG at Leibniz
Supercomputing Centre (www.lrz.de) with the allocations pn36ge and pn36jo.

## References

  * Abbott  _et al._ [2017a] B. P. Abbott _et al._ (Virgo, LIGO Scientific), GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral, Phys. Rev. Lett. 119, 161101 (2017a), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ [2017b] B. P. Abbott _et al._ (GROND, SALT Group, OzGrav, DFN, INTEGRAL, Virgo, Insight-Hxmt, MAXI Team, Fermi-LAT, J-GEM, RATIR, IceCube, CAASTRO, LWA, ePESSTO, GRAWITA, RIMAS, SKA South Africa/MeerKAT, H.E.S.S., 1M2H Team, IKI-GW Follow-up, Fermi GBM, Pi of Sky, DWF (Deeper Wider Faster Program), Dark Energy Survey, MASTER, AstroSat Cadmium Zinc Telluride Imager Team, Swift, Pierre Auger, ASKAP, VINROUGE, JAGWAR, Chandra Team at McGill University, TTU-NRAO, GROWTH, AGILE Team, MWA, ATCA, AST3, TOROS, Pan-STARRS, NuSTAR, ATLAS Telescopes, BOOTES, CaltechNRAO, LIGO Scientific, High Time Resolution Universe Survey, Nordic Optical Telescope, Las Cumbres Observatory Group, TZAC Consortium, LOFAR, IPN, DLT40, Texas Tech University, HAWC, ANTARES, KU, Dark Energy Camera GW-EM, CALET, Euro VLBI Team, ALMA), Multi-messenger Observations of a Binary Neutron Star Merger, Astrophys. J. 848, L12 (2017b), arXiv:1710.05833 [astro-ph.HE] .
  * Toro [1999] E. F. Toro,  _Riemann Solvers and Numerical Methods for Fluid Dynamics_ , 2nd ed. (Springer-Verlag, 1999).
  * Bernuzzi  _et al._ [2012] S. Bernuzzi, A. Nagar, M. Thierfelder, and B. Brügmann, Tidal effects in binary neutron star coalescence, Phys.Rev. D86, 044030 (2012), arXiv:1205.3403 [gr-qc] .
  * Radice  _et al._ [2014a] D. Radice, L. Rezzolla, and F. Galeazzi, Beyond second-order convergence in simulations of binary neutron stars in full general-relativity, Mon.Not.Roy.Astron.Soc. 437, L46 (2014a), arXiv:1306.6052 [gr-qc] .
  * Radice  _et al._ [2014b] D. Radice, L. Rezzolla, and F. Galeazzi, High-Order Fully General-Relativistic Hydrodynamics: new Approaches and Tests, Class.Quant.Grav. 31, 075012 (2014b), arXiv:1312.5004 [gr-qc] .
  * Bernuzzi and Dietrich [2016] S. Bernuzzi and T. Dietrich, Gravitational waveforms from binary neutron star mergers with high-order weighted-essentially-nonoscillatory schemes in numerical relativity, Phys. Rev. D94, 064062 (2016), arXiv:1604.07999 [gr-qc] .
  * Sweby [1984] P. K. Sweby, High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws, SIAM Journal on Numerical Analysis 21, 995 (1984).
  * Doulis  _et al._ [2022] G. Doulis, F. Atteneder, S. Bernuzzi, and B. Brügmann, Entropy-limited higher-order central scheme for neutron star merger simulations, (2022), arXiv:2202.08839 [gr-qc] .
  * Guermond and Pasquetti [2008] J.-L. Guermond and R. Pasquetti, Entropy-based nonlinear viscosity for fourier approximations of conservation laws, Comptes Rendus Mathematique 346, 801 (2008).
  * Guermond  _et al._ [2011] J.-L. Guermond, R. Pasquetti, and B. Popov, Entropy viscosity method for nonlinear conservation laws, Journal of Computational Physics 230, 4248 (2011), special issue High Order Methods for {CFD} Problems.
  * Guercilena  _et al._ [2017] F. Guercilena, D. Radice, and L. Rezzolla, Entropy-limited hydrodynamics: a novel approach to relativistic hydrodynamics, Comput. Astrophys. Cosmol. 4, 3 (2017), arXiv:1612.06251 [gr-qc] .
  * Gourgoulhon  _et al._ [2001] E. Gourgoulhon, P. Grandclement, K. Taniguchi, J.-A. Marck, and S. Bonazzola, Quasiequilibrium sequences of synchronized and irrotational binary neutron stars in general relativity: 1. Method and tests, Phys.Rev. D63, 064029 (2001), arXiv:gr-qc/0007028 [gr-qc] .
  * Tichy [2012] W. Tichy, Constructing quasi-equilibrium initial data for binary neutron stars with arbitrary spins, Phys. Rev. D 86, 064024 (2012), arXiv:1209.5336 [gr-qc] .
  * Banyuls  _et al._ [1997] F. Banyuls, J. A. Font, J. M. A. Ibanez, J. M. A. Marti, and J. A. Miralles, Numerical 3+1 General Relativistic Hydrodynamics: A Local Characteristic Approach, Astrophys. J. 476, 221 (1997).
  * Mignone  _et al._ [2010] A. Mignone, P. Tzeferacos, and G. Bodo, High-order conservative finite difference GLM-MHD schemes for cell-centered MHD, J.Comput.Phys. 229, 5896 (2010), arXiv:1001.2832 [astro-ph.HE] .
  * Thierfelder  _et al._ [2011] M. Thierfelder, S. Bernuzzi, and B. Brügmann, Numerical relativity simulations of binary neutron stars, Phys.Rev. D84, 044012 (2011), arXiv:1104.4751 [gr-qc] .
  * Jiang [1996] G. Jiang, Efficient Implementation of Weighted ENO Schemes, J. Comp. Phys. 126, 202 (1996).
  * Borges  _et al._ [2008] R. Borges, M. Carmona, B. Costa, and W. S. Don, An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws, Journal of Computational Physics 227, 3191 (2008).
  * Liu and Osher [1998] X. Liu and S. Osher, Convex eno high order multi-dimensional schemes without field by field decomposition or staggered grids, J. Comput. Phys. 142, 304 (1998).
  * Del Zanna, L. and Bucciantini, N. [2002] Del Zanna, L. and Bucciantini, N., An efficient shock-capturing central-type scheme for multidimensional relativistic flows - i. hydrodynamics, Astron. Astrophys. 390, 1177 (2002), astro-ph/0205290 .
  * Hesthaven [2018] J. S. Hesthaven,  _Numerical Methods for Conservation Laws_ (Society for Industrial and Applied Mathematics, Philadelphia, PA, 2018) https://epubs.siam.org/doi/pdf/10.1137/1.9781611975109 .
  * Brügmann  _et al._ [2008] B. Brügmann, J. A. Gonzalez, M. Hannam, S. Husa, U. Sperhake, _et al._ , Calibration of Moving Puncture Simulations, Phys.Rev. D77, 024027 (2008), arXiv:gr-qc/0610128 [gr-qc] .
  * Dietrich  _et al._ [2015a] T. Dietrich, S. Bernuzzi, M. Ujevic, and B. Brügmann, Numerical relativity simulations of neutron star merger remnants using conservative mesh refinement, Phys. Rev. D91, 124041 (2015a), arXiv:1504.01266 [gr-qc] .
  * Berger and Oliger [1984] M. J. Berger and J. Oliger, Adaptive Mesh Refinement for Hyperbolic Partial Differential Equations, J.Comput.Phys. 53, 484 (1984).
  * Nakamura  _et al._ [1987] T. Nakamura, K. Oohara, and Y. Kojima, General Relativistic Collapse to Black Holes and Gravitational Waves from Black Holes, Prog. Theor. Phys. Suppl. 90, 1 (1987).
  * Shibata and Nakamura [1995] M. Shibata and T. Nakamura, Evolution of three-dimensional gravitational waves: Harmonic slicing case, Phys. Rev. D52, 5428 (1995).
  * Baumgarte and Shapiro [1999] T. W. Baumgarte and S. L. Shapiro, On the numerical integration of Einstein’s field equations, Phys. Rev. D59, 024007 (1999), arXiv:gr-qc/9810065 .
  * Bernuzzi and Hilditch [2010] S. Bernuzzi and D. Hilditch, Constraint violation in free evolution schemes: comparing BSSNOK with a conformal decomposition of Z4, Phys. Rev. D81, 084003 (2010), arXiv:0912.2920 [gr-qc] .
  * Hilditch  _et al._ [2013] D. Hilditch, S. Bernuzzi, M. Thierfelder, Z. Cao, W. Tichy, and B. Bruegmann, Compact binary evolutions with the Z4c formulation, Phys. Rev. D88, 084057 (2013), arXiv:1212.2901 [gr-qc] .
  * Douchin and Haensel [2001] F. Douchin and P. Haensel, A unified equation of state of dense matter and neutron star structure, Astron. Astrophys. 380, 151 (2001), astro-ph/0111092 .
  * Read  _et al._ [2009] J. S. Read, B. D. Lackey, B. J. Owen, and J. L. Friedman, Constraints on a phenomenologically parameterized neutron- star equation of state, Phys. Rev. D79, 124032 (2009), arXiv:0812.2163 [astro-ph] .
  * Shibata  _et al._ [2005] M. Shibata, K. Taniguchi, and K. Uryu, Merger of binary neutron stars with realistic equations of state in full general relativity, Phys. Rev. D71, 084021 (2005), arXiv:gr-qc/0503119 .
  * Bauswein  _et al._ [2010] A. Bauswein, H.-T. Janka, and R. Oechslin, Testing Approximations of Thermal Effects in Neutron Star Merger Simulations, Phys.Rev. D82, 084043 (2010), arXiv:1006.3315 [astro-ph.SR] .
  * Liang [1977] E. P. T. Liang, Relativistic simple waves - Shock damping and entropy production, Astrophys. J. 211, 361 (1977).
  * Anile [1990] A. M. Anile,  _Relativistic Fluids and Magneto-fluids, by A. M. Anile, pp. 348. ISBN 0521304067. Cambridge, UK: Cambridge University Press, February 1990._ (1990).
  * Sod [1978] G. Sod, A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws, J. Comput. Phys. 27, 1 (1978).
  * Marti and Müller [1999] J. M. Marti and E. Müller, Numerical hydrodynamics in special relativity, Living Rev. Rel. 2, 3 (1999), arXiv:astro-ph/9906333 .
  * Font  _et al._ [2000] J. A. Font, M. A. Miller, W.-M. Suen, and M. Tobias, Three Dimensional Numerical General Relativistic Hydrodynamics I: Formulations, Methods, and Code Tests, Phys. Rev. D61, 044011 (2000), gr-qc/9811015 .
  * Font  _et al._ [2002] J. A. Font _et al._ , Three-dimensional general relativistic hydrodynamics. II: Long-term dynamics of single relativistic stars, Phys. Rev. D65, 084024 (2002), arXiv:gr-qc/0110047 .
  * Cordero-Carrion  _et al._ [2009] I. Cordero-Carrion _et al._ , An improved constrained scheme for the Einstein equations: an approach to the uniqueness issue, Phys. Rev. D79, 024017 (2009), arXiv:0809.2325 [gr-qc] .
  * Dietrich  _et al._ [2017] T. Dietrich, S. Bernuzzi, and W. Tichy, Closed-form tidal approximants for binary neutron star gravitational waveforms constructed from high-resolution numerical relativity simulations, Phys. Rev. D96, 121501 (2017), arXiv:1706.02969 [gr-qc] .
  * Dietrich  _et al._ [2018] T. Dietrich, D. Radice, S. Bernuzzi, F. Zappa, A. Perego, B. Brügmann, S. V. Chaurasia, R. Dudi, W. Tichy, and M. Ujevic, CoRe database of binary neutron star merger waveforms, Class. Quant. Grav. 35, 24LT01 (2018), arXiv:1806.01625 [gr-qc] .
  * Dietrich  _et al._ [2015b] T. Dietrich, N. Moldenhauer, N. K. Johnson-McDaniel, S. Bernuzzi, C. M. Markakis, B. Brügmann, and W. Tichy, Binary Neutron Stars with Generic Spin, Eccentricity, Mass ratio, and Compactness - Quasi-equilibrium Sequences and First Evolutions, Phys. Rev. D92, 124007 (2015b), arXiv:1507.07100 [gr-qc] .
  * Tichy [2006] W. Tichy, Black hole evolution with the BSSN system by pseudo-spectral methods, Phys.Rev. D74, 084005 (2006), arXiv:gr-qc/0609087 [gr-qc] .
  * Tichy [2009a] W. Tichy, A New numerical method to construct binary neutron star initial data, Class.Quant.Grav. 26, 175018 (2009a), arXiv:0908.0620 [gr-qc] .
  * Tichy [2009b] W. Tichy, Long term black hole evolution with the BSSN system by pseudo-spectral methods, Phys.Rev. D80, 104034 (2009b), arXiv:0911.0973 [gr-qc] .
  * Tichy [2011] W. Tichy, Initial data for binary neutron stars with arbitrary spins, Phys.Rev. D84, 024041 (2011), arXiv:1107.1440 [gr-qc] .
  * Kyutoku  _et al._ [2014] K. Kyutoku, M. Shibata, and K. Taniguchi, Reducing orbital eccentricity in initial data of binary neutron stars, Phys. Rev. D90, 064006 (2014), arXiv:1405.6207 [gr-qc] .
  * Reisswig and Pollney [2011] C. Reisswig and D. Pollney, Notes on the integration of numerical relativity waveforms, Class.Quant.Grav. 28, 195015 (2011), arXiv:1006.1632 [gr-qc] .
  * Baumgarte and Shapiro [2010] T. Baumgarte and S. Shapiro,  _Numerical Relativity_ (Cambridge University Press, Cambridge, 2010).