Towards an accurate description of an accretion induced collapse and the associated ejected mass

TL;DR

This paper models the accretion-induced collapse of white dwarfs into neutron stars, focusing on accurately simulating the ejected mass and its potential observational signatures, with implications for fast radio bursts and optical transients.

Contribution

It provides detailed hydrodynamic simulations of AIC with improved treatment of the equation of state and nuclear burning, estimating ejected mass across various parameters.

Findings

Ejected mass ranges from 0.01 to 0.1 solar masses.

Simulation results are explained by a simple theoretical model.

Focus on accurate hydrodynamics without neutrino or GR effects.

Abstract

We revisit the accretion-induced collapse (AIC) process, in which a white dwarf collapses into a neutron star. We are motivated by the persistent radio source associated with the fast radio burst FRB 121102, which was explained by Waxman as a weak stellar explosion with a small ($\sim 10^{-5}M_{\odot}$) mildly relativistic mass ejection that may be consistent with AIC. Additionally, the interaction of the relatively low ejected mass with a pre-collapse wind might be related to fast optical transients. The AIC is simulated with a one-dimensional, Lagrangian, Newtonian hydrodynamic code. We put an emphasis on accurately treating the equation of state and the nuclear burning, which is required for any study that attempts to accurately simulate AIC. We leave subjects such as neutrino physics and general relativity corrections for future work. Using an existing initial profile and our own initial profiles, we find that the ejected mass is $\sim 10^{-2}$ to $10^{-1}M_{\odot}$ over a wide range of parameters, and we construct a simple model to explain our results.

Equations

$$R\approx {\left(\frac{3}{4\pi }\frac{M_{\text{Ch}}}{{\rho }_{nuc}}\right)}^{1/3}\cong 20\text{km},$$

$$ϵ\sim \frac{3G{M}_{\text{Ch}}}{5R}\sim 50\text{MeV/baryon},$$

$$F=\xi (\rho )F_{\text{SNA}}+(1-\xi (\rho ))F_{\text{NSE}},$$

$$\xi (\rho )=\frac{1}{2}\left\{1+\tanh \left[\frac{{\log }_{10}\left(\rho \right)-{\log }_{10}\left({\rho }_t\right)}{{\rho }_{\delta }}\right]\right\},$$

$$\frac{\partial p}{\partial r}=-\frac{Gm\rho }{r^2},$$

$$u^2=\left(p_1-p_0\right)\left(V_0-V_1\right)$$

$${ϵ}_1-{ϵ}_0=\frac{1}{2}\left(p_1+p_0\right)\left(V_0-V_1\right),$$

$$u={\kappa }^{\frac{1}{2}}{G}^{\frac{1-\gamma }{2}}{\left(t_c-t_{\text{bounce}}\right)}^{1-\gamma }U(X),$$

$$D(0)=G{\left(t_c-t_{\text{bounce}}\right)}^2{\rho }_{\text{bounce}}.$$

$$Tds=\delta Q,$$

$$\delta Q_{nuc}=de-\sum _i{\mu }_{i}d{y}_i,$$

$$\frac{du}{dt}=-\frac{1}{\rho }\frac{\partial p}{\partial r}-\frac{mG}{r^2},$$

$$Tds=-\delta y_e\left(\frac{1}{6}{\mu }_e-\left({\mu }_n-{\mu }_p\right)\right),$$

$$p=\kappa {\rho }^{\gamma },$$

$$X=r/r_0(t),r_0(t)={\kappa }^{\frac{1}{2}}{G}^{\frac{1-\gamma }{2}}{\left(t_c-t\right)}^{2-\gamma }.$$

Full text

hydrodynamics – white dwarfs – stars: neutron

††journal: ApJ

.

1 Introduction

accretion-induced collapse (AIC) is a theorized process in which a white dwarf (WD) collapses into a neutron star (NS), followed by an explosion that ejects a fraction of the star’s mass at mildly relativistic velocities. AIC has been proposed as the outcome of an accreting WD and as a possible NS formation channel (Canal et al., 1980; Nomoto, 1986). Existence of young pulsars in globular clusters (Lyne et al., 1996) suggests that some fraction of NSs indeed were created in this way. A renewed interest in this process has arisen recently, following new astronomical discoveries that may be related (Waxman, 2017; Lyutikov & Toonen, 2019; Moriya, 2019).

One discovery that may be related to AIC is the persistent radio source associated with FRB 121102 (Scholz et al., 2016). The fast radio burst (FRB) source resides in a dwarf galaxy at a distance of $\sim 970\,\text{Mpc}$, with a persistent radio source located in that direction (Chatterjee et al., 2017). Waxman (2017) suggested that the persistent radio source was created by a weak stellar explosion, with a small ($\sim 10^{-5}M_{\odot}$) mildly relativistic mass ejection, which may be consistent with AIC, while the resulting NS acted as the source of the FRBs (see also Margalit et al., 2019; Kashiyama & Murase, 2017).

Another window for observing AIC might be fast-rising blue optical transients (Drout et al., 2014). The interaction of the ejected mass with a pre-collapse wind could be related to these events (Lyutikov & Toonen, 2019), which are characterized by a short optical rise time, $<10\,\text{days}$, and a peak luminosity comparable to supernovae. For analysis of the expected light curve of AIC, see Metzger et al. (2009) and Darbha et al. (2010).

The AIC process is the consequence of a WD, which is being held mostly by electron degeneracy pressure, that accumulates enough mass to surpass the Chandrasekhar mass limit, $M_{\text{Ch}}$, at which point it can no longer support its own mass and collapses. The collapse should happen without a prior thermonuclear runaway that will explode the star before it is able to collapse. The conditions that determine the fate of the WD were studied, for example, by Nomoto (1986) and Tauris, T. M. et al. (2013). In this work we assume that a nuclear runaway did not occur and that the WD was able to collapse. In this case, the energies gained by nuclear reactions ($\sim\text{MeV/baryon}$) are not sufficient for a nuclear runaway, since the gravitational binding energy is much deeper, as shown below. The WD collapses over its own freefall time, $t_{ff}\sim 0.1\,\text{s}$, until the core reaches nuclear densities ($\rho_{nuc}\sim 10^{14}\,\text{g}\>\text{cm}^{-3}$) and repelling strong interactions halt the collapse. The radius of the core at bounce is

[TABLE]

roughly a factor of $5$ larger than the Schwarzschild radius, $r_{s}=2GM_{\text{Ch}}/c^{2}\cong 4\,\text{km}$. The characteristic velocity of the ejected material is comparable to the escape velocity, $v_{esc}=\sqrt{r_{s}/R}c\cong 0.5c$, and is mildly relativistic. The gravitational energy per baryon around the core is

[TABLE]

which is converted to internal and kinetic energy, corresponding to a temperature of $T\lesssim 5\times 10^{11}\,\text{K}$. During the collapse, the outer layers bounce off the dense core, creating a shock wave that propagates outward, ejecting some mass in the process. Since most of the gravitational energy is converted to internal energy, only some fraction of the initial mass will be ejected, while the remaining core becomes an NS.

Simulating AIC raises some challenges that should be carefully addressed. The range of densities and temperatures throughout the process varies greatly, from nuclear densities ($\rho_{nuc}\sim 10^{14}\,\text{g}\>\text{cm}^{-3}$) at core bounce to very low densities ($<10^{5}\,\text{g}\>\text{cm}^{-3}$) at the edge of the star. The entire range should be described accurately and smoothly by the equation of state (EOS). The high temperatures and densities cause nuclear reactions to be very rapid, and these reactions should be treated in a self-consistent way. Neutrinos play an important part in the process (Woosley & Baron, 1992), and general relativity (GR) corrections should also be accounted for, as the radius of the core after bounce is comparable to the Schwarzschild radius, leading to alterations of the proto-neutron star (PNS Liebendörfer et al., 2001).

Several previous studies have estimated the amount of ejected mass from an AIC. Woosley & Baron (1992) estimated an ejected mass of $M_{\text{ej}}\sim 0.01M_{\odot}$, using a progenitor made by Nomoto (1986), and one-dimensional (1D) hydrodynamic simulations that include neutrino transport. They used the EOS that was developed by Baron, Cooperstein, and Kahana (BCK; Baron et al., 1985), and no GR corrections were included. Matter was flash burned to nuclear statistical equilibrium (NSE) when the temperature exceeded $5\times 10^{9}\,\text{K}$. In their simulations, the prompt bounce shock failed to explode the star, and the mass loss occurred through neutrino-driven wind. Fryer et al. (1999) used the same initial profile and the EOS of Herant et al. (1994), which couples the Lattimer-Swesty (LS) nuclear EOS (Lattimer & Swesty, 1991) with a low-density EOS, and found that the prompt bounce shock ejects most of the mass. They obtained ejecta mass of $M_{\text{ej}}\sim 0.1-0.3\,M_{\odot}$ for different variants of the LS EOS, nuclear burning assumptions, neutrino physics assumptions, and GR effects. They showed that the main difference with the work of Woosley & Baron (1992) was the different EOS, where the results from the BCK EOS never agree with the results from the LS EOS, despite a wide range of attempted physical parameters. This was also the case in Swesty et al. (1994), who compared the two EOSs and argued that the BCK EOS is less accurate. Fryer et al. (1999) further compared the obtained low-$Y_{e}$ ejected mass with the abundance of heavy elements in the Galaxy to derive an upper limit on the AIC rate in the Galaxy of $\sim 10^{-5}\,\text{yr}^{-1}$. Dessart et al. (2006) performed multidimensional simulations and found an ejected mass of $M_{\text{ej}}\sim 10^{-3}M_{\odot}$. They used the HShen EOS (Shen et al., 1998), extended to low densities and temperatures, and included neutrino emission and transport, but they did not include nuclear burning or GR corrections. The ejection was through a successful shock followed by a neutrino-driven wind. Dessart et al. (2006) suggested that the discrepancy with the previous works is because of the multidimensional treatment and the high-rotating progenitor.

1.1 Objectives

A careful treatment of the EOS is an essential part in core-collapse supernova simulations, since the dynamics is determined by the behavior of matter at nuclear densities (Bethe, 1990; Schneider et al., 2017). The same holds also for the AIC process, although not studied in a comparable detail to the core-collapse supernova problem. A significant complication for a calculation of the AIC problem is that the edge of the star must be included in the simulation, in which the density drops significantly, and thus an accurate description of both the high- and low-density regimes, as well as the intermediate regime that connect them, is required. For this purpose, the preparation of a relevant EOS for several regimes (of density, temperature, composition) is necessary, and the different regimes should be smoothly connected. Since only a small fraction of the star is being ejected, the intermediate regimes, in which matter can change its density and temperature by several orders of magnitudes during the collapse, play an important role in determining the amount of ejected mass. In this work we introduced and applied a scheme that smoothly and accurately connects the different regions, in a manner that enables every part of the star to be described by the relevant EOS at any time. There are two types of transitions that should be taken into account. The first occurs when the matter starts from having some composition of isotopes, and as it heats up during the collapse and undergoes nuclear reactions, the composition can reach NSE. The second is the transition from ideal gas to nuclear density EOS, occurring when the density rises above $\sim 10^{11}\,\text{g}\>\text{cm}^{-3}$. Previous works did not treat self-consistently both transitions. It was either that nuclear reactions were neglected (or matter instantly converted to NSE state when the temperature crossed some threshold) or that the transition to nuclear density EOS was not smooth (Herant et al., 1994).

Nuclear burning is also important in order to have a star with appropriate initial conditions that will successfully collapse in a feasible manner. Having the entire star at NSE prior to collapse is unrealistic, since the density and temperature near the edge are far below NSE conditions. Nuclear burning is required for matter originally not in NSE to have a smooth transition to NSE when the temperature is high enough. If nuclear burning is not enabled and the composition for some part of the star is fixed, the EOS will eventually fail to describe the matter when its density becomes high enough for nuclear interactions to be important ($\rho\gtrsim 10^{11}\,\text{g}\>\text{cm}^{-3}$). Additionally, electron degeneracy causes the temperature to be very sensitive to the internal energy, resulting in large fluctuations during the collapse. This issue was also handled in our work (see Appendix A.2), in a manner that has not been done in previous AIC studies.

We aim in this paper to provide a numerical calculation of an AIC (for given initial conditions and input physics) for the spherical case, where neutrinos and GR corrections are neglected. We demonstrate that the amount of ejected mass following AIC is largely affected by the EOS and the methods in which the transitions between regimes are handled. To show that, we use a 1D, Lagrangian, Newtonian hydrodynamics scheme and focus on an accurate treatment of the EOS and the nuclear reaction network in order to resolve the uncertainties associated with these parameters. These steps are necessary in order to simulate AIC, but they are not sufficient. In order to completely describe AIC, neutrino emission and transport, as well as GR effects, must be taken into account. These processes have a large effect on the collapse dynamics, considering the high densities and temperatures that take place during the process, and are being omitted at this point. Consequently, we do not aim to fully solve AIC in this work, but to provide the necessary first step before solving the full problem. Our solution for this case may serve as the foundation for adding more physical processes, such as the effect of deleptonization and neutrinos, GR corrections, rotation, and additional dimensions. Calculations were done with the VULCAN 1D hydrodynamic code (Livne, 1993), together with MESA routines for the EOS and for the nuclear reaction network (Paxton et al., 2010, 2015), modified for our purposes.

1.2 Paper Structure

The structure of the paper is as follows: In section 2, we describe the hydrodynamic scheme we use for the simulations. Section 3 describes our attempts to reproduce the results of Fryer et al. (1999). In section 4, we describe the results of simulations with our own initial profiles. In section 5, we propose a simple model to describe the outcome of the collapse for initial profiles with isentropic cores. In appendix A, we provide additional numerical details regarding the scheme, while appendix B contains a comparison with an analytic solution of the collapse for an ideal gas polytrope (Yahil, 1983).

2 Methods

We use the 1D Lagrangian hydrodynamics code VULCAN 1D (Livne, 1993) along with our modifications. One such modification is replacing the energy with the entropy as an EOS variable. The reason is that the electrons, which contribute most of the pressure, are highly degenerate, and therefore small deviations in the energy can lead to large fluctuations in the temperature. For more details, see Appendix A.2.

The EOS consists of a few terms: electron-positron plasma, radiation, nuclei, nuclear level excitations, and Coulomb corrections, where full ionization is assumed at all times. The electrons and positrons are treated with the Timmes EOS (Timmes & Arnett, 1999) at all regimes. At low densities ($\rho\ll 10^{11}\text{g}\>\text{cm}^{-3}$), nuclei are treated as an ideal gas mixture, described by the density, $\rho$, the specific entropy, $s$, and the mass fraction of each isotope, ${X_{i}}$. We use the MESA routines (Paxton et al., 2010) for the EOS of the ions at low densities and for the Coulomb corrections, with some modifications, to accurately describe the entropy (for more details, see Kushnir (2018)). Coulomb corrections are based on Chabrier & Potekhin (1998). It is important to include all isotopes with nonnegligible mass fractions at the densities and temperatures where the MESA EOS is active. We have found that a list of 183 isotopes is sufficient (for more details about the isotope selection process, see A.1).

The ions are assumed to be in NSE at high temperatures and low densities, with the composition being determined by the density, $\rho$, entropy, $s$ and electron fraction, $Y_{e}$. At this regime, and at high densities, a tabulated EOS, which describes matter in NSE and at nuclear densities, is used (see appendix A.1 for more details).

As densities approach $\sim 10^{11}\,\text{g}\,\text{cm}^{-3}$, nuclear interactions become important and the ions can no longer be treated as an ideal gas. For these regimes, we use the tables and routines provided in stellarcollapse.org/microphysics (O’Connor & Ott, 2010; Schneider et al., 2017). These tables use the nuclear EOS (Lattimer & Swesty, 1991; Shen et al., 1998, 2011b, 2011a) and, as in the NSE case, receive the triplet $(\rho,s,Y_{e})$ and return the rest of the variables. As the density approaches nuclear densities ($\sim 10^{14}\,\text{g}\,\text{cm}^{-3}$), these EOSs become stiff, although there is some variability in modeling the transition to nuclear densities (Lattimer, 2012). The routines provided by Schneider et al. (2017) allow the creation of an EOS table for matter in NSE and a table for the nuclear EOS using the single nucleus approximation (SNA; Lattimer & Swesty, 1991) and their merger in a thermodynamically consistent manner. This is done by summing the free energy density $F$ of both EOSs at a crossover region:

[TABLE]

where $\xi(\rho)$ is the volume fraction of the SNA EOS, given by

[TABLE]

the transitional density is $\rho_{t}=1.7\times 10^{11}\text{g}\>\text{cm}^{-3}$ and the dimensionless thickness is $\rho_{\delta}=0.33$. We have made some modifications in these routines to fit smoothly to the modified MESA EOS, such as the Coulomb correction term. In what follows, we add "/NSE" to the name of the high-density EOS in order to describe the resulting EOS. When the EOSs of O’Connor & Ott (2010) are used, the entire table is provided in advance with no editing option, so the transition between the MESA and the tabulated nuclear EOS is not guaranteed to be smooth.

Figure 1 shows the density-temperature trajectory of a typical mass element during AIC (taken from the profile described in Section 4.1). At the early stages of the collapse, matter travels along a low-density and low-temperature isentrope, described by the EOS where the composition is evolved using a nuclear reaction network (full EOS; solid line). When the temperature and density are high enough, nuclear burning takes place, quickly increasing the entropy and changing the composition to NSE. At this stage, the tabulated EOS is being used (dashed line). The cell then crosses the crossover region, centered at $\rho_{t}=1.7\times 10^{11}\,\text{g}\,\text{cm}^{-3}$ (vertical dotted line). Later on, the collapse halts by a shock wave, and the element expands. It is evident from the figure that a large span of densities and temperatures and various transitions between different types of EOS are required to describe AIC.

At the low-temperature and low-density regimes, nuclear burning takes place and is implemented using the MESA routines (Paxton et al., 2015), with raw reaction rates taken from the JINA reaclib database (Cyburt et al., 2010). We keep $Y_{e}$ constant in each cell during the simulations, and we ignore neutrino emission, except in Section 4.3. Neutrino physics of Section 4.3 is based on the deleptonization scheme of Liebendörfer (2005) and the thermal neutrino creation of Itoh et al. (1996). For more details, see appendix A.1 and A.3.

3 Reproducing previous results

In order to test our scheme, we compared our results with the analytic solution of Yahil (1983), which describes the collapse of a star with a polytrope EOS. This solution is quite relevant, since a complicated EOS can be approximated as polytropes under some conditions (mainly due to the degeneracy of the electrons). The simulation results deviate by less than $6\%$ from the analytic solution over a range of $10$ orders of magnitude. More details regarding the comparison to the analytic solution are given in appendix B.

We next tried to reproduce the results of Fryer et al. (1999) for their initial profile (taken from Woosley & Baron (1992)). This initial profile is taken in the midst of the collapse, where all parts of the star have already begun falling toward the center. The inner part of the star, up to an enclosed mass of $m=0.54M_{\text{tot}}$, where $M_{\text{tot}}$ is the total mass of the profile, is in NSE and has gone through some deleptonization, with $Y_{e}$ starting at $\sim 0.4$ at the center and linearly increasing with the mass until it reaches $0.5$ at the edge of the inner part. For $m>0.54M_{\text{tot}}$, the star has not gone through nuclear burning or deleptonization. There, the composition is divided equally between carbon and oxygen (CO) until $m=0.82M_{\text{tot}}$, and the outer layers of the star are composed of helium. The temperature at the boundary of the NSE and non-NSE regions sharply drops from $\sim 9\times 10^{9}\,\text{K}$ to $\sim 2\times 10^{7}\,\text{K}$, since the nuclear reactions that had taken place at the inner parts of the star contribute energy to these parts and increase their temperature. Fryer et al. (1999) ran this initial profile for $0.2\,\text{s}$, testing the sensitivity to several parameters, such as the EOS and neutrino physics.

We tried to reproduce a run where the EOS used for the high densities was the LS EOS, with an incompressibility of the bulk nuclear matter parameter of $K_{s}=180\,\text{MeV}$ (LS180 of O’Connor & Ott, 2010), and without any neutrino treatment, for which Fryer et al. (1999) had obtained an ejected mass of $0.2M_{\odot}$. There were several issues we encountered while trying to reproduce these results. The original profile used in Woosley & Baron (1992) and in Fryer et al. (1999) could not be found 111We thank Chris Fryer and Stan Woosley for their effort to locate this profile.. Nevertheless, a similar profile was kindly provided to us by Eddie Baron. This profile did not exactly match the original profile, as the central density, radius, and maximum collapse speed of the profile we were given are $4.05\times 10^{10}\,\text{g}\,\text{cm}^{-3}$, $1.1\times 10^{8}\,\text{cm}$, and $3.2\times 10^{8}\,\text{cm}\,\text{s}^{-1}$ compared to $6.53\times 10^{10}\,\text{g}\,\text{cm}^{-3}$, $0.96\times 10^{8}\,\text{cm}$ and $3.0\times 10^{8}\,\text{cm}\,\text{s}^{-1}$, respectively. Additionally, the innermost region of the part of the star that is not in NSE, from $m=0.55M_{\text{tot}}$ up to $m=0.58M_{\text{tot}}$, has negative entropy according to our MESA EOS. This is due to the degeneracy of the ions, which the MESA EOS does not handle, as it assumes that the ions are composed of an ideal gas. We have not modified this region, as the negative entropy does not interrupt the execution of the EOS. For these low densities, Fryer et al. (1999) assumed that the composition is in NSE at all times (it is unclear to us how this was implemented with an initial profile starting at a certain non-NSE composition).

We ran this profile with our scheme, using the LS EOS for the high-density regime and NSE for high temperatures and low densities, with nuclear burning taking place, and no neutrinos. The entropy in the region that is initially negative quickly rises to positive values owing to nuclear burning and small shock waves. The obtained ejected mass in our calculations is $\simeq 0.173M_{\odot}$, with this number converging for a resolution of 400 cells or higher, with the cells divided such that the resolution around the mass cut (the Lagrangian mass coordinate separating the star and the ejecta) is increased. This is within $~{}15\%$ of the results of Fryer et al. (1999), despite all the issues we encountered.

We reran this simulation with LS220 and LS375 EOSs from O’Connor & Ott (2010). While the results for the LS220 case were similar to the LS180 case, the LS375 case resulted in a slightly lower ejecta mass, $\simeq 0.16M_{\odot}$. We also reran with the LS220/NSE EOS and we obtained an ejecta mass of $\simeq 0.095M_{\odot}$, about half the amount obtained with the EOS of O’Connor & Ott (2010). These results demonstrate the large impact of the low density regime of the EOS on the ejected mass (compared with the low impact of the high-density regimes), and the need to describe it accurately.

4 A parameter study

The initial profile of a collapsing WD has many parameters that can affect the ejected mass, such as the mass, composition, temperature profile, and electron fraction. We do not aim here to determine what is the initial profile, but instead we parameterize a few simple initial profiles, similar to the profile from the previous section, in order to provide an estimate for the range of possible ejected mass during AIC.

Our initial profiles consist of a Chandrasekhar-mass star with an isentropic core in a hydrostatic equilibrium with most of its mass having an adiabatic index below $4/3$. While the entropy profile from the previous section increases from $\approx 0.8\,k_{B}/$baryon at the center up to $\approx 1.5\,k_{B}/$baryon at the edge of the core, we choose to keep the core entropy constant for simplicity. It is assumed that the inner parts of the star have gone through nuclear burning and some deleptonization and are in NSE, with the electron fraction, $Y_{e}$, starting at some value $Y_{e}(0)$ at the center and rising linearly with the mass until it is $0.5$ at about $m=0.65M_{\text{tot}}$. Constructing the profile is done inside-out, where the density, pressure, temperature, and other thermodynamic quantities are determined from hydrostatic equilibrium and the predetermined entropy and electron fraction. As the distance from the center increases, the temperature decreases until it crosses a threshold $T_{\text{thr}}=9\times 10^{9}\,\text{K}$, after which it is assumed that nuclear burning is negligible. In this region, the temperature drops to $2\times 10^{8}\,\text{K}$ and the mass fraction is equally divided between carbon and oxygen. The entropy of the outer parts also drops by a constant factor, determined by the ratio between the entropy before the temperature drops and afterward, taken with the density at the transition. The initial entropy is chosen such that the total mass of the star is the Chandrasekhar mass.

The star is driven to collapse by giving its layers initial infall velocities or reducing the pressure in the center. The total kinetic energy from the infall velocities is orders of magnitude smaller than the gravitational energy, and from the kinetic energy at bounce. Since the adiabatic index is lower than $4/3$, the collapse continues until the EOS stiffens. An external pressure is applied to the star to maintain the hydrostatic condition at the beginning of the run and is kept constant until bounce. After bounce and expansion, as the pressure of the ejecta becomes smaller as it expands, the external pressure is modified so that the expansion will continue smoothly. We verified that the application of the external pressure has no effect on the ejected mass.

4.1 A specific example

The collapse of one specific profile is shown in Figure 2 and Figure 3. The nuclear EOS used in the run was the LS220/NSE EOS, created from the routines of Schneider et al. (2017). The initial density and electron fraction at the center were $5\times 10^{10}\,\text{g}\,\text{cm}^{-3}$ and $0.4$, respectively. The uniform initial core entropy was $1.755\,k_{B}/\text{baryon}$ and the resolution was 400 cells, divided with equal radial spacing. A file containing the initial structure is provided in Zenodo doi:10.5281/zenodo.3740458. We did not include deleptonization, and the whole $Y_{e}$ profile (as a function of mass) was kept constant throughout the collapse. Panel (a) of Figure 2 shows the trajectory of chosen mass elements that represent the different trajectories, and panel (b) of Figure 2 shows the total specific energy (gravitational, kinetic, and internal) of the mass elements, taken at several snapshots (before the time of bounce, at the time of bounce, during the shock propagation, and after the shock breakout). Figure 3 shows the density, entropy, and velocity profiles at various times (the earliest time corresponds to the initial configuration while later epochs take place before the time of bounce, during the shock propagation and after the shock breakout). As the star collapses, the density increases in most of the star (the reason for the decrease at the edge is discussed below). After bounce, the densities reach nuclear densities of $\sim 10^{14}\,\text{g}\,\text{cm}^{-3},$ and the formation of the shock wave is clearly visible. The ejected mass of this run is $2.84\times 10^{-2}M_{\odot}$, which can be read from the point where the energy curve intersects with the x-axis, showing the amount of mass with positive energy. The behavior of the simulations from the previous section is qualitatively similar.

The typical energies of the ejecta are tens of MeV/baryon, with the outermost shells rising above $100\,\text{MeV/baryon}$. These energies correspond to mildly relativistic velocities of about $0.15c-0.3c$ for most of the ejected mass. The ejecta is composed mainly of ${}^{56}\text{Ni}$, since it is the most abundant isotope of matter in NSE at low temperatures and $Y_{e}=0.5.$ As seen in the latest epoch in Figure 3, the ejecta density profile is exponential with respect to the radius (and to the velocity, since the velocity profile is homologous) to a very good approximation.

Some mass elements started moving outward before bounce, seen at $t\simeq-10\,\text{ms}$ in Figure 2a and from the decrease in the density near the edge of the star at $t\simeq-4\,\text{ms}$ in Figure 3. This is due to the energy gained by nuclear burning, causing these elements to reverse their motion before bounce. After bounce, they are quickly caught up by the emerging shock, which moves at much higher velocities (note that the energy gained by the difference of the gravitational binding energy is tens of MeV/baryon, an order of magnitude higher than the typical nuclear binding energy of $\sim\text{MeV/baryon}$).

In order to check convergence, we examined additional scenarios, where resolution was altered while the rest of the parameters remained the same. The resolution increased up to 988 cells that were divided in the same manner. Figure 4 shows the convergence of the run, by plotting, as a function of the resolution, the divergence of the ejected mass from that of the run with the 988 cells (where the ejected mass was $0.0285M_{\odot}$). For resolutions higher than 200 cells, the the deviation of the calculated ejected mass from that of the highest resolution does not exceed $3\times 10^{-4}M_{\odot}$, reflecting a $1\%$ error.

4.2 Varying the EOS and initial conditions

In order to assess the effects of the EOS and initial conditions, we calculated the ejected mass for several nuclear EOSs and initial uniform core entropy values. The ejected mass in these runs is provided in Figure 5. The LS220/NSE and SKRA/NSE EOSs are created from the routines of Schneider et al. (2017) with different nuclear parameters (Lattimer & Swesty, 1991; Rashdan, 2000), and the EOSs of HShen (Shen et al., 1998) and GShen (Shen et al., 2011a, with the NL3 parameterization) are taken from O’Connor & Ott (2010). The LS180 EOS used in the previous section is not part of the list of EOSs used here, because its parameters, and specifically the incompressibility of the bulk nuclear matter being $K_{s}=180\,\text{MeV}$, are currently not favored (Lattimer, 2012). The LS220 of O’Connor & Ott (2010) is also shown, for comparison, and has larger ejecta masses than the LS220/NSE, as in Section 3.

The ejected masses in Figure 5 are plotted as a function of the uniform initial core entropy, all starting with a central density of $5\times 10^{10}\,\text{g}\text{cm}^{-3}$. Note that for each value of the initial core entropy there is a value for the central electron fraction $Y_{e,c}$ that ensures that the mass of the star remains the Chandrasekhar mass. Since the pressure rises with the entropy and with $Y_{e,c}$, the electron fraction should be decreased for a larger initial core entropy in order to keep the mass of the star unchanged. As in the previous section, deleptonization is not included. The ejected mass in all runs is a few$\times 10^{-2}M_{\odot}$. The EOS and initial core entropy alter the ejected mass by a factor of a few but do not change it by orders of magnitude. The EOSs with the low-density treatment of Schneider et al. (2017) all result in lower ejected mass, showing the importance of treating this region correctly. Yasin et al. (2020) argued that the nucleon effective mass in the high-density part of the EOS, $\rho\simeq 3\times 10^{14}\,\text{g}\,\text{cm}^{-3}$, plays an important role in collapsing stars, with larger effective masses leading to a faster shock evolution and explosion. This might also be important when determining the ejecta mass in AIC studies (note that neutrinos, which are neglected here, play a large role in the work of Yasin et al., 2020). Comparing EOSs with the same low-density treatment, we find that the large effective mass of the LS220/NSE EOS has a higher ejected mass than the low nucleon effective mass EOS of SKRA/NSE. Comparing the LS220 to the HShen EOS and the GShen EOS, which also have a low effective mass, seem to agree with this only at high specific entropies.

The initial profile also has a considerable effect, as the results of this section are smaller by a factor of a few than the results in section 3, where the ejected mass is $\simeq 0.095M_{\odot}$ for the LS220/NSE EOS. To explain this discrepancy, we attempted to build a hydrostatic initial profile with the same entropy structure as in the initial profile in Section 3 (with non-isentropic core), where the negative entropy region was replaced by a low entropy value. The total mass in the initial profile we built, using the LS220/NSE EOS, was $\sim 0.05M_{\odot}$ less than the Chandrasekhar mass of the original profile, with this mass deficiency affecting the density profile at the outer parts of the star. The cause for this discrepancy is the hydrostatic equilibrium requirement that does not hold for the original profile. The ejected mass of our own hydrostatic profile was about $0.025M_{\odot}$, similar to the results in Figure 5 for the LS220/NSE EOS.

We have also created a Chandrasekhar mass version of the same profile, by slightly altering the entropy structure at the core. The extra mass was mainly added to the core and not to the outer layers. The ejected mass in this case is $\approx 0.032M_{\odot}$, compared to $\approx 0.095M_{\odot}$ from the original profile (see Section 3), despite the same progenitor mass of both profiles. This illustrates that the original profile’s excess mass of $0.05M_{\odot}$ in its outer layers, which cannot be in hydrostatic equilibrium in our profiles, is the main cause of the large ejected mass of this profile.

Calculations without nuclear burning, where the entire star is assumed to begin in NSE, are described in Section 5 and the obtained ejecta masses are shown in Figure 7. In this case, much higher ejecta masses were obtained, further demonstrating the importance of accurately describing the low-density regime and the outer layers of the progenitor star.

4.3 Effects of neutrinos

Neutrino creation and emission, which are not taken into account in the previous simulations, will lead to energy and entropy losses, making it more difficult for the mass shells to escape from the star. On the other hand, the neutrinos can deposit some of the energy back when interacting with the outer layers, before escaping the star. To grasp the effect of the neutrinos, we tested the case where the neutrinos only carry energy away without any additional interaction (see appendix A.3). We stress that since this is not the main focus of the paper, we only use a crude approximation for the effect of neutrinos. Detailed calculations are required to fully investigate the effects of neutrinos. The addition of neutrino emission, mainly due to deleptonization, leads to a reduction in the ejected mass by roughly an order of magnitude, to $\sim 10^{-3}M_{\odot}$, for the same initial profiles of Section 4.2. The simulations include $\sim 300$ cells, where the outer cell mass resolution is $\sim 10^{-4}\,M_{\odot}$. This can be considered a lower limit for these profiles, EOSs and input physics.

5 A simple estimate of the ejected mass

We now construct a simple model to describe the structure of the PNS at the end of the collapse and to evaluate the ejected mass. We assume that the tabulated EOS, which describes nuclear matter for high densities and ideal gas in NSE for low densities, can be used for the entire star. We therefore assume that the matter is in NSE and that there is no burning, and we neglect neutrinos. Another assumption is, as in the previous section, that the initial profile is isentropic (throughout the entire star in this case). We assume that the flow is isentropic throughout the collapse, until the matter is shock heated by the shock that was created from the bounce of the infalling matter, which alters the entropy profile, resulting in hydrostatic equilibrium. The model requires two unknown parameters - the central density of the final hydrostatic star, $\rho_{\text{static}}$, and the (larger) central density at the time of bounce, $\rho_{\text{bounce}}$, which is used in order to estimate the shock strength. The two parameters, together with the hydrostatic equilibrium and isentropic flow conditions above, are sufficient to construct the PNS. The two parameters are retrieved from the numerical simulations in the previous section. Figure 6 shows the central density as a function of time for a typical run, and how we extracted the density of the hydrostatic star and at bounce that we used in the model. The ejected mass can be calculated from the difference in mass between the initial WD and the PNS.

To begin constructing a hydrostatic star, the central density, $\rho_{\text{static}}$, and entropy (or any other two thermodynamic variables, given an EOS) should be predetermined. The latter is determined from the initial entropy, while the former, as mentioned above, is extracted from the simulation result. The final central density can be approximately estimated as the nuclear density, but the ejected mass was found to be quite sensitive to this value, changing by $10-30\%$ as a result of a $5\%$ change in the density. The density-pressure profile is then determined from the hydrostatic equilibrium,

[TABLE]

and the shock heating. To estimate the shock heating, we assume that the infalling matter bounced and immediately stopped; hence, the strength of the shock will be determined out of the infall velocity, $u$, through the Hugoniot conditions,

[TABLE]

where $V=1/\rho$ and the 0 and 1 subscripts refer to the upstream and downstream, respectively. Once the infall velocity is known, it is possible to retrieve all the variables using the conditions mentioned above.

To estimate the infall velocity, we make use of the self-similar solution discussed in Appendix B. The velocity is given by

[TABLE]

where $\kappa=p/\rho^{\gamma}$; $X$ and $U$ are the known self-similar coordinate and velocity, respectively; and $t_{c}-t_{\text{bounce}}$ is the time from bounce (for the self-similar solution, the density diverges at $t_{c}$). The adiabatic index $\gamma$ and the constant $\kappa$ are determined from the EOS at the given initial entropy and an averaged value of $Y_{e}$. The quantity $t_{c}-t_{\text{bounce}}$ is related to the second unknown parameter - the central density at bounce – by using the (known) self-similar density:

[TABLE]

In such a way, the velocity can be obtained given the radius, $r$, and the PNS can be constructed.

We note that the quantity $t_{c}-t_{\text{bounce}}$ can be roughly estimated by comparing the final central density (the first parameter) to the self-similar density using the same relation, Equation (9), but this gives an overestimate of $t_{c}-t$, since the density of the end state is not as high as it was at the time of the bounce and shock formation. This will result in a weaker shock and lower temperatures and pressures, which can support less mass, leading to a higher ejected mass. The estimation of $t_{c}-t_{\text{bounce}}$ from $\rho_{\text{bounce}}$ is not entirely optimal but gives the correct ejected mass to within $15\%$ for the LS220 and the LS220/NSE EOSs, where in some cases the difference is less than $5\%$. When the initial entropy decreases and the ejected mass becomes low ($\lesssim 0.035M_{\odot}$), some of the matter being ejected by the shock falls back into the star. When it meets the star, another shock wave is created, which further increases the entropy of this matter. This effect is not captured in our model, so the estimates for the runs with low initial entropy give a lower mass compared to the simulations.

Figure 7 presents the simulation results and the model estimations for profiles with the LS220/NSE and HShen EOSs and an initial central density of $5\times 10^{10}\,\text{g}\,\text{cm}^{-3}$. The model results for the HShen EOS are not as accurate as for the LS220 EOS but provide the correct results to within $50\%$. Since the model calculates the ejecta from the mass difference between the two stars, it is valid only because the ejecta is not a very small fraction of the star, about $2-5\%$. For an ejecta of $5\%$ of the star, an error of $1\%$ of the mass of the PNS would result in an error of $20\%$ and will increase as the ejecta becomes a smaller fraction. If the range of ejected mass would have been an order of magnitude smaller, the accuracy of the model would have to be high enough to calculate the mass of the PNS within $<0.1\%$ in order to have reasonable results.

Figure 8 compares the energy and entropy profiles of the simulation to the ones derived with the simple model for a run with the LS220/NSE EOS, starting with an initial central density of $5\times 10^{10}\,\text{g}\,\text{cm}^{-3}$, uniform entropy of $1.45\,k_{B}/\text{baryon}$ and $Y_{e}(0)=0.4$. The mass with positive energy is being ejected. The entropy change shows how well the shock strength is described by the model, since it is independent of the hydrostatic equilibrium constraint, unlike the other parameters, such as density and pressure. The ejected mass of this run is $0.0517\,M_{\odot}$, while the model gives a mass of $0.0537M_{\odot}$, a difference of $\sim 4\%$.

The deviation in the density and entropy profiles between the model and the simulation is less than $2\%$ in density and less than $20\%$ in entropy, except for the model’s last part, $10^{-3}M_{\odot}$ from the edge, where the entropy quickly rises.

6 Conclusions

In this work, we implemented a numerical scheme to calculate the amount of ejected mass following AIC. We assumed a spherical collapse and neglected neutrinos and GR corrections, while accurately treating the EOS and nuclear burning. We found the amount of ejected mass to be a few $\times 10^{-2}M_{\odot}$ for a large range of initial conditions and EOSs, and, as expected, to always move at mildly relativistic velocities. The low-density regime of the EOS was found to have a nonnegligible impact on the EOS, as seen in figures 5 and 7, which show that all regimes in the EOS should be accurately treated. We have suggested a simple model that estimates the ejected mass at the end of the AIC, for the LS EOS, to within $15\%$ given two free parameters, while describing the PNS structure to reasonable accuracy ($<2\%$ in density and $<20\%$ in entropy).

Our numerical scheme can serve as a basis for future studies, where more physical aspects can be taken into account and integrated into the scheme. Such processes were already studied in the context of AIC, and include GR (Fryer et al., 1999), neutrino physics (Woosley & Baron, 1992; Fryer et al., 1999; Dessart et al., 2006) and rotation (Dessart et al., 2006). Marek et al. (2006) developed a method to include GR corrected potential for a Newtonian hydrodynamic code, which could be integrated in an AIC simulation as well. We believe that our work is essential for any study that aims to accurately describe the AIC process.

Waxman (2017) suggested that the persistent radio source, associated with FRB 121102, can be explained by the propagation of a spherical shell of plasma into a surrounding medium. The properties of the persistent source imply that the mass of the shell is $\sim 10^{-5}M_{\odot}$ moving at mildly relativistic velocities, for a duration of $t<10^{2.5}\,\text{yr}$. This suggests that the source was created by a weak stellar explosion, such as AIC, since the mass of the shell is much smaller than the ejecta mass from a typical supernova. The obtained ejected mass of $\sim 10^{-2}M_{\odot}$ is larger than the estimated ejected mass shell of the FRB 121102 persistent source ($\sim 10^{-5}M_{\odot}$). This inconsistency could be a result of the neglected physical processes, such as neutrino physics, GR corrections and rotation, or nonrepresentative initial conditions. A preliminary neutrino emission calculation did reduce the ejecta mass by an order of magnitude, but not down to $\sim 10^{-5}M_{\odot}$. Future work could determine more reliably the connection between the persistent source of FRB 121102 and the AIC process.

Acknowledgements

We thank Eddie Baron, Boaz Katz, Eran Ofek, and Eli Waxman for useful discussions. D.K. is supported by the Israel Atomic Energy Commission, the Council for Higher Education, the Pazi Foundation, and by a research grant from The Abramson Family Center for Young Scientists.

Appendix A Numerical details

A.1 Nuclear reaction network and NSE

Nuclear reactions were calculated with the MESA routines (Paxton et al., 2015). The forward reaction rates are taken from the JINA reaclib database (Cyburt et al., 2010). The rates in JINA are only valid up to $10^{10}\,\text{K}$, which does not cover the temperature range in our simulations, even for matter that is not in NSE. We therefore used the rate at $10^{10}\,\text{K}$ for higher temperatures. Backward rates are calculated from detailed balance. We allow the EOS of each numerical cell to change from the MESA routines for ideal gas to NSE/nuclear, and vice versa. When the cell is at the MESA EOS, isotopes approach NSE owing to nuclear burning. When the temperature is higher than $5\times 10^{9}\,\text{K}$ and the mass fraction $X_{i}$ of every isotope does not differ by more than 0.01 from the NSE composition, the EOS is replaced. Since the EOSs are not always identical, we kept the pressure and density of the previous EOS, while the other quantities (temperature, entropy, etc..) are set by the new EOS. The pressure and density are chosen so that the force applied by the cell and its mass remain unchanged. The NSE composition is calculated for 3335 isotopes in order to accurately describe the EOS over a wide range of $T,\rho$, and $Y_{e}$. In order to have a smooth transition when changing the EOS, it is required that the MESA EOS will contain all isotopes whose NSE composition at transition has nonnegligible mass fractions. The list of these isotopes was prepared in the same method as in Kushnir (2018), by observing the isotopes with a molar mass fraction, $Y_{i}=X_{i}/A_{i}$ larger than $10^{-5}$ at the range of parameters for which the MESA EOS is relevant, $T\in[2\times 10^{9},3\times 10^{10}]\,\text{K}$, $\rho\in[100,10\times 10^{10}]\,\text{g}\>\text{cm}^{-3}$ and $Y_{e}\in[0.495,0.5]$. We have found that list NSE5 (with 179) from Kushnir (2018) together with four additional isotopes (34P, 35P, 62Ga, 65Ge), is sufficient. As the temperature of a cell in NSE drops below $5\times 10^{9}\,\text{K}$, the EOS changes to the MESA EOS, with the composition being the NSE composition at the time of the transition. The isotope list, chemical data, and reaction rates can be found on Zenodo doi: 10.5281/zenodo.3740458

A.2 Entropy as an EOS variable

At the large densities reached during the simulations, matter is highly degenerate and the temperature is a very sensitive function of the internal energy. Since the temperature is used in all our EOS as one of the input variables, it should be accurately obtained. The dependence of the specific entropy $s$ on the temperature is not a very sensitive function, and thus the entropy, instead of the energy, is used to evolve the EOS quantities, using the second law of thermodynamics:

[TABLE]

with $Q$ being the heat flow (per unit mass) in or out of the system. The heating can be a result of shock, neutrino emission, and nuclear burning. Note that when the composition changes owing to nuclear burning, the heat added is the sum of the energy increase due to rest mass difference and chemical potential times the composition difference:

[TABLE]

where $\delta Q_{\text{nuc}}$ is the added heat, $de$ is the added energy, and $\mu_{i}$ and $dy_{i}$ are the chemical potential and the change in the molar fraction of the i-th isotope, respectively. In NSE, matter reaches thermal equilibrium with respect to nuclear reactions, and these terms cancel each other out, such that the entropy no longer changes owing to nuclear reactions. This equation, together with the equation of motion and the EOS, is solved in the hydrodynamic code. The equation of motion is

[TABLE]

where $r$ is the radial coordinate, $u$ is the velocity, $p$ is the pressure, $\rho$ is the density, and $m$ is the mass up to $r$. This is done in two steps within the hydrodynamic code in a leapfrog integration scheme, where the velocities and positions are calculated at interleaved time points.

A.3 Neutrino physics

Although neutrino physics has an important role at the large densities and temperatures that take place in the AIC process, it is only partly considered in this work. In order to give a lower bound to the ejected mass, we omit the interaction of the neutrinos after creation and assume they freely escape the star, reducing the energy and entropy as a result. We consider two processes - deleptonization due to electron capture and creation of thermal neutrino-anti-neutrino pairs. To account for deleptonization, we use the prescription of Liebendörfer (2005), in which the electron fraction $Y_{e}$ is determined from the instantaneous density alone. This is only a crude approximation, as the prescription was calibrated from a full calculation that depends on an EOS and an initial profile that are not the ones used in our simulations. Also, this prescription is relevant only until bounce, so we let the $Y_{e}$ profile freeze afterward. Entropy loss is calculated by (Bethe, 1990)

[TABLE]

where $\mu_{\{e,n,p\}}$ is the chemical potential of the electrons, neutrons, and protons, respectively.

Creation and emission of thermal neutrinos are calculated from expressions $\{2.1,\,3.1,\,4.1,\,5.1,\,5.17,\,5.27,\,5.29,\,6.5\}$ of Itoh et al. (1996), taking into account all the processes that contribute to the specific energy loss rate, which are given as a function of the density, temperature, mean number of nucleons $\bar{A}$, and mean charge $\bar{Z}$. Since we assume that the composition does not change in this process, the energy loss is entirely converted to heat loss. This process, however, is quite negligible when deleptonization occurs.

Appendix B Comparison to analytic solution

In order to test our numerical scheme, we compared it with the self-similar analytic solution of Yahil (1983) for the collapse of a spherically symmetric star with a polytropic EOS:

[TABLE]

for the adiabatic index $\nicefrac{{6}}{{5}}\leq\gamma\leq\nicefrac{{4}}{{3}}$, where $p$ is the pressure, $\rho$ is the density, and $\kappa$ is a constant in space and time. The solution uses the self-similar dimensionless coordinate

[TABLE]

The interval from the catastrophe time $t_{c}-t$ determines the scaling of the physical variables (coordinate, velocity, etc.) to the self-similar ones. This interval is unknown during the simulation and is computed by best fitting to the analytic solution. Figure B.1 shows the velocity and density during the collapse of a star with an ideal gas EOS with $\gamma=1.3$, normalized to the self-similar variables, done with a resolution of 2000 cells, initially divided with equal radial spacing. The self-similar solution is plotted on top for comparison. The agreement is within $6\%$ over a range of 10 orders of magnitude. A convergence test we have performed shows that the deviation between the simulation and the analytic solution is (roughly) inversely proportional to the resolution.