A Consistent Modeling of Neutrino-driven Wind with Accretion Flow onto a Protoneutron Star and Its Implications for ⁵⁶Ni Production

In this study, we use the spherically symmetric and general relativistic semi-analytic wind model in Wanajo (2013). We followed Thompson et al. (2001) for the detailed calculation method. Previous works studied the physical state of the neutrino wind for a wide range of neutron star masses and neutrino luminosities (e.g., Otsuki et al. 2000; Bliss et al. 2018). Based on these results, they then studied the behavior of the neutrino-driven wind and the PNS evolution was studied by stitching this semi-analytic wind model (e.g., Wanajo et al. 2001; Wanajo 2013).

The basic equations to describe the spherically symmetric and steady-state winds in the Schwarzschild geometry are

$$\begin{eqnarray}&&\dot{M}=4\pi {r}^{2}\rho v,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&v\displaystyle \frac{{dv}}{{dr}}=-\displaystyle \frac{1+{(v/c)}^{2}-(2{GM}/{c}^{2}r)}{\rho (1+\epsilon /{c}^{2})+P/{c}^{2}}\displaystyle \frac{{dP}}{{dr}}-\displaystyle \frac{{{GM}}_{\mathrm{PNS}}}{{r}^{2}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\dot{Q}=v\left(\displaystyle \frac{d\epsilon }{{dr}}-\displaystyle \frac{P}{{\rho }^{2}}\displaystyle \frac{d\rho }{{dr}}\right)\,,\end{eqnarray} \tag{ 3 }$$

where $\dot{M}$ is the constant mass outflow rate, r is the distance from the center of the protoneutron star, $\dot{Q}$ is the heating rate, ρ is the baryon mass density, v is the radial velocity of the wind, P is the pressure, and is the specific internal energy. The system of Equations (1)–(3) is closed with the Helmholtz equation of state (Timmes & Swesty 2000), which describes the stellar plasma as a mixture of arbitrarily degenerate and relativistic electrons and positrons, blackbody radiation, and ideal Boltzmann gases of a defined set of fully ionized nuclei, taking into account corrections for the Coulomb effects. The source term $\dot{Q}$ includes both heating and cooling by neutrino interactions. Heating is due to the following three processes; i) neutrino and antineutrino captured by free nucleons, ii) neutrino scattering on electrons and positrons, and iii) neutrino–antineutrino pair annihilation into electron–positron pairs. Cooling is due to electron and positron capture by free nucleons, and annihilation of electron–positron pairs into neutrino–antineutrino pairs (for more details, see Equations (8)–(16) in Otsuki et al. 2000).

To determine the luminosity of each type of neutrino L_ν, we use the assumptions of $\dot{{Y}_{e}}=0$ (Bliss et al. 2018), in which we assume electron/positron captures are in equilibrium and an initial composition consists mainly of neutrons and protons. Then, we get

$$\begin{eqnarray}&&{Y}_{e}={\left[1+\displaystyle \frac{{L}_{{\overline{\nu }}_{e}}^{n}\langle {\sigma }_{{\overline{\nu }}_{e}p}\rangle }{{L}_{{\nu }_{e}}^{n}\langle {\sigma }_{{\nu }_{e}n}\rangle }\right]}^{-1}=0.5,\end{eqnarray} \tag{ 4 }$$

where ${L}_{\nu }^{n}={L}_{\nu }/\langle {E}_{\nu }\rangle$ is the number luminosity and is assumed to be the same for electron neutrinos and antineutrinos. The cross-sections of electron neutrino absorption at neutrons ($\langle {\sigma }_{{\nu }_{e}n}\rangle$) and electron antineutrino absorption at protons ($\langle {\sigma }_{{\overline{\nu }}_{e}p}\rangle$) depend on the average neutrino and antineutrino energies. Thus, given ${L}_{{\nu }_{e}}$, $\langle {E}_{{\nu }_{e}}\rangle$ and Y_e as parameters, Equation (4) and the assumption of ${L}_{{\nu }_{e}}^{n}={L}_{{\overline{\nu }}_{e}}^{n}$ leads to the antineutrino energy and luminosity. With respect to the choice of Y_e, we are interested in an environment that maximizes the production of ⁵⁶Ni. Previous studies have shown that the abundance of ⁵⁶Ni is extremely suppressed at Y_e < 0.5 in nuclear statistical equilibrium (NSE) (Seitenzahl et al. 2008), and the same behavior is known to occur in the neutrino-driven wind environment (Wanajo et al. 2018). In order to focus on maximizing ⁵⁶Ni production, we here fix Y_e value to be 0.5 throughout this paper.

Regarding the inner boundary condition, we assume, for simplicity, that the neutrino-driven wind starts at the gain radius R_gain, where the heating and cooling due to neutrino interactions are in equilibrium ($\dot{Q}\approx 0$). It is because the wind blows from the heating region outside the gain radius. It allows us to determine the temperature of the inner boundary, T₀. While the inner boundary is set to be ρ₀ = 10¹⁰ g cm⁻³ in the ordinary neutrino-driven wind models (e.g., Otsuki et al. 2000; Wanajo et al. 2001; Wanajo 2013), we set the density at the inner boundary in our model to ${\rho }_{0}={10}^{10}\,{\rm{g}}\,{\mathrm{cm}}^{-3}{L}_{\nu ,51}^{1/2}$, following the method discussed by Fujibayashi et al. (2015), in order to solve with the very high neutrino luminosity. Given the density of the inner boundary ρ₀, Equation (1) determines the initial velocity v₀ at r = R_gain for each $\dot{M}$.

The mass outflow rate $\dot{M}$ corresponds physically to determining how much material is ejected by the neutrino-driven wind, and algebraically to giving the behavior of the velocity in solving the equations. In other words, the solutions of Equations (1)–(3) depend on this mass outflow rate $\dot{M}$ (Qian & Woosley 1996).

For the case of small $\dot{M}$, the wind moves slowly, and dv/dt can become zero before the velocity v reaches the sonic speed v_s. When this occurs, the velocity decreases after reaching a maximum value and becomes always subsonic (breeze or subsonic solution). When $\dot{M}$ increases, it can happen that the acceleration term becomes zero at the same time as v reaches v_s. This case corresponds to a critical value $\dot{M}={\dot{M}}_{\mathrm{tran}}$ (transonic solutions). In this critical case, the velocity increases through the sound speed to supersonic values, eventually becoming a constant when almost all the internal energy is converted into mechanical kinetic energy. Mass outflow rates larger than ${\dot{M}}_{\mathrm{tran}}$ are unphysical because, for these values of $\dot{M}$, v has reached v_s when the acceleration term is still positive, resulting in an infinite acceleration. Thus, we only need to focus our attention on cases of $\dot{M}\,\leqq \,{\dot{M}}_{\mathrm{tran}}$. Here we briefly review two particular cases of the neutrino-driven wind with different mass outflow rates $\dot{M}$, which give transonic and subsonic solutions. Figure 2 shows the fluid velocity and temperature as a function of the radius from the center of the PNS for each $\dot{M}$, where the PNS mass M_PNS = 1.4 M_⊙, the gain radius R_gain = 40 km and the neutrino luminosity L_ν = 10⁵² erg s⁻¹. In these PNS parameters (L_ν, R_gain, M_PNS), the transonic and subsonic solutions correspond to the mass outflow rate $\dot{M}={\dot{M}}_{\mathrm{tran}}$ = 8.25 × 10⁻³ M_⊙ s⁻¹ and $\dot{M}=8.00\times {10}^{-3}\,{M}_{\odot }\,{{\rm{s}}}^{-1}$, respectively.

Zoom In Zoom Out Reset image size

Note that $\dot{M}={\dot{M}}_{\mathrm{tran}}$ gives the maximum mass ejection. In the following, we focus our discussion on transonic solutions which give the maximum ejected amount of the wind.

The system of Equations (1)–(3) provides a (transonic) wind solution ${\dot{M}}_{\mathrm{tran}}$ for each set of three PNS parameters (L_ν, R_gain, M_PNS). In the following, we consider modeling the relation between the maximum mass ejection of the wind ${\dot{M}}_{\mathrm{wind},\mathrm{iso}}={\dot{M}}_{\mathrm{tran}}$ and the three PNS parameters. Figure 3 illustrates mass ejection rates ${\dot{M}}_{\mathrm{wind},\mathrm{iso}}$ as a function of each PNS parameter. We confirmed that our results are in close agreement with the results of the previous 1D numerical (Sumiyoshi et al. 2000) and semianalytical studies (Otsuki et al. 2000; Wanajo et al. 2001) with the same set of parameters, except for the treatment of Y_e.

Zoom In Zoom Out Reset image size

It is found that ejection rates ${\dot{M}}_{\mathrm{wind},\mathrm{iso}}$ have a relation with the PNS parameters, which is approximated by a power-law function as

$$\begin{eqnarray}&&\begin{array}{l}{\dot{M}}_{\mathrm{wind},\mathrm{iso}}\approx 8.3\times {10}^{-3}\,{M}_{\odot }\,{{\rm{s}}}^{-1}\\ \ \times \,{\left(\displaystyle \frac{{L}_{{\nu }_{e}}}{{10}^{52}\mathrm{erg}{{\rm{s}}}^{-1}}\right)}^{\alpha }{\left(\displaystyle \frac{{R}_{\mathrm{gain}}}{4\times {10}^{6}\mathrm{cm}}\right)}^{\beta }{\left(\displaystyle \frac{{M}_{\mathrm{PNS}}}{1.4{M}_{\odot }}\right)}^{\gamma },\end{array}\end{eqnarray} \tag{ 5 }$$

where the index we adopted is as follows; α = 7/4, β = 5/2 and γ = − 7/2. To credit this modeling, we should mention two points in comparison to previous studies. First, our model focuses on a mainly larger region for the gain radius R_gain (10–60 km), compared to previous studies (10–30 km; Qian & Woosley 1996; Otsuki et al. 2000). We confirmed that almost the same relation can be found as in Qian & Woosley (1996) when we focus on the same parameter region of R_gain as the literature. Second, while the result of Wanajo et al. (2001) tends to deviate from the power-law relation at high luminosity, our solutions do not. It is due to our treatment of the density at the inner boundary (gain radius) that follows the method of Fujibayashi et al. (2015). This method takes into account the physical equilibrium conditions, especially at high neutrino luminosity, which is different from those of the Wanajo et al. (2001). We also confirmed the same trend as in Wanajo et al. (2001) when using the same boundary conditions.

Since we employ this power-law relation to construct our wind model in a later section, here we discuss the error between this equation and the numerical solution. Table 1 shows the wind solution from the numerical calculations (${\dot{M}}_{\mathrm{wind},\mathrm{cal}}$), the value estimated from the power-law relation (${\dot{M}}_{\mathrm{wind},\mathrm{model}}$), and the error of ${\dot{M}}_{\mathrm{wind},\mathrm{cal}}$ with respect to ${\dot{M}}_{\mathrm{wind},\mathrm{model}}$, for typical values of each PNS parameter. The error values in Table 1 indicate that our subsequent analytical discussion using the power-law relation keeps an error within 30% of the more precise numerical calculations.

Table 1.  The wind Solution from the Numerical Calculations (${\dot{M}}_{\mathrm{wind},\mathrm{cal}}$), the value Estimated from the Power-law Relation (${\dot{M}}_{\mathrm{wind},\mathrm{model}}$), and the Error of ${\dot{M}}_{\mathrm{wind},\mathrm{cal}}$ with Respect to ${\dot{M}}_{\mathrm{wind},\mathrm{model}}$, for Typical values of each PNS Parameter

${L}_{{\nu }_{e}}$	Rgain	MPNS	${\dot{M}}_{\mathrm{wind},\mathrm{cal}}$	${\dot{M}}_{\mathrm{wind},\mathrm{model}}$	Errora
(1051erg s−1)	(km)	(M⊙)	(M⊙ s−1)	(M⊙ s−1)	(%)
10	40	1.4	8.25 × 10−3	8.25 × 10−3	0.0
100	40	1.4	5.98 × 10−1	4.64 × 10−1	28.9
1	40	1.4	1.23 × 10−4	1.47 × 10−4	−16.3
10	50	1.4	1.79 × 10−2	1.44 × 10−3	24.4
10	10	1.4	3.00 × 10−4	2.58 × 10−3	16.2
10	40	2.0	2.40 × 10−3	2.37 × 10−3	1.3
10	40	1.2	1.54 × 10−2	1.42 × 10−2	8.5

Note.

^aWe denote the difference between the value of the numerical solution to the model value, normalized by the model value $\left(\tfrac{{\dot{M}}_{\mathrm{wind},\mathrm{cal}}-{\dot{M}}_{\mathrm{wind},\mathrm{model}}}{{\dot{M}}_{\mathrm{wind},\mathrm{model}}}\right)$, as a percentage.

Download table as:  ASCIITypeset image

In this section, we construct a phenomenological mass accretion model based on Müller et al. (2016). We assume that matter reaches the PNS with a freefall timescale. The isotropic mass accretion rate ${\dot{M}}_{\mathrm{acc},\mathrm{iso}}$ is thus related to the mass coordinate M of the infalling shell as (Woosley & Heger 2015),

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{acc},\mathrm{iso}}=\displaystyle \frac{2M}{{t}_{f}}\,\displaystyle \frac{\rho }{\overline{\rho }-\rho },\end{eqnarray} \tag{ 6 }$$

where ρ is the initial density of a given mass shell, $\overline{\rho }$ is the average density inside a given mass shell located at an initial radius r (i.e., $M=4\pi \overline{\rho }{r}^{3}/3$ ), and t_f is the infall time, which is defined as a function of the average density $\overline{\rho }$ inside a given mass shell by

$$\begin{eqnarray}&&{t}_{f}=\sqrt{\displaystyle \frac{\pi }{4G\overline{\rho }}}.\end{eqnarray} \tag{ 7 }$$

Figure 4 shows the mass accretion rates ${\dot{M}}_{\mathrm{acc},\mathrm{iso}}$ as a function of t, which is identified with t_f. In our model, the time origin is when the neutrino-driven wind starts to blow, corresponding to the time of the SN shock revival, which is given by the time at the mass shell of ${M}_{s=4{k}_{B}}$ accreting onto the PNS. Here ${M}_{s=4{k}_{B}}$ gives the mass coordinate with the entropy being 4k_B baryon⁻¹. This is because recent hydrodynamics simulations show that the shock launch takes place when a mass element with s = 4k_B baryon⁻¹ accretes onto the shock (Ertl et al. 2016; Suwa et al. 2016). Gray lines are the mass accretion rates of progenitor stars from Sukhbold et al. (2018) (on 0.1 M_⊙ steps over the range M_ZAMS = 12.0–20.0 M_⊙). We approximate our accretion model as the red line, which is given by

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{acc},\mathrm{iso}}(t)={\dot{M}}_{\mathrm{acc},0}{\left(\displaystyle \frac{t}{{t}_{0}}+1\right)}^{-2},\end{eqnarray} \tag{ 8 }$$

where ${\dot{M}}_{\mathrm{acc},0}$ and t₀ are free parameters.

Zoom In Zoom Out Reset image size

Our description above is based on a one-dimensional radial flow. In order to construct a consistent model of the wind with mass accretion, it is necessary to take into account the geometric structure. Hereafter, we present the collimation of the wind and the accretion flow due to the asymmetric structure of the supernova explosion with a geometrical factor f_Ω (≤1). f_Ω relates the wind intrinsic properties and isotropic equivalents as

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{wind}}={f}_{{\rm{\Omega }}}{\dot{M}}_{\mathrm{wind},\mathrm{iso}},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{acc}}=(1-{f}_{{\rm{\Omega }}}){\dot{M}}_{\mathrm{acc},\mathrm{iso}}.\end{eqnarray} \tag{ 10 }$$

We next rewrite the outflow rate ${\dot{M}}_{\mathrm{wind}}$, which is obtained for a given set of the three parameters (L_ν, R_gain, M_PNS) in Section 2.1, to the equation which is approximately determined by a given set of the two parameters (${\dot{M}}_{\mathrm{acc}}$, M_PNS). We need the model of the gain radius R_gain and the neutrino luminosity L_ν as a function of ${\dot{M}}_{\mathrm{acc}}$ and M_PNS. Müller et al. (2016) found that, in the condition of ${\dot{M}}_{\mathrm{acc}}\gg {10}^{-3}\,{M}_{\odot }\,{{\rm{s}}}^{-1}$, the gain radius R_gain can be described as

$$\begin{eqnarray}&&{R}_{\mathrm{gain}}\approx 40\,\mathrm{km}\,{\left(\displaystyle \frac{{\dot{M}}_{\mathrm{acc}}}{0.1{M}_{\odot }{{\rm{s}}}^{-1}}\right)}^{1/3}{\left(\displaystyle \frac{{M}_{\mathrm{PNS}}}{1.4{M}_{\odot }}\right)}^{-1}\,.\end{eqnarray} \tag{ 11 }$$

This approximation has been confirmed to show reasonably consistent results with the contraction of the PNS in hydrodynamics simulations.

In our system, the neutrino luminosity L_ν is assumed to be dominated by accretion luminosity L_ν,acc. The accretion luminosity L_ν,acc is roughly given by the mass accretion rate ${\dot{M}}_{\mathrm{acc}}$ and the gravitational potential at the neutron star surface (Fischer et al. 2009),

$$\begin{eqnarray}&&{L}_{{\nu }_{e}}\approx {L}_{\nu ,\mathrm{acc}}=\eta \displaystyle \frac{{{GM}}_{\mathrm{PNS}}{\dot{M}}_{\mathrm{acc}}}{{R}_{\mathrm{PNS}}}\,,\end{eqnarray} \tag{ 12 }$$

where R_PNS = 5R_gain/7 (Müller et al. 2016) and η is an efficiency parameter, which specifies the conversion of accretion energy into neutrino luminosity. Note that the neutrino luminosity includes two components: accretion luminosity and diffusion luminosity (Fischer et al. 2009). It is difficult to completely separate these components. By taking into account the contribution of the diffusion luminosity L_ν,diff, η would exceed unity (see Müller & Janka 2014). In this study, we calibrate our model with the electron neutrino luminosity at ∼1 s after the SN shock revival, which is well studied in the first-principles calculations of CCSN explosions (Müller & Janka 2014), and use η = 1.

From Equations (5) to (10), we can write the ejection rate of the neutrino-driven wind with the accretion flow onto the PNS as

$$\begin{eqnarray}\begin{array}{rcl}{\dot{M}}_{\mathrm{wind}} & \approx & 1.3\times {10}^{-2}\,{M}_{\odot }\,{{\rm{s}}}^{-1}\\ & & \times \,{f}_{{\rm{\Omega }}}{\left(\displaystyle \frac{(1-{f}_{{\rm{\Omega }}}){\dot{M}}_{\mathrm{acc},\mathrm{iso}}}{0.1{M}_{\odot }{{\rm{s}}}^{-1}}\right)}^{\tfrac{2\alpha +\beta }{3}}{\left(\displaystyle \frac{{M}_{\mathrm{PNS}}}{1.4{M}_{\odot }}\right)}^{2\alpha -\beta +\gamma },\end{array}\end{eqnarray} \tag{ 13 }$$

where the value of the indices estimated by adopting our model are $\tfrac{2\alpha +\beta }{3}=2$ and 2α − β + γ = − 5/2.

In this section, we describe a potential of the neutrino-driven wind to solve the Ni problem. A summary of the claim of the Ni problem is that while on average ∼ 0.07 M_⊙ of ⁵⁶Ni should be synthesized in a canonical CCSNe, the first-principles simulations have been able to reproduce less than half of it. Our purpose is to investigate whether the wind can eject 0.07 M_⊙ of ⁵⁶Ni or not.

⁵⁶Ni is primarily synthesized where a material of an electron fraction Y_e ≈ 0.5 experiences NSE in a high temperature environment (≳5 × 10⁹ K). This is because in NSE it is dominated by Fe-peak nuclei due to the largest binding energy per nucleon and a proton-to-nucleon ratio being close to the electron fraction of the environment. The neutrino-driven wind of interest now experiences NSE, as can be seen in Figure 2, and is thus a sufficient environment for the synthesis of ⁵⁶Ni. In addition, detailed nucleosynthesis calculations confirmed the dominance of ⁵⁶Ni in an outflow of electron fraction Y_e ≈ 0.5 (Wanajo et al. 2018). Therefore, assuming that all the ejected material has Y_e = 0.5 as described in Section 2.1, it could be understood that the maximum ejected amount of Y_e = 0.5 material is available as a robust upper limit on the ejected ⁵⁶Ni mass. Namely, in this section, we estimate the ejectable maximum amount of Y_e = 0.5 material (which corresponds to the maximum amount of ⁵⁶Ni) by integrating the mass ejection rate of our wind model.

We integrate Equation (13) with Equation (8) from the shock revival time (t = 0) as follows,

$$\begin{eqnarray}\begin{array}{rcl}{M}_{\mathrm{ej},\infty } & = & {\displaystyle \int }_{0}^{\infty }{dt}\,\,{\dot{M}}_{\mathrm{wind}}(t)\\ & \approx & 1.3\times {10}^{-2}\,{M}_{\odot }\,{{\rm{s}}}^{-1}\,{f}_{{\rm{\Omega }}}{\left(\displaystyle \frac{(1-{f}_{{\rm{\Omega }}}){\dot{M}}_{\mathrm{acc},0}}{0.1{M}_{\odot }{{\rm{s}}}^{-1}}\right)}^{2}\\ & & \times \,{\left(\displaystyle \frac{{M}_{\mathrm{PNS},0}}{1.4{M}_{\odot }}\right)}^{-5/2}{\displaystyle \int }_{0}^{\infty }{dt}{\left(\displaystyle \frac{t}{{t}_{0}}+1\right)}^{-4}\\ & = & 4.3\times {10}^{-3}\,{M}_{\odot }\\ & & \times \,{f}_{{\rm{\Omega }}}\left(\displaystyle \frac{{t}_{0}}{1\,{\rm{s}}}\right){\left(\displaystyle \frac{(1-{f}_{{\rm{\Omega }}}){\dot{M}}_{\mathrm{acc},0}}{0.1{M}_{\odot }{{\rm{s}}}^{-1}}\right)}^{2}{\left(\displaystyle \frac{{M}_{\mathrm{PNS},0}}{1.4{M}_{\odot }}\right)}^{-5/2},\end{array}\end{eqnarray} \tag{ 14 }$$

where we neglect the mass evolution of the PNS within the integration, which decreases the mass ejection by ${ \mathcal O }(10) \%$. Moreover, when we adopt f_Ω = 1/3, which gives the geometric effect term ${f}_{{\rm{\Omega }}}{(1-{f}_{{\rm{\Omega }}})}^{2}$ maximum, Equation (14) is written as

$$\begin{eqnarray}\begin{array}{rcl}{M}_{\mathrm{ej},\infty }^{\max } & = & 6.4\times {10}^{-4}\,{M}_{\odot }\\ & & \times \,\left(\displaystyle \frac{{t}_{0}}{1\,{\rm{s}}}\right){\left(\displaystyle \frac{{\dot{M}}_{\mathrm{acc},0}}{0.1{M}_{\odot }{{\rm{s}}}^{-1}}\right)}^{2}{\left(\displaystyle \frac{{M}_{\mathrm{PNS},0}}{1.4{M}_{\odot }}\right)}^{-5/2}.\end{array}\end{eqnarray} \tag{ 15 }$$

The main goal of this paper is to find out the ejectable maximum mass of our wind model from Equation (15). When we adopt the initial mass of the PNS to be 1.4 M_⊙ (e.g., Müller & Janka 2014), then two free parameters, ${\dot{M}}_{\mathrm{acc},0}$ and t₀, remain. We first adopt ${\dot{M}}_{\mathrm{acc},0}=1\,{M}_{\odot }\,{{\rm{s}}}^{-1}$ for ${\dot{M}}_{\mathrm{acc},0}$ as a phenomenological upper limit from Figure 4. t₀ is constrained by the maximum PNS mass as follows. Taking into account the geometric effects (Equation (10)), and ignoring the mass decreases due to ejected wind because it is relatively small, the mass of the PNS can be written as

$$\begin{eqnarray}\begin{array}{rcl}{M}_{\mathrm{PNS}}(t) & \approx & {M}_{\mathrm{PNS},0}+{\displaystyle \int }_{0}^{t}{dt}^{\prime} \,{\dot{M}}_{\mathrm{acc}}(t^{\prime} )\\ & = & {M}_{\mathrm{PNS},0}+(1-{f}_{{\rm{\Omega }}}){\dot{M}}_{\mathrm{acc},0}{t}_{0}\left(\displaystyle \frac{t}{t+{t}_{0}}\right),\end{array}\end{eqnarray} \tag{ 16 }$$

where M_PNS,0 is the initial mass of the PNS (at the time of the SN shock revival). We assume that the wind ceases when the PNS mass exceeds a black hole mass. These assumptions of the total accretion mass (Equation (16) and the PNS mass up to 2.1 M_⊙) gives us $(1-{f}_{{\rm{\Omega }}}){\dot{M}}_{\mathrm{acc},0}\,{t}_{0}\leqslant 0.7\,{M}_{\odot }$. To conclude, the ejectable maximum mass of our wind model is given with ${\dot{M}}_{\mathrm{acc},0}=1\,{M}_{\odot }\,{{\rm{s}}}^{-1}$, t₀ = 1.05 s, and f_Ω = 1/3 in Equation (15) as follows,

$$\begin{eqnarray}&&{M}_{\mathrm{ej},\infty }^{\max }=0.067\,{M}_{\odot }.\end{eqnarray} \tag{ 17 }$$

If most of this compensation from the wind is added at late phase, which is later than the computation time of the first-principles simulations, this value is then sufficient to compensate for the lack of ⁵⁶Ni in the recent Ni problem. In the following, we investigate the time evolution of the cumulative ejected mass of the wind and discuss the nature of the explosion that could solve the Ni problem.

Figure 5 shows the time evolution of the cumulative ejected mass of the wind model. The cumulative ejected mass is given as

$$\begin{eqnarray}\begin{array}{rcl}{M}_{\mathrm{ej}}({t}_{e}) & = & {\displaystyle \int }_{0}^{{t}_{e}}{dt}\,\,{\dot{M}}_{\mathrm{wind}}(t)\\ & \approx & 6.4\times {10}^{-4}\,{M}_{\odot }\left[1-{\left(\displaystyle \frac{{t}_{0}}{{t}_{0}+{t}_{e}}\right)}^{3}\right]\\ & & \times \,\left(\displaystyle \frac{{t}_{0}}{1\,{\rm{s}}}\right){\left(\displaystyle \frac{{\dot{M}}_{\mathrm{acc},0}}{0.1{M}_{\odot }{{\rm{s}}}^{-1}}\right)}^{2}{\left(\displaystyle \frac{{M}_{\mathrm{PNS},0}}{1.4{M}_{\odot }}\right)}^{-5/2},\end{array}\end{eqnarray} \tag{ 18 }$$

where t_e is time after the neutrino-driven wind starts to blow, corresponding to the time after the SN shock revival. We adopt f_Ω = 1/3, which gives the geometric effect term ${f}_{{\rm{\Omega }}}{(1-{f}_{{\rm{\Omega }}})}^{2}$ maximum. When we fix M_PNS,0 = 1.4 M_⊙ for simplicity, then two free parameters, ${\dot{M}}_{\mathrm{acc},0}$ and t₀, remain to determine the trajectory of the time evolution. As with the condition for Equation (17), these two parameters are given by the conditions ${\dot{M}}_{\mathrm{acc},0}\leqslant 1.0\,{M}_{\odot }\,{{\rm{s}}}^{-1}$ and $(1-{f}_{{\rm{\Omega }}}){\dot{M}}_{\mathrm{acc},0}\,{t}_{0}\leqslant 0.7\,{M}_{\odot }$ , respectively, from the phenomenological accretion model (Figure 4) and the limits of the total accretion mass (Equation (16) and the PNS mass up to 2.1 M_⊙). A degenerate set of parameters that converge to the same M_ej,∞ are shown in the same color in Figure 5. We further compare the time evolution with the multi-D first-principles simulations, especially at t_e ≲ 1 s (Wanajo et al. 2018), which are shown in Figure 5 as rhombus points. It indicates that, within parameter ambiguities, the total ejectable amount of the neutrino-driven wind is roughly determined within 1 s from the onset of the blowing, which is reachable for first-principles simulations. Moreover, we also find that the supplementable amount from the wind at a later phase (t_e ≳ 1 s) remains ≲ 0.01 M_⊙. It also shows that, comparing with first-principles simulations, the expected total amount from the wind in a canonical CCSN explosion is about ≲ 0.03 M_⊙ at most.

Zoom In Zoom Out Reset image size

We should mention that, in fact, some models in first-principles simulations of CCSNe have produced a sufficient amount of ⁵⁶Ni in the total of explosive and the wind nucleosynthesis (Bruenn et al. 2016; Eichler et al. 2018; Bollig et al. 2020; Burrows et al. 2020). However, the claim of the Ni problem is that we are now struggling to reproduce sufficient ⁵⁶Ni amount as a canonical nature of the CCSNe explosion, and many previous studies place their hopes on the supplement from the neutrino-driven wind at a later phase (e.g., Wongwathanarat et al. 2017; Wanajo et al. 2018). Our conclusion on this issue is that, in order to compensate for sufficient ⁵⁶Ni by the neutrino-driven wind, it is preferred to have an active ejection in the early phase rather than a continuous ejection until the later phase. We note that while the wind driven by magneto-rotational explosions may be more energetic, it is not expected to solve the Ni problem due to low ⁵⁶Ni as it is ejected with a low electron fraction (Nishimura et al. 2015). It should also be emphasized that the total amount of synthesized ⁵⁶Ni can be estimated robustly if the first-principles simulations are performed up to ∼2 s.

To further expand the discussion in Equations (14) and (18) about the ejectable amount of the wind, here we discuss the effect of the time evolution of the PNS mass on the results above. We first introduce the ratio of the total accretion mass to the initial PNS mass as

$$\begin{eqnarray}&&\epsilon =\displaystyle \frac{1}{{M}_{\mathrm{PNS},0}}{\int }_{0}^{\infty }{dt}\,\,{\dot{M}}_{\mathrm{acc}}(t).\end{eqnarray} \tag{ 19 }$$

Using this ratio , the time evolution of the PNS mass expressed in Equation (16) is written as

$$\begin{eqnarray}&&{M}_{\mathrm{PNS}}(t)\approx {M}_{\mathrm{PNS},0}\left(1+\epsilon \,\displaystyle \frac{t}{t+{t}_{0}}\right),\end{eqnarray} \tag{ 20 }$$

where we note that this ignores the mass loss due to wind, as in Equation (16). No matter how light the initial PNS mass is assumed to be at SN shock revival (e.g., M_PNS,0 ∼ 1.2 M_⊙), this is less than unity, since this mass will never grow more than twice as large due to the constraint that the wind stops when it becomes a black hole. Then we can update the time evolution term of the integral in Equation (14) as follows

$$\begin{eqnarray}&&\begin{array}{l}{\displaystyle \int }_{0}^{\infty }{dt}{\left(\displaystyle \frac{t}{{t}_{0}}+1\right)}^{-4}{\left(1+\epsilon \displaystyle \frac{t}{t+{t}_{0}}\right)}^{-5/2}\\ \ =\,\displaystyle \frac{2\left({\epsilon }^{2}-4\epsilon +8\sqrt{\epsilon +1}-8\right)}{3{\epsilon }^{3}}{t}_{0}\,\\ \ \approx \,\left(\displaystyle \frac{1}{3}-\displaystyle \frac{5}{24}\epsilon +\displaystyle \frac{7}{48}{\epsilon }^{2}\right){t}_{0}.\end{array}\end{eqnarray} \tag{ 21 }$$

where = 0 corresponds to the result in Equation (14). We find that, for < 1, the effect of mass dependence only decreases the wind mass loss. The same can be argued for the time evolution of the cumulative mass, i.e., the discussion in Equation (18) and Figure 5.

Figure 6 shows a comparison of the results between analytical integration without the effect of the PNS masses and the numerical integration with the effect. This figure shows that the result of Equation (14) gives a robust upper limit on the ejectable mass of the wind. This comparison indicates that the results of Equations (14) and (18) give robust upper limits on the ejectable mass of the wind.

Zoom In Zoom Out Reset image size