Spritz: general relativistic magnetohydrodynamics with neutrinos

In order to properly include neutrino emission and absorption, we need to add one evolution equation for the electron fraction, which we define as

$$\begin{equation}{Y}_{\mathrm{e}}=\frac{{n}_{\mathrm{e}}}{{n}_{\mathrm{p}}+{n}_{\mathrm{n}}},\end{equation} \tag{ 5 }$$

being n_e, n_p, and n_n the electron, proton, and neutron number densities.

From the local conservation of the total baryon number, neglecting the mass difference between neutrons and protons, we obtain the following equation for the electron fraction, valid in absence of neutrino emission/absorption:

$$\begin{equation}{\nabla }_{\mu }\left({Y}_{\mathrm{e}}\rho {u}^{\mu }\right)=0,\end{equation} \tag{ 6 }$$

expressing the fact that Y_e is advected along the fluid lines. This equation is commonly referred to as the electron fraction advection and can be expressed in a hyperbolic conservative form as

$$\begin{equation}{\partial }_{t}\left(\sqrt{\gamma }D{Y}_{\mathrm{e}}\right)+{\partial }_{i}\left[\alpha \sqrt{\gamma }D{Y}_{\mathrm{e}}\left({v}^{i}-\frac{{\beta }^{i}}{\alpha }\right)\right]=0.\end{equation} \tag{ 7 }$$

In presence of reactions involving neutrinos, the local electron fraction obtained from the above equation is then modified according to equation (19) (see section 2.3).

The Spritz code can handle tabulated finite-temperature and composition dependent EOS via the EOS_Omni thorn included in the Einstein Toolkit. This is crucial since a proper description of the matter composition depending on temperature is necessary in order to estimate the emission and absorption rates associated with the different processes involving neutrinos (see the next section). Moreover, as a consequence of such processes, Y_e necessarily undergoes changes that must be estimated accurately when dealing with dynamical scenarios.

The exact matter composition at the typical densities reached in the core of an NS is still unknown and so is the correct EOS. A large number of proposed tabulated EOS inspired by nuclear physics calculations can be found in the literature (see, e.g. the database in [32, 37] for several examples). These EOS are usually three-dimensional tables where every hydrodynamical variable, such as the gas pressure P or the specific internal energy ɛ, can be related to the rest-mass density ρ, the temperature T, and the electron fraction Y_e.

When building initial data, however, a one-dimensional (i.e. barotropic) EOS is typically needed, where P is just a function of ρ. In this case, reducing the three-dimensional table P = P(ρ, T, Y_e) to a simpler one-dimensional relation P = P(ρ) becomes necessary, implying that two conditions on the NS matter should be imposed. The first and most common one is to assume the NS to be initially in β-equilibrium, which is a reasonable assumption for old NSs, such as those encountered in BNS or NSBH binary systems prior to merger. As a second assumption, one may decide to fix either a constant value for the entropy (S-slicing condition) or for the temperature (T-slicing condition). The latter is the one typically used in BNS or NSBH merger simulations since it is reasonable to expect NSs to be cold prior to merger. In this paper, along with the standard T-slicing condition, we have also used the S-slicing condition to test the ability of our code in dealing with 'hot' NSs.

All the computations presented in this paper are performed adopting the LS220 EOS [38], that has been already used in a number of papers dealing with the evolution of BNS systems (e.g. [16, 39, 40]).

During the merger of BNS or NSBH systems, temperatures as high as T ∼ 10 MeV ∼ 10¹¹ K can be produced and also the electron fraction Y_e may change considerably. In this scenario, neutrinos play a key role in both the transport of energy and in determining the evolution of Y_e and temperature, which are in turn crucial parameters for the r-process nucleosynthesis taking place in the ejected matter and the subsequent production of heavy elements. A proper estimate of the rates of the different reactions involving neutrinos is thus necessary in order to compute the nucleosythesis yields and to model the radioactively-powered kilonova signals accompanying such mergers (as the one already observed after GW170817; e.g. [14, 15]).

The typical timescale for weak processes producing neutrinos can be estimated from the changing electron fraction as

$$\begin{equation}{t}_{\text{WP}}\sim \left\vert \frac{{Y}_{\mathrm{e}}}{{\dot {Y}}_{\mathrm{e}}}\right\vert \ll {t}_{\text{dyn}},\end{equation} \tag{ 8 }$$

being t_dyn the dynamical timescale of the simulated astrophysical event [41]. By carrying away energy, neutrinos can significantly cool down the (meta)stable NS remnant of a BNS merger or the accretion disk around the spinning BH resulting from either a BNS or an NSBH merger (e.g. [16]). Moreover, a fraction of the emitted neutrinos may be reabsorbed by the outer material, inducing heating and leptonization of the material itself. The surface where the neutrino optical depth is τ = 2/3 conventionally defines the 'neutrinosphere' (e.g. [42]), which separates the diffusive regime of the high-density interiors (≳10¹² g cm⁻³; e.g. [43]) and the nearly free streaming regime of the exterior. The intermediate region between τ ≪ 1 and τ ≫ 1 (i.e. where neutrinos are neither free to escape nor fully trapped) is the challenging one for neutrino transport. In its energy averaged version, the optical depth along each path ξ followed by neutrinos can be defined as [44]

$$\begin{equation}{\tau }_{\xi }={\int }_{\xi }\rho \left(\mathbf{x}\right)k\left(\mathbf{x}\right)\sqrt{{\gamma }_{ij}\enspace \mathrm{d}{x}^{i}\enspace \mathrm{d}{x}^{j}},\end{equation} \tag{ 9 }$$

being k(x) the energy averaged opacity at position x. The path giving the minimum optical depth is the favoured one for neutrino escape and allows us to define a single optical depth for each given location

$$\begin{equation}\tau \left(\mathbf{x}\right)=\underset{\xi \in {\Xi}}{\mathrm{min}}\enspace {\tau }_{\xi }=\underset{\xi \in {\Xi}}{\mathrm{min}}{\int }_{\xi }\rho \left(\mathbf{x}\right)k\left(\mathbf{x}\right)\sqrt{{\gamma }_{ij}\enspace \mathrm{d}{x}^{i}\enspace \mathrm{d}{x}^{j}},\end{equation} \tag{ 10 }$$

where Ξ is the set of all possible paths including position x.

The complexity and extremely high computational cost of the full neutrino transport problem solved via the Boltzmann radiation transport equations forced the introduction of approximate schemes and simplifying assumptions (e.g. [33] and references therein). We consider here a so-called neutrino leakage scheme, already employed successfully in BNS and NSBH simulations (e.g. [25, 41, 45–47]). In particular, we adopt the leakage method presented in [44, 48], which has been implemented in the publicly available ZelmaniLeak code [32]. In what follows, we introduce the leakage scheme and the basic physical assumptions. The numerical implementation is instead discussed in the next section (and in particular in section 3.4).

In the neutrino leakage scheme adopted in this work, we consider three neutrino species, electron neutrino ν_e, electron antineutrino ${\bar{\nu }}_{\mathrm{e}}$, and heavy-lepton neutrinos ν_x (including ν_μ , ${\bar{\nu }}_{\mu }$, ν_τ , ${\bar{\nu }}_{\tau }$), and for each one we compute the local number and energy emission rates according to the following steps.

The neutrino optical depths, which are crucial to determine the emission rates (see below), are computed under the assumption that neutrinos escape along radial paths from the centre (ray-by-ray approach). For each species, we compute the local spectral averaged opacity as the sum of the opacities due to the scattering off nucleons, neutrino-nucleus scattering, and neutrino absorption by free nucleons (see [49] for details). Then, we use these mean opacities to compute the optical depths along each radial path (equation (9)).

In the diffusive regime, the number and energy rates (i.e. number and energy per unit volume, per unit time) can be written as [49]

$$\begin{equation}{R}_{{\nu }_{i}}^{\text{diff}}=\frac{4\pi c{g}_{{\nu }_{i}}}{{\left(hc\right)}^{3}}\frac{{\zeta }_{{\nu }_{i}}}{3{\chi }_{{\nu }_{i}}^{2}}T{F}_{0}\left({\eta }_{{\nu }_{i}}\right),\end{equation} \tag{ 11 }$$

$$\begin{equation}{Q}_{{\nu }_{i}}^{\text{diff}}=\frac{4\pi c{g}_{{\nu }_{i}}}{{\left(hc\right)}^{3}}\frac{{\zeta }_{{\nu }_{i}}}{3{\chi }_{{\nu }_{i}}^{2}}{T}^{2}{F}_{1}\left({\eta }_{{\nu }_{i}}\right),\end{equation} \tag{ 12 }$$

where i = 1, 2, 3 and ν₁ = ν_e, ${\nu }_{2}={\bar{\nu }}_{\mathrm{e}}$, ν₃ = ν_x , while ${g}_{{\nu }_{1}}={g}_{{\nu }_{2}}=1$ and ${g}_{{\nu }_{3}}=4$. Moreover, $\zeta ={\left({E}^{2}\lambda \right)}^{-1}$, χ = τ/E², with E the average neutrino energy (computed assuming a Fermi–Dirac distribution at the local temperature T) and λ the mean free path, and F₀(η), F₁(η) are the Fermi integrals defined in [50] as function of the neutrino chemical potential η. Energy and number rates are also computed for the free neutrino emission regime (${Q}_{{\nu }_{i}}^{\text{free}}$ and ${R}_{{\nu }_{i}}^{\text{free}}$), taking into account capture processes, electron–positron pair annihilation, plasmon decay, and nucleon–nucleon bremsstrahlung (see [48, 49]). Finally, the actual emission rates are found by combining the free emission and diffusive ones as follows

$$\begin{equation}{R}_{{\nu }_{i}}^{\text{eff}}={R}_{{\nu }_{i}}^{\text{free}}\left(1+\frac{{R}_{{\nu }_{i}}^{\text{free}}}{{R}_{{\nu }_{i}}^{\text{diff}}}\right),\end{equation} \tag{ 13 }$$

$$\begin{equation}{Q}_{{\nu }_{i}}^{\text{eff}}={Q}_{{\nu }_{i}}^{\text{free}}\left(1+\frac{{Q}_{{\nu }_{i}}^{\text{free}}}{{Q}_{{\nu }_{i}}^{\text{diff}}}\right).\end{equation} \tag{ 14 }$$

For a given radial direction (θ, ϕ), the isotropic-equivalent neutrino luminosity incoming from below at a distance r can be computed (in the coordinate frame) as

$$\begin{equation}\begin{aligned}\hfill {L}_{{\nu }_{i}}^{\text{iso}}\left(r,\theta ,\phi \right)& =4\pi {\int }_{0}^{r}\left[\frac{\alpha \left({r}^{\prime },\theta ,\phi \right)}{\alpha \left(r,\theta ,\phi \right)}\right]{Q}_{{\nu }_{i}}^{\text{eff}}\left({r}^{\prime },\theta ,\phi \right)\alpha \left({r}^{\prime },\theta ,\phi \right)W\left({r}^{\prime },\theta ,\phi \right)\hfill \\ \hfill & \quad {\times}\left[1+{v}^{r}\left({r}^{\prime },\theta ,\phi \right)\right]\sqrt{{g}_{rr}\left({r}^{\prime },\theta ,\phi \right)}{{r}^{\prime }}^{2}\enspace \mathrm{d}{r}^{\prime },\hfill \end{aligned}\end{equation} \tag{ 15 }$$

being v^r the radial velocity ¹² . We can also define a fluid rest frame (FRF) luminosity as

$$\begin{equation}{L}_{{\nu }_{i}}^{\text{iso,FRF}}\left(r\right)=\frac{{L}_{{\nu }_{i}}^{\text{iso}}\left(r\right)}{\alpha \left(r\right)W\left(r\right)\left[1+{v}^{r}\left(r\right)\right]}\enspace .\end{equation} \tag{ 16 }$$

The heating and leptonization due to the reabsorption of a fraction of neutrinos by the material along their path (i.e. ν_e and ${\bar{\nu }}_{\mathrm{e}}$ reabsorption on neutrons and protons, respectively) is taken into account via the local heating rate [48]

$$\begin{align}\hfill {Q}_{\left({\nu }_{\mathrm{e}},{\bar{\nu }}_{\mathrm{e}}\right)}^{\text{heat}}& ={f}_{\text{heat}}\frac{{L}_{\left({\nu }_{\mathrm{e}},{\bar{\nu }}_{\mathrm{e}}\right)}^{\text{iso,FRF}}}{4\pi {r}^{2}}{\sigma }_{\left({\nu }_{\mathrm{e}},{\bar{\nu }}_{\mathrm{e}}\right)}^{\mathrm{h}\mathrm{e}\mathrm{a}\mathrm{t}}\frac{\rho }{{m}_{\left(n,p\right)}}{X}_{\left(\mathrm{n},\mathrm{p}\right)}\left(4.275{\tau }_{\left({\nu }_{\mathrm{e}},{\bar{\nu }}_{\mathrm{e}}\right)}+1.15\right){\mathrm{e}}^{-2{\tau }_{\left({\nu }_{\mathrm{e}},{\bar{\nu }}_{\mathrm{e}}\right)}},\hfill \end{align} \tag{ 17 }$$

where f_heat is a scaling factor of order one (we set f_heat = 1), ${\sigma }_{\left({\nu }_{\mathrm{e}},{\bar{\nu }}_{\mathrm{e}}\right)}^{\mathrm{h}\mathrm{e}\mathrm{a}\mathrm{t}}$ is the reabsorption cross-section (see below), m_(n,p) and X_(n,p) are the neutron or proton masses and mass fractions, and the factor ${\mathrm{e}}^{-2{\tau }_{\left({\nu }_{\mathrm{e}},{\bar{\nu }}_{\mathrm{e}}\right)}}$ is added to suppress heating at very large optical depths. For the reabsorption cross-section, we adopt the following expression [44]

$$\begin{equation}{\sigma }_{\left({\nu }_{\mathrm{e}},{\bar{\nu }}_{\mathrm{e}}\right)}^{\text{heat}}=\frac{1+3{\alpha }_{\text{EC}}^{2}}{4}{\sigma }_{0}\frac{{\langle {E}^{2}\rangle }_{\left({\nu }_{\mathrm{e}},{\bar{\nu }}_{\mathrm{e}}\right)}^{\mathrm{N}\mathrm{S}}}{{\left({m}_{\mathrm{e}}{c}^{2}\right)}^{2}}\langle 1-{f}_{\left({\mathrm{e}}^{-},{\mathrm{e}}^{+}\right)}\rangle ,\end{equation} \tag{ 18 }$$

where α_EC = −1.25, σ₀ = 1.76 × 10⁻⁴⁴ cm², ${\langle {E}^{2}\rangle }^{\mathrm{N}\mathrm{S}}$ is the mean squared neutrino energy at the neutrinosphere, and $\langle 1-{f}_{\left({\mathrm{e}}^{-},{\mathrm{e}}^{+}\right)}\rangle$ are the blocking factors defined in [51].

The full neutrino emission and reabsorption problem at a given time is solved along each radial direction by moving outwards from the centre and, at each radius, subtracting the heating rate from the emission rate, i.e. ${Q}_{{\nu }_{i}}^{\text{eff}}\to {Q}_{{\nu }_{i}}^{\text{eff}}-{Q}_{{\nu }_{i}}^{\text{heat}}$ and ${R}_{{\nu }_{i}}^{\text{eff}}\to {R}_{{\nu }_{i}}^{\text{eff}}-{Q}_{{\nu }_{i}}^{\text{heat}}/{\langle E\rangle }_{{\nu }_{i}}^{\mathrm{N}\mathrm{S}}$, with ${\langle E\rangle }_{{\nu }_{i}}^{\mathrm{N}\mathrm{S}}$ the average neutrino energy at the neutrinosphere and ${Q}_{{\nu }_{x}}^{\text{heat}}=0$ ¹³ .

In order to couple the result to the GRMHD evolution, the Y_e and ɛ are then modified as follows:

$$\begin{equation}{Y}_{\mathrm{e}}\to {Y}_{\mathrm{e}}+{\Delta}t\frac{\partial {Y}_{\mathrm{e}}}{\partial t},\end{equation} \tag{ 19 }$$

being Δt the local time step, and where

$$\begin{equation}\frac{\partial {Y}_{\mathrm{e}}}{\partial t}=\frac{{R}_{\bar{{\nu }_{\mathrm{e}}}}^{\text{eff}}-{R}_{{\nu }_{\mathrm{e}}}^{\text{eff}}}{\rho }{m}_{\mathrm{n}},\end{equation} \tag{ 20 }$$

being m_n the rest-mass of the neutron, and

$$\begin{equation}\varepsilon \to \varepsilon +{\Delta}t\frac{\partial \varepsilon }{\partial t},\end{equation} \tag{ 21 }$$

where

$$\begin{equation}\frac{\partial \varepsilon }{\partial t}=-\frac{{{\Sigma}}_{i}{Q}_{{\nu }_{i}}^{\text{eff}}}{\rho }.\end{equation} \tag{ 22 }$$

As already stated in section 2.2, the Spritz code adopts the EOS_Omni thorn of the Einstein Toolkit software infrastructure. This thorn is able to handle a large variety of EOS, including ideal fluid, polytropic, and tabulated ones.

During our first tests with the EOS_Omni thorn and tabulated EOS, we noticed that the EOS_Omni thorn presented some limitations in dealing with such EOS type. In particular, to compute the temperature T from the specific internal energy ɛ, the thorn adopts a Newton–Raphson routine with a fall-back to a bisection routine in case of too many iterations, after verifying that the root is bracketed. We found this algorithm to be not robust enough in cases when T weakly depends on ɛ, which may lead T to go out of the bounds present in the chosen table (see [52]). This problem was present in particular when dealing with NS initial data using the lowest T available in the EOS table. Such initial data undergo a sharp temperature increase in the core of the NS due to numerical readjustment of the initial data given by the solution of the TOV equations. We proposed a modification of the EOS_Omni thorn to the Einstein Toolkit developers that consisted in preferring the fall-back to the more robust bisection method in such cases. In this way, we verified the temperature T to be always contained in the range available in the table. This modification was accepted and it is now included in the publicly available Einstein Toolkit since May 2020 [53].

We performed all the simulations discussed in section 4 using this new version of the EOS_Omni thorn. We therefore caution the reader that the Spritz code should be used with the May 2020 release of the Einstein Toolkit (or later versions) when using tabulated EOS.

In order to compute the initial data, one needs to reduce the 3D EOS table to a 1D EOS, in which the pressure P is only a function of the rest-mass density ρ. To do this, we assume β-equilibrium and then apply the S-slicing or T-slicing condition mentioned in section 2.2. We coded a python script for this purpose (available with the public version of Spritz) that produces a 1D tabulated EOS starting from a 3D tabulated EOS in .h5 format, such as the ones provided in [32]. The 1D EOS is saved in the CompOSE format [37] that can be easily used with LORENE [54]. The initial data used in this paper, reproducing a single non-rotating NS (TOV), were in particular produced with the code Nrotstar that can compute equilibrium solutions for non-rotating or uniformly rotating NSs. These solutions are non-magnetized, but a magnetic field can be easily added to the initial data as long as the field strength is ≲10¹⁷ G, such that no significant effects on the NS structure nor significant violations of the constraint equations are introduced.

To read the initial data in the Spritz code we developed the ID_Nrotstar thorn which is simply a reader that makes use of the LORENE library to read initial data produced with Nrotstar and import them in the Cartesian grid used by the code. Since the initial data were produced assuming β-equilibrium, we also developed an additional thorn, Spritz_SetBeta, that instead makes sure that, when computing the conservative variables from the primitive ones at iteration 0, the code uses the same 1D EOS used to compute the initial data. After the initial data are correctly imported and conserved variables computed, the evolution starts and the full 3D EOS table is used.

Both ID_Nrotstar and Spritz_SetBeta are part of version 1.1.0 of the Spritz code [28].

When using tabulated EOS we employ the 1D method for the conservative-to-primitive inversion presented by Palenzuela et al [25]. This method is a modification to the 1D method already used in GRHD [52]. It consists of rewriting the conserved variables in the following way

$$\begin{equation}q\equiv \frac{\tilde {\tau }}{D},\hspace{25.0pt}r\equiv \frac{{S}^{2}}{{D}^{2}},\hspace{25.0pt}s\equiv \frac{{B}^{2}}{D},\hspace{25.0pt}t\equiv \frac{{B}_{i}{S}^{i}}{{D}^{\frac{3}{2}}},\end{equation} \tag{ 23 }$$

and searching for the independent variable x ≡ hW. One then looks for the solution of f(x) = 0, where f(x) = x − hW. We point the reader to [25] for more details about the algorithm. Here, it is only important to note that the Brent's method [55] is used for the root finding, where the independent variable x should be properly bracketed, thus $\left.x\in \right]{x}_{\text{L}},{x}_{\text{R}}\left[\right.$, with $f\left({x}_{\text{L}}\right)\cdot f\left({x}_{\text{R}}\right){< }0$. The left and right bounds can be defined in the following way (see [56]):

$$\begin{equation}\begin{aligned}\hfill {x}_{\text{L}}& =1+q-s,\hfill \\ \hfill {x}_{\text{R}}& =2+2q-s.\hfill \end{aligned}\end{equation} \tag{ 24 }$$

If no consistent bound is found, then the point is set to atmosphere.

As we will show in section 4, we are also interested in performing simulations where the initial temperature T is forced to be constant. This may be useful in order to avoid spurious neutrino production in particular scenarios, e.g. during BNS inspiral (for some examples, see [57–59]) or when evolving a single cold NS (that may undergo a sharp initial rise of temperature as already mentioned in section 3.1). However, the aforementioned 1D conservative-to-primitive scheme cannot be used in such cases and we adopt a modification of the 3eqs method that was already implemented in the Spritz code (see [27, 60] for details), where the $\tilde {\tau }$ variable is not used in computing the primitive variables. In particular, we refer to equation (45) of [60], defining the function

$$\begin{equation}f\left({W}_{\text{guess}}\right)\equiv {S}^{2}-\left[{\left(\hat{Z}+{B}^{2}\right)}^{2}\frac{{W}_{\text{guess}}^{2}-1}{{W}_{\text{guess}}^{2}}-\frac{2\hat{Z}+{B}^{2}}{{\hat{Z}}^{2}}{\left({B}^{i}{S}_{i}\right)}^{2}\right],\end{equation} \tag{ 25 }$$

where

$$\begin{equation}\hat{Z}={W}_{\text{guess}}^{2}\left(\hat{\rho }+\hat{\rho }\hat{\varepsilon }+\hat{P}\right),\end{equation} \tag{ 26 }$$

$$\begin{equation}\hat{\rho }=\frac{D}{{W}_{\text{guess}}},\end{equation} \tag{ 27 }$$

and $\hat{P}$ and $\hat{\varepsilon }$ can be computed via the EOS using $\hat{\rho }$ and the constrained value of T. The algorithm proceeds as follows:

• (a)  
  The initial guess for the solution is assumed to be W_guess ∈ [1.0, 1.5]; ¹⁴
• (b)  
  $\hat{\rho }$, $\hat{P}$, $\hat{\varepsilon }$, and $\hat{Z}$ are computed using the EOS with the constrained value of T and the conserved variables;
• (c)  
  If, using equation (25), f(1) ⋅ f(1.5) > 0, the point is actually set to the atmosphere;
• (d)  
  The Brent's method [55] is applied to the function f defined in equation (25).

We note that we use equation (45) but not equation (46) of [60] when forcing the temperature to be constant. Therefore, we also need to update the value of $\tilde {\tau }$ after each conservative-to-primitive calculation in order to guarantee consistency between primitive and conservative variables. This is similar to what is done in other codes when using a cold EOS during the evolution. As we will show in section 4, the code is able to easily switch from a constrained to a free temperature evolution without particular problems.

Our implementation of the neutrino leakage scheme described in section 2.3 is based on the thorn ZelmaniLeak available at the stellarcollapse website [32] and firstly presented in [44]. In particular, we employ version 20161117 of such thorn. The thorn ZelmaniLeak uses all the cross-sections and heating rate described in section 2.3 and these cannot be modified by the user unless the code itself is modified. Nevertheless, the user can choose whether to activate neutrino heating or not as well as to include or not neutrino emission since the beginning of the simulation or after some time. Moreover, the user can freely set the number of radii across which the optical depth is computed.

In order to test the implementation of the tabulated EOS treatment, we here report the results of all the simulations performed without enabling the leakage scheme, starting from both S-slicing and T-slicing initial data.

The results for the evolution of the maximum of ρ and T for the S-slicing initial condition are shown in figure 2. In these models the maximum of the temperature is located at the NS centre and it shows an increase of less than 1% by the end of the simulation (likely due to shocks produced by the NS oscillations). In particular, the figure shows exact match for simulations 01, 02, 03, and 13, as expected (see tables 1 and 2). Noticeably, adopting octant symmetry in pure-hydro simulations 02 and 03, performed with the GRHydro and the Spritz codes respectively, produces the same results as adopting full-3D in simulations 01 and 13. Moreover, the magnetic field of simulation 13 is correctly handled during the evolution and does not significantly alter the hydrodynamic quantities as expected (we remind that, even if large, a magnetic field of ∼10¹⁶ G provides a magnetic energy which is still ∼2 orders of magnitude below equipartition).

Zoom In Zoom Out Reset image size

The same comparison for T-slicing initial condition is shown in figures 3 and 4. In this case the maximum of the temperature is located instead on the NS surface. Simulation 09 is the most delicate in the pure-hydro setting, since it forces the temperature T to be constant for the first ∼2 ms and then allows it to evolve (see section 3.3). When the temperature is free to evolve, an artificial shock is produced at the surface of the NS (as expected), but, after this initial transient, the maximum of ρ follows closely the results given by the simulations 07 and 08, where T is evolved since the beginning. Also the temperature, after the initial transient, tends to a constant value. In addition, figure 4 shows perfect match between simulation 09, performed in pure-hydro, and the magnetized simulation 15.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Based on the above results, we conclude that the tabulated EOS treatment is correctly handled by our implementation and we can then proceed in testing the neutrino leakage scheme.

Here we report the results of simulations involving neutrino leakage with constant-S and constant-T initial data, including the evolution of the total neutrino luminosity for each neutrino species, computed according to equation (29). We first present the results of tests performed without the heating contribution of equation (17) and then including it.

4.2.1. Tests without heating

Figure 5 shows the comparison of tests evolving S-slicing initial data with neutrino leakage, but without the contribution of neutrino absorption and heating: the maxima of ρ and T normalized to their initial values are shown in the top panels, while the bottom panels show the results for the luminosity of each neutrino species (electron neutrinos, electron antineutrinos, and the μ and τ species going from left to right) as computed in equation (29). In particular, the luminosity plots show that the scenario is clearly dominated by electron capture. Also in this case we can see that the maximum temperature, which for the S-slicing initial data is located at the NS centre, shows an increase of less than 1%. Neutrino cooling at the centre of the star is not effective due to the high density (and thus high optical depths) and therefore it does not significantly affect the temperature evolution in that region.

Zoom In Zoom Out Reset image size

A similar comparison for T-slicing initial data is shown in figure 6 (we remind that in this case the maximum of the temperature is located on the NS surface and it is strongly affected by the artificial shocks that develop there). Despite minor differences due to the different implementations in the GRHydro and Spritz codes, the results appear in good agreement.

Zoom In Zoom Out Reset image size

4.2.2. Tests including heating

We now turn to consider how the heating contribution alters the results of simulations. In figures 7 and 8 we compare the results respectively of one S-slicing and one T-slicing ID performed with and without such contribution. As already seen in figure 3, starting from cold NS initial data produces a sharp transient for the maximum of T (located at the NS surface for the T-slicing ID) in the first few time steps, where the NS internal temperature undergoes a re-adjustment (due also to the expected production of shocks at the NS surface). This transition may be an issue when considering neutrino leakage since it may produce luminosities much larger than expected. Moreover, we recall that the heating given by equation (17) is not self-consistent in terms of energy balance (see also section 3.4). Therefore, when considering the heating contribution (figure 8), we activated the leakage 1 ms later, i.e. after the initial transient.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 9 collect results for S-slicing ID and neutrino leakage including heating. In this case, without an initial temperature readjustment, the heating contribution does not need to be activated after 1 ms. We also show the maximum magnetic field evolution for the magnetized cases 13 (without leakage) and 14 (with leakage and the heating contribution) in figure 10. We found an exact match.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figure 11, we compare the cases with cold NS initial data (T-slicing) and neutrino leakage including heating. For simulation 11, which evolves the temperature since the beginning, we enable the leakage after only 1 ms. For simulations 12 and 16, evolving the temperature only after 2 ms, we enable the leakage at 3 ms. Despite the difference in the activation times of T evolution and leakage, and in the presence or absence of magnetic fields, all the results show a very good agreement in the maximum rest-mass density and the late-time electron neutrino luminosities. Finally, looking again at the maximum magnetic field evolution, figure 12 shows that also simulations 15 and 16 are perfectly matching each other.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

All the test results presented in this section are indicative of a correct implementation of the neutrino leakage scheme and that the code is ready to be used in more complex astrophysical scenarios, e.g. BNS mergers including tabulated EOS, magnetic fields, and neutrino emission and absorption (with the intrinsic limitations of the leakage scheme itself; see discussion below).

The first step in the development of a high-order scheme is the choice of the reconstruction method. Here, we consider the fifth-order WENOZ algorithm [68]. In the following, we will consider only one dimension without loss of generality: the multidimensional scheme is simply retrieved by considering the fluxes in each direction separately.

The fifth-order WENO scheme employs a 5-points stencil, S⁵, which is subdivided into three 3-points substencils, {S⁰, S¹, S²}. The polynomial approximation f_i+1/2, which is the reconstruction of the grid function f_i on the left side of the interface ¹⁶ , is built through the following convex combination of the interpolated values ${f}_{i+1/2}^{k}$, that are third degree polynomials defined on each substencil S^k , k = 0, 1, 2:

$$\begin{equation}{f}_{i+1/2}=\sum\limits _{k=0}^{2}{\omega }_{k}{f}_{i+1/2}^{k}.\end{equation} \tag{ A.1 }$$

The polynomial on each substencil is given by the quadratic interpolations

$$\begin{equation}{f}_{i+1/2}^{0}=\frac{1}{8}\left(3{f}_{i-2}-10{f}_{i-1}+15{f}_{i}\right),\end{equation} \tag{ A.2 }$$

$$\begin{equation}{f}_{i+1/2}^{1}=\frac{1}{8}\left(-{f}_{i-1}+6{f}_{i}+3{f}_{i+1}\right),\end{equation} \tag{ A.3 }$$

$$\begin{equation}{f}_{i+1/2}^{2}=\frac{1}{8}\left(3{f}_{i}+6{f}_{i+1}-{f}_{i+2}\right).\end{equation} \tag{ A.4 }$$

The weights ω_k are defined as

$$\begin{equation}{\omega }_{k}=\frac{{\alpha }_{k}}{{\sum }_{j=0}^{2}{\alpha }_{j}}.\end{equation} \tag{ A.5 }$$

For WENOZ, the unnormalized weights α_k are defined as

$$\begin{equation}{\alpha }_{k}={d}_{k}\left(1+\frac{\vert {\beta }_{0}-{\beta }_{2}\vert }{{\beta }_{k}+\varepsilon }\right),\end{equation} \tag{ A.6 }$$

with ɛ = 10⁻²⁶ (which avoids a possible division by zero), optimal weights d_k = (1/16, 10/16, 5/16), corresponding to the weights obtained for smooth fields, and smoothness indicators

$$\begin{equation}{\beta }_{0}=\frac{13}{12}{\left({f}_{i-2}-2{f}_{i-1}+{f}_{i}\right)}^{2}+\frac{1}{4}{\left({f}_{i-2}-4{f}_{i-1}+3{f}_{i}\right)}^{2},\end{equation} \tag{ A.7 }$$

$$\begin{equation}{\beta }_{1}=\frac{13}{12}{\left({f}_{i-1}-2{f}_{i}+{f}_{i+1}\right)}^{2}+\frac{1}{4}{\left({f}_{i-1}-{f}_{i+1}\right)}^{2},\end{equation} \tag{ A.8 }$$

and

$$\begin{equation}{\beta }_{2}=\frac{13}{12}{\left({f}_{i}-2{f}_{i+1}+{f}_{i+2}\right)}^{2}+\frac{1}{4}{\left(3{f}_{i}-4{f}_{i+1}+{f}_{i+2}\right)}^{2},\end{equation} \tag{ A.9 }$$

that measure the regularity of the kth polynomial approximation ${f}_{i}^{k}$ at the stencil S^k .

Note that the choices of the coefficients in (A.2)–(A.4) and of the optimal weights d_k follow the one in [69], which differ from the one in the original paper, because it has been noted that these values suit better the high order scheme in combination with the derivation operation.

The derivation operation is a high-order procedure which allows one to obtain a high order approximation from the point value quantities calculated at the intercell location.

This step has to be performed right after the computation of the fluxes via an approximate Riemann solver and it is necessary to preserve the accuracy in the calculation of spatial derivatives for schemes with order n > 2. As we did before, we will restrict the discussion to one dimension. The procedure described here follows the one outlined in the ECHO paper [69]. Using this procedure, we will provide the numerical flux function ${\hat{f}}_{i+1/2}$, given a stencil of intercell fluxes {f_i+1/2}.

The finite difference approximation of the first derivative in the point x_i can be written as

$$\begin{align}\hfill h{f}^{\prime }\left({x}_{i}\right)& \approx {\hat{f}}_{i+1/2}-{\hat{f}}_{i-1/2}\hfill \\ \hfill & =a\left({f}_{i+1/2}-{f}_{i-1/2}\right)+b\left({f}_{i+3/2}-{f}_{i-3/2}\right)+c\left({f}_{i+5/2}-{f}_{i-5/2}\right),\hfill \end{align} \tag{ A.10 }$$

where the approximation has been truncated at sixth order and h is the constant grid spacing.

If we now expand both sides of the equation in Taylor series around x_i we find

$$\begin{equation}h{f}_{i}^{\left(1\right)}=\sum\limits _{k=0}^{+\infty }{f}_{i}^{\left(k\right)}\frac{{h}^{k}}{k!{2}^{k}}\left[1-{\left(-1\right)}^{k}\right]\left[a+{3}^{k}b+{5}^{k}c\right],\end{equation} \tag{ A.11 }$$

where the exponents indicate the corresponding order of derivation, and the first derivative has been rewritten as ${f}_{i}^{\left(1\right)}\equiv {f}^{\prime }\left({x}_{i}\right)$. It is clear that all terms with even k vanish. For n = 2, where b = c = 0, we find a = 1. For n = 4, where c = 0, we have a = 9/8 and b = −1/24. Finally, for n = 6, the solution is a = 75/64, b = −25/384, c = 3/640. The next step is to write

$$\begin{equation}{\hat{f}}_{i+1/2}={d}_{0}{f}_{i+1/2}+{d}_{2}\left({f}_{i-1/2}+{f}_{i+3/2}\right)+{d}_{4}\left({f}_{i-3/2}+{f}_{i+5/2}\right),\end{equation} \tag{ A.12 }$$

and the comparison with (A.11) gives the relations d₀ = a + b + c, d₂ = b + c, d₄ = c. The numerical values of d₀, d₂, and d₄ for the different order of approximation are provided in table A1. Note that for n = 2 one gets ${\hat{f}}_{ j+1/2}={f}_{j+1/2}$ as expected.

In order to highlight the nature of this procedure as a correction for higher than second order approximation, it is convenient to rewrite equation (A.12) as

$$\begin{equation}{\hat{f}}_{i+1/2}={f}_{i+1/2}-\frac{1}{24}{{\Delta}}^{\left(2\right)}{f}_{i+1/2}+\frac{3}{640}{{\Delta}}^{\left(4\right)}{f}_{i+1/2},\end{equation} \tag{ A.13 }$$

where only the first term is used in the case n = 2, the second is added for n = 4 and the complete expression is used for n = 6. For a generic index i the second and fourth order numerical derivative are given by

$$\begin{equation}{{\Delta}}^{\left(2\right)}{f}_{i}={f}_{i-1}-2{f}_{i}+{f}_{i+1},\end{equation} \tag{ A.14 }$$

and

$$\begin{align}\hfill {{\Delta}}^{\left(4\right)}{f}_{i}& ={{\Delta}}^{\left(2\right)}{f}_{i-1}-2{{\Delta}}^{\left(2\right)}{f}_{i}+{{\Delta}}^{\left(2\right)}{f}_{i+1}\hfill \\ \hfill & ={f}_{i-2}-4{f}_{i-1}+6{f}_{i}-4{f}_{i+1}+{f}_{i+2},\hfill \end{align} \tag{ A.15 }$$

respectively.

The first test performed to check the convergence of the total procedure is the evolution of a relativistic simple wave [70, 71]. We have run this test using WENOZ as reconstruction method along with n = 2, 4, 6 correction to the HLLE Riemann solver (in the following, they will be addressed as HLLE2, HLLE4, and HLLE6, respectively).

The initial data are set up by choosing a reference state: following [35], we chose a right-propagating simple wave with ρ₀ = 1 and v₀ = 0. Assuming a polytropic EOS with Γ = 5/3 and K = 100, one can compute the sound speed in the reference frame via

$$\begin{equation}{c}_{0}^{2}=\frac{K{\Gamma}\left({\Gamma}-1\right){\rho }_{0}^{{\Gamma}}}{\left({\Gamma}-1\right){\rho }_{0}+K{\Gamma}{\rho }_{0}^{{\Gamma}}},\end{equation} \tag{ A.16 }$$

obtaining, in the specific case, c₀ ≈ 0.815. After the reference state has been defined, the velocity is perturbed with a sin-like function, so that its profile becomes (dashed line in the left panel of figure A1)

$$\begin{equation}v=a{\Theta}\left(X-\vert x\vert \right){\mathrm{sin}}^{6}\left[\frac{\pi }{2}\left(\frac{x}{X}-1\right)\right],\end{equation} \tag{ A.17 }$$

where Θ(x) is the Heaviside function, a = 0.5, and X = 0.3. Finally, the new sound speed is computed according to the Riemann invariant [71]

$$\begin{equation}{c}_{\mathrm{s}}=\sqrt{{\Gamma}-1}\frac{\frac{\sqrt{{\Gamma}-1}+{c}_{0}}{\sqrt{{\Gamma}-1}-{c}_{0}}{\left(\frac{1+v}{1-v}\right)}^{\sqrt{{\Gamma}-1}/2}-1}{\frac{\sqrt{{\Gamma}-1}+{c}_{0}}{\sqrt{{\Gamma}-1}-{c}_{0}}{\left(\frac{1+v}{1-v}\right)}^{\sqrt{{\Gamma}-1}/2}+1},\end{equation} \tag{ A.18 }$$

so that c_s = c₀ at v = 0 and ${c}_{\mathrm{s}}\to \sqrt{{\Gamma}-1}$ as v → 1. The other quantities follow from the EOS:

$$\begin{equation}\hat{\varepsilon }=\frac{{c}_{\mathrm{s}}^{2}}{{\Gamma}\left({\Gamma}-1-{c}_{\mathrm{s}}^{2}\right)},\end{equation} \tag{ A.19 }$$

$$\begin{equation}\hat{\rho }={\varepsilon }^{1/\left({\Gamma}-1\right)},\end{equation} \tag{ A.20 }$$

$$\begin{equation}\hat{p}={\varepsilon }^{{\Gamma}/\left({\Gamma}-1\right)},\end{equation} \tag{ A.21 }$$

where $\hat{\varepsilon }$, $\hat{\rho }$, and $\hat{p}$ are, respectively, the specific internal energy, the density, and the pressure normalized over the corresponding quantities in the reference state. The solutions are computed on a one-dimensional domain [−1.5, 1.5], employing RK4 integrator for HLLE2 and HLLE4, and RK65 for HLLE6 ¹⁷ , with a CFL factor of 0.125.

Zoom In Zoom Out Reset image size

During the evolution, the profile of the wave begins to steepen until a shock is formed at t ≈ 0.63 (see [70]). In order to quantify the convergence properties of the various methods, we computed the self convergence factor defined as

$$\begin{equation}p\equiv {\mathrm{log}}_{2}\left(\frac{{\Vert}f\left(4{\Delta}x\right)-f\left(2{\Delta}x\right){\Vert}}{{\Vert}f\left(2{\Delta}x\right)-f\left({\Delta}x\right){\Vert}}\right).\end{equation} \tag{ A.22 }$$

The functions f(Δx), f(2Δx), and f(4Δx) represent the numerical solutions calculated on uniform grids with corresponding grid spacing, and the norm employed is the L2-norm. In this test, the three different resolutions are Δx = 0.0075, 0.003 75, 0.001 875, corresponding to 400, 800, 1600 points.

As it can be seen from the right panel of figure A1, the nominal convergence order is reached until the appearance of the shock. For both WENOZ + HLLE2 and WENOZ + HLLE4, the convergence is dominated by the order of the derivation operation; otherwise, for WENOZ + HLLE6, the convergence is dominated by the order of the reconstruction method and this is why we cannot get an order of convergence higher than fifth. As expected, the convergence order goes down for all methods when the shock is formed.

A second test that has been performed is the evolution of a non-magnetized TOV star; the setup is the same used in the first paper of Spritz [27]. In particular, the initial configuration is generated using a polytropic EOS with Γ = 2.0 and K = 100, and initial rest-mass density ρ = 1.28 × 10⁻³. The evolution of the system is then carried out adopting an ideal fluid EOS with the same value of Γ. The physical domain is [−20, 20] for x-, y-, and z-coordinates, with low, medium, and high resolution having 32³, 64³, and 128³ cells, respectively. All the tests lasted for 5 ms using the WENOZ reconstruction method and the three approximation for the Riemann solver (HLLE2, HLLE4, and HLLE6). In the cases of HLLE2 and HLLE4, RK4 method is employed for time stepping, while RK65 is used in HLLE6 case, with a CFL factor of 0.25.

In the continuum limit, the evolution of this kind of system is trivial; however, the discretization of the problem brings errors (due to the discretization itself) that cause radial oscillations, which are observable, for example, in the central rest-mass density (see figure A2). The amplitude of these oscillations becomes smaller as the number of points increases. In the right panels of figure A2 it is possible to note that the density has a peculiar behaviour for low and medium resolution at late times. This fact can be traced back to the choice of the ideal fluids EOS in the evolution of the system: it is known that truncation errors with this EOS are very large, because significant unphysical shock-heating is observed at low densities [35].

Zoom In Zoom Out Reset image size

In order to verify the convergence of the high-order methods, we compute the self-convergence factor (based on deviations of central rest-mass density with respect to the initial value), which oscillates around the value p = 3 for both HLLE4 and HLLE6. Such order of convergence is maintained until the aforementioned truncation errors become significant, i.e. until ∼4 ms.

In the end, figure A3 reports the power spectrum of the evolution of the rest-mass densities of the different runs. The power spectrum is computed via a fast Fourier transform (FFT) in order to extract the amplitudes and the frequencies of the oscillations, and then the amplitudes are normalized to the maximum for each simulation. Figure A3 also shows the peaks' frequencies of the oscillations taken from the literature [72], that were obtained with independent codes. All the simulations show a good agreement with each other and the independent results. In particular, it is worth noting that the high-order reconstruction coupled to high-order Riemann solvers (black-dotted and green-dashed curves in the figure) is evidently capable of better resolving the overtones (i.e. the higher frequency peaks in the spectrum) with respect to the lower-order methods (red-solid and blue-dash-dotted curves in the figure).

Zoom In Zoom Out Reset image size

Finally, we performed a test evolving a magnetized TOV star using a tabulated EOS (LS220), with an S-slicing initial condition, and employing WENOZ as reconstruction method and the 4th order approximation for the HLLE Riemann solver (this case has been called WENOZ + HLLE4). This case is then compared with the same case evolved using PPM reconstruction method and the 2nd order approximation to HLLE, denoted with PPM + HLLE2.

The upper panels of figure A4 show the evolution of the central rest-mass density ρ_c and of the maximum of the temperature T_max, both normalized over their initial values, respectively ρ_c,0 and T_max,0. The results obtained with the use of the high-order scheme, in particular the setup WENOZ + HLLE4, are more precise than the ones obtained with the older version of the Spritz code; using high-order methods helps reducing the oscillations around the real value. Moreover, enabling WENOZ and the fourth-order correction to HLLE softens the slight increasing behaviour of T_max, as shown in the upper right panel of figure A4.

Zoom In Zoom Out Reset image size

The gain in accuracy is particularly evident in the plot for the power spectrum of the evolution of the central rest-mass density, shown in the lower panel of figure A4. Each power spectrum is computed, as before, via the FFT and the amplitudes are normalized over their maximum for each simulation. It can be easily seen that, while the lower-order version of Spritz shows a noticeable peak only for the fundamental frequency, the high-order upgrade can resolve very well also the first overtone, which results to be more prominent than the one of the PPM + HLLE2 case.