Long-term simulations of dynamical ejecta: Homologous expansion and kilonova properties

Accurate numerical-relativity simulations are essential to study the rich phenomenology of binary neutron star systems. In this work, we focus on the material that is dynamically ejected during the merger process and on the kilonova transient it produces. Typically, radiative transfer simulations of kilonova light curves from ejecta make the assumption of homologous expansion, but this condition might not always be met at the end of usually very short numerical-relativity simulations. In this article, we adjust the infrastructure of the bam code to enable longer simulations of the dynamical ejecta with the aim of investigating when the condition of homologous expansion is satisfied. In fact, we observe that the deviations from a perfect homologous expansion are about $\lesssim 30\%$ at roughly 100 ms after the merger. While the calculation of the kilonova light curves is affected by the resolution as well as our method of simplifying the ejecta simulation, these deviations from the homologous expansion also bias the results. We determine this influence by extracting the ejecta data for different reference times and use them as input to radiative transfer simulations. Our results show that the light curves for extraction times later than 80 ms after the merger deviate by $\lesssim 0.4\mathrm{mag}$ and are mostly consistent with numerical noise. Accordingly, deviations from the homologous expansion for the dynamical ejecta component are negligible from $\sim 80\mathrm{ms}$ for the purpose of kilonova modeling.

Figure 1

The timeline of a simulation. The orange arrow and orange frame represents the “normal” simulation and the blue arrow and blue frame the “modified” simulation with the Cowling approximation. The red triangle illustrates the checkpoint used for the grid modification. As an example we show snapshots of the rest-mass density of the H4-128 and the H4-128–30 ms simulation in the $x\text{−}y$ (the orbital) plane, see Tables 1 and 2. The given time of the snapshots is relative to the merger time.

Figure 2

Modification of the grid structure to enable long-term simulations. In the top row, the grid itself is visualized: in yellow the data of the initial simulation and in blue the Schwarzschild spacetime for the extended region. In the bottom row, the rest-mass density $D$ for the grid change of one simulation is shown in the $x\text{−}y$ (the orbital) plane. The bound matter is removed.

Figure 3

Evolution of the total rest mass when bound matter is removed (solid lines) and the ejecta mass (dashed lines) in the simulation for low (red), medium (blue), and high (green) resolutions. As a comparison, the ejecta masses of the normal simulations before the grid modification are shown in faint lines. The top panel shows the results for the simulations with H4; the bottom panel shows the results for the simulations with SLy. The masses are extracted from refinement level $l=0$.

Figure 4

Snapshot of the rest-mass density $D$ in the orbital plane at $t=59\mathrm{ms}$ after the merger. From left to right: H4-128–30 ms, H4-128–34 ms, and H4-128–39 ms.

Figure 5

The evolution of the outflowing material. The snapshots show the rest-mass density $D$ in $x\text{−}y$ and $x\text{−}z$ planes on refinement level $l=2$ for the H4-144–35 ms and SLy-144–25 ms simulation at six different times after the merger. The radial velocities $v_r$ for (0.1, 0.2, 0.3, 0.4, 0.5, 0.6)$c$ are plotted as contour lines from white to dark green.

Figure 6

Analyzing the expansion of the ejecta by considering the evolution of the masses $m_r$ within a sphere with radius $r$. The quantities are extracted from refinement level $l=2$. The contour lines in white to green indicate the mean radial velocities $\overline{v_r}$ for $(0.1,0.2,0.3,0.4,0.5,0.6)c$.

Figure 7

The homology parameter $\chi$, Eq. (5), as a function of time after the merger for all simulations. The left panels show the parameter for the normal simulations extracted from refinement level $l=0$, and the right panels show the parameter for the modified simulations extracted from refinement level $l=1$.

Figure 8

The deviation of the radial velocity from homologous expansion $⟨|\mathrm{\Delta }{v}_r|⟩$, Eq. (7), as a function of time after the merger for the modified simulations extracted from refinement level $l=2$.

Figure 9

Bolometric light curves for the H4-144–35 ms (top panel) and SLy-144–25 ms (bottom panel) simulations. The input data are extracted from refinement level $l=2$. The labels give the reference time $t_0$ in milliseconds after the merger for modeling the light curves. We compute the light curves for the azimuth $\mathrm{\Phi }=0$. While the results are shown for all viewing angles ${\theta }_{\mathrm{obs}}$, the light curves for ${\theta }_{\mathrm{obs}}\ne 0$ are shown as faint lines. The deposition curve, based on the amount of energy available, is shown in dashed lines for each model.

Figure 10

Top panels: light curves for the H4-144–35 ms simulation. The input data are extracted from refinement level $l=2$. The labels give the reference time $t_0$ in milliseconds after the merger for modeling the light curves. We compute the light curves for the azimuth $\mathrm{\Phi }=0$. While the results are shown for all viewing angles ${\theta }_{\mathrm{obs}}$, the light curves for ${\theta }_{\mathrm{obs}}\ne 0$ are shown as faint lines. Bottom panels: differences of the light curves relative to the curve of $t_0=133.4\mathrm{ms}$. To compare with the Monte Carlo noise, we have calculated the light curve for $t_0=133.4\mathrm{ms}$ a total of 6 times and plotted twice the maximum difference from the mean value as a gray shadow in the background.

Figure 11

Top panels: light curves for the SLy-144–25 ms simulation. The input data are extracted from refinement level $l=2$. The labels give the reference time $t_0$ in milliseconds after the merger for modeling the light curves. We compute the light curves for the azimuth $\mathrm{\Phi }=0$. While the results are shown for all viewing angles ${\theta }_{\mathrm{obs}}$, the light curves for ${\theta }_{\mathrm{obs}}\ne 0$ are shown as faint lines. Bottom panels: differences of the light curves relative to the curve of $t_0=114.5\mathrm{ms}$. To compare with the Monte Carlo noise, we have calculated the light curve for $t_0=114.5\mathrm{ms}$ a total of 6 times and plotted twice the maximum difference from the mean value as a gray shadow in the background.

Figure 12

Maps in the $v_y-v_z$ velocity plane for the density $\rho$, electron fraction $Y_e$, temperature $T$, and “infrared” (i.e., $(1.1–2.3)\mathrm{\mu }\mathrm{m}$, corresponding to the JHK bands of Figs. 10 and 11) Planck-mean opacity ${\kappa }_{\text{planck}}^{\mathrm{IR}}$. The maps are computed at 1 day after the merger for the H4-144–35 ms model extracted at 133.4 ms (top panels) and the SLy-144–25 ms model extracted at 114.5 ms (bottom panels). Note that the velocity range is different between the two models. The black lines show the location of the infrared photosphere, calculated as the surface from which the (radial) optical depth to the grid boundary equals unity.