Formulating bulk viscosity for neutron star simulations

In order to extract the precise physical information encoded in the gravitational and electromagnetic signals from powerful neutron-star merger events, we need to include as much of the relevant physics as possible in our numerical simulations. This presents a severe challenge, given that many of the involved parameters are poorly constrained. In this paper we focus on the role of nuclear reactions. Combining a theoretical discussion with an analysis connecting to state-of-the-art simulations, we outline multiple arguments that lead to a reactive system being described in terms of a bulk viscosity. The results demonstrate that in order to properly account for nuclear reactions, future simulations must be able to handle different regimes where rather different assumptions/approximations are appropriate. We also touch upon the link to models based on the large-eddy-strategy required to capture turbulence.

Figure 1

Illustrating the behavior of $\beta (\omega )$ in two cases: the solution to the Cattaneo equation (i.e., the full extended irreversible thermodynamics (EIT) relaxation toward the Navier-Stokes limit, blue solid curve) and the parabolic case (the Navier-Stokes limit, orange dashed line). The shaded region indicates frequencies we assume we may not have “access” to numerically (this region moves toward higher frequencies when the numerical resolution is increased).

Figure 2

The maximum (for each point we assume that $\omega =\mathcal{A}$) potential relative contribution of the bulk viscous approximation ${\mathrm{\Pi }}^{\max }/p^{\mathrm{eq}}$ at each point in phase space using the APR [33, 49] equation of state. We see that the bulk viscous pressure contribution can be large for most temperatures when $n\gtrsim {10}^{-3}{n}_{\mathrm{sat}}$. However, the bulk viscous approximation should only be used where the reaction rate cannot be resolved by the numerical simulation, which is where the grid frequency is greater than $\mathcal{A}$. Also shown are contours at $\mathcal{A}=\{{10}^3,{10}^4,{10}^7,{10}^9\}{\mathrm{s}}^{-1}$ (solid, dashed, dot-dash, dot). For current simulations, frequencies of $\sim {10}^6{\mathrm{s}}^{-1}$ are resolvable. This shows that the bulk viscous approximation should only be used for $T\gtrsim 10\mathrm{MeV}$, and as the grid resolution improves becomes less necessary.

Figure 3

Plot of $\mathcal{A}$ for the APR equation of state used in [29]. The restoring term $\gamma$ is calculated assuming the Fermi surface approximation remains valid. Contours are at $\mathcal{A}=\{{10}^3,{10}^4,{10}^7,{10}^9\}{\mathrm{s}}^{-1}$ (solid, dash, dot-dash, dot).

Figure 4

Plot of the inverse of the neutrino equilibration timescale $1/{\tau }_{\mathrm{eq}}$. Contours are at $1/{\tau }_{\mathrm{eq}}=\{{10}^3,{10}^4,{10}^7,{10}^9\}{\mathrm{s}}^{-1}$ (solid, dash, dot-dash, dot).

Figure 5

Copy of Fig. 2 with contours from Fig. 4.