Tabulated NSE and aprox19

The aprox19 network can be run in a manner where we blends the standard aprox19 network with a table for nuclear statistic equilibrium resulting from a much larger network at high density and temperatures. This option is enabled by building with:

NETWORK_DIR=aprox19 USE_NSE_TABLE=TRUE

Interface

The main interface for the NSE table is defined in nse_table.H:

AMREX_GPU_HOST_DEVICE AMREX_INLINE
void nse_interp(nse_table_t& nse_state, bool skip_X_fill=false)

Here, nse_table_t is a simple struct that holds the thermodynamic and NSE state and the optional parameter skip_X_fill can be used to skip interpolating the mass fractions if desired (since there are a lot of them, and the interpolation can be expensive).

Composition and EOS

The NSE table was generated using pynucastro [10, 11], using 96 nuclei and electron/positron capture/decay rates from [37]. The table takes $Y_e$ as the primary composition variable and provides a set of mass fractions that is mapped into those used by aprox19. Using the value allows us to attain a lower $Y_e$ than aprox19 can represent.

Note

The full details of the NSE table are provided in [36]. The table can be regenerated using the script nse_tabular/make_nse_table.py.

When we are using the NSE network, we always take the composition quantities in the EOS directly from eos_state.aux[] instead of from eos_state.xn[]. The AUX_THERMO preprocessor variable is enabled in this case, and the equations of state interpret this to use the auxiliary data for the composition. This is described in Auxiliary composition.

NSE Table Outputs

The NSE table provides values for the auxiliary composition, $\overline{A}$, and $⟨B/A⟩$ resulting from the full 96 nuclei network. It also provides a set of 19 $X_k$ that map into the isotopes carried by aprox19. These three quantities are stored as aux data in the network and are indexed as iye, iabar, and ibea. Additionally, when coupling to hydrodynamics, we need to advect these auxiliary quantities.

The evolution equations for the auxiliary variables are:

$$\begin{array}{r}\begin{array}{rl}\frac{D{Y}_e}{Dt}& =\sum _k\frac{Z_k}{A_k}{\dot{\omega }}_k\\ \frac{D\overline{A}}{Dt}& =-{\overline{A}}^2\sum _k\frac{1}{A_k}{\dot{\omega }}_k\\ \frac{D}{Dt}\left(\frac{B}{A}\right)& =\sum _k\frac{B_k}{A_k}{\dot{\omega }}_k\end{array}\end{array}$$

Therefore each of these auxiliary equations obeys an advection equation in the hydro part of the advancement.

The table also provides $d{Y}_e/dt$, $(d⟨B/A⟩/dt{)}_{\mathrm{weak}}$, and ${ϵ}_{\nu ,\mathrm{react}}$, the weak rate neutrino losses. These quantities are used to update the thermodynamic state as we integrate.

NSE Flow

The time integration algorithm is described in detail in [36]. Here we provide an outline:

• initialize the problem, including $X_k$

• fill the initial aux data with $Y_e$, $\overline{A}$, and $(B/A)$

• in hydro, we will update these quantities simply via advection (for Strang-split evolution)

• for the reactive update:

  • check if NSE applies (see below)

  • if we are in an NSE region:

    • Compute the initial temperature given $\rho$, $e$, and $Y_e$, using an EOS inversion algorithm that understands NSE (in particular that the composition depends on $T$ in NSE)

    • Compute the plasma neutrino losses, ${ϵ}_{\nu ,\mathrm{thermal}}$

    • Use $\rho$, $T$, and $Y_e$ to evaluate the NSE state and construct $[\mathbf{R}({\mathcal{U}}^{\prime }){]}^n$, the source term from reactions to the reduced conserved state ${\mathcal{U}}^{\prime }$ (this is the state used by the SDC algorithm and includes the internal energy density, mass fractions, and auxiliary variables).

      This is done via finite differencing in time (through a step $\tau \ll \mathrm{\Delta }t$), and the reactive sources are constructed to exclude the advective contributions. The size of $\tau$ is controlled via integrator.nse_deriv_dt_factor.

      In particular, the energy source is constructed as:

      $$R(\rho e)=N_A\frac{\mathrm{\Delta }(\rho ⟨B/A⟩)}{\tau }+N_A\mathrm{\Delta }{m}_{np}{c}^2\rho \frac{d{Y}_e}{dt}-\rho ({ϵ}_{\nu ,\mathrm{thermal}}+{ϵ}_{\nu ,\mathrm{react}})$$

      where $\mathrm{\Delta }{m}_{np}$ is the difference between the neutron and H atom mass.

      Important

      It only makes sense to include the weak rate neutrino losses, ${ϵ}_{\nu ,\mathrm{react}}$, if the initial model that you are using in your simulation also included those losses. Otherwise, the energy loss from our NSE table will likely be too great and that simulation will not be in equilibrium. This is an issue, for example, when using a MESA model constructed with aprox21, which does not have all of the weak rates we model here.

      The weak rate neutrino losses can be disabled by integrator.nse_include_enu_weak=0.

    • Predict ${\mathcal{U}}^{\prime }$ to the midpoint in time, $n+1/2$ and construct $[\mathbf{R}({\mathcal{U}}^{\prime }){]}^{n+1/2}$.

    • Do the final update to time $n$ as:

      $${\mathcal{U}}^{\prime ,n+1/2}={\mathcal{U}}^{\prime ,n}+\frac{\mathrm{\Delta }t}{2}[\mathcal{A}\left({\mathcal{U}}^{\prime }\right){]}^{n+1/2}+\frac{\mathrm{\Delta }t}{2}[\mathbf{R}({\mathcal{U}}^{\prime }){]}^{n+1/2}$$

      where $[\mathcal{A}\left({\mathcal{U}}^{\prime }\right){]}^{n+1/2}$ are the advective updates carried by the SDC algorithm.

    • Compute the energy generation rate from the change in internal energy from ${\mathcal{U}}^{\prime ,n}$ to ${\mathcal{U}}^{\prime ,n+1}$, excluding advection.

    • Update the total energy.

    • Set the mass fractions carried on the grid from the NSE table (with the new temperature and $Y_e$).

  • if we are not in NSE:

    • integrate the aprox19 network as usual

    • update the aux quantities at the end of the burn

NSE check

For a zone to be consider in NSE, we require $\rho$ > network.rho_nse and either

• $T$ > network.T_nse together with the composition check

• $T$ > network.T_always_nse

where we assume that T_always_nse > T_nse.

The composition check considers the following nuclei groups:

• He-group: atomic numbers 1 to 2 (H to He)

• C-group: atomic numbers 6 to 7 (C to N)

• O-group: atomic number 8 (O)

• Si-group: atomic number 14 (Si)

• Fe-group: atomic numbers 24 to 30 (Cr to Zn)

and we then say that a composition supports NSE if:

• $X({\mathrm{C}}_{\mathrm{group}})$ < network.C_nse

• $X({\mathrm{O}}_{\mathrm{group}})$ < network.O_nse

• $X({\mathrm{Si}}_{\mathrm{group}})$ < network.Si_nse

• $X({\mathrm{Fe}}_{\mathrm{group}})+X({\mathrm{He}}_{\mathrm{group}})$ > network.He_Fe_nse

NSE table ranges

The NSE table was created for:

• $9.4<{\log }_{10}(T)<10.4$

• $7<{\log }_{10}(\rho )<10$

• $0.43<Y_e<0.5$