Time integration for neutrino radiation transport using minimally implicit Runge-Kutta methods

TL;DR

This paper introduces a new minimally implicit Runge-Kutta method for stable, efficient neutrino-radiation hydrodynamics simulations, suitable for stiff equations in astrophysics.

Contribution

The paper presents a novel minimally implicit Runge-Kutta approach that allows analytical inversion, reducing computational cost while handling stiffness in neutrino-matter interaction equations.

Findings

Method effectively handles stiff neutrino-matter reactions.

Application to supernova simulations demonstrates stability and efficiency.

Results show accurate modeling of neutrino transport in astrophysical environments.

Abstract

The evolution of many astrophysical systems is dominated by the interaction between matter and radiation such as photons or neutrinos. The dynamics can be described by the evolution equations of radiation hydrodynamics in which reactions between matter particles and radiation quanta couples the hydrodynamic equations to those of radiative transfer (see Munier & Weaver (1986a) and Munier & Weaver (1986b)). The numerical treatment has to account for their potential stiffness (e.g., in optically thick environments). In this article, we will present a new method to numerically integrate these equations in a stable way by using minimally implicit Runge-Kutta methods. With these methods, the inversion of the implicit operator can be done analytically, so the computational cost is equivalent to that of an explicit method. We strongly take into account the physical behavior of the evolved variables in the limit of the stiff regime in the derivation of the methods. We will show the results of applying these methods to the reactions between neutrinos and matter in some tests and also in realistic core-collapse supernovae simulations.

Tables

Table 1. Table 1: List of first order MIRK simulations. The first three columns give the name of the simulation, and the values of the parameters a 𝑎 a and b 𝑏 b . In the last column, the symbols ✓ ✓ \checkmark , △ △ \bigtriangleup , and × \times indicate simulations that ran stably and with correct results, simulations that ran stably into the post-bounce phase, but gave wrong results, and simulations that turned unstable when the core reached optically thick conditions, respectively.

model	$a$	$b$	result
M1-1	$0$	$0$	$✓$
M1-2	$1/2$	$1/2$	$\times$
M1-3	$0$	$1/2$	$\times$
M1-4	$1/2$	$0$	$△$

Table 2. Table 2: List of second order MIRK simulations. The first five columns give the name of the simulation, the values of the parameters a 𝑎 a , a ′ superscript 𝑎 ′ a^{\prime} , b 𝑏 b , and b ′ superscript 𝑏 ′ b^{\prime} . The symbols in the last column have the same meaning as in Tab.  1 .

model	$a$	$a^{\prime }$	$b$	$b^{\prime }$	result
M2-11	$+1/2$	$-1/4$	$+1/2$	$-1/4$	$\times$
M2-12	$+1/2$	$-1/4$	$+1/2$	$+1/4$	$\times$
M2-13	$+1/2$	$-1/4$	$-1/2$	$-3/4$	$\times$
M2-14	$+1/2$	$-1/4$	$-1/2$	$-9/4$	$△$
M2-21	$+1/2$	$+1/4$	$+1/2$	$-1/4$	$\times$
M2-22	$+1/2$	$+1/4$	$+1/2$	$+1/4$	$\times$
M2-23	$+1/2$	$+1/4$	$-1/2$	$-3/4$	$\times$
M2-24	$+1/2$	$+1/4$	$-1/2$	$-9/4$	$△$
M2-31	$-1/2$	$-3/4$	$+1/2$	$-1/4$	$\times$
M2-32	$-1/2$	$-3/4$	$+1/2$	$+1/4$	$\times$
M2-33	$-1/2$	$-3/4$	$-1/2$	$-3/4$	$\times$
M2-34	$-1/2$	$-3/4$	$-1/2$	$-9/4$	$△$
M2-41	$-1/2$	$-9/4$	$+1/2$	$-1/4$	$\times$
M2-42	$-1/2$	$-9/4$	$+1/2$	$+1/4$	$\times$
M2-43	$-1/2$	$-9/4$	$-1/2$	$-3/4$	$\times$
M2-44	$-1/2$	$-9/4$	$-1/2$	$-9/4$	$✓$
M2-44-1	$-1/4$	$-25/8$	$-1/4$	$-25/8$	$✓$
M2-44-2	$-1/16$	$-289/32$	$-1/16$	$-289/32$	$✓$
M2-44-3	$-1/2$	$-9/4$	$-1/4$	$-25/8$	$✓$
M2-44-4	$-1/4$	$-25/8$	$-1/2$	$-9/4$	$✓$
M2-51	$-1/2$	$-9/4$	$3/4$	$1/24$	$\times$
M2-52	$-1/2$	$-9/4$	$3/4$	$-1/8$	$\times$
M2-53	$3/4$	$1/24$	$3/4$	$1/24$	$\times$
M2-54	$3/4$	$-1/8$	$-1/2$	$-9/4$	$△$
M2-55	$3/4$	$1/24$	$-1/2$	$-9/4$	$✓$

Table 3. Table 3: Values of x = ( a , b ) 𝑥 𝑎 𝑏 x=(a,b) (first column) and the corresponding values of the functions used to compute x ′ = ( a ′ , b ′ ) superscript 𝑥 ′ superscript 𝑎 ′ superscript 𝑏 ′ x^{\prime}=(a^{\prime},b^{\prime}) (second and third columns).

$x$	$\frac{x-1}{2}$	$\frac{{\left(1-x\right)}^2}{2x}$
$-1/2$	$-3/4$	$-9/4$
$-1/4$	$-5/8$	$-25/8$
$-1/16$	$-17/32$	$-289/32$
$1/2$	$-1/4$	$1/4$
$3/4$	$-1/8$	$1/24$

Equations

$${\partial }_{t}E$$

$${\partial }_t{F}^i$$

$$C^{(0)}$$

$$C^{(1),i}$$

$$(u^{n+1}-u^n)/\Delta t$$

$$(w^{n+1}-w^n)/\Delta t$$

$$E^{n+1}$$

$$(F^i{)}^{n+1}$$

$${\left(\begin{array}{c}E\\ F^i\end{array}\right)}^{n+1}={\left(\begin{array}{c}E\\ F^i\end{array}\right)}^n+\left(\begin{array}{cc}\frac{\Delta t}{1+\Delta t{\kappa }^n(1-a)}& 0\\ 0& {\delta }^{ij}\frac{\Delta t}{1+\Delta t{\kappa }^{\prime n}(1-b)}\end{array}\right){\left(\begin{array}{c}c{\kappa }_a(E_{\text{eq}}-E)+S_E\\ -c{\kappa }_{\text{tra}}{F}^i+S_F^i\end{array}\right)}^n,$$

$${\left(\begin{array}{c}E\\ F^i\end{array}\right)}^{n+1}=\left(\begin{array}{cc}\frac{-a}{1-a}& 0\\ 0& {\delta }^{ij}\frac{-b}{1-b}\end{array}\right){\left(\begin{array}{c}E\\ F^i\end{array}\right)}^n+\left(\begin{array}{c}\frac{E_{\text{eq}}^n}{1-a}\\ 0\end{array}\right).$$

$$a<1/2,b<1/2,$$

$$E^{n+1}=E_{\text{eq}}^{n+1}+\mathcal{O}(\Delta t),$$

$$(F^i{)}^{n+1}=0+\mathcal{O}(\Delta t).$$

$$E^{n+1}$$

$$(F^i{)}^{n+1}$$

$$E^{(1)}$$

$$(F^i{)}^{(1)}$$

$$E^{n+1}$$

$$(F^i{)}^{n+1}$$

$$E^{(1)}$$

$$(F^i{)}^{(1)}$$

$$\begin{array}{rl}{E}^{n+1}=& {\left[1+\Delta t{\kappa }^{(1)}\left(\frac{a}{2}-a^{\prime }\right)\right]}^{-1}\left\{\left[\frac{1}{2}-\Delta t{\kappa }^{(1)}\left(\frac{1-a}{2}\right)\right]E^n+\left[\frac{1}{2}-\Delta t{\kappa }^{(1)}{a}^{\prime }\right]E^{(1)}+\frac{\Delta t{S}_E^{(1)}}{2}+\frac{\Delta t{\kappa }^{(1)}{E}_{\text{eq}}^{(1)}}{2}\right\}\\ =& {\left[1+\Delta t{\kappa }^{(1)}\left(\frac{a}{2}-a^{\prime }\right)\right]}^{-1}\{\left[\frac{1}{2}-\Delta t{\kappa }^{(1)}\left(\frac{1-a}{2}\right)\right]E^n+\left[\frac{1}{2}+\Delta t{\kappa }^{(1)}\left(\frac{1}{2}-a^{\prime }\right)\right]E^{(1)}\\ & +\frac{\Delta t}{2}\left[S_E^{(1)}+{\kappa }^{(1)}(E_{\text{eq}}^{(1)}-E^{(1)})\right]\},\end{array}$$

$$\begin{array}{rl}(F^i{)}^{n+1}=& {\left[1+\Delta t{\kappa }^{\prime (1)}\left(\frac{b}{2}-b^{\prime }\right)\right]}^{-1}\left\{\left[\frac{1}{2}-\Delta t{\kappa }^{\prime (1)}\left(\frac{1-b}{2}\right)\right](F^i{)}^n+\left[\frac{1}{2}-\Delta t{\kappa }^{\prime (1)}{b}^{\prime }\right](F^i{)}^{(1)}+\frac{\Delta t(S_F^i{)}^{(1)}}{2}\right\}\\ =& {\left[1+\Delta t{\kappa }^{\prime (1)}\left(\frac{b}{2}-b^{\prime }\right)\right]}^{-1}\{\left[\frac{1}{2}-\Delta t{\kappa }^{\prime (1)}\left(\frac{1-b}{2}\right)\right](F^i{)}^n+\left[\frac{1}{2}+\Delta t{\kappa }^{\prime (1)}\left(\frac{1}{2}-b^{\prime }\right)\right](F^i{)}^{(1)}\\ & +\frac{\Delta t}{2}\left[(S_F^i{)}^{(1)}-{\kappa }^{\prime (1)}(F^i{)}^{(1)}\right]\}.\end{array}$$

$$\frac{a}{2}-a^{\prime }>0,1-a>0;\frac{b}{2}-b^{\prime }>0,1-b>0.$$

$$E^{n+1}$$

$$(F^i{)}^{n+1}$$

$${\lambda }_E=\frac{a^{\prime }a-\left(\frac{1-a}{2}\right)(1-a)}{\left(\frac{a}{2}-a^{\prime }\right)(1-a)},$$

$${\lambda }_F=\frac{b^{\prime }b-\left(\frac{1-b}{2}\right)(1-b)}{\left(\frac{b}{2}-b^{\prime }\right)(1-b)}.$$

$$E^n$$

$$(F^i{)}^n$$

$$E^{n+1}$$

$$(F^i{)}^{n+1}$$

$$a^{\prime }=\frac{a-1}{2},b^{\prime }=\frac{b-1}{2}.$$

$$a^{\prime }=\frac{(1-a{)}^2}{2a},a\in (-\infty ,0)\cup (1/2,1),$$

$$b^{\prime }=\frac{(1-b{)}^2}{2b},b\in (-\infty ,0)\cup (1/2,1),$$

$$E^{n+1}$$

Full text

Radiative transfer — Methods: numerical — Supernovae: general

1 Introduction

Radiation plays a crucial role in many fields of astrophysics. Besides representing the most important channel for observations of astronomical objects, its interaction with matter shapes their structure and drives their dynamics in quiescent as well as highly dynamic phases. Various types of reactions such as emission, absorption, and scattering exchange energy and momentum between radiation quanta and matter particles. The corresponding cross sections and interaction rates can depend in very complicated ways on the thermodynamic properties of the gas and the radiation spectra. The resulting mean free paths of radiation can be very small compared to characteristic structural length scales (e.g., density scale heights) in optically thick regions with, for example, high gas densities, and at the same time exceeding the dimensions of the system in transparent regions such as atmospheric layers. While the former limit is commonly reached for photons in the interior of ordinary stars, densities close to that of nuclear matter are required to turn matter optically thick to neutrino radiation. Thus, neutrinos are extremely important in core-collapse supernovae (CCSNe) and neutron stars, where this condition is met. Their importance is such that the explosion mechanism of CCSNe cannot be understand without a detailed account of the generation and transport of neutrinos.

Much of the complexity of theoretically modelling the aforementioned systems comes from the equation underlying radiative transfer. At a basic level, the Boltzmann equation describes the evolution of the distribution function, $f$, of radiation quanta a phase-space of seven dimensions, viz. time, position, and neutrino or photon momentum. It takes into account transport terms involving divergence operators in position or momentum space as well as collision terms representing the interaction of one quantum with others or with matter that take the form of integrals over momentum space. As a rigorous treatment of the Boltzmann equation is feasible only in special cases, numerical methods used to model, e.g., CCSNe rely on approximations. In a very common approach, a momentum-space integration of the distribution function multiplied by the tensorial product of $n=0,1,...$ unit vectors yields a series of moments, now only functions of space and time. The first ones have a direct physical interpretation: $n=0,1,$ and $2$ correspond to the radiation energy and momentum densities and the radiation pressure, respectively. The resulting infinite series of evolution equations takes the form of conservative form, in which the moment of order $n+1$ appears as a (spatial) flux of the moment of degree $n$ and source terms accounting for the reactions follow from moments of the collision integrals. We note that, when coupling radiative transfer to the dynamics of the gas to form the system of radiation hydrodynamics, the latter appear with an opposite sign in the fluid equations. Truncating the series at a finite degree $n$ and closing the system with a local algebraic relation for the higher moment(s) defines the family of $M_{n}$ methods.

$M_{0}$ or (flux-limited) diffusion and $M_{1}$ or algebraic Eddington tensor methods offer a good compromise between accuracy and numerical costs, and are thus widely used in relativistic astrophysics. $M_{1}$ methods are very good at modelling radiation in the optically thick and the transparent regimes and also work well in the intermediate, semi-transparent regime. Nonetheless, a few difficulties remain. A particularly important one pertains to the time integration in the optically thick regime, in which the typical time scales of interactions between radiation and matter (the inverse of the reaction rates) can be many orders of magnitude smaller than the time scales associated to the radiation propagation or the dynamical time scales: the equations are stiff. A similar behaviour can be found in other conservation laws in many physical problems such as Resistive Relativistic Magnetohydrodinamics (Cordero-Carrión et al., 2023), general relativistic force-free electrodynamics (Mahlmann et al., 2021), rarefied gases problems (Koellermeier & Samaey, 2022) or shallow water equations (Koellermeier & Pimentel-García, 2022).

Designing methods for stiff equations requires specific considerations. Explicit time integration is only stable if the numerical time step is reduced to the characteristic time scales of the fastest evolving term, which in this case would be the radiation-matter interaction ones. Implicit methods, on the other hand, allow for a stable evolution even when using the–much larger–time steps set by, e.g., radiation propagation or hydrodynamics. They, however, can be very complicated to implement due to the inversion of the operators involved, in particular for parallel execution, and suffer from low computational efficiency. As a compromise, Implicit-Explicit (IMEX) Runge-Kutta methods (Pareschi & Russo, 2005) combine an implicit integration of only the stiff terms with an explicit integration of the rest of the equations. This strategy has been used very recently by Izquierdo et al. (2022). Semi-implicit numerical schemes (Foucart et al., 2016) have also been used very recently in neutron star mergers (see, for example, Radice et al., 2022). We present in this manuscript a new numerical scheme to solve the $M_{1}$ neutrino-hydrodynamics equations to first and second order in time that preserves stability properties of implicit methods and, at the same time, has a computing speed similar to that of an explicit method. The scheme proposed is based on implicit evaluations of only evolved variables and requires a slight modification from an explicit method. We implement the new solver in the neutrino-hydrodynamics code of Just et al. (2015) as an higher order alternative to the implicit-explicit scheme that is already implemented.

Other very different numerical approaches can be considered. For example, Monte Carlo methods can be used to include neutrino transport in the context of core-collapse supernovae and this strategy has been implemented in the general relativistic code SpEC (Foucart, 2018; Foucart et al., 2021).

This article is structured as follows. In Sect. 2 we describe the main features of the equations governing the dynamics of the neutrino transport within the $M_{1}$ radiative transfer scheme. In Sect. 3 we explain in detail how the MIRK methods, up to second order, can be used to numerically integrate these equations, including the stability analysis in the stiff limit. In Sect. 4 we show some numerical successful simulations of stellar core-collapse in spherical symmetry using our proposed MIRK scheme, first and second order, and we compare the results with some other numerical approaches previously used. Finally, in Sect. 5 we present our conclusions.

2 Equations for neutrino transport in M1

The basic variables of $M_{1}$ radiative transfer, the energy density, $E(t,\vec{r},\omega)$, and the momentum density, $\vec{F}(t,\vec{r},\omega)$, of the radiation field, are functions of time $t$, position $\vec{r}$ and particle energy $\omega$ (Munier & Weaver, 1986a, b). Owing to their conservative character, the corresponding evolution equations take the form of balance laws including the spatial transport and the redistribution across particle energies by differential operators. Exchange of energy and momentum with matter enter the equations via source terms that typically depend only on the local state of radiation field and the matter, but not on their derivatives. Since we will deal with the latter terms, we write the $M_{1}$ system in the following short-hand notation (see also Just et al., 2015):

[TABLE]

where the terms with spatial or energy derivatives, $S_{E}$ and $S_{F}^{i}$, are split off from the interaction source terms, $C^{(0)}$ and $C^{(1),i}$. The form of the interaction source terms depends on the choice of interactions and possible approximations used to describe them. We specialise to the important case of thermal emission and absorption and isotropic scattering. Then, the rates of energy and momentum exchange are proportional to the absorption and transport opacities, $\kappa_{\mathrm{a}}$ and $\kappa_{\mathrm{tra}}$, respectively:

[TABLE]

Eq. (3) describes how matter emits radiation thermally with an equilibrium distribution $E_{\text{eq}}$ (e.g., Maxwell-Boltzmann for photons, Fermi-Dirac for neutrinos), and how it absorbs the local radiation energy. Eq. (4) accounts for the transfer of momentum to the gas by means of absorption and scattering reactions. We note that the same terms appear with the opposite sign (and integrated over particle energy) as sources in the hydrodynamic equations for the gas.

We do not delve deeper into the detailed, potentially very complicated, dependence of the opacities on $E$ and $\vec{F}$ as well as on the composition and thermodynamic state of the gas because our method is valid for general opacity laws. Our main focus lies on the stiff, optically thick limit, in which the opacities are very high, $\kappa_{\mathrm{a,tra}}\gg 1$, and the interaction terms dominate over $S_{E}$ and $S_{F}^{i}$ in Eqs. (1) and (2). Under these conditions, numerical difficulties arise due to the need to simultaneously follow all the terms with characteristic time scales that can differ by many orders of magnitude.

The physically correct stiff limit consists of $E$ approaching the equilibrium energy density $E_{eq}$. Furthermore, Eq. (4) indicates that high opacities will reduce $\vec{F}$ to zero. However, the precise manner in which $\vec{F}$ vanishes matters a lot for getting the correct solution. In a non-uniform radiation field, $\vec{F}$ has to approach the diffusion limit satisfying $\displaystyle\vec{F}\to\vec{F}_{\mathrm{diff}}=\frac{1}{3\,c\,\kappa_{\mathrm{tra}}}\vec{\nabla}E$. While some $M_{1}$ methods (Jin & Levermore, 1996; Pons et al., 2000; Jin et al., 2000; Audit et al., 2002) deal with this requirement by explicitly enforcing the diffusion flux for high optical thickness, others (Just et al., 2015) found that an appropriate treatment of the flux terms in $S_{F}$ is sufficient to reproduce the correct limit. In practise, approaches such as the one of Just et al. (2015) allow us to offload the issue of the correct diffusion limit to the solution of $S_{F}$. As long as our method for $C^{(1),i}$ ensures that $\vec{F}$ vanishes in the optically thick limit in the absence of $S_{F}$, the coupled solution of $S_{F}$ and $C^{(1),i}$ will behave correctly.

We focus now on the numerical schemes that can be used to solve Eqs. (1) and (2).

The time-integration strategy for the transport equations, (1) and (2), is usually chosen based on a trade-off between stability, accuracy, and numerical costs. These goals are somewhat at odds with each other: the most stable schemes, implicit time integrators, and the most accurate ones, high-order methods, are also the most expensive ones; furthermore, high-order implicit methods tend to be particularly complex. The difficulties are exacerbated when applying the integrators to terms involving spatial as well as temporal derivatives. For this reason, an operator-splitting approach is common in which the transport terms, $S_{E,F}$, the interaction terms, $C^{0,1}$, and, in the case of coupled radiation hydrodynamics, the flux and source terms of the hydrodynamics equations not connected to neutrino interactions, are treated separately using suitable methods. In the applications we are mostly interested in, CCSNe, we follow the evolution of the system on the hydrodynamical time scales, which leads us to select an explicit time integrator for the latter group of terms. Furthermore, the maximum hydrodynamic flow and sound speeds are similar to the characteristic velocities of the neutrino transport terms, which allows us to use an explicit time integrator for them with roughly the same stability constraint on the time step. On the other hand, their stiffness makes an implicit time integration scheme the only feasible option for the interaction terms.

The IMEX strategy described above is commonly employed in neutrino-hydrodynamics codes in high-energy astrophysics. Among the proposed methods, we follow the one implemented by Just et al. (2015), whose discretised schematics we briefly summarise in the following. We denote the conserved variables of hydrodynamics (the densities of mass, momentum, energy) and of the neutrino radiation ($E,\vec{F}$), collectively as $u$ and $w$, respectively, and use superscripts n,n+1 to indicate the states at discrete time steps $t^{n}$ and $t^{n+1}=t^{n}+\Delta t$. Then our prescription to update the variables to the next time step is given by

[TABLE]

where we the symbols $\mathcal{L}_{\mathrm{hydro}}$ and $\mathcal{L}_{\mathrm{tr}}$ stand for the discretised operators including the fluxes and sources of hydrodynamics and the fluxes of neutrino transport, respectively. Without entering into further details, we note that they are evaluated explicitly with data of the previous time step, $t^{n}$. The neutrino-matter interactions, represented by the operator $\mathcal{L}_{\mathrm{int}}$, i.e., the discretised version of Eqs. (3) and (4), depends on both $w$ and $u$. Its dependence on the hydrodynamic variables is a result of both the opacities and the equilibrium energy density, $E_{\mathrm{eq}}$, being functions of the thermodynamic state of the gas. We note the its counterpart in the hydrodynamic equations, $\bar{\mathcal{L}}_{\mathrm{int}}$, can be computed once $\mathcal{L}_{\mathrm{int}}$, and thus presents no further complication.

A fully implicit treatment of $\mathcal{L}_{\mathrm{int}}$ in (6) would entail evaluating all the variables it depends on at the new time step, i.e., setting $\mathcal{L}_{\mathrm{int}}(u^{n+1},w^{n+1})$. The intricate dependence of $\kappa_{\mathrm{a,tra}}$ and $E_{\text{eq}}$ on $u$ makes this task computationally costly, which burden the numerical solution. This step would require at multiple times the recovery of the primitive (thermodynamic) variables, in particular the temperature, from $u$, i.e., the inversion of non-linear relations.

We propose a Minimally Implicit (MIRK) Runge-Kutta method that minimize computation cost of the process of recovery of variables. Our alternative approach differs in that we evaluate only the conserved neutrino variables at the new time step, $w^{n+1}$, but treat the hydrodynamic variables and the variables derived from them, opacities and equilibrium energy density, explicitly by using $u^{n}$. This simple change permits preserving the stability properties and simultaneously reducing the computational cost to that of an explicit method, as there is no need to apply the recovery multiple times. In the following we explain the method in detail. Such an approach has been implemented by Just et al. (2015) without exploring the mathematical framework presented in the next section. Here we go beyond their method, and this mathematical framework also allows for a higher order extension.

3 Numerical methods

This section present the equations of a general MIRK method at first and second order. The general expressions contain undetermined coefficients that we will choose adequately in order to guarantee a correct behaviour in the stiff limit regime.

3.1 First order method

The equations of a first order MIRK method for Eqs. (1) and (2) take the form

[TABLE]

where $a,b$ are arbitrary real coefficients that we will select later according to stability criteria. From previous equations, the explicit expressions for $E^{n+1}$ and $(F^{i})^{n+1}$ can be derived easily; they can be cast in matrix form as: {widetext}

[TABLE]

where $\kappa:=c\,\kappa_{a}$ and $\kappa^{\prime}:=c\,\kappa_{\text{tra}}$. The conditions $a,b<1$ must be satisfied to force non-zero (and positive) denominators always. Notice that the equations in this form resemble a pure explicit method with effective time steps $\displaystyle\Delta t_{E}=\frac{\Delta t}{1+\Delta t\,\kappa^{n}(1-a)}$ and $\displaystyle\Delta t_{F}=\frac{\Delta t}{1+\Delta t\,\kappa^{\prime n}(1-b)}$ for the $E$ and $F^{i}$ evolution equations, respectively. The previous matrix expression has been easily and analytically derived thanks to the fully explicit evaluation of the non conserved variables (e.g., all the variables different from $E$ and $F^{i}$). Due to this reason, one would expect to have a computational cost similar to that of applying a fully explicit method. We now analyze the behaviour in the stiff limit regime.

Mathematically speaking, the stiff limit refers to $\kappa_{a},\kappa_{\text{tra}}\rightarrow\infty$. In that limit, Eq. (9) reads

[TABLE]

Conditions

[TABLE]

must be fulfilled for the spectral radius of the updated matrix to be strictly bounded by 1, thus having a stable numerical method. This is a more restrictive condition in comparison with previous conditions $a,b<1$ (needed to avoid zero values in the denominators). In order to guarantee a correct behaviour of the numerical solution at the stiff limit, and assuming a well-behaved and smooth data for the previous time step, $E^{n}=E_{\text{eq}}^{n}+\mathcal{O}(\Delta t)$ and $(F^{i})^{n}=0+\mathcal{O}(\Delta t)$, we have that, for all $a,b$,

[TABLE]

[TABLE]

So, independently on the values of the coefficients $a,b$ we get well-behaved and smooth data in the next time step at first order. This give us, in principle, full freedom for choosing $a$ and $b$. However, the behaviour of the evolved variables are far from been smooth in supernovae simulations. Therefore, we should guarantee their correct behaviour at the stiff limit even when we are dealing with non-smooth data, and regardless the possible presence of numerical errors in the previous time steps. The choice $b=0$ guarantees the correct behaviour for $F$ at the stiff limit, i.e., $(F^{i})^{n+1}=0$. It remains to choose a value for $a$. By analogy with $b$, and taking into account the particular case $E_{\text{eq}}=0$, we will simply consider $a=0$. This means that the behaviour of $E$ at the stiff limit is not controlled by previous values of this quantity, but only by evaluations of $E_{\text{eq}}=E_{\text{eq}}(u)$, which only depends on the hydrodynamic variables $u$. With this choice, it is satisfied that $E^{n+1}=E_{\text{eq}}^{n}=E_{\text{eq}}^{n+1}+\mathcal{O}(\Delta t)$. Finally, in the case $a=b=0$ the method reads:

[TABLE]

3.2 Second order method

Hereafter we follow the same strategy as in the first order case. Two stages are needed for the second order method. We denote the intermediate step by a $(1)$ superindex and the final step by $n+1$. In general, we have four coefficients, $a,a^{\prime},b,b^{\prime}$, to be determined based on stability arguments. The first stage reads

[TABLE]

and the second stage can be written as {widetext}

[TABLE]

Isolating $E^{(1)}$ and $(F^{i})^{(1)}$, we get similar expressions to those of first order, just substituting the superindex $n$ by $(1)$:

[TABLE]

Then, $E^{n+1}$ and $F^{i,n+1}$ can be expressed explicitly in terms of previous evaluations of these quantities as:

[TABLE]

[TABLE]

The following conditions are necessary for forcing non-zero (positive) denominators always:

[TABLE]

We finally determine the coefficients of the method taking into account previous conditions and the behaviour of the numerical solution in the stiff limit.

The stiff limit refers to $\kappa^{n},\kappa^{(1)},\kappa^{\prime n},\kappa^{\prime(1)}\to\infty$. In that limit, Eqs. (20) and (21) read

[TABLE]

where

[TABLE]

$|\lambda_{E}|\leq 1$ and $|\lambda_{F}|\leq 1$ must be satisfied to guarantee stability of the numerical method.

If we assume well-behaved and smooth data at second order in time at $t=t^{n}$,

[TABLE]

from Eqs. (23) and (24) we get

[TABLE]

where we have used a Taylor expansion over $E^{(1)}_{\text{eq}}$, $X$ represents the vector of all the evolved variables (i.e. $X=(u,w)$), $\nabla_{X}=(\partial_{X_{1}},\partial_{X_{2}},...)$ and $\tilde{S}_{X}$ is the source term in the evolution equation of the form $\partial_{t}X=\tilde{S}_{X}$. So in order to satisfy Eqs. (27) and (28) in the next time step $t^{n+1}$, we need to impose $a^{\prime}=\frac{a-1}{2}$. We could choose $b^{\prime}$ in resemblance with $a^{\prime}$,

[TABLE]

This choice for $a^{\prime}$ preserves the second order behaviour of our numerical solution in the next time step and at the stiff limit when smooth data are involved, $a/2-a^{\prime}>0$ is trivially satisfied and $\lambda_{E}=-1$; note that having $\lambda_{E}=-1$ means consider a value for $\lambda_{E}$ at the border of the stability region, which seems to be not very convenient. In addition, as commented for the first order method, we must take into account that the evolved variables have a non-smooth behaviour in supernovae simulations.

Our proposal is then to consider $b\neq 0$ and $b^{\prime}=(b-1)^{2}/(2b)$ to guarantee $\lambda_{F}=0$ (as we did choosing $b=0$ for the first order method), preserving the correct behaviour of our numerical solution regardless the possible presence of numerical errors or non-smooth data. By analogy with $b^{\prime}$, and taking into account the particular case $E_{\text{eq}}=0$, we will consider $a\neq 0$ and $a^{\prime}=(a-1)^{2}/(2a)$, or $\lambda_{E}=0$ equivalently. For these choices, the conditions (22) result in

[TABLE]

[TABLE]

and Eq. (29) stays

[TABLE]

Second order for smooth data at stiff limit would be guaranteed if $a=1$, but this is incompatible with conditions (32). We only get first order in time at the stiff limit for $E$, a price to pay when non-smooth data is considered. Within these constraints, we still have some freedom for the election of the coefficients $a,b$.

We can write our numerical method in such a way it resembles a pure explicit scheme of the form

[TABLE]

For the choices (31), Eqs. (18) and (19) keep the same form, while Eqs. (20) and (21) can be written as: {widetext}

[TABLE]

where $a,b$ have not been chosen yet. For $\kappa,\kappa^{\prime}\to 0$, we recover the structure of the previous second order pure explicit method.

Finally, with conditions (32) and (33), Eqs. (18) and (19) keep the same form, while Eqs. (20) and (21) can be written as: {widetext}

[TABLE]

with $a,b$ still to be chosen. For $\kappa,\kappa^{\prime}\to 0$, we recover the structure of the previous second order pure explicit method.

4 Numerical simulations and results

4.1 Input physics

To assess its properties, we apply our method to the radiation-hydrodynamics of stellar core-collapse in spherical symmetry. This setup tests the scheme in a highly dynamic system including both the optically thin and the optically thick, stiff regimes of neutrino-matter interactions. As such, it represents a demanding problem for numerical codes. While the neglect of non-spherical flows limits the degree of realism, it makes the problem more standardised and controllable. Therefore, our tests follow in the footsteps of many previous studies of new schemes that used similar setups (e.g., Rampp & Janka, 2002; Liebendörfer et al., 2004; Sekiguchi, 2010; O’Connor & Ott, 2010; Müller et al., 2010; O’Connor, 2015; Just et al., 2015; Kuroda et al., 2016; Perego et al., 2016; O’Connor et al., 2018; Just et al., 2018; Laiu et al., 2021).

All simulations presented in the remainder of this section use the neutrino-(magneto-)hydrodynamics code Alcar (Just et al., 2015) and, except where explicitly stated, the same input physics, initial conditions, and, except for the time integration, numerical methods and parameters. We solve the equations of special relativistic hydrodynamics including a balance law for the electron fraction of the gas, $Y_{e}$. We account for the self-gravity of the star using a pseudo-relativistic gravitational potential (potential A of Marek et al., 2006). The spectral $M_{1}$ transport modules evolve the radiation energy and momentum density in a reference frame comoving with the fluid. The coupling between neutrino particle energies via velocity and gravitational terms, e.g., Doppler or gravitational red-/blue-shifts, are included up to first order in $v/c$. We describe the thermodynamic properties of the gas using the nuclear equation of state (EoS) SFHo (Steiner et al., 2013). Strictly speaking, an EoS of this type, assuming that the composition of the gas is given by nuclear statistical equilibrium, is not valid for low temperatures and densities. Nonetheless, we simplify our setup by not including a transition to a sub-nuclear EoS below a threshold density. This choice has no implication for the tests at hand because the neutrino-matter interaction rates are very small at the densities where the transition between EoS regimes would take place.

We employ the spectral $M_{1}$ transport methods for three species of neutrinos: electron neutrinos, $\nu_{e}$, and antineutrinos, $\bar{\nu}_{e}$, and one species, $\nu_{X}$, including the mu and tau neutrinos and their antiparticles. Our set of neutrino-matter reactions contains the important interactions that dominate the dynamics of core collapse (see Just et al., 2015, 2018, for implementation details):

• •

absorption and emission of $\nu_{e}$ and $\bar{\nu_{e}}$ by $\beta$ processes of free neutrons and protons and nuclei,

• •

iso-energetic scattering of neutrinos of all flavours off nucleons and nuclei,

• •

pair creation of neutrinos of all flavours by electron-positron annihilation and nucleonic bremsstrahlung,

• •

non-iso-energetic scattering of neutrinos of all flavours off electrons and positrons.

We note that the last process is not written in terms of an opacity and thus our MIRK method does not apply. We treat it in an operator-split manner in the same way as described in Just et al. (2018). In principle, the same holds for the pair processes. However, the approximate treatment of O’Connor (2015) reformulates the interaction in terms of opacities, which allows us to include them in the MIRK scheme.

4.2 Initial data and reference simulation

As a test case, we used the same model as in the comparison of neutrino-hydrodynamics codes of Rampp & Janka (2002), i.e., the core of a star of a zero-age main-sequence mass of $M_{\mathrm{ZAMS}}=15\,M_{\odot}$. Before presenting results of the new MIRK method implementation, we describe the dynamics of a reference simulation (denoted RK2) computed with the traditional scheme used in the Alcar (Just et al., 2015). It uses a method similar to our first-order MIRK scheme as a building block in a second-order Runge-Kutta time integrator.

As the central density increases during collapse, electron captures deleptonize the matter and drive the electron fraction at the center to values $Y_{e,\mathrm{c}}\approx 0.28$ and the lepton fraction, including the net lepton number corresponding to $\nu_{e}$ and $\bar{\nu_{e}}$, to $Y_{l,\mathrm{c}}\approx 0.34$. These quantities assume a roughly constant level once the neutrinos are trapped inside the inner core as densities $\rho\gtrsim 10^{12}\,\textrm{g\,cm}^{-3}$ render the gas optically thick. The shock wave launched at core bounce (time $t_{\mathrm{b}}$) and the formation of the proto-neutron star (PNS) stalls about 70 ms later after having reached a maximum radius of $r_{\mathrm{sh;max}}\approx 145\,\mathrm{km}$, i.e., still inside the collapsing iron core (see Fig. 1, second panel). It recedes slowly for another 90 ms to $r_{\mathrm{sh}}\approx 125\,\mathrm{km}$. Matter continues to fall through the shock wave and settles onto the PNS, which gradually contracts from a maximum radius of up to $r_{\mathrm{PNS}}\lesssim 80\,\mathrm{km}$ immediately after bounce to $r_{\mathrm{PNS}}\gtrsim 20\,\mathrm{km}$ at $t_{\mathrm{pb}}=1\,\mathrm{s}$ (here we use the radius of the $\nu_{e}$-sphere as a proxy for the PNS radius). By $t_{\mathrm{pb}}=t-t_{\mathrm{b}}\approx 150\,\mathrm{ms}$, the entire iron core has been accreted. Consequently, the density and ram pressure of the accreting matter drops, which causes a brief expansion of the shock by about 10 km. Neutrinos heat the post-shock gas, but, as is typical for spherically symmetric core collapse, the conditions for shock revival and an explosion are never met. Thus, the shock wave gradually contracts to a radius below 50 km over the course of 1 s after bounce. The early neutrino emission is characterized by the intense burst of $\nu_{e}$ emitted in the first few tens of ms after bounce (third panel of Fig. 1). After the burst, the $\nu_{e}$ luminosity, $L_{\nu_{e}}$, and those of the other two flavors reach slowly varying values of several $10^{52}\,\mathrm{erg}/\mathrm{s}$. We find the typical ordering with almost equal luminosities of the electronic flavors and a lower emission of the heavy-lepton neutrinos as well as the dependence of the luminosities on the mass accretion rate that leads to the lower levels of $L_{nu}$ after $t_{\mathrm{pb}}\sim 200\,\mathrm{ms}$. The mean energies (bottom panel) with values in the range of 10–25 MeV reflect the rising temperatures near three neutrinospheres of the three flavors with the sequence $e_{\nu_{e}}<e_{\bar{\nu_{e}}}<e_{\nu_{X}}$ following from the hierarchy of neutrino-matter cross sections.

4.3 First order MIRK numerical simulations

The first order MIRK scheme has two free parameters, $a$ and $b$. We compare the four combinations of setting them to zero and to a non-zero value of $1/2$ (see Tab. 1). Simulation M1-1 with $a=b=0$ satisfies the correct optically thick limit. It produces a stable simulation whose results are very close to those of the reference simulation, both in terms of the global evolution shown in Fig. 1 and in terms of the radial profiles of Fig. 2. The density profiles at representative epochs during the evolution (top panel) are almost identical to the ones of model RK2. The PNS at the center as well as the surrounding region of decreasing density do not show any notable difference between the two simulations. The only small discrepancies appear right at the shock wave where $\rho$ falls by about an order of magnitude over a few km. The entropy and the electron/lepton fractions (second and third panels) are more sensitive than the density to the details of the neutrino treatment. Nevertheless, there are only minor deviations of model M1-1 from RK2. Apart from the shock wave, we only find small differences in the precise pattern of oscillations in the entropy behind the shock wave at late times ($t_{\mathrm{pb}}=400\,\mathrm{ms}$). If anything, the MIRK simulations might be able to resolve the shock wave more sharply. Furthermore, there is a minor offset of $Y_{l}$ outside of the shock wave. This deviation turns out to be connected to the neutrinos, not the matter, as we find a similar offset in the profiles of the neutrino luminosities (bottom panel) exterior to the shock. Among the neutrinos, we point out the relatively pronounced temporal fluctuations of the heavy lepton species, $\nu_{X}$, in particular of its mean energy, which we attribute to the fact that these neutrinos are generated and absorbed only via the relatively subdominant pair processes. Thus, they tend to bear the imprint of fluctuations at their production site at larger radii to a higher degree than the electron type neutrinos. In any case, the differences between the first order MIRK run and the reference solution are entirely within the margins of uncertainty of the latter alone.

The correct limit in the optically thick limit of the momentum equation is a crucial requisite for the stability of the simulations. Models M1-2 and M1-3 with $(a,b)=(1/2,1/2)$ and $(0,1/2)$, respectively, which do not satisfy the asymptotically correct behavior (13) for all, smooth and non-smooth, initial data, turn unstable once the core becomes optically thick at a central density $\rho_{\mathrm{c}}\gtrsim 3\times 10^{12}\,\textrm{g\,cm}^{-3}$ (see Fig. 3). The instability appears first in the form of strong fluctuations near the origin that spread outward and lead to a termination of the simulation before the bounce can occur.

Model M1-4 with $(a,b)=(1/2,0)$ evolves stably and correctly through collapse and until immediately before bounce (Fig. 4). The evolution of the central electron and lepton fractions agrees well with the reference model. After bounce, however, differences between the two models appear. Most notably, the central values of $Y_{e,l}$ do not stabilize at the levels they reached during neutrino trapping, but decrease further (note the drop of the two green lines in the top panel of Fig. 4 for $\rho_{\mathrm{c}}\gtrsim 2\times 10^{14}\,\textrm{g\,cm}^{-3}$). After about 30 ms more, they reach a minimum around $Y_{e,l}\approx 0.065$, i.e., far below the correct values. Unlike in the reference case, the two variables are almost equal at $r=0$ throughout the entire post-bounce evolution (see third panel of Fig. 5), i.e., the net neutrino lepton number is close to zero. Additionally, the center of the PNS is hotter at almost twice the entropy of that of RK2 (second panel).

Discrepancies between the models are present in the $\nu_{e}$ burst, mostly in the form of larger fluctuations in all three flavors. Afterwards model M1-4 emits a considerably lower $L_{\nu_{e}}$ and higher $L_{\bar{\nu}_{e}}$ than RK2. The mean neutrino energies lie below the ones of RK2.

Additional differences appear in the structure of the core. Until $t_{\mathrm{pb}}\approx 20\,\mathrm{ms}$, the shock wave transiently expands faster than in RK2. This phase is characterized by the appearance of a bump in density (see top panel of Fig. 5, $t_{\mathrm{pb}}$, at $r\approx 60\,\mathrm{km}$) absent from the reference case and marked differences in the entropy and $Y_{e,l}$ profiles (second and third panels). Whereas the shock wave starts to recede in the reference model after it reached its furthest expansion at $t_{\mathrm{pb}}\approx 70\,\mathrm{ms}$, it stays in M1-4 at the same radius for another 90 ms before expanding up to $r_{\mathrm{sh}}\approx 158\,\mathrm{km}$ and retreating only thereafter. The main differences between the models are found in the PNS, which is less dense with a shallower $\rho$ profile, hotter and neutron-richer in M1-4 than in RK2. In the surrounding hot bubble, the differences are less pronounced, but still far larger than the ones between the reference model and M1-1. Finally, a numerical instability develops in the PNS after almost 400 ms of post-bounce evolution.

To summarize, our MIRK scheme is able to reproduce the results of the reference simulations stably and correctly if the parameters $a$ and $b$ for the energy and momentum equations, respectively, are chosen such that they satisfy the correct limit in the optically thick regime. Using a different parameter in the momentum equation, i.e., $b\neq 0$ makes the simulations unstable once neutrinos are trapped by scattering reactions. The choice $b=0$, but $a\neq 0$, i.e., obeying the constraints in the momentum, but not the energy equation, cures this instability, but results in incorrect results once, near bounce, the emission/absorption reactions become stiff as well. The simulation can continue for several 100 ms thereafter, but the PNS properties are wrong and eventually a numerical instability ensues.

4.4 Second order MIRK numerical simulations

We performed a series of simulations using the second order MIRK scheme and explore the evolution for various combinations of the four parameters $a,a^{\prime},b,b^{\prime}$ (see Tab. 2). The basic set of simulations consists of the 16 models M2-11, …, M2-44, in which we set, following the possible choices introduced in Sect. 3, $a\in\{1/2,-1/2\}$ and $a^{\prime}=\frac{a-1}{2}\in\{-1/4,-3/4\}$ or $a^{\prime}=\frac{(1-a)^{2}}{2a}\in\{1/4,-9/4\}$, and analogously for $b$ and $b^{\prime}$ (see Tab. 3 for the values of these two functions for all the parameters we have used). The nomenclature of the models is given by the following systematic scheme: the last two digits of the generic model name M2-AB indicate the values of the parameters $a$ and $b$. Indices $\mathrm{A}=1,2,3,4$, stand for values $(a,a^{\prime})=(a,\frac{a-1}{2})=(1/2,-1/4)$, $(a,a^{\prime})=(a,\frac{(1-a)^{2}}{2a})=(1/2,1/4)$, $(a,a^{\prime})=(a,\frac{a-1}{2})=(-1/2,-3/4)$, and $(a,a^{\prime})=(a,\frac{(1-a)^{2}}{2a})=(-1/2,-9/4)$, respectively, and analogously for index $\mathrm{B}$ and parameters $(b,b^{\prime})$.

We find the same three evolutionary paths as in the first order case. Most combinations result in a numerical instability at the onset of neutrino trapping, as shown for the example of model M2-11 in Fig. 3. As in the unstable first order runs, the instability develops in the optically thick core and causes catastrophic oscillations in the electron and lepton fractions which quickly lead to a termination of the simulations.

All simulations with $(b,b^{\prime})=(b,\frac{(1-b)^{2}}{2b})=(-1/2,-9/4)$ avoid this instability, irrespective of the values of $a$ and $a^{\prime}$. However, within the basic set of the 16 models M2-11 – M2-44 only the choice $a<0$ and $a^{\prime}=\frac{(1-a)^{2}}{2a}$ produces stable and correct results that, like for M1-1, agree very well with the reference simulation both in the evolution of global quantities (Fig. 1) and in the profiles at specific times (Fig. 2). The differences with model RK2 are limited to minor details such as the width of the shock wave (second panel of Fig. 2) or a small offset in the neutrino luminosity outside the shock wave.

Models M2-41, M2-42, and M2-43 with $(b,b^{\prime})=(-1/2,-9/4)$ and $(a,a^{\prime})\neq(-1/2,-9/4)$ show the same behavior as model M1-4 (see Fig. 4 and Fig. 5). Until close to the point at which the source terms in the energy equation become stiff, they follow the reference simulation. At that point, however, they yield an incorrect PNS with too low $Y_{e,l}^{\mathrm{cnt}}$, too shallow density profiles, and too low entropy. The luminosities and mean energies show the same deviations from RK2 as found in M1-4, and all models suffer the same numerical instabilities after a time of $t_{\mathrm{pb}}\sim 300-400\,\mathrm{ms}$.

We added models M2-44-1 – M2-44-4 similar to M2-44. Their results agree well with those of model M2-44, indicating that stability and accuracy do not depend on the specific values as long as $(a^{\prime},b^{\prime})=\left(\frac{(1-a)^{2}}{2a},\frac{(1-b)^{2}}{2b}\right)$ and $a,b<0$.

Another group of simulations, models M2-51 – M2-55, probe positive values of $a$ and $b$ between $1/2$ and $1$, which according to Eq. (32) and (33) should also lead to a stable evolution. However, we find that all simulations with $b=3/4$ are unstable. If we set, as in M2-44, $(b,b^{\prime})=(-1/2,-9/4)$, we obtain a stable and correct simulation with $(a,a^{\prime})=(3/4,1/24)$ and a stable, but incorrect one with $(a,a^{\prime})=(3/4,-1/8)$.

Hence, we find that, similarly to the first order schemes, the stability is set by the parameters for integrating the momentum equation: only $b^{\prime}=\frac{(1-b)^{2}}{2b}$ and $b<0$ are stable. Among these, the ones for which the parameters for the energy equation fulfil the constraint $a^{\prime}=\frac{(1-a)^{2}}{2a}$ and $a<0$ or $1/2<a<1$ are also correct.

5 Conclusions

We have derived a Minimally Implicit Runge-Kutta method for $M_{1}$ equations for neutrino transport. We use it to treat the neutrino matter interaction terms describing reactions such as absorption, emission, and scattering in an operator-split manner separately from the (hyperbolic) transport terms of an $M_{1}$ method. In general, the stiffness of the interaction terms in the optically thick regime poses a stability problem for their time integration. The problem can be overcome by fully implicit methods, but these can be very costly because of the complex dependence of the reaction rates on the neutrino fields and the thermodynamic state of the matter. We propose a simplified approach that reduces the use of implicit terms to the minimum required for stability by evaluating the opacities and the thermodynamics at the original time step. This choice makes the resulting scheme take a form similar to that of an explicit method. The first order method is a straightforward modification of an explicit scheme with an effective, reduced time step that guarantees the stability. This method is similar to the one that has already been used in this context in the neutrino-hydrodynamics code Alcar (Just et al., 2015). Here we give a mathematical framework and a generalization to a second order method that retains the simplicity of the first order one. We already applied similar schemes in (Cordero-Carrión et al., 2023) for the resistive relativistic magnetohydrodynamic equations.

The second order method depends on two numerical parameters. We demonstrate that these parameters can be chosen in such a way that they satisfy an algebraic condition that guarantees the correct optically thick limit for the source terms for neutrino energy and momentum. In this case, the new time integrator gives stable and accurate results in simulations of the collapse of stellar cores, as we show by implementing it in the code Alcar and comparing it to its traditional solver. If, on the other hand, these conditions are violated, the simulations become unstable once the core turns optically thick.

Our scheme is simple and efficient and can be used in a wide range of similar applications. Examples can be found in contexts with rarefied gases (Koellermeier & Samaey, 2022), shallow water equations, (Koellermeier & Pimentel-García, 2022) or force-free electrodynamics in General Relativity (Mahlmann et al., 2021). There, stiff terms appear in the conservative laws coming from the corresponding scenarios. Regarding future works in supernovae simulations, we will explore more general sets of parameters, the incorporation of more complex reactions such as pair processes, and a combination with a well-balanced method (see (Castro & Parés, 2020)) to manage the fluxes in order to preserve stationary solutions.

Acknowledgments

The authors acknowledge support by the Spanish Agencia Estatal de Investigación / Ministerio de Ciencia, Innovación y Universidades through the Grants No. PID2021-125458NB-C21, PID2021-127495NB-I00 funded by MCIN/AEI/10.13039/501100011033 and by the European Union, by the Generalitat Valenciana through the Grants No. PROMETEO/2019/071 and ACIF/2019/169 - European Social Fund and the Astrophysics and High Energy Physics programme of the Generalitat Valenciana ASFAE/2022/026 funded by MCIN and the European Union NextGenerationEU (PRTR-C17.I1). MO was supported by the Ramón y Cajal programme of the Agencia Estatal de Investigación (RYC2018-024938-I). This research was partially supported by the Perimeter Institute for Theoretical Physics through the Simons Emmy Noether program. Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development and by the Province of Ontario through the Ministry of Research and Innovation.