Late-time post-merger modeling of a compact binary: effects of relativity, r-process heating, and treatment of transport effects

The initial state for simulations is the final state of a 3D black hole-neutron star simulation carried out in Foucart et al [27]. The simulation consists of a black hole with mass 7 $M_{\odot}$ and dimensionless spin $\chi=0.8$ (aligned with the orbital angular momentum of the binary) merging with an irrotational neutron star with mass 1.4 $M_{\odot}$. During merger, the neutron star matter is modeled using the Lattimer-Swesty equation of state [28] with nuclear incompressibility $K_{0}=220\,\mathrm{MeV}$, but by the end of the simulation, when only an accretion disk of density $\rho_{0}\lesssim 10^{12}$g cm^-3 remains outside the black hole, nuclear forces outside nuclei have become unimportant (except for providing binding energy for nuclei). This is one of the older black hole-neutron star simulations carried out using the Spectral Einstein Code, but we have deliberately chosen the same system as that evolved with a Newtonian code by Fernandez et al [17]. The post-merger simulation begins about 14 ms after merger, at which time the disk has baryonic mass of about $M_{0}=0.07M_{\odot}$.

We evolve on a 2D spherical-polar grid. In the asymptotically Minkowski coordinates inherited from the merger simulation, the grid extends in radius from $r_{\rm min}=26$ km (close to the horizon) to $r_{\rm max}=120,000$ km. This outer radius, unusual for studies of this kind, allows us to study the evolution of ejecta velocity over a wide spatial range. The radial grid has 260 points. To cover so many decades in radius, we use as our radial coordinate (in which the grid spacing is uniform) $u=\log(r)$. The colatitude angular variable covers $0<\theta<\pi$ with 168 grid points. We concentrate grid points near the equator by making the standard coordinate transformation

and evolving on a uniform grid in $\zeta$. We choose $h=0.4$.

Our resulting grid is comparable to that of Fernandez et al [17], with somewhat higher equatorial angular resolution but a bit lower radial resolution (in order to have a farther outer boundary with about the same number of radial points). This grid size is fairly modest (smaller in radial resolution than Fernandez et al’s more recent studies) [18] and may well have errors of order 10% in some quantities, but since the grid is constant across runs, we can readily identify the magnitudes and directions of effects of different physical inputs.

Converting from the 3D merger grid to the 2D post-merger grid requires azimuthally averaging. The X-weighed azimuthal average of 3D function $f$ is

We average all metric functions with $X=1$: the lapse $\alpha$, shift $\beta^{i}$, 3-metric $\gamma_{ij}$, and extrinsic curvature $K_{ij}$. (Naturally, polar rather than Cartesian components must be averaged.) From $\gamma_{ij}$ we know the determinant $\gamma=\det(\gamma_{ij})$. We also average the density variable $\rho_{\star}$ with $X=1$, where $\rho_{\star}=\alpha\sqrt{\gamma}u^{t}\rho_{0}$, and $\rho_{0}$ is the rest (baryonic) density. We compute averages of the electron fraction $Y_{e}$, temperature $T$, and covariant spatial components of the 4-velocity $u_{i}$ with weight $X=\rho$. We have also tried weighting these averages by $\rho_{0}$ and averaging the pressure $P$ rather than the temperature [and then recovering $T(r,\theta)$ from the averaged $P$, $\rho_{0}$, and $Y_{e}$]; these sorts of choices have very little effect on the disk’s initial global quantities and early evolution.

We do not add a fallback matter component, as did some of the runs in [17]. The fallback is far from the black hole, making a Newtonian vs. relativistic comparison of its effects perhaps less interesting. Also, the fallback component of these early runs was not tracked far and is not accurately known. In future work, we will use more recent simulations where the bounded matter never leaves the merger grid. (Fallback interaction with the disk makes the 10-20 ms post-merger state of these newer runs less axisymmetric, so that azimuthal averaging might cause large initial disk oscillations. One could avoid this by computing the azimuthal Reynolds stress$\sim\overline{\rho v_{i}v_{j}}-\overline{\rho}\overline{v}_{i}\overline{v}_{j}$ and adding this to the Euler equation. However, one would need a model for how the azimuthal Reynolds stress decreases with time as the system axisymmetrizes.)

The post-merger grid must extend much farther than the merger grid. The matter fields can be assumed to be vacuum outside of the original merger grid, but the metric must be extrapolated. In particular, we extrapolate $\alpha-1$, $\gamma_{ij}-\delta_{ij}$, $\beta^{i}$, $K_{ij}$, which all fall off to zero as $r\rightarrow\infty$. We find that simply setting $\alpha=M/r$ (even for a best fit of $M$) and setting the others to zero is insufficiently smooth and leads to observable artifacts when outflow crosses the merger grid outer boundary. Instead, for each $\theta$, we use the last two radial points on the azimuthally averaged grid to fit each of the above metric coefficients to functions of the form $f_{1}(r)=A_{1}/r$ and $f_{2}(r)=A_{2}/r+B_{2}/r^{2}$. We then use the one that falls off most quickly in $r$. The resulting extrapolations appear at least visually smooth and reasonable, and we do not notice any artificial behavior when matter crosses into the extrapolated region.

To avoid unphysical disturbances from suddenly turning on viscous stresses, we activate the viscosity gradually and continuously, starting 7.5 ms after the simulation begins and reaching its full strength about a millisecond later.

We use the Helmholtz equation of state [36] with its standard treatment of lepton and photon contributions. This equation of state is coded to use derivatives of the free energy so as to guarantee thermodynamic consistency, which could be important when the advective disk becomes convectively unstable, and modeling the resulting convection requires being able to simulate adiabatic flows at constant composition accurately. At very high densities ($\rho>10^{1}1$ g cm^-3) and temperatures ($T>10$ MeV), we switch to analytic expressions for the leptons corresponding to the relativistic limit [37].

The baryon component is a sum of free proton, free neutron, and alpha particle ideal gases in nuclear statistical equilibrium (NSE). We implement this exactly as in [14], using the Saha-like equation

where $x_{n}$, $x_{p}$, and $x_{\alpha}$ are the mass fraction in free neutrons, free protons, and alpha particles, respectively, $n$ is the number density of nucleons, and $Q=28.29566$ MeV is the binding energy of a helium-4 nucleus. This nuclear binding energy contributes $-x_{\alpha}Q/(4m_{n})$ to the total specific internal energy.

The above prescription automatically incorporates the energy released by helium recombination, and at all densities considered, the ideal gas approximation for baryons is completely adequate. However, more sophisticated NSE models find that the equilibrium nuclear state at the low densities of our simulation is actually heavy nuclei, which will affect the binding energy released (albeit not by a large amount) and the number density. Therefore, we have run two simulations using the low-density baryon component taken from the DD2 equation of state [38]. We use the isolated baryonic component, available in tabulated form at M. Hempel’s website [39]. Similar to the Helmholtz handling of the leptonic component, we use bicubic interpolation (in $\rho_{0}$ and $T$) of the free energy and its first derivatives. We interpolate linearly in $Y_{e}$. Second derivatives (needed, e.g. for the sound speed), are found by finite differencing tabular values and smoothing in density.

We must extrapolate to densities and temperatures far below the DD2 table minimum (1660 g cm^-3 and 0.1 MeV). For this, we presume an ideal gas, taking the mean molecular weight (equivalently, the mean mass number of baryonic free particles) and nuclear binding energy per baryon to be equal to the nearest corresponding point on the table–e.g. the minimum density point of the given $T$ and $Y_{e}$ for points with $\rho_{0}$ below the table minimum.

Even the above prescription for the EoS is questionable, because it assumes the continued validity of NSE at densities and temperatures far below where this is a safe assumption. In fact, neither an assumption of fast nuclear reactions (NSE) nor an assumption of very slow reactions (constant isotope abundances) is justified, and the only adequate treatment would be to explicitly include a nuclear reaction network with many isotope abundances as new evolution variables. In particular, it is expected that r-process nucleosynthesis may have important effects on the thermal evolution of outflows [25].

In one run, we incorporate the effect of r-process heating in an approximate but physically motivated way, following the prescription in [40]. We add the following driving term to the evolution of $\rho_{\star}Y_{e}$, which automatically acts as a heating agent, and an expected partially-compensating non-NSE neutrino energy loss term.

where $t_{\rm rp}=1$ second and $Y_{e,f}=$0.38.

In our simulation, these terms are applied everywhere in the fluid. This is not realistic in the inner disk, since above $T\gtrsim 0.5$ MeV, NSE is applicable. However, the accretion time is less than $t_{\rm rp}$ for most of the matter in this region, which automatically limits the effect of the r-process on it, and including it everywhere errors to the side of the greatest possible r-process effect ¹¹1Actually, our code includes density and temperature cutoffs on the r-process terms, but due to a bug they were not applied in the run reported in this paper..

For our equation of state, we find that $\left.\frac{\partial h}{\partial Y_{e}}\right|_{\rho_{0},T}=-0.026$ and is indeed very nearly constant for $\rho_{0}\lesssim 10^{9}$ g cm^-3, $T\lesssim 0.1$ MeV, and $Y_{e}$ between 0.1 and 0.35 and predicts an energy release per rest energy of 0.007 when $Y_{e}$ changes from 0.1 to 0.38. [At $Y_{e}>0.35$, the slope of $h(Y_{e})$ flattens, but the effect of the $0.35<Y_{e}<0.38$ region on the overall energy release is very small.] The cooling term in Equation 13 removes an energy per rest energy of 0.002 for the same $Y_{e}$ change, so the overall effect is indeed a heating with efficiency 0.005. Furthermore, there is an amplification effect in the disk, in that r-process heating will increase pressure, which increases the viscous heating rate.

This r-process model should only be regarded as a crude approximation, and our implementation has several foreseeable limitations. The model has the matter move in NSE to higher $Y_{e}$, while physically it would not be in NSE. However, Foucart et al [40] point out that the difference in binding energy from NSE is typically $\sim$10% of the energy released by the r-process, so the NSE approximation is probably acceptable.

The assumption of NSE is also invoked in the computation of neutrino emission rates in the M1 code. These will, then, also be only approximately accurate in non-NSE regions. We note that the neutrino emissions leading to the cooling processes of Equations 12 and 13 come from different reactions than those modeled by the M1 emission, so it is appropriate for both to operate simultaneously.

Furthermore, note that [40] only envisioned these approximate r-process terms acting on the outflow. It is unclear how r-process nucleosynthesis would proceed for material orbiting at low density and low temperature in an accretion disk. For material in tidal tails, Metzger et al [25] find that energy deposition due to r-process nucleosynthesis has a strong dependence on the trajectory of the outflow, with the r-process being strongly suppressed after apocenter, when the density of matter starts increasing. Future studies with a full nuclear reaction network will be required to determine what heating rate is appropriate for matter evolving for multiple seconds in an accretion disk. R-process heating has also been included in a couple of other studies of post-merger dynamics [22, 41].

Viscosity transports angular momentum outward, leading to accretion of gas into the black hole and expansion of the disk outward. Observing the matter flow at a fixed radii in what are initially the outer disk, we see first an outflow and then at later times an inflow, which is the expected behavior for a viscosity-driven accretion disk. (See Section 5.2 of Frank, King, and Raine [42].) Viscous heating in the equator drives vertical convection, visible especially in density snapshots. We note that convection was not seen before restricting viscosity to $T_{r\phi}$ and $T_{\theta\phi}$ and implementing the Helmholtz equation of state.

Figure 1 shows the rate of advection of baryonic mass into the black hole and out of a large observation radius, i.e. the accretion rate $\dot{M}_{\rm acc}$ and outflow rate $\dot{M}_{\rm out}$. The initial state is not in equilibrium, and the initial $\dot{M}_{\rm acc}$ is already non-zero even before viscosity is turned on. There is a dip at early times as this initial accretion subsides and before viscosity-driven evolution begins. Once viscosity is activated, the accretion rate jumps to $\sim M_{\odot}$ s^-1 and subsequently follows an approximate power law in time, $\dot{M}_{\rm acc}\propto t^{-n}$ where $n=2.1\pm 0.1$ for run Std. This is similar to the $n=1.9$ found for viscous Newtonian hydrodynamics by Fernandez et al [8] for a slightly different system. It is observably steeper than the $n=1.5$ found for viscous relativistic hydrodynamics by Fujibayashi et al [21] for slightly different systems. (This paper used isentropic equilibrium initial data with analytically specified $Y_{e}(\rho_{0})$ and $hu_{\phi}(\Omega)$ initial functions, with $\Omega$ the angular velocity. Their black hole masses closest to ours are 6 and 10 $M_{\odot}$.)

The other runs have the same $\dot{M}_{\rm acc}$ for the first second but show noticeably slower accretion thereafter. The reason appears to be that the inner disk is more depleted at $t>1$s for other runs as compared to Std. Simulations with heavy nuclei (Nuc and rProc) unbind more material in outflow after about a second, which can clearly be seen in the plot of $\dot{M}_{\rm out}$. It can be seen even more directly in Figure 2, which plots the total bound and unbound mass on the grid as a function of time (with unbound defined as $u_{t}<-1$). The accretion rate for Diff drops well below that of Std at later times (about 5 seconds post-merger). Diff does not have greatly more unbound mass than Std, so the explanation is different. The last 5 seconds of Diff was run at a lower resolution, so its late-time evolution does have larger error than the other runs. However, the divergent evolution from Std is quite believable when one looks at disk profiles which show density, entropy, and $Y_{e}$ growing distinct for the Diff run well before the switch to low resolution.

These profiles are shown in Figure 3 for a time 2 seconds after merger. We plot the inner 9000 km. While profiles appear continuous, the fraction of matter that is unbound raises from 0.1 at $r\sim 4000$ km to 0.9 at around $r\sim 9000$ km. There is no unbound material inside $r<2300$ km. The density of the disk is indeed seen to be higher for Std than for other runs.

Most of the profiles change significantly with time, the exception being the temperature profile, which is roughly constant because it is set by the requirement of dynamical equilibrium for a thick disk: $k_{B}T\sim Mm_{p}/r$. Note, therefore, that since the density of the disk is continually decreasing, its specific entropy is continually increasing. It is important to be clear about what processes do and do not affect the specific entropy. Adiabatic compression and recombination of nuclei in NSE are reversible effects which increase temperature but not entropy. The source of entropy growth in these simulations is viscous heating. (The only other conceivable sources, shock heating and neutrino absorption, are comparatively unimportant.) During the first half-second or so, this heating is partly offset by neutrino cooling, but afterwards the disk enters and advective state, and the observed continual entropy growth is expected. It should be recognized, however, that recombination in NSE can indirectly affect the entropy, because increasing the temperature increases the viscosity $\eta$ (cf. Equation 4). This provides one means whereby Std and Nuc, which differ primarily in the latter’s inclusion of recombination to heavy nuclei, will have distinct entropy evolutions.

By 2 seconds, the $Y_{e}$ source term in rProc has increased $Y_{e}$ to near its target value everywhere. The r-processing of the innermost disk is mostly the result of the unphysical inclusion of the r-process in those high temperature regions, but it should also be remembered that the gas near the black hole at $t=2$ seconds is not the gas that was initially there, but gas that was initially farther away at lower temperature that has accreted inward. The viscous timescale $\tau_{\rm visc}\sim\alpha^{-1}(r/H)^{2}\sqrt{r^{3}/M}$ is 0.33 s for $r\approx 100$ km and 2 s for $r\approx 340$ km, around where $T=0.5$ MeV (although the gas near the black hole at 2 s would have presumably spent most of its time inside this radius).

For Diff, the diffusion effects have succeeded in partially leveling $Y_{e}$ (increasing it in the inner region) and $S$ (decreasing it in the very hot inner region). Some care must be taken with this interpretation, though, since the effective heat conduction alters the vertical entropy gradient and partially suppresses convection [26]. Thus, Std and Diff are really two models of mixing.

Neutrino luminosities are shown in Figure 4. We plot three neutrino species: electron neutrinos, electron anti-neutrinos, and all other flavors and their anti-particles (labeled “x”). There are no noticeable differences between cases. Very early in the simulation, the disk becomes completely optically thin. Neutrino cooling remains important for the first few hundred milliseconds. Its effect on the disk is most easily seen in the mass-averaged $Y_{e}$, which climbs during this time due to greater emission of $\nu_{\overline{e}}$ than $\nu_{e}$. After $t=0.5$ s, the luminosities begin to fall off more steeply, and the disk becomes radiatively inefficient. After one second, neutrino emission is negligible, the average $Y_{e}$ on the grid is nearly constant, and evolution of the neutrino fields is discontinued.

In Figure 5 we plot the mass-weighted average of $Y_{e}$ of bound matter on the grid, denoted $\langle Y_{e}\rangle$. During the first second, evolution of $\langle Y_{e}\rangle$ is dominated by neutrino processes (mostly emission). After the first second, neutrino effects become negligible, and $Y_{e}$ of a fluid element can only change by the r-process (for run rProc). $\langle Y_{e}\rangle$ can continue to change by accretion and outflow, e.g. it will decrease if gas of above average $Y_{e}$ falls into the black hole or becomes unbound. In run Nuc, $\langle Y_{e}\rangle$ drops and then rises mainly because first higher $Y_{e}$ matter and then lower $Y_{e}$ matter is marked as unbound. The average $Y_{e}$ for all matter on the grid is more nearly constant during the fall and rise times of $\langle Y_{e}\rangle$ for bound matter, although the average $Y_{e}$ of all matter on the grid drops between about 3 and 5 seconds when higher-$Y_{e}$ outflow leaves the grid.

Disks with efficient neutrino cooling are subject to a regulating mechanism that keeps $Y_{e}\approx 0.1$ [43, 44], but this will not operate for this disk for most of its evolution. This figure also shows the mass-weighted average entropy of bound matter, with increases as the disk density decreases while the temperature at each radius changes much less.

Given an observation radius $r_{\rm obs}$, we compute outflows as follows. First, we only consider unbound matter. We note that matter at $r>4\times 10^{4}$ km is all unbound by the definition $u_{t}<-1$ and hence also by the definition $hu_{t}/h(T=0)<-1$. The definition is more important in the interior, and we use the stricter condition $u_{t}<-1$. The ejected mass measured at radius $r_{\rm obs}$ at time $t$ is defined to be the sum of the unbound mass that has passed through $r_{\rm obs}$ by time $t$ and the unbound mass interior to $r_{\rm obs}$ at time $t$.

where $f_{\rm ub}=1$ if the fluid at a point is unbound and $=0$ if it is bound. By $t=10$ s, the second term is negligible for $r_{\rm obs}<10^{4}$ km.

The mass-averaged value of intensive variable $X$ is defined as

One important intensive variable is the asymptotic velocity $V_{\infty}$, a proxy for the asymptotic Lorentz factor, which, continuing to use $u_{t}$ as an energy proxy, we compute as

Finally, we compute the mass distribution of $V_{\infty}$ and $Y_{e}$ in the outflow. This can be defined in two ways. Consider intensive quantity $X$, whose range is divided into bins with equal-width $\Delta X$. Ejecta histograms often display the amount of outflow mass in each bin, i.e. for the $i$-th bin, $\Delta M_{X}(X_{i}+\Delta X/2)$ is defined to be the amount of outflow mass with $X$ between $X_{i}$ and $X_{i}+\Delta X$. A plot of $\Delta M_{X}$ shows at a glance the amount of mass in each bin, but of course this number depends on the chosen $\Delta X$. This dependence can be removed by defining the distribution $M_{X}=\Delta M_{X}/\Delta X$. Below, we plot $\Delta M_{X}$ but provide $\Delta X$ so that readers can convert to $M_{X}$ if they wish.

Ejecta masses and mass-averaged quantities are given in Table 1. Considering first the ejecta mass itself, we do see noticeable effects of our modeling choices. In particular, 60% more matter is unbound from the disk if one uses the nuclear equation of state component of Nuc and rProc. R-process heating itself seems to have a much more modest effect. This is consistent with Just et al [~\cite[cite]{[\@@bibref{}{10.1093/mnras/stv009}{}{}]}], whose simulations with an r-process heating prescription show little difference in ejected mass. For comparison, Fernandez et al [17] evolving this same case (their model Fdisk), found an ejected mass of $0.0048M_{\odot}$. Our model with the same equation of state (run Std) has $M_{\rm out}=0.0054M_{\odot}$, which is quite close given the differences in the treatment of gravity and neutrinos.

We see a greater difference from the Newtonian leakage simulation [17] in $\langle V_{\infty}\rangle$. The average $V_{\infty}$ measured at $r_{\rm ob}=10^{4}$ km is 0.039 $c$ in the Newtonian simulation but 0.054 $c$ in Std. Might this be a difference of the relativistic treatment? Relativistic MHD simulations (including 2D MHD dynamo simulations) do produce much faster average outflow velocities, of around 0.1 $c$ [8, 12], but this is largely because of an early-time high-velocity MHD outflow, which will also be absent in the relativistic viscous simulations reported here. The late-time “thermal” disk outflow in MHD simulations is closer in speed to that of their viscous counterparts. 2D relativistic viscous hydrodynamic simulations of tori resembling post-merger configurations have been carried out by Fujibayashi et al [20, 21]. Their initial conditions are artificial (exact equilibria, constant entropy), but they vary black hole and disk mass, as well as initial angular velocity profile. They find ejecta velocities in the range 0.05-0.1 $c$. The closest configuration in that paper to the one studied here was perhaps model M06L05, with $M_{\rm BH}=6M_{\odot}$, $M_{\rm disk}=0.1M_{\odot}$, $a_{\rm BH}/M_{\rm BH}=0.8$. At simulation end, this model has $M_{\rm out}=0.0087M_{\odot}$, $\langle Y_{e}\rangle=0.379$, $\langle S\rangle=27.5k_{B}$, $\langle V_{\infty}\rangle\approx 0.07c$. Their equation of state is DD2 at high densities, so run Nuc is presumably the best comparison. Given that our initial states are still somewhat different, their $M_{\rm out}$, $\langle S\rangle$, and $\langle V_{\infty}\rangle$ are reasonably close to ours. Overall, Fujibayashi et al’s results do indicate that ejecta velocities $\gtrsim 0.05c$ are common in relativistic viscous evolutions of post-merger scenarios, and our having simulated an exact case with Newtonian counterpart confirms that the difference is probably due to the relativistic treatment. In fact, our average asymptotic velocities measured at 10,000 km underestimate the true asymptotic velocity, because measurement at further radii demonstrates that the outflow is still gaining kinetic energy beyond that point, a point to which we will return.

Our $\langle Y_{e}\rangle\approx 0.3$ is somewhat smaller than that found by Fernandez et al [17] for this case, $\langle Y_{e}\rangle=0.35$, and that Fujibayashi et al [21] found for their model M06L05, $\langle Y_{e}\rangle=0.379$. It is close to the values of $\langle Y_{e}\rangle$ that Fujibayashi et al [20] and Just et al [19] found for outflows from systems with lower mass black holes. It is difficult to identify the cause of this difference. All simulations except that of Fernandez et al involve somewhat different initial states, and neutrino transport is handled at least slightly differently by all codes. One concern about the simulations in this paper is the use of lower density cutoffs on neutrino interactions and on turning off neutrino effects altogether after one second. To check the effect of these simplifications, we reran Nuc and Diff for 2 seconds with density cutoffs a factor of 100 lower. We note that the mass-averaged $Y_{e}$ for unbound mass on the grid settles by about one second. We find that the change in outflow composition is completely negligible. It remains possible that the cause of the differences is in the transport and leakage algorithms themselves, in which case the correct outflow composition for viscous hydrodynamic evolution will only be determined by a full neutrino transport scheme.

Distributions of $Y_{e}$ and $V_{\infty}$ are shown in Figure 6. Particle diffusion effects seem to be able to modestly reduce the spread of $Y_{e}$ in the Diff run, as might be expected. In fact, diffusion terms exert a less obvious but stronger effect on the $V_{\infty}$ distribution, which is bimodal but has a very sharp peak at around 0.09 $c$.

The $Y_{e}$ distributions are often found to exhibit a double peak structure. In binary neutron star mergers, there is often a low-$Y_{e}$ equatorial ejecta and a higher-$Y_{e}$, more isotropic ejecta, but in our case all components of the ejecta with appreciable mass are equatorially concentrated. Instead, the two peaks represent ejecta from different times. As mentioned in the discussion of evolution of average $Y_{e}$ for bound matter, the ejecta from the first half of the evolution has a higher average $Y_{e}$, while the later ejecta has a narrower $Y_{e}$ distribution peaked at around 0.28. An indication of this can be seen in the left-hand panel of Figure 7, which shows $Y_{e}$ ejecta distributions for Std at times of 5 and 10 seconds, and the effect of the time interval from 5 to 10 seconds is to accumulate outflow at the lower $Y_{e}$ peak (although even at 5 seconds, the late-time peak has begun to appear).

The r-process source terms have the anticipated effect of driving all the outflow in the rProc run to the target $Y_{e}$ i.e. the expected $Y_{e}$ at the end of the r-process). We note that this should not be understood as the rProc simulation predicting different values of $Y_{e}$ at the beginning of the r-process, and thus different r-process yields. Instead, the rProc simulation shows $Y_{e}$ at the end of the r-process, while all other simulations show $Y_{e}$ before r-process nucleosynthesis; the distributions of $Y_{e}$ at the end of these simulations are thus not directly comparable.

The simulated r-process does not increase the total mass of ejecta or its entropy, despite having deposited thermal energy into the gas. In fact, these numbers are a bit lower for rProc than they are for Nuc. Just et al [19] also attempted to incorporate r-process effects by adding a heating term to outgoing fluid, and also found rather small effects on global outflow quantities. That the effect in our case is to lower $M_{\rm out}$ and $\langle S\rangle$, however, remains counter-intuitive. Possibly this is because the r-process timescale is similar to the timescale during which neutrinos are important, so during the first second, gas could undergo r-process alteration of $Y_{e}$ and radiated the resulting heat both from the $\dot{Q}$ of our model (Eq. 12) and through the M1 neutrino transport active then, and this might greatly reduce the dominance of heating over cooling compared to what was anticipated from the analysis in Section 2.4. To check this, we re-ran Nuc from 1 to 5 seconds with r-process effects starting only at 1 second, after neutrino effects become negligible. We find that adding the r-process source terms increases the outflow mass (defined as in Eq. 14 but with the $t$ integral from 1 to 5) by 30%. This confirms that our r-process terms in a simpler setting do energize and unbind matter. However, a reliable model of these effects will require a consistent treatment of neutrino losses. Another expected effect of $r$-process heating on ejecta is a smoothing of spatial inhomogeneities, but this would not be visible on the length and timescales of this simulation [45, 22, 46].

In Figure 7, we investigate the sensitivity of our distributions to the observation radius and time. $M_{Y}$ is insensitive to observation radius, since $Y_{e}$ of the outflow just advects at these late times. However, we see that a significant fraction of the outflow is generated in the last half of the simulation, the last five seconds, so that evolution to at least ten seconds was indeed necessary. The $V_{\infty}$ profiles are more sensitive to extraction radius. This could already be seen in the higher $\langle V_{\infty}\rangle$ in Table 1. The question arises whether this acceleration is a physical effect that can be explained by pressure forces. One way to estimate the possible importance of pressure forces is to compute the specific energy, and hence asymptotic velocity, from $hu_{t}/h(T=0)$ rather than from $u_{t}$. Because the outflow is not steady, $hu_{t}$ is not a constant for fluid elements–it could not be, since for outflow of finite width, pressure forces must push fluid at the inner edge inward rather than outward–so this does not provide a more realistic asymptotic velocity, but it does indicate the magnitude of the effect pressure forces can have. We see that pressure forces can quite easily push the asymptotic velocity above 0.1 $c$. Another concern is whether pressure could be anomalously high in our simulations. However, the average specific entropy reported in Table 1 is well in line with entropy reported in other studies [18, 21]. Thus, it is likely that the observation radius of $r_{\rm ob}=$10,000 km, which is fairly common in prior studies, underestimates the asymptotic velocities in earlier works, even apart from nuclear heating effects which were known to be missing. Confirmation of this interpretation comes from the study of Kasen et al [47], in which disk outflow at $10^{4}$ km was injected onto a much larger grid, and it was found that homology was not achieved until $\sim 100$ s.

Strong acceleration from $10^{4}$ km and $10^{5}$ km is present in Nuc as well as in Std, as clearly seen from Table 1 and Figure 7. For Nuc, there seems to be less internal energy available for acceleration beyond $r>10^{5}$ km than there is for Std. The asymptotic velocity using $hu_{t}$ is similar in both cases, but the outflow does more of its acceleration inside $10^{5}$ km in the Nuc case.