SIMULATIONS OF ACCRETION POWERED SUPERNOVAE IN THE PROGENITORS OF GAMMA-RAY BURSTS

To include rotation and angular momentum transport in our one-dimensional model, we track the specific angular momentum ℓ ≡ rv_ϕ, where v_ϕ is the azimuthal velocity, which we interpret as the mass-weighted spherical average of an underlying polar-angle-dependent angular momentum ℓ(r, θ). If, e.g., spherical shells rotate rigidly, ℓ(r, θ)∝sin ²θ, and the fluid density is spherically symmetric, then the one-dimensional specific angular momentum is two-thirds of the midplane value, ℓ = 2/3ℓ_mid. The azimuthal Navier-Stokes equation, combined with the equation of mass continuity, then implies the one-dimensional angular momentum transport equation (see, e.g., Thompson et al. 2005)

$$\begin{eqnarray} & &\frac{\partial (\rho \ell)}{\partial t} + \frac{1}{r^2}\frac{\partial (r^2 v_r \rho \ell)}{\partial r} - \frac{1}{r^2}\frac{\partial }{\partial r}(r^3 \nu \rho \sigma _{r \phi }) = 0, \end{eqnarray} \tag{ 1 }$$

where ν is a shear viscosity and

$$\begin{equation} \sigma _{r \phi } = r \frac{\partial }{\partial r} \left(\frac{\ell }{r^2}\right) \end{equation} \tag{ 2 }$$

is the r − ϕ component of the shear tensor. The energy dissipated through shear viscosity was accounted for by including the specific heating rate (see, e.g., Landau & Lifshitz 1959)

$$\begin{equation} \dot{\epsilon }_{\rm visc} \equiv \frac{Q_{\rm visc}}{\rho }= \nu \sigma _{r \phi }^2, \end{equation} \tag{ 3 }$$

where Q_visc denotes the volumetric viscous heating rate. The dimensional reduction in Equation (1) is inaccurate in regions where the disk is geometrically thin. There the mass-weighted spherical average closely approximates the midplane value, ℓ ∼ ℓ_mid. We ignore this effect, but we do incorporate thermodynamic corrections addressing the transition to a thin disk in Section 2.6.

Our treatment of shear viscosity is similar to our methodology in Lindner et al. (2010), and for completeness we reproduce our methodology here. Since we do not simulate the magnetic field of the fluid, we utilize a local definition of the shear viscosity to emulate the magnetic stress arising from the nonlinear development of the magnetorotational instability (MRI; Balbus & Hawley 1998 and references therein). It should be kept in mind, however, that the effects of MRI are in some respects very different from those of the viscous stress. For example, the thick disk surrounding our collapsar black hole is convective; in unmagnetized accretion flows convection transports angular momentum inward, toward the center of rotation (Ryu & Goodman 1992; Stone & Balbus 1996; Igumenshchev et al. 2000), whereas in magnetized flows, convection can also transport angular momentum outward (Balbus & Hawley 2002; Igumenshchev 2002; Igumenshchev et al. 2003; Christodoulou et al. 2003). Although we include treatment for convective energy flux and compositional mixing (see Section 2.5), we do not include angular momentum transport by convection.

Our definition of the local viscous stress emulating the MRI must be valid under rotationally supported, pressure supported, and freely falling conditions, and we proceed as in Lindner et al. (2010). Thompson et al. (2005) suggest that since the wavenumber of the fastest growing MRI mode, which is given by the dispersion relation v_Ak ∼ Ω where v_A is the Alfvén velocity and Ω = v_ϕ/r is the angular velocity, should be about the gas pressure scale height in the saturated quasi-state state, k∝H⁻¹, the Maxwell ρv²_A and viscous νρΩ stresses (up to factors in |dln Ω/dln r| that we neglect) can be equated if the viscosity is given by

$$\begin{equation} \nu _{\rm MRI} = \alpha H^2 \Omega, \end{equation} \tag{ 4 }$$

where α is a dimensionless parameter. If the pressure scale height is defined locally,

$$\begin{equation} H=|\mathbf {\nabla }\ln P|^{-1}, \end{equation} \tag{ 5 }$$

the viscosity defined in Equation (4) suffers from divergences at pressure extrema. To alleviate this problem, as in Lindner et al. (2010), we define a second viscosity according to the Shakura & Sunyaev (1973) prescription

$$\begin{equation} \nu _{\rm SS} = \alpha \frac{P}{\rho } \Omega ^{-1}. \end{equation} \tag{ 6 }$$

Shakura–Sunyaev viscosity overestimates the magnetic stress in stratified hydrostatic atmospheres. We thus set the viscosity in Equations (1) and (3) to equal the harmonic mean of the above two viscosities

$$\begin{equation} \nu = \frac{2\,\nu _{\rm MRI}\,\nu _{\rm SS} }{\nu _{\rm MRI}+\nu _{\rm SS}}, \end{equation} \tag{ 7 }$$

where the pressure gradient in Equation (5) is calculated by the finite differencing of pressure in neighboring fluid cells. Additionally, we have applied a Gaussian kernel smoothing to the radial dependence of H to help filter short-wavelength numerical instabilities. We describe this procedure in Section 2.5.

In FLASH, we treat specific angular momentum as a conserved "mass scalar" that is being advected with the fluid, which makes ρℓ a conserved variable; the corresponding centrifugal force is then incorporated in the calculation of the gravitational acceleration as we explain in Section 2.2 below. Then the third parabolic term in Equation (1) is computed explicitly through the inclusion of the radial ρℓ-flux −rνρσ_rϕ in the advection of ℓ.

Numerical stability of an explicit treatment of a parabolic term places an upper limit on the time step

$$\begin{equation} \Delta t < \frac{\Delta r^2}{2\nu }, \end{equation} \tag{ 8 }$$

where Δr is the grid resolution. For α ≫ 0.01, the viscous time step in our simulations becomes significantly shorter than the Courant time step. In our test integrations with a γ-law equation of state (EOS; Lindner et al. 2010), we find that, while not implying an outright instability, a choice of Δt that saturates the limit in Equation (8) results in weak stationary staggered perturbations in the fluid variables. We ignore this complication and allow our time step to be set by the limit in Equation (8) of the cell with the smallest viscous diffusion time across the cell.

We calculate contributions to the gravitational potential from a central point mass and a spherically symmetric extended envelope. General relativistic effects become important at the innermost radius, which in some simulations is as small as r_min = 25 km. At radii r ∼ r_min, the black hole dominates the enclosed mass after about 0.5 s. Thus, we describe the gravity of the black hole using the approximate, pseudo-Newtonian gravitational force for a rotating black hole proposed by Artemova et al. (1996), which is a generalization of the Paczyński & Wiita (1980) pseudopotential to rotating black holes. However, we continue to calculate the gravity of the fluid in the Newtonian limit. The Artemova et al. gravitational acceleration in the equatorial plane of a rotating black hole is given by

$$\begin{equation} \mathbf {g}_{\rm BH} (r,\theta =\pi /2)= -\frac{{\it G M}_{\rm B H}}{r^{2-\beta }(r - r_{\rm H})^\beta } \mathbf {\hat{r}}, \end{equation} \tag{ 9 }$$

where r_H = [1 + (1 − a²)^1/2]GM_BH/c² is the radius of the event horizon expressed in terms of the dimensionless spin parameter a, and β = r_ISCO/r_H − 1 is a dimensionless exponent with r_ISCO denoting radius of the innermost stable prograde equatorial circular orbit. We assume a dimensionless spin parameter of a = 0.9 in these calculations. Our treatment does not incorporate general relativistic corrections to the viscous stress and momentum equations (see, e.g., Beloborodov 1999).

We adopt the form of the gravitational acceleration in Equation (9), which was derived for the equatorial plane of the black hole, to represent the mass-weighted spherical average of the gravitational acceleration, by setting g_BH(r) = g_BH(r, θ = π/2). This approximation is appropriate when the accreting mass is concentrated in the equatorial plane, especially when the innermost disk is geometrically thin, and is probably rather inaccurate for an accretion flow that is geometrically thick down to r_ISCO. Our simulations predict a geometrically thin disk at r ≲ 100 km or greater radii after material has circularized in our simulation, so this assumption seems adequate.

For each zone, the gravitational acceleration due to fluid self-gravity is calculated from

$$\begin{eqnarray} \mathbf {g}_{\rm self} (r_i) &=& -\frac{4 \pi }{3} \frac{G}{r^2} \left\lbrace \rho _i \left[r_i^3 - \left(r_i - \frac{\Delta r_i}{2}\right)^3 \right] \right. \nonumber \\ & & + \left. \sum _{r_k<r_i} \rho _k \left[ \left(r_k + \frac{\Delta r_k}{2}\right)^3 - \left(r_k - \frac{\Delta r_k}{2}\right)^3 \right] \right\rbrace \mathbf {\hat{r}},\nonumber\\ \end{eqnarray} \tag{ 10 }$$

where Δr_i and Δr_k are the radial widths of the grid cells. The net gravitational and inertial acceleration in our calculation is then given by

$$\begin{equation} \mathbf {a}_{\rm {tot}} = \mathbf {g}_{\rm BH} + \mathbf {g}_{\rm self} + \mathbf {a}_{\rm cent}, \end{equation} \tag{ 11 }$$

where

$$\begin{equation} \mathbf {a}_{\rm cent} = \frac{\ell ^2}{r^3} \mathbf {\hat{r}} \end{equation} \tag{ 12 }$$

is the centrifugal acceleration.

To calculate the internal energy of the fluid, we use the Helmholtz EOS of Timmes & Swesty (2000) included with the FLASH distribution, which accounts for the contributions to pressure and other thermodynamic quantities from radiation, ions, electrons, positrons, and Coulomb corrections. We track the abundances of 47 nuclear isotopes treated in the nuclear statistical equilibrium (NSE) calculations of Seitenzahl et al. (2008) and pass the local nuclear composition to the EOS as input. Given density, temperature, and nuclear composition, the Helmholtz EOS provides the internal energy, density, pressure, entropy, specific heats, adiabatic indices, electron chemical potential, and various derivative thermodynamic quantities. During the course of the thermodynamic update and the cooling update which is operator split from the thermodynamic update, the temperature must be derived from the internal energy, and in the Helmholtz EOS this is achieved by numerically solving for the implicit relation

$$\begin{equation} \epsilon _{\rm EOS}(\rho,T,{\mathbf X}) = \epsilon \end{equation} \tag{ 13 }$$

for the temperature, where is the specific internal energy and $\mathbf {X}\equiv (X_1, \ldots, X_{47})$ is the vector of isotopic mass fractions X_i.

The fluid heats and cools in response to nuclear compositional transformation. We do not integrate a nuclear reaction network, but instead model the change of the nuclear composition as a gradual convergence to NSE in the part of the flow where the convergence timescale τ_NSE is comparable to or shorter than the age of the system. In this model, as we explain below, the nuclear composition responds instantaneously to a change of the temperature, implying that the dependence of the composition on the temperature must be taken into account, in a manner that conserves the combined specific internal and nuclear energy + _nuc when solving the EOS for temperature. Here, _nuc is the specific (negative) nuclear binding energy of the fluid

$$\begin{equation} \epsilon _{\rm nuc} = \sum _{i} \frac{X_i E_{{\rm B},i}}{A_i m_p}, \end{equation} \tag{ 14 }$$

while E_{B, i} is the negative nuclear binding energy of the isotope and A_i is the atomic mass of the isotope.

The timescale for convergence to NSE can be approximated via (Khokhlov 1991; see also Calder et al. 2007)

$$\begin{equation} \tau _{\rm NSE}=\rho ^{0.2} \exp \left(\frac{179.7}{T_9} - 40.5\right) \,{\rm s}, \end{equation} \tag{ 15 }$$

where T = 10⁹ T₉ K and ρ is the density in g cm⁻³. At relevant densities, this timescale is of the order of 1 s for T_NSE ≈ 4 × 10⁹ K. We calculate the nuclear mass fractions using the publicly available solver of Seitenzahl et al. (2008) which solves for the NSE mass fractions X_{NSE, i} of 47 nuclear isotopes as a function of density ρ, temperature T, and proton-to-nucleon ratio Y_e = ∑_iZ_iX_i/A_i, where Z_i is the atomic number an isotope. At temperatures T > 3 × 10⁹ K we model convergence to NSE via

$$\begin{equation} \left(\frac{\partial X_i}{\partial t}\right)_{\rm nuc} = \frac{X_{i,{\rm NSE}} (\rho,T_{\rm NSE},Y_e)- X_i}{\tau _{\rm NSE}(\rho,T)}, \end{equation} \tag{ 16 }$$

where T_NSE is the temperature that the fluid element would have given enough time to relax into NSE while keeping the total specific energy + _NSE and proton-to-nucleon ratio Y_e fixed. The temperature T_NSE is implicitly defined by the condition (cf. Equation (13))

$$\begin{eqnarray} & &\epsilon _{\rm EOS}[\rho,T_{\rm NSE},{\mathbf X}_{{\rm NSE}}(\rho,T_{\rm NSE},Y_e)] + \epsilon _{\rm nuc}[{\mathbf X}_{{\rm NSE}}(\rho,T_{\rm NSE},Y_e)] \nonumber \\ &&\quad= \epsilon (\rho,T,\mathbf {X})+\epsilon _{\rm nuc}(\mathbf {X}). \end{eqnarray} \tag{ 17 }$$

This condition ensures that the sum of the internal and nuclear energy densities in NSE would equal the sum of the two energy densities in the model. We solve Equation (17) for T_NSE(ρ, Y_e, , _nuc) iteratively and then update the abundances by discretizing Equation (16) with

$$\begin{eqnarray} X_i (t&+&\Delta t) = X_{i,\rm {NSE}}(\rho,T_{\rm NSE},Y_e)\nonumber \\ &+& [ X_i(t) - X_{i,\rm {NSE}} (\rho,T_{\rm NSE},Y_e) ]\exp \left[-\frac{\Delta t}{ \tau _{\rm {NSE}}(\rho,T)}\right].\nonumber \\ \end{eqnarray} \tag{ 18 }$$

Following the update of the nuclear mass fractions, we update the specific internal energy to account for heating or cooling due to any change in specific nuclear binding energy

$$\begin{equation} \epsilon (t+\Delta t) = \epsilon (t) + \epsilon _{\rm nuc}[\mathbf {X}(t)] - \epsilon _{\rm nuc}[\mathbf {X}(t+\Delta t)], \end{equation} \tag{ 19 }$$

and finally update the temperature from Equation (13).

This prescription does not affect the proton-to-nucleon ratio Y_e; that latter is a conserved mass scalar in our simulations. Thus, the expected partial neutronization in the mildly degenerate innermost segment of the accretion flow is not calculated and our prescription cannot be used to accurately estimate the ⁵⁶Ni fraction within the Fe-group elements synthesized in the simulation.

The hot innermost accretion flow cools via neutrino emission. At the densities observed in our simulation, the disk and stellar atmosphere are transparent to neutrinos. The two most significant neutrino-emission channels (e.g., Di Matteo et al. 2002, and references therein) are as follows.

• 1.  
  Pair capture on free nucleons (the Urca process). p + e⁻ → n + ν and $n + e^+ \rightarrow p + \bar{\nu }$. The cooling rate is

  $$\begin{equation} Q_{eN} = 9\times 10^{33} \rho _{10} T_{11}^6 X_{\rm nuc} \,{\rm erg \,cm^{-3}\,s^{-1}}, \end{equation} \tag{ 20 }$$

  where ρ = 10¹⁰ρ₁₀ g cm⁻³, T = 10¹¹ T₁₁ K, and X_nuc = X_p + X_n is the mass fraction in free nucleons.
• 2.  
  Pair annihilation ($e^- + e^+ \longrightarrow \nu + \bar{\nu }$ ). The cooling rate is

  $$\begin{equation} \vspace*{-1.5pt} Q_{e^+e^-} = 1.5\times 10^{33} T_{11}^9 \,{\rm erg \,cm^{-3} \,s^{-1}}. \end{equation} \tag{ 21 }$$

  All three flavors, e, μ, and τ, of neutrinos are included.

We have included the above neutrino cooling rates in our calculations, where losses are computed via

$$\begin{equation} \vspace*{-1.5pt} \epsilon (t+\Delta t)= \epsilon (t) - \frac{Q_\nu }{\rho }\Delta t, \end{equation} \tag{ 22 }$$

where $Q_\nu =Q_{eN} \,{+}\, Q_{e^+e^-}$ is the total volumetric neutrino cooling rate. The update of the internal energy due to cooling is operator split from the update due to nuclear compositional change.

We introduce convective energy transport and compositional mixing within the framework of mixing length theory (e.g., Kuhfuß 1986). In the calculation of the convective transport fluxes, we ignore the radial variation of the mean molecular weight as well as rotation, and the condition for instability is simply the Schwarzschild criterion, ∂s/∂r < 0. Then, in unstable zones, the convective energy flux is

$$\begin{equation} \vspace*{-1.5pt} F_{\rm conv} = -\frac{1}{2} c_P \rho v_{\rm conv} \lambda _{\rm conv}\left(\frac{\partial T}{\partial s} \right)_P \, \frac{\partial s}{\partial r}, \end{equation} \tag{ 23 }$$

where c_P is the specific heat at constant pressure, λ_conv is the length over which convection occurs, s is specific entropy, and v_conv is the convective velocity. The convective velocity can be approximated by

$$\begin{equation} \vspace*{-1.5pt} v_{\rm conv}\sim \frac{1}{2}\lambda _{\rm conv} \left[-\frac{g}{\rho } \left(\frac{\partial \rho }{\partial T}\right)_P \frac{T}{c_P} \frac{ds}{dr}\right]^{1/2}, \end{equation} \tag{ 24 }$$

where g < 0 is the gravitational acceleration in the local rest frame of the convectively unstable fluid

$$\begin{equation} \vspace*{-1.5pt} g = g_{\rm grav} - \frac{dv}{dt} = \frac{1}{\rho } \frac{\partial P}{\partial r} = -\frac{P}{\rho H}, \end{equation} \tag{ 25 }$$

and g_grav = g_BH + g_self is the net gravitational acceleration in the inertial frame, v in the second step denotes the mass-weighted spherical average of the fluid velocity at radius r. To parameterize our uncertainty regarding the value of the convective mixing length, we introduce a dimensionless parameter $\xi _{\rm conv}\sim {\cal O}(1)$ defined as

$$\begin{equation} \vspace*{-1.5pt} \xi _{\rm conv}\equiv \left(\frac{\lambda _{\rm conv}}{H}\right)^2. \end{equation} \tag{ 26 }$$

Then, combining Equations (23)–(26), we obtain the standard expression

$$\begin{equation} \vspace*{-1.5pt} F_{\rm conv} = \frac{1}{4} \xi _{\rm conv} H^2 c_P \left[ -\frac{P}{H} \left(\frac{\partial \rho }{\partial T} \right)_P \right]^{1/2} \left(-\frac{T}{c_P}\frac{\partial s}{\partial r} \right)^{3/2},\hspace{4pt} \end{equation} \tag{ 27 }$$

which is appropriate even when the fluid is not in hydrostatic equilibrium and v_r ≠ 0.

In evaluating the convective energy flux at a boundary (face) of a computational cell, we use face-centered linear interpolation of the density, temperature, and pressure. The zone-centered values of the specific heat c_P, specific entropy s, and thermodynamic derivatives (∂P/∂T)_ρ and (∂P/∂ρ)_T are returned by the EOS routine, and the face-centered values are again computed by linear interpolation. Then, (∂ρ/∂T)_P is calculated from

$$\begin{equation} \left(\frac{\partial \rho }{\partial T} \right)_P = -\left(\frac{\partial P}{\partial T} \right)_\rho \bigg/ \left(\frac{\partial P}{\partial \rho } \right)_T . \end{equation} \tag{ 28 }$$

The convective energy flux never exceeds

$$\begin{equation} F_{\rm conv} \le \rho \epsilon c_{\rm s}, \end{equation} \tag{ 29 }$$

where c_s = (γ_cP/ρ)^1/2 is the adiabatic sound speed and γ_c is the adiabatic index.

We anticipate that a local application of MLT, in which the expression for the convective energy flux contains a pressure derivative in the denominator, may contain an instability. The instability is an artifact of modeling the intrinsically nonlocal convective energy transport with a local nonlinear differential operator. To control—if not entirely prevent—undesirable outcomes of the instability, we filter short wavelength perturbations in the calculation of the pressure scale height H that enters our estimates of the viscosity and the energy flux transported by convection by applying a Gaussian smoothing

$$\begin{equation} P_{\rm smooth} (r) = \frac{\sum _i k_i (r) \,P_i}{\sum _i k_i(r)}, \end{equation} \tag{ 30 }$$

where the summations are over all of the cells in the simulation, and the spherically averaged smoothing kernel k_i is given by

$$\begin{eqnarray} &&k_i (r)\!=\! \frac{1}{2\sqrt{2\pi }}\frac{\Delta r_i r_i}{r\sigma } \left\lbrace \exp \left[ -\frac{(r - r_i)^2}{2 \sigma ^2} \right] -\exp \left[ -\frac{(r + r_i)^2}{2 \sigma ^2} \right] \right\rbrace . \nonumber\\ \end{eqnarray} \tag{ 31 }$$

Here, σ is a radius-dependent smoothing length that we set to (1/2)r. Similarly, in the evaluation of the specific entropy derivative in Equation (27), we smooth the specific entropy s via

$$\begin{equation} s_{\rm smooth} (r) = \frac{\sum _i k_i (r) \,\rho _i \,s_i}{\sum _i k_i(r) \, \rho _i }. \end{equation} \tag{ 32 }$$

The filtering affects only the evaluation of F_conv and helps avoid breakdown of our transport scheme, but residual artificial non-propagating waves do develop, and saturate, on wavelengths comparable to the smoothing length.

The accretion shock formally presents a negative entropy gradient but physically does not give rise to convection. The upstream of the shockwave is marginally convectively stable as the shockwave traverses the progenitor's convective core, and becomes absolutely stable in the radiative envelope. To prevent spurious convection across the shock transition, we modify the convective flux to decline to zero linearly near the shock

$$\begin{equation} F_{{\rm conv},{\rm mod}}(r) =\left\lbrace \begin{array}{@{}l@{\quad }l@{}}(1-r/r_{\rm shock}) F_{\rm conv}(r), & r<r_{\rm shock}, \\ \ms 0, & r\ge r_{\rm shock},\end{array}\right. \end{equation} \tag{ 33 }$$

where r_shock is the radius of the accretion shock front which we track during the simulation.

Murphy & Meakin (2011) argue that on physical grounds, in quasi-stationary "stalled" shocks in the standard core-collapse context, the distance from the shock r_shock–r is the appropriate convective length scale near the shock, as convective eddies can grow to the largest size available to them. If we had set the convective mixing length λ_conv proportional to the distance from the shock, which is the adaptation of MLT that Murphy & Meakin suggest, Equation (33) would have contained a quadratic factor (1 − r/r_shock)², instead of the linear factor (1 − r/r_shock) that we employ. The physically motivated modification of λ_conv of Murphy & Meakin, which we became aware of after the completion of this work, and our ad hoc version should give rise to similar dynamics, especially when the shock travels outward as in our simulations.

Convection also gives rise to compositional mixing in the convective region. We model the mixing of nuclear species in the diffusion approximation (e.g., Cloutman & Eoll 1976; Kuhfuß 1986)

$$\begin{equation} \left[\frac{\partial (\rho X_i)}{\partial t}\right]_{\rm mix} = - \frac{1}{r^2} \frac{\partial }{\partial r} (r^2 {\cal F}_{{\rm mix},i}), \end{equation} \tag{ 34 }$$

where

$$\begin{equation} {\cal F}_{{\rm mix},i} = - \frac{1}{3} \nu _{\rm conv} \rho \,\frac{\partial X_i}{\partial r} \end{equation} \tag{ 35 }$$

is the mass flux of species i transported by convection, while ν_conv is the compositional diffusivity which we take to be proportional to the convective velocity multiplied by the pressure scale height

$$\begin{equation} \nu _{\rm conv} = \xi _{\rm mix} v_{\rm conv} \lambda _{\rm conv}, \end{equation} \tag{ 36 }$$

and $\xi _{\rm mix}\sim {\cal O}(1)$ is a dimensionless parameter. We again apply the flux limitation behind the shock front in the form of the linear factor in Equation (33). The compositional diffusion is also subject to the time step limitation imposed in Equation (8). It is worth noting that compositional diffusion implies a flux of nuclear energy given by

$$\begin{equation} F_{\rm nuc,mix} = \sum _i \frac{E_{{\rm B},i}\,{\cal F}_{{\rm mix},i} }{A_i\, m_p}. \end{equation} \tag{ 37 }$$

The entropy transport equation implied by our algorithm is

$$\begin{eqnarray} &&\rho T \frac{d s}{d t}+ \frac{1}{r^2} \frac{\partial }{\partial r} (r^2 F_{\rm conv, mod}) = Q_{\rm visc} - Q_\nu + Q_{\rm nuc}, \end{eqnarray} \tag{ 38 }$$

where d/dt = ∂/∂t + v_r∂/∂r and

$$\begin{equation} Q_{\rm nuc} = - \rho \sum _i \frac{E_{{\rm B},i}}{A_i m_p} \frac{\partial X_i}{\partial t} \end{equation} \tag{ 39 }$$

is the rate of heating or cooling associated with nuclear compositional transformation (see Equations (14) and (16)).

We have thus far assumed a rotating, quasi-spherical accretion flow. However, near the black hole, where cooling by neutrino emission and nuclear photodisintegration into nucleons is significant, the flow can become geometrically thin (e.g., MacFadyen & Woosley 1999; Popham et al. 1999; Kohri et al. 2005; Chen & Beloborodov 2007). In this case, the quasi-spherical treatment underestimates the density, pressure, and temperature near the midplane of thin disk, where the bulk of the neutrino emission takes place. We introduce a correction that adjusts the temperature of the flow to be closer to the physical, thin-disk value.

In what follows, the quantities applying to the thin disk will be marked with tilde. Let $\tilde{H}_z$ denote the vertical half-thickness of the thin disk. We assume that the vertical half-thickness of the quasi-spherical flow is H_z ∼ (1/2)r. Since the thin and the quasi-spherical flow must contain the same column density, $\tilde{H}_z \tilde{\rho }\sim ({1}/{2}) r \rho$. Ignoring differences in nuclear composition between the thin and quasi-spherical flows, the same correspondence must apply to the total internal energies integrated along the vertical column, $\tilde{H}_z \tilde{\rho }\tilde{\epsilon }\sim ({1}/{2}) r \rho \epsilon$, and thus, $\tilde{\epsilon }\sim \epsilon$ while $\tilde{H}_z \tilde{P} \sim ({1}/{2}) r P$. Vertical force balance in the thin disk requires $\tilde{P}/\tilde{H}_z = \tilde{\rho }\tilde{|}g_z|$, where $\tilde{g}_z = -{\rm sgn}(z)\,g\,\tilde{H}_z/r$ is the gravitational acceleration in the z-direction and $g\equiv (\mathbf {g}_{\rm BH}+\mathbf {g}_{\rm self})\cdot \hat{\mathbf {r}}$ is the radial gravitational acceleration (see Section 2.2). Enforcing that $\tilde{H}_z \le H_z$, we obtain

$$\begin{equation} \tilde{H}_z = {\rm min} \left[\left(-\frac{rP}{\rho g}\right)^{1/2},\frac{r}{2}\right]. \end{equation} \tag{ 40 }$$

To account for the higher density in the thin disk, we could pass $\tilde{\rho }$ to the EOS. This, however, would result in a modified pressure $\tilde{P}$. Out of a possibly unfounded concern that a modification of the pressure would introduce spurious dynamics in the spherically averaged flow, we opted to modify our estimate of the disk midplane temperature in a manner not directly affecting the fluid pressure. This corrected temperature then enters the calculation of the neutrino cooling rate and the NSE composition, both of which are highly sensitive to the midplane temperature.

We estimate the midplane temperature $\tilde{T}$ from the following extrapolation:

$$\begin{equation} \ln \left(\frac{\tilde{T}}{T}\right)\approx \left(\frac{\partial \ln T}{\partial \ln \rho }\right)_\epsilon \ln \left(\frac{\tilde{\rho }}{\rho }\right), \end{equation} \tag{ 41 }$$

where the partial derivative, which we denote with χ, is evaluated at constant specific internal energy and can be expressed as

$$\begin{eqnarray} \chi &\equiv & \left(\frac{\partial \ln T}{\partial \ln \rho } \right)_\epsilon \nonumber \\ &=& -\left(\frac{\partial s}{\partial \ln \rho } \right)_T \bigg/ \left(\frac{\partial s}{\partial \ln T} \right)_\rho - \frac{P}{\rho c_V T}, \end{eqnarray} \tag{ 42 }$$

where c_V is the specific heat at constant volume. The quantities on the right-hand side of Equation (42) are all provided by the Helmholtz EOS.

Since $\tilde{\rho }/\rho \sim r / (2\tilde{H}_z)$, the midplane temperature of the disk can be approximated via

$$\begin{equation} \tilde{T} = \left(\frac{r}{2\tilde{H}_z}\right)^\chi T \equiv \Xi \,T, \end{equation} \tag{ 43 }$$

where the last equality defines the dimensionless temperature correction factor Ξ. To ensure continuity near the shock transition, we modify the correction factor to linearly approach unity at the shock transition by defining

$$\begin{equation} \Xi _{\rm mod} = 1 + \left(1-\frac{r}{r_{\rm shock}}\right) (\Xi - 1). \end{equation} \tag{ 44 }$$

For clarity of notation, we drop the subscript in Ξ_mod in what follows.

The correction introduced in Equation (43) affects both the temperature calculated from internal energy via the EOS and the NSE temperature calculated from the total internal and nuclear energy as described in Section 2.3 above. Equation (13) is corrected to become

$$\begin{equation} \epsilon _{\rm EOS}(\rho,\Xi ^{-1}\tilde{T},{\mathbf X}) = \epsilon, \end{equation} \tag{ 45 }$$

while Equation (17) is corrected to become

$$\begin{eqnarray} & &\epsilon _{\rm EOS}[\rho,\Xi ^{-1}\tilde{T}_{\rm NSE},{\mathbf X}_{{\rm NSE}}(\tilde{\rho },\tilde{T}_{\rm NSE},Y_e)] \nonumber \\ & &\quad + \,\epsilon _{\rm nuc}[{\mathbf X}_{{\rm NSE}}(\tilde{\rho },\tilde{T}_{\rm NSE},Y_e)] \nonumber \\ & &\quad= \epsilon (\rho,\Xi ^{-1}\tilde{T},{\bf X}) + \epsilon _{\rm nuc}(\mathbf {X}). \end{eqnarray} \tag{ 46 }$$

Note that since Ξ ⩾ 1, the estimated midplane disk temperature is higher than the temperature calculated without this correction, but this allows the disk to cool faster than it would otherwise. The rate of cooling by neutrino emission is then calculated from Equations (20) and (21) but at density $\tilde{\rho }$ and temperature $\tilde{T}$.

The initial model is the rotating M_star ≈ 14 M_☉ W-R star 16TI of Woosley & Heger (2006), evolved to pre-core-collapse from a 16 M_☉ main-sequence progenitor.³ To prepare the model 16TI, Woosley & Heger assumed that the rapidly rotating progenitor, which is near breakup at its surface at r_star ≈ 4 × 10⁵ km, had low initial metallicity, 0.01 Z_☉, and became a W-R star shortly after central hydrogen depletion, which implied an unusually small amount of mass loss. For illustration, the specific angular momentum at the three-quarters mass radius was ℓ_3/4 ∼ 8 × 10¹⁷ cm² s⁻¹, implying circularization around a 5 M_☉ black hole at r ∼ 2500 km, much larger than ISCO. The circularization radii of the outermost layers of the star are in the range 10⁴–10⁵ km. Woosley & Heger provide a radius-dependent angular momentum profile ℓ_16TI(r). We introduce the dimensionless parameter ξ_ℓ to scale the specific angular momentum ℓ(r) of our initial model relative to that of 16TI

$$\begin{equation} \ell (r) = \xi _{\ell } \,\ell _{\rm 16TI} (r). \end{equation} \tag{ 47 }$$

The plots of density, temperature, angular momentum, and composition in Section 3 show the initial conditions. The angular momentum profile is specific to our fiducial Run 1 with ξ_ℓ = 0.5, half of the rotation rate of 16TI.

The iron core of the model 16TI, with a mass ∼1 M_☉, has mass too low to collapse directly into a black hole, but should instead first collapse into a neutron star. The latter could, but need not, be driven to a successful explosion by the delayed neutrino mechanism. A black hole can form by fallback. We do not in any way account for the core bounce and its consequences, nor for the heating by the neutrinos emitted from the proto neutron star. Our central compact object is a point mass from the outset equipped with, as we clarify below, an absorbing boundary condition.

Pseudo-logarithmic gridding is achieved by capping the adaptive resolution at radius r with Δr > (1/8)ηr where η is a dimensionless parameter. We choose η = 0.15 for all but Run 2, where η = 0.075. Beyond the outer edge of the star we place a cold (10⁴ K), low-density, stellar-wind-like medium with density profile ρ(r) = 3 × 10⁻⁷ (r/r_star)⁻² g cm⁻³.

The simulation was carried out in the spherical domain r_min < r < r_max. We placed the inner boundary at r_min ∼ 25 km and the outer boundary well outside the star at r_max = 10⁷ km. In Table 1, we summarize the main parameters of our simulations, and also present some of the key measurements, defined in Section 3, characterizing the outcome of each simulation. Each simulation was run for ∼10⁷ hydrodynamic time steps and required ∼5000 CPU hours to complete.

Table 1. Summary of Simulation Parameters

Run	Δrmin	αb	ξℓc	ξconvd	ξconv, mixe
	(km)a
1	10.0	0.1	0.5	2.0	6.0
2f	$\mathbf {0.5}$	0.1	0.5	2.0	6.0
3	10.0	$\mathbf {0.2}$	0.5	2.0	6.0
4	10.0	$\mathbf {0.025}$	0.5	2.0	6.0
5	10.0	0.1	0.5	$\mathbf {5.0}$	$\mathbf {15.0}$
6	10.0	0.1	0.5	$\mathbf {0.5}$	6.0
7	10.0	0.1	0.5	$\mathbf {1.0}$	6.0
8	10.0	0.1	$\mathbf {0.25}$	2.0	6.0
9	10.0	0.1	0.5	2.0	$\mathbf {3.0}$

Notes. ^aThe minimum resolution element size. ^bThe dimensionless viscous stress-to-pressure ratio. ^cRotational profile parameter (see Equation (47)). ^dConvective efficiency parameter (see Equation (27)). ^eConvective compositional mixing efficiency parameter (see Equation (34)). ^fThis run also had additional angular resolution (see Section 2.7).

Download table as:  ASCIITypeset image

The boundary condition at r_min was unidirectional "outflow" that allowed free flow from larger to smaller radii (v_r < 0) and disallowed flow from smaller to larger radii (v_r > 0) by imposing a reflecting boundary condition. We imposed the torque-free boundary condition via (see, e.g., Zimmerman et al. 2005)

$$\begin{equation} \frac{\partial }{\partial r}\left(\frac{\ell }{r^2}\right)_{r=r_{\rm min}} = 0. \end{equation} \tag{ 48 }$$

As in other Eulerian codes, the boundary conditions in FLASH are set by assigning values to fluid variables in rows of "guard" cells just outside the boundary of the simulated domain. Let r_1/2 denote the leftmost cell within the simulated domain, and let $r_{\mathcal G}$ where ${\mathcal G}=(-{7}/{2},-{5}/{2},-{3}/{2},-{1}/{2})$ be the four guard cells to the left of r_1/2 such that the grid separation corresponds to $\Delta {\mathcal G}=1$. The torque-free boundary condition, if assumed to apply for r ⩽ r_min, implies $\ell _{\mathcal G}/r_{\mathcal G}^2=\ell _{1/2}/r_{1/2}^2$. All other fluid variables X were simply copied into the guard cells, $X_{\mathcal G}=X_{1/2}$, and were subsequently rendered thermodynamically consistent. This simple prescription approximates free inflow (toward smaller r) across r_min. The guard cell values for other fluid variables are assigned ignoring curvature of the coordinate mesh and formally violate conservation laws at r < r_min.

The mass of the black hole M_BH was initialized with the mass of the initial stellar model contained within r_min. The black hole mass was evolved by integrating the mass crossing the boundary at r = r_min,

$$\begin{equation} \frac{dM_{\rm BH}}{dt}= (- 4 \pi r^2 \rho v_r)_{r=r_{\rm min}}. \end{equation} \tag{ 49 }$$

The sum of the mass of the black hole and the mass contained on the computational grid remains constant to a high level of precision throughout each simulation.

We conducted tests of internal energy conservation, angular momentum transport, and spatial resolution convergence. The time-integrated equation for the conservation of internal energy in absence of nuclear and thermal energy interconversion in spherical coordinates reads

$$\begin{eqnarray} &&E_{\rm int, tot}(t_{\rm max}) - E_{\rm int, tot}(t_{\rm min}) = 4 \pi r_{\rm min,test}^2 \int _{t_{\rm min}}^{t_{\rm max}} (v_r \rho \epsilon + F_{\rm conv}) dt \nonumber \\ &&\quad-4\pi \int _{t_{\rm min}}^{t_{\rm max}} \int _{r_{\rm min,test}}^{r_{\rm max}} \left[ \frac{P}{r^2}\frac{\partial (r^2v_r) }{\partial r}- Q_{\rm visc}+ Q_{\nu } \right] r^2 dr dt =0,\nonumber\\ \end{eqnarray} \tag{ 50 }$$

where

$$\begin{equation} E_{\rm int, tot} = \int _{r_{\rm min,test}}^{r_{\rm max}}{4 \pi r^2 \rho \epsilon dr} \end{equation} \tag{ 51 }$$

and r_{min, test} ⩾ r_min is a reference radius defining the inner boundary of the spherical annulus in which we test energy conservation. We have ignored any flow of energy through r_max, since stellar material does not reach this radius in the course of any simulation.

In Figure 1, we utilize Equation (50) to test the global conservation of internal energy in a run with identical test parameters to Run 2, except that we disabled nuclear compositional change and used a Newtonian gravitational potential without relativistic corrections. In the legend, the apparent error ΔE is defined as the absolute value of the difference between the left- and right-hand sides of Equation (50). The evaluation of the various terms in Equation (50) was carried out in post-processing from cell-centered data recorded in Δt = 0.01 s intervals, which, in retrospect, is prone to the introduction of various spatial and temporal discretization artifacts not present in the actual simulation.

Zoom In Zoom Out Reset image size

We find that the apparent error is <1% of the largest term in Equation (50) when calculated for r_{min, test} = 50 km for the time interval 0 s ⩽ t ⩽ 70 s. The apparent error is most significant, ΔE ∼ 4 × 10⁵⁰ erg, prior to and during the first few seconds after shock passage. The apparent error that accrues after the first few seconds following shock passage is less than 10⁵⁰ erg. This can be compared to the total binding energy change on the simulation grid, which if sufficiently large can imply a supernova. We calculate this energy in Section 3.5 below and find that it increases by ∼(1.5–2) × 10⁵¹ erg following shock formation and a significant fraction (∼5 × 10⁵⁰ erg) of the increase is accrued later than a few seconds after shock formation, when the change in the cumulative apparent error is very small. Therefore, it does not seem that the apparent energy conservation error at the levels seen in the simulations should significantly impact the prospects for explosion. We note that in our calculations we explicitly transport specific internal energy rather than the total energy, by setting the parameter eintSwitch to a very large value. We would like to reiterate that it is likely that the apparent error is an artifact of post-processing and the true energy conservation is better. To demonstrate the latter, however, one would have to reconstruct the diagnostic energy fluxes using the very same interpolation procedure as is performed within the PPM in FLASH. We anticipate carrying out such a test in an extension of this work.

In steady state accretion, mass accretion associated with viscous angular momentum transport should occur at the rate

$$\begin{equation} \dot{M}_{\rm s.s.} = - 4 \pi \left(\frac{\partial \ell }{\partial r} \right)^{-1} \frac{\partial }{\partial r} \left(r^4 \nu \rho \frac{\partial \Omega }{\partial r} \right). \end{equation} \tag{ 52 }$$

In Figure 2 we show $\dot{M}(r)$ for Run 2 at t = 50 s along with the analytic steady-state estimate of Equation (52). In the rotationally supported region 100 km ≲ r ≲ 1000 km, the deviation from the analytical value is ≲ 5%, which lends credence to the accuracy of our angular momentum transport scheme.

Zoom In Zoom Out Reset image size

Our spatial resolution was chosen such that we resolve the innermost, neutrino-cooled region of the disk over several zones. One caveat is that we do not resolve the sonic radius of the flow, an issue discussed in McKinney & Gammie (2002). Because we use a torque-free boundary condition, the calculation is not subject to spurious viscous dissipation at the inner boundary. However, our boundary condition may still influence the fluid flow, and it therefore may be more apt to consider values of $\dot{M}$, as opposed to, for example, α or ℓ, when comparing our work to other simulations.

Run 1 and Run 2 contained identical hydrodynamic parameters and differed only in spatial resolution. Run 2 was capable of one additional level of resolution refinement over Run 1 and the parameter η described in Section 2.7 was set to one-half the value in Run 1, allowing for significantly higher resolution as a function of radius. Figures 3–6 show the density, temperature, specific entropy, and the mean atomic weight in Run 1, and also in the higher-resolution Run 2 at different times. Substantial agreement is seen between the low- and high-resolution simulations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The primary limitations of the model of collapsar accretion that we have presented here include: (1) a very approximate one-dimensional treatment of the intrinsically two- and three-dimensional flow structures; (2) limited adequacy of the Navier–Stokes viscous stress as a model for the magnetic stress arising in the nonlinear development of the MRI; (3) no a priori knowledge of the expected efficiency of convective energy transport ξ_conv and of compositional transport ξ_{conv, mix}; (4) treatment of nuclear compositional transformation through relaxation to NSE rather than by integrating the nuclear reaction network that would have allowed us to make predictions about the nucleosynthetic output; (5) neglect of ambient compositional stratification and nuclear compositional transformation inside convective cells in the calculation of the convective heat flux; (6) the lack of modeling of the axial relativistic jet and its enveloping cocoon that are thought to be present in LRGB sources; and (7) the use of a relatively low mass progenitor star, which may or may not be able to yield a black hole and an explosion with an energy as high as has been inferred in supernovae associated with LGRBs. Overcoming limitations (1) through (6) will require much more computationally expensive multidimensional hydrodynamic and MHD simulations. Limitation (7) can be addressed by applying our current method to other, more massive stellar models; here, we speculate what collapsars in higher mass progenitors may behave like in Section 4 below.

It would be tempting in view of limitation (5) to try to incorporate the effects of compositional stratification ambient to convective cells in the convective energy flux, which is normally achieved by multiplying the energy flux in Equation (27) with a factor

$$\begin{equation} \left[1-\left(\frac{dT}{d\mu }\right)_{P,\rho }\frac{d\mu }{dr} \bigg/ \left(\frac{T}{c_p} \frac{ds}{dr}\right) \right]^{1/2}, \end{equation} \tag{ 53 }$$

where μ is the mean nuclear mass, and utilizing the Ledoux instead of the Schwarzschild criterion (see, e.g., Bisnovatyi-Kogan 2001). This would be meaningful as long as the convective eddy turnover time τ_conv ∼ λ_conv/v_conv were shorter than the nuclear timescale τ_nuc ≲ τ_NSE, so that the convective cells can be treated as adiabatic before they mix. However, at radii where Ledoux convection would differ most from Schwarzschild convection, namely, where the photodissociation into helium nuclei and free nucleons is substantial, the convective timescale is much longer than the nuclear timescale, τ_conv ≫ τ_nuc. The internal composition of a convective cell evolves as it rises, and the associated entropy change is a much stronger effect than the variation of the ambient composition treated in Ledoux convection. Magnetization of the medium may play a role in this regime but its effects are poorly understood. Research into the interplay of convection and nuclear burning is ongoing (see, e.g., Arnett & Meakin 2011). Cognizant of these and other limitations, we adopt the Schwarzschild model and consider it but a parameterization of complex, still-to-be-explored physics.

Each simulation exhibits unshocked radial accretion of the inner, low-angular momentum mass shells of the progenitor star through the inner boundary lasting ∼(20–30) s at relatively steady accretion rates of ∼(0.1–0.2) M_☉ s⁻¹. The central mass accretion rate, black hole mass, and total mass on the computational grid as a function of time in each simulation are shown in Figure 7. The abrupt drop of the central accretion rate at ∼(20–30) s is associated with the appearance of an accretion shock precipitated by the arrival of the mass shells with specific angular momentum sufficient to lead to circularization around the black hole.

Zoom In Zoom Out Reset image size

In Figure 8, we show the location of the shock r_shock and its velocity v_shock ≡ dr_shock/dt as a function of time. For each of the runs, we identify the time when the shock first reaches radius 10 r_min = 250 km as the shock formation time t_shock and list the shock formation times in Column 2 of Table 2. We also provide the mass of the black hole at this point, M_BH(t_shock), in Column 8. The black hole mass at the time of shock formation was M_BH(t_shock) ∼ (5.2–5.5) M_☉ in Runs 1–7 and 9. In Run 8, which was initiated with reduced initial angular momentum, the accretion shock appeared later and the black hole mass is correspondingly larger.

Zoom In Zoom Out Reset image size

After the formation of the shock, the fluid nearest the inner boundary is rotationally supported and accretes as a result of angular momentum transport driven by the viscous shear stress. Subsequent to shock formation, the accretion rate declines rapidly either promptly or following a short delay. The typical rapid drop of the accretion rate is by a factor ∼5–10 (Runs 1, 2, 3, 5, 7, and 9), and this is followed by a continued power-law-like decline. By the end of each simulation at 100 s, the accretion rate has typically declined to ∼(10⁻³ − 10⁻⁴) M_☉ s⁻¹ (Runs 1, 2, 3, 4, 7, 8, and 9) or a factor of 100–1000 of the pre-shock value. Final black hole masses were ∼(6–7) M_☉ in the simulations with ξ_ℓ = 0.5 and ∼10 M_☉ in Run 8 with reduced initial angular momentum ξ_ℓ = 0.25.

In Run 4, with a low value of the viscosity parameter α = 0.025, the shock first made a very slow progress from 300 km to 2000 km during the first 10 s from its appearance. Then, at 30 s, the shock suddenly accelerated to v_shock ∼ 5000 kms⁻¹. The near-stagnation of the shock can be understood by noticing that during the 10 s, the neutrino cooling rate matches the viscous heating rate; the rapid cooling prevents the central entropy rise and convection seen in all other runs (see Section 3.3 below). In Section 4.2 of Lindner et al. (2010), we discussed the scenario in which the shock stagnation brought about by efficient neutrino cooling prolongs the LGRB central engine activity resulting in a longer prompt emission.

In Run 6, which had convective efficiency ξ_conv = 0.5, the shock stalled at the radius ∼10⁴ km for ∼5 s before proceeding outward. Note that the reinvigoration of the shock is solely driven by the convective energy transport, as we do not simulate the negligible neutrino energy and momentum deposition. The stalling and restarting of the shock was reflected in a strong variability of the central accretion rate.

Shock passage leaves a shock- and convection-heated, pressure supported envelope which contains much more mass than the disk, consistent with what we saw in Lindner et al. (2010). Figure 3 shows that the density in the envelope is an approximate power law of radius ρ∝r^−0.9. Figure 4 indicates that the temperature is also a power-law T∝r^−0.4. The pressure (not shown) is an approximate power law P∝r^−1.8. The profiles extend inward into the regime in which rotational support dominates pressure support. The mass of the rotationally supported material in the grid, where a_cent > (1/2)|g_self + g_BH|, promptly following disk formation was typically ≲ 5% of the total mass on the grid. Most of the mass on the grid was in the pressure supported atmosphere seamlessly connecting to the disk. The mass of the disk in each simulation is shown in Figure 7. In some of the runs, certain variability is seen in the disk mass over the first few seconds of disk formation. Afterward, the disk mass in each simulation declines monotonically. In most runs (1–4 and 7–9) the disk mass declines to M_disk ≲ 10⁻⁵ M_☉ by the end of the simulation, while in Run 6, the mass at the end of the simulation is somewhat larger but still very small, M_disk ∼ 3 × 10⁻⁴ M_☉.

Specific angular momentum as a function of radius is shown in Figure 9. In the initial angular momentum profile of the model, compositional boundaries coincide with discontinuities in the profile, but in 16TI these occur only at mass coordinates that are accreted directly onto the black hole, prior to the initial circularization. The angular momentum profile of the mass shells remaining at initial circularization is monotonically increasing and most of the remaining mass has nearly the same angular momentum, ∼(1–2) × 10¹⁷ cm² s⁻¹.⁴ This implies that the shocked atmosphere has nearly uniform specific angular momentum everywhere except at the radii where the timescale on which the viscous torque transport angular momentum is shorter than the time since circularization. At (25–50) s, there is a mild, sub-Keplerian inward downturn in ℓ(r) at r ≲ 1000 km. Angular momentum transport is too slow within the initial ∼100 s to affect the radii ≳ 10⁴ km.

Zoom In Zoom Out Reset image size

The mass accretion rate as a function of radius in Run 1 at t = 50 s is shown in Figure 2. The accretion rate is independent of radius for r ≲ 2000 km, which is the radii where the angular momentum profile has relaxed to a viscous quasi-equilibrium. The analytic expectation, given in Equation (52), is shown as well. Figure 10 shows the radial velocity v_r, angular velocity v_ϕ, and Keplerian velocity as a function of radius at t = 18 s, just as material begins to circularize outside of the black hole, and at t = 30 s, after an accretion disk has formed. At t = 18 s, the velocity v_ϕ reaches the Keplerian value at the innermost radii.

Zoom In Zoom Out Reset image size

Throughout the simulations, we tracked the value of our estimate of the vertical (z-directed) pressure scale height-to-radius ratio $\tilde{H}_z/r$, as described in Section 2.6. When the estimated ratio is below one-half, this indicates that in two dimensions, the flow should be disk-like, and when $\tilde{H}_z/r\ll 0.5$, the flow is a geometrically thin disk. We found that $\tilde{H}_z/r$ is below 0.5 but is still always above a minimum of 0.3 everywhere, except in Run 4, which had the lowest viscosity. The disk-like radii where the vertical pressure scale height-to-radius ratio is below one-half are r ≲ 200 km immediately following circularization and shrink to r ≲ 100 km by the end of the simulation. In Run 4 with a reduced viscous stress-to-pressure ratio α, neutrino cooling drove the disk to be geometrically thin, where H_z/r ≲ 0.3 in the inner r ≲ 500 km. In the innermost zone in Run 4, H_z/r = 0.1 at t = 20 s, the lowest seen in any simulation. By t = 35 s, no thin disk is present. In Figure 11 we show the value of Ξ_mod defined in Equation (44) throughout the simulation in Runs 1, 4, and 6; it does not drop below ∼0.77. Only in Run 4 is a genuinely thin accretion disk present, and there it is limited to small radii. The outer radius of the thin disk decreases as the neutrino luminosity drops (see Section 3.3). We attribute the observed moderate thinning of the accretion flow to the cooling of the flow by the photodisintegration of helium nuclei into free nucleons, and in Run 4, the additional contribution of neutrino cooling is also significant.

Zoom In Zoom Out Reset image size

To understand the energetics of the accretion flow in a collapsar, we need to consider the transport of mechanical, thermal, and nuclear binding energy, as well as the loss to neutrino emission. Before turning to energy transport, we discuss the neutrino losses, which turn out to be not significant in the regime we consider.

The integrated neutrino luminosity is dominated by the emission from the inner ∼100 km. The luminosity as a function of time in the representative Runs 1, 4, and 6 is shown in Figure 12. In simulations with α = 0.1, neutrino luminosities integrated over the entire computational domain peaked immediately following shock formation at ∼(1–200) × 10⁵¹ erg s⁻¹. The peak luminosity lasted anywhere from less than a second in the runs with high peak luminosities to a few seconds in the runs with low peak luminosities. After the peak, the luminosity decays first very rapidly until it has dropped to ∼10⁵⁰ erg s⁻¹, and then continues to decay approximately exponentially by several orders of magnitude to settle at ∼(10⁻⁶ to 10⁻⁵) erg s⁻¹ after ∼50 s. The sharp luminosity peak is an artifact of the abrupt nature of shock formation in our 1.5D dimensional treatment and is probably not physically significant. The total energy that would be deposited by an absorption of the emitted neutrinos, which we do not calculate, is negligible.

Zoom In Zoom Out Reset image size

Now turning to energy transport, we examine the radial transport of all forms of energy, the thermal and kinetic energies, the nuclear binding energy, and the gravitational potential energy. The gravitational potential energy is a nonlocal functional of the mass distribution. However, ignoring relativistic effects, one can define the gravitational potential energy per unit volume to be ρ(Φ_BH + (1/2)Φ_s), where Φ_BH is the gravitational potential of the black hole which we define via Φ_BH(r) ≡ ∫^∞_rg_BH(r')dr' with g_BH given in Equation (9) and Φ_self is that of the self-gravity of the star. Then, ρv_rΦ, where Φ = Φ_BH + Φ_self, can be interpreted as the flux of gravitational energy advected by the fluid, but one must additionally include the flux of gravitational energy transported by self-gravity (see, e.g., Binney & Tremaine 2008, their Appendix F), which equals

$$\begin{equation} F_{\rm grav,self}=\frac{1}{8\pi G} \left(\Phi _{\rm self}\nabla \frac{\partial \Phi _{\rm self}}{\partial t}-\frac{\partial \Phi _{\rm self}}{\partial t}\nabla \Phi _{\rm self}\right). \end{equation} \tag{ 54 }$$

This term is significant only in the outer envelope of the star. The rate with which the sum total of these energies is transported radially is given by

$$\begin{eqnarray} \dot{E} &=& 4 \pi r^2\left[\rho v_r \left(\epsilon + \epsilon _{\rm nuc}+ \frac{P}{\rho } + \frac{1}{2} v_r^2 + \frac{1}{2}\frac{\ell ^2}{r^2} + \Phi \right) + F_{\rm grav,self}\right. \nonumber \\ &&-\left. \rho \nu \ell \frac{\partial \Omega }{\partial r} + F_{\rm conv}+F_{\rm nuc,mix}\right], \end{eqnarray} \tag{ 55 }$$

where the convention is such that $\dot{E}>0$ implies the transport of positive energy outward, opposite from the convection employed in the definition of the mass accretion rate $\dot{M}$. Here, −ρνℓ∂Ω/∂r is flux of energy transported by the viscous stress.

Specific entropy as a function of radius at several times in Runs 1 and 2 is shown in Figure 5. After the formation of the accretion shock and the rotationally supported flow, viscous dissipation heats the fluid, thus producing a negative entropy gradient in the shock downstream. The negative specific entropy gradient extends almost to the shock front, and thus the energy injected at small radii can travel to raise the entropy of the entire post-shock region. Figure 13 shows the specific entropy in Run 4. Here, the high neutrino luminosity after the accretion shock has formed keeps the entropy in the post shock region relatively low. For the first ∼10 s after shock formation, no specific entropy inversion is seen, and the fluid is stable against convection. When the neutrino luminosity begins to drop around t ≈ 33 s, the entropy rises, a negative specific entropy gradient appears in the post-shock region, and convection starts transporting the viscously dissipated energy outward.

Zoom In Zoom Out Reset image size

Figure 14 plots the net transport rate $\dot{E}$ and the various constituent terms in Run 2 at t = 30 s; the radii and other observables quoted in the remainder of this section will be specific to this particular simulation snapshot and will vary across different simulations and different times within a simulation. Approximate radial independence of the energy transport rate, $\partial \dot{E}/\partial r\approx 0$ for 200 km ≲ r ≲ 4000 km, where the transport rate is positive $\dot{E}\approx 10^{50}\,{\rm erg}\,{\rm s}^{-1}>0$, is indicative of quasi-steady-state accretion. At larger radii, r ≳ 5000 km, where the inner inflow gives way to an outer outflow—a precursor of the brewing explosion—no quasi-steady state is present and the fluid variables evolve on the dynamical time in the wake of the expanding shock. At small radii, r ≲ 100 km, where one expects a steady state, the curve $\dot{E}(r)$ exhibits a small positive gradient, as well as a sawtooth consistent with that seen in the accretion rate $\dot{M}(r)$. The constancy of the plotted energy transport rate is contingent on an accurate cancellation of the other transport terms. We suspect that the observed nonconstancy is arising from relatively small inconsistencies in the discretization or gravitational source terms in FLASH and in the calculation of the gravitational energy during post-processing.

Zoom In Zoom Out Reset image size

At r ≲ 1000 km, the inward advection of thermal and kinetic energy dominates over the outward transport by convection and the viscous stress. Therefore, the innermost flow is an advection-dominated accretion flow (ADAF; Narayan & Yi 1994, 1995; Blandford & Begelman 1999). At r ≳ 1000 km, the outward transport of thermal energy by convection dominates the inward transport by advection and this region is thus a convection-dominated accretion flow (CDAF, see, e.g., Stone et al. 1999; Igumenshchev et al. 2000; Blandford & Begelman 2004). Inward nuclear binding energy advection and nuclear compositional mixing both act to transport the total energy outward if one counts the negative nuclear binding energy in the total energy budget. Convection transports energy from the ADAF–CDAF transition radius where the magnitude of the enthalpy advection flux ∼|v_r(ρ + P)| equals the convection flux F_conv to the shock radius r_shock. In Figure 14 the former is at r_ADAF ≈ 1000 km and the latter is at r_shock ≈ 3.5 × 10⁴ km. In Run 1 and similar runs, r_ADAF increases very slowly from ∼1000 km to ∼2000 km from shock formation until t = 100 s. In Run 4, the radius is ∼200 km throughout the simulation, and in Run 6, the radius grows from ∼5000 km to over ∼10⁴ km in the course of the simulation.

We suspect that the location of r_ADAF determines the amount of energy that can be carried to the shock front, and that in turn, the ADAF–CDAF transition radius is primarily a function of the convective efficiency ξ_conv and the viscous stress-to-pressure ratio α. Simulations with larger values of ξ_conv resulted in stronger shocks and larger amounts of unbound material. The convective compositional mixing parameter ξ_{conv, mix} had little effect on the final outcome of our simulations. Turning to the viscosity parameter, Run 3 with a large values of α produced a somewhat less energetic shock with less unbound material at the end of the simulation. Run 4, the simulation with the lowest viscosity, produced the most energetic explosion, even in the presence of more pervasive cooling by neutrino emission and by photodisintegration in the low-α regime; these flows are denser and hotter (see, e.g., Popham et al. 1999; Chen & Beloborodov 2007). In Milosavljević et al. (2012), we show that r_ADAF is expected to be smaller under low viscosity conditions because the advective luminosity is proportional to v_r, which is proportional to α. This trend is reproduced in the present simulations.

Our simplified treatment of nuclear compositional transformation, which entails relaxation to NSE on a temperature-dependent timescale, is designed to track the impact of nuclear photodisintegration and recombination on the thermodynamics of the flow. However, it does not allow the computation of the ultimate nucleosynthetic output in the presence of out-of-NSE burning. Thus, the results presented here can only be understood in light of the limitations of the method.⁵ It is also worth recalling that we do not calculate neutronization that could modify the proton-to-nucleon ratio Y_e.

In the hottest, innermost accretion flow, photodisintegration of heavy nuclei saps energy that could otherwise be transported by convection to larger radii to energize the shocked envelope. However, once the nuclei are broken down, convective mixing can dredge up free nucleons to larger and cooler radii, where they can recombine and heat the fluid locally. Figure 6 shows the mean atomic mass $\bar{A}$ as a function of radius at various times in Run 1 and Run 2. The mean atomic mass drops below 4 in the innermost (200–300) km. The positive gradient in $\bar{A}$ seen in portions of the convective region would in the Ledoux picture enhance the convective energy flux, but our Schwarzchild treatment of convection does not capture this effect. We argue in Section 2.9 that since the nuclear timescale is shorter than the convective timescale at radii where Ledoux convection implies an enhanced energy transport, the nuclear compositional transformation inside convective cells, not considered in the Ledoux treatment, should dominate. Lacking a theory of convection in this regime we adhere to the simpler Schwarzschild parameterization.

In Figure 15, we show the mass-weighted abundances of the most common elements in our simulations in Run 1 at various times. By t = 30 s, again, the inner 200 km is made up almost entirely of free nucleons in nearly equal portions, as Y_e ≈ 0.5 everywhere. The effect of convective compositional transport of the reprocessed nuclear species—the free nucleons, helium, and iron—from the hot innermost accretion flow is seen in the power-law tails extending to near the location of the shock in the right panels.

Zoom In Zoom Out Reset image size

Although in our calculations we do not allow the evolution of Y_e, we can still speculate about the effects of neutronization. Chen & Beloborodov (2007) computed the structure of time-independent accretion disks around Kerr black holes including the effects of pair capture and neutronization. In their models with α = 0.1, the same as our fiducial viscous stress-to-pressure ratio, at the radii where ρ ∼ 10⁷ g cm⁻³, corresponding to the density in the innermost disk in our simulations, they find Y_e ≈ 0.5, the same as in our non-neutronizing treatment. In their models with α = 0.01, however, Chen & Beloborodov (2007) find that at densities corresponding to the innermost disk in our simulations, significant neutronization was in effect and Y_e dropped well below neutron–proton equality. It is therefore possible that in the very innermost regions of the disks in our Run 4 with a low viscosity α = 0.025, the true value of Y_e should be lower than we assume. This would modify the abundances and thermodynamics of the portion of the flow in NSE. The key question of consequence for the viability of the mechanism we propose for the production of luminous supernovae is, will the neutron-rich material pollute larger radii and drive a tendency toward the synthesis of iron instead of ⁵⁶Ni? Because this neutronization only seems to be most significant in the hottest innermost regions, where neutrino cooling is efficient and the flow is predominantly rotationally supported, it is possible that most of the neutron-enhanced material would be advected into the black hole. This is a quantitative question that can be answered only with multidimensional simulations.

In Figure 16, we show the location of r_NSE, the largest radius at which nuclear compositional transformation, in the form of gradual convergence toward NSE, is taking place in our Run 1. This is the only region in which our calculations will capture nucleosynthesis and photodisintegration. Effectively, it is the radius at which T > 3 × 10⁹ K. Early in the simulation, the hot stellar core accretes into the black hole, and r_NSE quickly drops from ≈3000 km to ≈800 km. After the shock forms, additional heating in the shock and by viscous dissipation rapidly increases r_NSE and it peaks at ≈4000 km. As the density and temperature drop and the viscous heating rate decreases, the inner regions of the star begin to cool, and r_NSE declines again.

Zoom In Zoom Out Reset image size

In Figure 17, we show the total mass of various nuclear species in the entire simulation domain as a function of time. The most notable change is the dip in the mass of iron-group elements. The initial decrease in the iron-group mass is the accretion of the core of the star onto the black hole. After shock formation, a rapid increase in the amount of free nucleons is seen, in addition to production of additional helium and iron-group elements. In Table 2, we show the total amount of newly fused Fe-group elements present at the end of the simulation, which fall in the range of 0.02 M_☉–0.09 M_☉. Since we do not calculate the out-of-NSE burning in convectively dredged up material, at least a fraction of the extended helium tail seen in Figure 15 could be expected to burn into iron and thus the iron-group mass in Figure 17 and Table 2 can be interpreted as a lower limit.

Zoom In Zoom Out Reset image size

The total thermal and mechanical energy, which we also refer to as the total binding energy, present on the grid was computed for each simulation via

$$\begin{equation} E_{\rm bind} = \int _{r_{\rm min}}^{r_{\rm max}}{ \rho \left(\epsilon + \frac{1}{2}v_r^2 + \frac{1}{2}\frac{\ell ^2}{r^2} +\Phi _{\rm BH} +\frac{1}{2}\Phi _{\rm self}\right) 4\pi r^2dr} \end{equation} \tag{ 56 }$$

and is shown in Figure 18. A positive binding energy indicates a potential for explosion. Runs 1, 2, 3, 4, 5, 8, and 9 acquired a positive total binding energy E_bind ∼ (3.5–6.2) × 10⁵⁰ erg by the end of the simulation, indicating the potential for explosion. These runs have ξ_conv ⩾ 2 in common. Run 7 with ξ_conv = 1 reached marginally unbound condition with E_bind ∼ (0.5–1) × 10⁵⁰ erg. Run 6 with ξ_conv = 0.5 remained gravitationally bound throughout the entire simulation.

Zoom In Zoom Out Reset image size

The total unbound mass in each simulation was calculated by summing the masses of any fluid element with a positive value of E_bind and a positive radial velocity. The total unbound masses in each simulation defined by this diagnostic are shown in Table 2. This criterion does not take into account any interaction that unbound and bound material might have subsequent to the measurement. This criterion also neglects future energy gain or loss from nuclear processes. The location and velocity of the outward moving accretion shock is shown in Figure 8. In exploding models, the typical shock velocities were (2000–4000) km s⁻¹. In Run 4 and Run 6, the shock stalled or slowed, and later restarted once or several times. In Figure 19, we show the evolution of various Lagrangian mass coordinates in Run 1 throughout the simulation. Once the shock moves beyond ∼4000 km, the infalling material obtains a positive velocity once it reaches r_shock.

Zoom In Zoom Out Reset image size

These results suggest that a high convective efficiency is required for sufficient transfer of energy from the inner accretion flow to the envelope to unbind the envelope. Simulations with higher values of ξ_conv had relatively larger shock velocities and amounts of unbound mass at the ends of their simulations, and Run 4 with a lower α had the largest value of E_bind. Run 8 with reduced initial specific angular momentum produced an explosion comparable to that seen in the fiducial model.