# Numerical Simulations of the Random Angular Momentum in Convection:
Implications for Supergiant Collapse to Form Black Holes

Andrea Antoni1 and Eliot Quataert1,2
1Astronomy Department and Theoretical Astrophysics Center, University of
California, Berkeley, CA 94720, USA
2Department of Astrophysical Sciences, Princeton University, Princeton, NJ
08544, USA
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

During the core collapse of massive stars that do not undergo a canonical
energetic explosion, some of the hydrogen envelope of a red supergiant (RSG)
progenitor may infall onto the newborn black hole (BH). Within the Athena++
framework, we perform three-dimensional, hydrodynamical simulations of
idealized models of supergiant convection and collapse in order to assess
whether the infall of the convective envelope can give rise to rotationally-
supported material, even if the star has zero angular momentum overall. Our
dimensionless, polytropic models are applicable to the optically-thick
hydrogen envelope of non-rotating RSGs and cover a factor of 20 in stellar
radius. At all radii, the specific angular momentum due to random convective
flows implies associated circularization radii of 10 - 1500 times the
innermost stable circular orbit of the BH. During collapse, the angular
momentum vector of the convective flows is approximately conserved and is
slowly varying on the timescale relevant to forming disks at small radii. Our
results indicate that otherwise failed explosions of RSGs lead to the
formation of rotationally-supported flows that are capable of driving outflows
to large radii and powering observable transients. When the BH is able to
accrete most of the hydrogen envelope, the final BH spin parameter is $\sim$
0.5, even though the star is non-rotating. For fractional accretion of the
envelope, the spin parameter is generally lower and never exceeds 0.8. We
discuss the implications of our results for transients produced by RSG
collapse to a black hole.

###### keywords:

black hole physics – convection – stars: massive – supernovae: general.

††pubyear: 2021††pagerange: Numerical Simulations of the Random Angular
Momentum in Convection: Implications for Supergiant Collapse to Form Black
Holes–A

##  1 Introduction

It is plausible that in some fraction of massive stellar deaths, iron core
collapse does not lead to a successful neutrino-driven supernova (SN; Zhang et
al., 2008; O’Connor & Ott, 2011; Ugliano et al., 2012; Sukhbold et al., 2016).
Instead, core collapse results in a failed supernova (FSN) in which the proto-
neutron star (PNS) quickly collapses to a black hole (BH) and a large fraction
of the star accretes onto the newborn BH. The outcome of stellar core collapse
depends on the final pre-SN structure of the star in a complex way (O’Connor &
Ott, 2011; Ugliano et al., 2012; Ertl et al., 2016; Sukhbold & Adams, 2020)
which probably cannot be disentangled without a statistically significant
sample of (4$\pi$) three-dimensional (3D) simulations of both the pre-SN
(Couch et al., 2015; Yadav et al., 2020; Fields & Couch, 2020) and post-bounce
(Burrows et al., 2020; Powell & Müller, 2020) phases of core collapse. The
result of the complex relationship between the explodability of a massive star
and its pre-SN evolution is that explodability is likely neither a monotonic
nor a singled-valued function of zero-age main sequence (ZAMS) mass (Sukhbold
et al., 2018).

The angular momentum content of a collapsing star has important consequences
for the outcome of core collapse, including the properties of the remnant, the
amount of material returned to the stellar environment, and the transient
signal produced. Collapse of rotating Wolf Rayet (WR) stars in successful or
failed SN (i.e. Type II or Type I collapsars, respectively) may power some
long gamma-ray bursts (Woosley, 1993; MacFadyen & Woosley, 1999). Accretion of
high angular momentum material in more extended rotating stars, e.g. red and
blue supergiants (RSGs and BSGs, respectively), in weak or failed SN may be
responsible for ultra-long gamma-ray transients (Woosley & Heger, 2012;
Quataert & Kasen, 2012; Perna et al., 2018).

Even in the absence of net angular momentum, FSN can still give rise to
detectable transients. During the PNS phase, neutrino losses reduce the
gravitational mass of the PNS by $\sim 0.1-0.3$ $M_{\odot}$ over a few
seconds. The envelope of the star, which is suddenly over-pressured due to the
nearly instantaneous reduction of the gravitational potential, reacts on a
dynamical time. A sound pulse forms and steepens into a weak shock that may
unbind a portion of the envelope of the star (Nadyozhin, 1980; Lovegrove &
Woosley, 2013; Fernández et al., 2018). For RSG progenitors, the ejection of
several $M_{\odot}$ of the hydrogen envelope results in a $\approx$ few day
shock breakout followed by a $\approx 300-400$ day transient of $\approx
10^{39}-10^{40}$ erg s-1 powered by hydrogen recombination in the expanding
ejecta (Lovegrove & Woosley, 2013; Piro, 2013; Lovegrove et al., 2017;
Fernández et al., 2018). The brightening and subsequent disappearance of a RSG
progenitor in the FSN candidate event N6946-BH1 is broadly consistent with
this model (Gerke et al., 2015; Adams et al., 2017; Basinger et al., 2020;
Neustadt et al., 2021). The neutrino-mass-loss mechanism can also launch weak
shocks in yellow supergiants (YSGs), BSGs, and WRs, but the ejected masses are
much lower and some progenitors eject no unbound mass at all (Fernández et
al., 2018). For semi-analytic work related to mass ejection from weak shocks,
see Coughlin et al. (2018a); Coughlin et al. (2018b, 2019); Ro et al. (2019);
Linial et al. (2021).

A known limitation of Lovegrove & Woosley (2013) and Fernández et al. (2018)’s
calculations of how much neutrino losses occur prior to BH formation was the
use of analytic functions to model the neutrino emission. Ivanov & Fernández
(2021) improves on this earlier work by using general-relativistic neutrino
radiation-hydrodynamic simulations to model the evolution of the inner core of
each progenitor to BH-formation using three different equations-of-state
(EOSs) for the PNS. Self-consistently modeling the core collapse for a given
progenitor is important because the neutrino-loss function and time to BH
formation set the amount of neutrino energy radiated and thus the energy of
the outgoing weak shock. For their RSG progenitor Ivanov & Fernández find that
this more careful treatment of BH-formation and a softer equation of state
reduces the ejected mass by a factor of a few and means that $1-3$ $M_{\odot}$
of the hydrogen envelope remains bound and will accrete onto the newborn BH.
For comparison, the characteristic mass loss from their YSG and lower-mass BSG
progenitors is $\approx 10^{-2}M_{\odot}$ while the remaining BSG and WR
progenitors eject $\approx 10^{-4}M_{\odot}$ in unbound mass. This work
suggests that the infall of a large fraction of the hydrogen envelope is
common in FSN of RSGs and YSGs. The accretion of this material is the focus of
our study.

Rotation is not the only source of angular momentum in FSN progenitors. The
random velocity fields in convective zones have an associated angular momentum
because in each radial shell there is a mean ‘horizontal’ velocity
perpendicular to the radial axis even if the net angular momentum of the star
is zero. For each convective zone, the key question is whether the mean
specific angular momentum arising due to convection, $j_{\rm rand}$, is larger
than $j_{\rm ISCO}$, the Keplerian specific angular momentum associated with
the innermost stable circular orbit (ISCO) of the BH. Gilkis & Soker (2014)
first explored the angular momentum present in the convective regions of pre-
SN stars in the context of jet-driven SN (Soker, 2010; Papish & Soker, 2011).
They derived an analytical estimate for $j_{\rm rand}$ contained in a shell of
material with randomly oriented velocities of magnitude $v_{\rm c}$, where
$v_{\rm c}$ is the convective velocity. For a sample of pre-SN stars computed
with MESA (Paxton et al., 2011; Paxton et al., 2013, 2015, 2018, 2019), their
analytical estimate predicts significant angular momentum content in the
helium and hydrogen layers of their supergiant models (specifically $j_{\rm
rand}>10^{5}j_{\rm ISCO}$ in the hydrogen envelopes of their RSG and YSG
models). They argue that the fallback of the helium layer is sufficient to
drive a canonical ($>10^{51}$ erg) SN, thus removing the possibility of a FSN.

Quataert, Lecoanet & Coughlin (2019), on the other hand, derived an analytical
scaling for $j_{\rm rand}$ as a function of $v_{\rm c}$ that is at least two
orders-of-magnitude smaller than the result of Gilkis & Soker (2014). To
validate their expression for $j_{\rm rand}$, Quataert et al. performed
Boussinesq simulations of convection in a Cartesian slab of material, which
confirmed their analytic scaling within a factor of 2. Gilkis & Soker (2016)
carried out 3D simulations of the convective helium zones in massive stars;
their simulations are also consistent with Quataert et al.’s estimate (to
within a factor of a few). When applied to MESA models of a RSG and a YSG,
Quataert et al. found that $j_{\rm rand}$ is significant only in the hydrogen
envelope. That is, material interior to the base of the hydrogen envelope
could accrete spherically onto the BH, but the hydrogen envelope had
sufficient angular momentum due to random convective flows that the material
would circularize beyond the ISCO of the BH.

The aim of our work is to determine whether the collapse of a non-rotating
supergiant envelope in a FSN results in rotationally-supported material
outside the ISCO of the newborn BH. We perform two numerical experiments. In
the first, we simulate convection in a spherical polytropic envelope in order
to relate the convective velocity, $v_{\rm c}$, to the magnitude of the
specific angular momentum, $j_{\rm rand}$, arising from the random motion of
the convective fluid. In the second experiment, we follow the collapse of the
envelope in order to study the evolution of $j_{\rm rand}$ during infall and
to measure the time-dependent angular momentum vector of the material accreted
onto the BH. Our study thus extends the results of Quataert et al. (2019) to
the spherical geometry of a star and determines how well $j_{\rm rand}$ of the
convective shells is conserved during infall.

This paper is organized as follows. Section 2 derives the analytical scaling
for $j_{\rm rand}$ of Quataert et al. (2019). We describe our supergiant model
and simulation methods in Section 3. Sections 4 and 5, respectively, present
the results of our convection and collapse simulations. We place these results
in the context of supergiants and FSN in Section 6. Our summary and
conclusions are given in Section 7. Appendix A describes how accretion rates
in the collapse simulations can be predicted from snapshots of the flow prior
to collapse using a simple test problem (we apply this technique to convective
flows in Section 5.2).

##  2 Background

Quataert et al. (2019)’s analytical estimate for $j_{\rm rand}$ will be useful
in what follows, so we review the derivation here. Consider a spherical shell
of the convective envelope at radius $r$ and with thickness $H$, where $H$ is
the local pressure scale height. The velocities of the convective material
have components along $r$, which do not give rise to angular momentum, and
‘horizontal’ components (perpendicular to $r$) that do. These velocities are
randomly oriented so a sum over an infinite number of eddies should result in
zero net angular momentum. But the volume of the shell is finite so there is a
finite number of eddies at radius $r$, given approximately by

$N_{\rm edd}\sim\frac{4\pi r^{2}H}{H^{3}}\sim 4\pi\frac{r^{2}}{H^{2}}.$

The convective velocity, $v_{\rm c}$, and the magnitude of the mean velocity
vector are related by $\|\langle\vec{v}\rangle\|\sim v_{\rm c}/\sqrt{N_{\rm
edd}}$. So to a factor of order unity, the mean horizontal velocity is related
to the convective velocity by,

$v_{h}\approx\frac{v_{\rm c}}{\sqrt{4\pi}}\frac{H}{r}.$ (1)

This implies a specific angular momentum content of

$j_{\rm rand}\sim rv_{h}\sim\frac{v_{\rm c}H}{\sqrt{4\pi}}.$ (2)

Or, normalizing to typical values for the convective regions of supergiants,

$j_{\rm rand}\sim 6\times
10^{18}\,\bigg{(}\frac{r}{1000R_{\odot}}\bigg{)}\bigg{(}\frac{v_{\rm
c}}{10\,{\rm km}/{\rm s}}\bigg{)}\bigg{(}\frac{H/r}{0.3}\bigg{)}\,{\rm
cm}^{2}\,{\rm s}^{-1}$ (3)

For comparison,

$j_{\rm ISCO}=1.15-3.46\frac{GM_{\bullet}}{c}\approx 0.5-1.5\times
10^{17}\bigg{(}\frac{M_{\bullet}}{10M_{\odot}}\bigg{)}\,{\rm cm}^{2}\,{\rm
s}^{-1}$ (4)

where $M_{\bullet}$ is the mass of the BH. The lower limit corresponds to
maximum BH spin (and prograde orbits) and the upper limit corresponds to zero
BH spin.

The upper panel of fig. 3 of Quataert et al. (2019) plots eq. (2) for MESA
models of a $15M_{\odot}$ RSG and a $22M_{\odot}$ YSG. Material interior to
the hydrogen envelope has $j_{\rm rand}\ll j_{\rm ISCO}$, implying spherical
accretion of the helium shell. The hydrogen envelope of both the RSG and YSG
models have $1\lesssim j_{\rm rand}/j_{\rm ISCO}\lesssim 10$ and therefore may
have sufficient angular momentum to feed a rotationally-supported structure
outside the BH horizon.

We emphasize that, in the preceding picture, the total angular momentum of the
star is zero. Random velocities give rise to non-zero mean angular momentum in
a given shell even though the sum over the star is zero. The key question we
study in this paper is the random angular momentum flows in spherical
convective zones and how, and whether, the envelope material restructures
itself during infall.

##  3 Numerical Methods

We model the collapse of non-rotating, convective supergiant envelopes in 3D
using Athena++ 111Version 19.0, https://princetonuniversity.github.io/athena
(Stone et al., 2020), an Eulerian hydrodynamic code based on Athena (Stone et
al., 2008). As the star is non-rotating, there is no symmetry axis to motivate
the use of spherical-polar coordinates. We instead perform our simulations in
a Cartesian grid. All simulations use third-order Runge-Kutta time integration
(integrator = rk3), piecewise parabolic spatial reconstruction (xorder = 3),
and the Harten-Lax-van Leer contact Riemann solver (`--flux hllc`). Our
supergiant model initially has radial symmetry and we use nested static mesh
refinement to allow resolution to fall off approximately logarithmically with
radial distance from the center of the star.

We use a gamma-law EOS and neglect the self-gravity of the gas. We do not
include radiation transport, but we do include a simplified cooling function
at the surface of the star to remove heat from the rising convective fluid. We
include a heating source term at small radii to drive convection in the
envelope and we run our simulations until thermal equilibrium is reached. Our
simulations solve a model problem that is similar to RSGs in many key ways. We
have a convective polytrope with a density and pressure profile that is
representative of the hydrogen envelope of RSGs with roughly the right
convective velocities.

One of the main goals of this work is to study the collapse phase of FSN. In
particular, we aim to understand whether there is reshuffling of the angular
momentum during the infall or whether the angular momentum content of the
convective zone prior to collapse is a good proxy for the accreted angular
momentum during collapse. We use our simulations of convection as the initial
conditions for the collapse calculations and our main focus is on the change
(or not) in angular momentum content of the material from the convection zone
through the collapse.

###  3.1 Equations Solved

We solve the equations of inviscid hydrodynamics

$\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\mn@boldsymbol{v})=0$ (5)
$\frac{\partial(\rho{\mn@boldsymbol{v}})}{\partial
t}+\nabla\cdot(\rho\mn@boldsymbol{v}\mn@boldsymbol{v}+P\textbf{{I}})=-\rho\nabla\Phi+\rho\mn@boldsymbol{a}_{\rm
damp}$ (6) $\frac{\partial E}{\partial
t}+\nabla\cdot\big{[}(E+P)\mn@boldsymbol{v}\big{]}=-\rho\mn@boldsymbol{v}\cdot\nabla\Phi-
Q_{\rm cool}+Q_{\rm heat}+\rho\mn@boldsymbol{v}\cdot\mn@boldsymbol{a}_{\rm
damp}$ (7)

where $\rho$ is the mass density, $\rho\mn@boldsymbol{v}$ is the momentum
density, $E=\epsilon+\rho\mn@boldsymbol{v}\cdot\mn@boldsymbol{v}/2$ is the
total energy density, $\epsilon$ is the thermal energy density, $P$ is the gas
pressure, and I is the identity tensor. The source terms for gravity, cooling,
heating, and damping ($\Phi$, $Q_{\rm cool}$, $Q_{\rm heat}$, and
$\mn@boldsymbol{a}_{\rm damp}$, respectively) are defined in Sec. 3.3. We
adopt an ideal gas EOS

$\epsilon=\frac{P}{\gamma-1}$ (8)

with adiabatic index $\gamma$.

###  3.2 Supergiant Model

In this section we introduce our semi-analytic model for the supergiant
envelope and atmosphere. Section 3.5 will describe how this model is mapped to
Athena++ to initialize the grid.

As discussed in Section 2, we expect everything interior to the base of the
hydrogen envelope to have $j_{\rm rand}\ll j_{\rm ISCO}$ and to be accreted
onto the core at the start of collapse. We represent all of this mass, $M$,
with a Plummer potential

$\Phi(r)=-\frac{GM}{(r^{n}+a^{n})^{\frac{1}{n}}}$ (9)

with softening length $a$ and index $n$. We take $n=8$, so that $\Phi(r)$
converges to a point mass potential for $r\gtrsim 1.4a$.

Supergiant envelopes have roughly power-law density profiles with nearly
constant entropy, so for $r\gtrsim a$ we will model the hydrogen envelope with
a density profile of the form $\rho\propto r^{-b}$. For RSGs and YSGs,
$1.5\lesssim b\lesssim 2.5$ (see Coughlin et al., 2018b, figs. 7 and 12). For
an ideal gas, a power-law density profile gives the envelope a temperature
profile of $T\propto r^{-1}$. To mimic the stellar photosphere where cooling
is rapid, the temperature is smoothly dropped to an isothermal value for
$r\gtrsim r_{\rm cool}$. The resultant temperature profile has the functional
form

$T_{\rm init}(r)=\frac{T_{\rm env}(r)}{2}\big{[}1-F(r)\big{]}+\frac{T_{\rm
iso}}{2}\big{[}1+F(r)\big{]}$ (10)

where

$T_{\rm env}(r)=\frac{\mu
m_{p}}{k_{b}(b+1)}\frac{GM}{(r^{n}+a^{n})^{\frac{1}{n}}},$ (11) $T_{\rm
iso}=f_{\rm iso}T_{\rm env}(r_{\rm cool}+0.5r_{0}),$ (12)

and

$F(r)=\tanh\bigg{[}\frac{r-r_{\rm cool}-0.5r_{0}}{0.9r_{0}}\bigg{]}.$ (13)

In the last expression, $r_{0}$ is the characteristic length scale in our
model. In our code units, $r_{0}=1$ (Sec. 3.6) and our simulations adopt
$r_{\rm cool}=5r_{0}$, $b=2.1$, $f_{\rm iso}=0.3$, and $a\leq 0.16r_{0}$ (Sec.
3.7).

For $r\lesssim r_{\rm cool}$, $T_{\rm init}\sim T_{\rm env}$, which is the
power-law temperature profile in the envelope that results from assuming an
ideal gas EOS and choosing the density profile for $r\lesssim r_{\rm cool}$ to
take the form

$\rho_{\rm
env}(r)=\rho_{0}\big{[}(r_{0}^{n}+a^{n})/(r^{n}+a^{n})\big{]}^{\frac{b}{n}},$
(14)

where $\rho_{0}$ is the density at $r_{0}$ and $\rho\propto r^{-b}$ for
$r\gtrsim a$. By contrast, for $r\gtrsim r_{\rm cool}$, $T_{\rm init}\sim
T_{\rm iso}=$ constant.

The constants in eqs. (12) and (13) set the location and width of the
transition region between the envelope temperature profile, $T_{\rm env}(r)$,
and the isothermal temperature, $T_{\rm iso}$. The shift of $0.5r_{0}$ ensures
that the transition region does not begin until $\sim r_{\rm cool}$ (otherwise
the temperature profile would depart from a power law interior to $r_{\rm
cool}$). With $r_{\rm cool}=5r_{0}$, the transition occurs between $5r_{0}$
and $7r_{0}$. When scaling our supergiant model to the physical parameters of
a star, we will associate this transition region with the photosphere in the
stellar model (see Sec. 3.6). In our simulations, we also assume that
convective material that reaches $r>r_{\rm cool}$ can cool. See Section 3.3.3
for a description of our cooling function implementation.

Given $T_{\rm init}(r)$, the ideal gas law, and $\Phi(r)$, the equation of
hydrostatic equilibrium takes the form

$\frac{dP_{\rm init}}{dr}=-\frac{\mu m_{p}}{k_{b}}\frac{P_{\rm init}}{T_{\rm
init}}\frac{d\Phi}{dr}.$ (15)

We integrate eq. (15) numerically to obtain the initial pressure profile,
$P_{\rm init}(r)$. The initial density profile, $\rho_{\rm init}(r)$ is
obtained from the ideal gas law.

The functions $T_{\rm init}(r)$, $P_{\rm init}(r)$, and $\rho_{\rm init}(r)$
are plotted versus radius in panels (a)-(c) of Fig. 1 (black, solid curves).
The upper $x$-axes give $r$ in code units of $r_{0}$. The curves adopt $r_{\rm
cool}=5r_{0}$ and $b=2.1$, as we do in all of our simulations, and
$a=0.16r_{0}$ as in our fiducial model. See Sec. 3.6 for a complete
description of the figure.

Figure 1: Panels (a)-(c): Comparison of our Athena++ model (with Plummer
softening length $a=0.16$) to an 18$M_{\odot}$ (at ZAMS) supergiant model,
computed with MESA, near the end of oxygen burning. The pink, dashed curve in
each panel shows the MESA model from the base of the convective hydrogen
envelope at $R\approx 3.5R_{\odot}$ to the photosphere at $r_{\rm
ph}=840R_{\odot}$ (this is the outer $10.49M_{\odot}$ of the now
$16.5M_{\odot}$ star). The black, solid lines are radial profiles of our
Athena++ model showing temperature (a), pressure (b), and mass density (c) at
initialization ($t=0$), scaled to the MESA model as described in Sec. 3.6. The
black, dotted lines for $t=4000t_{0}$ show the same quantities after the
simulation achieves thermal equilibrium (see Sec. 4.1). The bottom axes show
radii in physical units and the upper $x$-axes give the original code units
for reference. The vertical dotted line and vertical solid line, respectively,
show $r_{\rm ph}$ and the $\tau=c/c_{s}$ radius ($R=818R_{\odot}$) in the MESA
model. To the left of these lines, the envelope is optically thick and roughly
adiabatic. Grey shading indicates where $r<r_{\rm heat,out}=1.5a$ or $r>r_{\rm
cool}$. The physically important part of our simulation domain, $r_{\rm
heat,out}<r<r_{\rm cool}$, falls between these two regions. The Athena++ model
shown thus covers a range of $33R_{\odot}$ to $700R_{\odot}$ in radius for
this MESA model. Our polytropic envelope with $b=2.1$ reproduces the power-law
hydrogen envelope of the MESA model well. Panel (d): The dashed line shows the
binding energy of the envelope exterior to $R$ (left-hand $y$-axis) and the
dotted line shows the mass enclosed at each $R$ (right-hand $y$-axis) for the
MESA model. The open and filled arrows on the mass axis mark the mass
coordinates outside of which the envelope would be unbound due to a weak shock
with the kinetic energy marked by a corresponding arrow on the binding energy
axis. These shock energies bound the range of values found for RSGs by Ivanov
& Fernández (2021). If the weak shock has energy of $3.6\times 10^{47}$ erg
(filled arrows), then roughly $3.6M_{\odot}$ could be ejected, leaving $\sim
6.9M_{\odot}$ of the hydrogen envelope to collapse onto the BH.

###  3.3 Source Terms

In this section, we define the gravity, heating, cooling, and damping source
terms appearing in the momentum and energy equations (eqs. 6 and 7,
respectively). All of these source terms are active during our convection
simulations. Only gravity is active during our collapse simulations.

####  3.3.1 Gravity

We neglect the self-gravity of the envelope gas, so $\Phi$ is simply eq. (9),
implemented in Cartesian coordinates.

####  3.3.2 Heating

We drive convection by including a heating source term in eq. (7) that
operates in the region $r_{\rm heat,in}\lesssim r\lesssim r_{\rm heat,out}$.
The constant, volume-integrated heating rate is

$L_{\rm heat}=f_{\rm heat}\times 4\pi r_{0}^{2}\rho_{0}c_{s,0}^{3}$ (16)

where $c_{s,0}$ is the adiabatic sound speed at $r_{0}$ and $f_{\rm heat}$ is
a constant. The heating source term is then

$Q_{\rm heat}=\frac{L_{\rm heat}}{4V}\bigg{[}1+\tanh\bigg{(}\frac{r-r_{\rm
heat,in}}{0.03r_{0}}\bigg{)}\bigg{]}\bigg{[}1-\tanh\bigg{(}\frac{r-r_{\rm
heat,out}}{0.03r_{0}}\bigg{)}\bigg{]}$ (17)

where $V$ is the volume of the heating region and the $\tanh$ factors smooth
the boundaries of the heating annulus. Our convection simulations adopt
$r_{\rm heat,in}=0.9a$ and $r_{\rm heat,out}=1.5a$. We choose $f_{\rm
heat}=6\times 10^{-4}$ to roughly reproduce the convective Mach numbers of RSG
envelopes (see Section 4.1.1).

####  3.3.3 Cooling

We mimic the stellar photosphere by including a cooling term in eq. (7) given
by

$Q_{\rm cool}=\frac{\rho k_{b}\big{[}T-T_{\rm init}\big{]}}{2\mu
m_{p}(\gamma-1)\tau_{\rm cool}}\bigg{[}1+\tanh\bigg{(}\frac{r-r_{\rm
cool}}{0.3r_{0}}\bigg{)}\bigg{]}$ (18)

if $T>T_{\rm init}$ and $Q_{\rm cool}=0$ if $T\leq T_{\rm init}$, where
$\tau_{\rm cool}=f_{\rm cool}\sqrt{{r_{\rm cool}^{3}}/{GM}}$ and the spatial
dependence contained in the $\frac{1}{2}[1+\tanh{x}]$ factor ensures that the
cooling function only operates at $r\gtrsim r_{\rm cool}$. As hot convective
material rises beyond $r_{\rm cool}$, its temperature is brought back to
$T_{\rm init}(r)$, as is the case for convective material that reaches the
photosphere in real stars. The timescale over which this cooling occurs,
$\tau_{\rm cool}$, is a fraction, $f_{\rm cool}$, of the dynamical time at
$r_{\rm cool}$. We set $f_{\rm cool}=0.01$ in all simulations.

We include cooling and place $r_{\rm cool}$ well inside the simulation domain
because the rising convective parcels carry finite angular momentum. Cooled
convective material sinks back into the envelope region rather than exiting
the domain and increasing the net angular momentum in the domain, as would
occur if we instead extended the convective envelope to the grid boundaries.

####  3.3.4 Damping

While cooling allows most of the convective material to return to $r\lesssim
r_{\rm cool}$, convection in this region still drives waves and some matter to
larger radii, causing mass and angular momentum to leave the domain. This
results in non-zero total angular momentum in the computation domain over the
very long simulation times required to achieve thermal equilibrium. For
numerical reasons, we are limited in how steep we can make the density
gradient at $r_{\rm cool}$ (parameterized by $f_{\rm iso}$). In real stars,
the density contrast near the photosphere is much larger, limiting mass and
momentum loss relative to our simple model. Instead, we damp motions at large
radii in order to limit the mass and angular momentum that leaves the domain.
We apply damping beyond a radius of $r_{\rm damp}$ with $r_{\rm damp}$
sufficiently larger than $r_{\rm cool}$ to avoid damping in the region where
nearly all of the cooling happens.

The damping source term in eqs. (6) and (7) is implemented as

$\mn@boldsymbol{a}_{\rm
damp}=-\frac{1}{2}\bigg{(}\frac{\mn@boldsymbol{v}\cdot\hat{r}}{\tau_{\rm
damp}}\bigg{)}\bigg{[}1+\tanh\bigg{(}\frac{r-r_{\rm
damp}}{0.3r_{0}}\bigg{)}\bigg{]}\hat{r}$ (19)

where $\tau_{\rm damp}=f_{\rm damp}\sqrt{r_{\rm damp}^{3}/GM}$ is the damping
timescale and the $\frac{1}{2}[1+\tanh x]$ factor limits damping to $r\gtrsim
r_{\rm damp}$.

In all simulations, we set $r_{\rm damp}=8.0r_{0}$ and $f_{\rm damp}=0.05$, so
that $\tau_{\rm cool}=0.16\tau_{\rm damp}$ and $\tau_{\rm damp}$ is similar
enough to the dynamical time at $r_{\rm damp}$ to avoid strong reflection of
waves back into the region of interest in our simulations.

In our tests, inclusion of the damping term did not effect the simulation
results for $r<r_{\rm cool}$. In particular, it did not modify the specific
angular momentum profiles in this region which are the key measurement for
this study. The damping simply prevents our star-in-a-box from expanding too
much over the simulation time and limits the loss of mass from the simulation
domain.

By including the region of significant cooling well within the domain
boundaries and by damping motions at large $r$, we keep the total specific
angular momentum in the box, which starts at zero, below $3\times 10^{-4}$ (in
our code units of $j_{0}$, see below). This is $1-2$ orders-of-magnitude
smaller than the mean specific angular momentum in any one shell (see Fig. 8),
so the change in total angular momentum in the box is negligible for our
purposes.

###  3.4 Boundary Conditions

####  3.4.1 Outer Boundaries

In all simulations, we use a box size $(40r_{0})^{3}$ centered at $x=y=z=0$.
At the outer boundaries, $x=y=z=\pm 20r_{0}$, we use zero-gradient, diode
boundary conditions that allow material to leave but not enter the domain.
That is, the density, pressure, and velocity components of the ghost zones are
copied from the last grid cell of the computational domain, but with the
following exception. The velocity component perpendicular to the boundary is
set to zero if the direction represents inflow.

####  3.4.2 Absorbing Sink Inner Boundary

In our simulations of collapse, we activate an absorbing sink centered at
$\mathbf{r}=\mathbf{0}$. The radius of the sink, $r_{s}$, is a runtime
parameter. When the sink is active, the density and pressure inside the sink
are set to $10^{-2}$ times the average values of the density and pressure just
outside the sink and the fluid velocities inside the sink are set to zero.
This allows material to flow into the sink unimpeded, representing perfect
accretion without feedback. The sink is not active during our convection
simulations, which do not have any inner boundary condition, since we are
using a Cartesian grid. Instead, the Plummer potential ensures that our
converging flows are well-behaved as $r\to 0$.

###  3.5 Initialization

Grid variables are initialized in Athena++ by specifying $\rho$, $E$ and the
Cartesian velocity components $\mn@boldsymbol{v}=(v_{x},v_{y},v_{z})$. The
initial density is $\rho_{\rm init}(\mathbf{r})$. The model is in hydrostatic
equilibrium, so $\mn@boldsymbol{v}(\mathbf{r},t=0)=\mathbf{0}$. Therefore, the
initial total energy density is $E(\mathbf{r},t=0)={P_{\rm
init}(\mathbf{r})}/{(\gamma-1)}$.

###  3.6 Code Units and Comparison to MESA Model

In our simulations, we set $GM=(r_{0}^{n}+a^{n})^{\frac{1}{n}}$ with
$r_{0}=1$. This yields a characteristic velocity of

$v_{0}=\bigg{(}\frac{GM}{r_{0}}\bigg{)}^{1/2}=\bigg{[}\frac{(r_{0}^{n}+a^{n})^{\frac{1}{n}}}{r_{0}}\bigg{]}^{1/2}\approx
1$ (20)

at $r_{0}$ and a time unit of

$t_{0}=\bigg{(}\frac{r_{0}^{3}}{GM}\bigg{)}^{1/2}=\Bigg{[}\frac{r_{0}^{3}}{(r_{0}^{n}+a^{n})^{\frac{1}{n}}}\Bigg{]}^{1/2}\approx
1$ (21)

where in each expression the last equality assumes $a^{n}\ll 1$. In these
units, the adiabatic sound speed at $r_{0}$ is

$c_{s,0}=\sqrt{\gamma/(b+1)}$ (22)

and the Keplerian specific angular momentum at $r_{0}$ is

$j_{0}=v_{0}r_{0}=\sqrt{GMr_{0}}=1.$ (23)

The unit of gas mass is $m_{0}=\rho_{0}r_{0}^{3}$, so mass accretion rates are
reported in units of $\dot{m}_{0}\equiv
m_{0}/t_{0}={\rho_{0}r_{0}^{3}}/{t_{0}}=\rho_{0}r_{0}^{3/2}(GM)^{1/2}$.

Our dimensionless setup can be scaled to a wide range of RSG and YSG envelopes
whose convective Mach number profiles are similar to those in Fig. 6 by
identifying the photosphere radius as the radius of our isothermal transition
at $r\approx 6r_{0}$. Rescaling is done most simply by computing the
quantities $r_{\rm ph}$, $M_{\bullet}$, and $\rho_{\star}$ from the stellar
model. These are, respectively, the photosphere radius, the mass interior to
the hydrogen envelope, and the density at $r_{\rm ph}/6$. Code units are then
scaled to the star by setting $M=M_{\bullet}$, $r_{0}=r_{\rm ph}/6$, and
$\rho_{0}=\rho_{\star}$ giving,

$r_{0}=1.16\times 10^{13}\bigg{(}\frac{r_{\rm ph}}{1000R_{\odot}}\bigg{)}{\rm
cm}$ (24)

$t_{0}=0.034\bigg{(}\frac{r_{\rm
ph}}{1000R_{\odot}}\bigg{)}^{{3}/{2}}\bigg{(}\frac{10M_{\odot}}{M_{\bullet}}\bigg{)}^{{1}/{2}}\text{yr}$
(25)

$\displaystyle j_{0}=1.24\times 10^{20}\bigg{(}\frac{r_{\rm
ph}}{1000R_{\odot}}\bigg{)}^{1/2}\bigg{(}\frac{M_{\bullet}}{10M_{\odot}}\bigg{)}^{1/2}{\rm
cm}^{2}\,{\rm s}^{-1}$ (26)

$\dot{m}_{0}=2.27\bigg{(}\frac{\rho_{\star}}{10^{-7}\text{g}/\text{cm}^{3}}\bigg{)}\bigg{(}\frac{r_{\rm
ph}}{1000R_{\odot}}\bigg{)}^{{3}/{2}}\bigg{(}\frac{M_{\bullet}}{10M_{\odot}}\bigg{)}^{{1}/{2}}M_{\odot}\,{\rm
yr}^{-1}.$ (27)

Panels (a)-(c) of Fig. 1 show the initial profiles for our Athena++ setup
(black, solid curves) scaled to an $18M_{\odot}$ (at ZAMS) star computed with
MESA (pink, dashed curves). The MESA model, which adopts Solar metallicity and
a mixing-length $\alpha$ of 1.5, is shown at the end of oxygen burning, at
which time the star is a $16.5M_{\odot}$ RSG. The Athena++ curves have been
scaled from dimensionless code units using the values $r_{\rm
ph}=840R_{\odot}$, $M_{\bullet}=6M_{\odot}$, and $\rho_{\star}=4.0\times
10^{-7}$ g cm-3 obtained from the MESA model. The upper $x$-axis shows $r$ in
the original code units, for reference. The vertical dotted line in each panel
marks $r_{\rm ph}$. The $\tau=c/c_{s}$ radius (vertical solid line) is
$818R_{\odot}$. To the left of this line, the MESA model is optically thick
and fluid motions in the envelope are roughly adiabatic as we assume in our
simulated envelope.

Panel (d) of Fig. 1 shows, for reference, the binding energy of the envelope
exterior to $R$ (dashed line, left-hand $y$-axis) and the mass enclosed at
each $R$ (dotted line, right-hand $y$-axis) for the MESA model. The open and
filled arrows on the left, respectively, mark $0.55\times 10^{47}$ and
$3.6\times 10^{47}$ erg, which are roughly the lowest and highest of the shock
kinetic energies obtained for RSG and YSG progenitors by Ivanov & Fernández
(2021) in their studies of mass-ejection in FSN (this range of energies covers
three different NS equations-of-state, one quite stiff and one quite soft).
The mass that is lost is roughly the material exterior to where the shock
kinetic energy is equal to the binding energy. The arrows on the right show
the mass coordinates that correspond to the two shock energies. If the weak
shock has energy of $3.6\times 10^{47}$ erg (filled arrow at left), then
material exterior to the mass coordinate of $12.9M_{\odot}$ (filled arrow at
right) would be ejected and roughly $\sim 6.9M_{\odot}$ of the hydrogen
envelope would remain bound and collapse onto the BH. For the lower shock
energy in Fig. 1, even more of the hydrogen envelope remains bound. Most of
the hydrogen envelope that we simulate is thus likely to remain bound and
accrete onto the BH in a FSN.

###  3.7 Summary of Runtime Parameters and Simulations Performed

In all simulations, we adopt a Plummer index of $n=8$ and a density power-law
index of $b=2.1$. Cooling occurs beyond a radius of $r_{\rm cool}=5r_{0}$ and
damping occurs outside of $r_{\rm damp}=8r_{0}$. The cooling timescale
fraction is $f_{\rm cool}=0.01$ and the damping timescale fraction is $f_{\rm
damp}=0.05$. We adopt $f_{\rm heat}=0.0006$ for the heating-rate parameter.
The adiabatic index in all simulations is $\gamma=1.4762$, which gives a flat
entropy profile in the envelope at initialization for our choice of $b$.

Table 1: Table of convection simulations.

Convect |  $\,\,a$ |  Basea |  Refineb |  $\Delta x_{\rm min}$c |  $t_{\rm max}$d |
---|---|---|---|---|---|---
Model |  $[r_{0}]$ |  Cells |  Levels | $[r_{0}]$ | $[t_{0}]$ | $a/\Delta x_{\rm min}$
F |  0.16 |  $128^{3}$ |  5 | 0.0098 | 5000 | 16.3
A |  0.08 |  $128^{3}$ |  6 | 0.0049 | 5000 | 16.3
R |  0.16 |  $64^{3}$ |  5 | 0.0195 | 4000 | 8.2

a Number of grid cells in the base, unrefined grid.
b Number of SMR refinement levels on top of the base resolution.
c Width of each grid cell in highest refinement region.
d Simulation run time.

Table 2: Table of collapse simulations.

Collapse | Restarta | Restart Timeb | Sink Sizec |
---|---|---|---|---
Run | Model | $t_{s}\,[t_{0}]$ | $r_{s}\,[r_{0}]$ | $r_{s}/\Delta x_{\rm min}$
1 | F | 4560 | 0.08 | 8.2
2 | F | 4900 | 0.08 | 8.2
3 | F | 4270 | 0.08 | 8.2
4 | F | 4800 | 0.08 | 8.2
1s | F | 4560 | 0.04 | 4.1

a Convection model used to initialize the collapse run.
b Convection model snapshot used to initialize the collapse run.
c Radius of low-pressure sink activated at start of collapse run.

Our convection simulations, summarized in Table 2, vary the softening length,
$a$, and the grid resolution. Model F is our fiducial model. Model A reduces
$a$ by a factor of 2 and adds one level of refinement in order to preserve
$a/\Delta x_{\rm min}$. Model R reduces the grid resolution by a factor of 2.
We note that the scale height $H\sim r$ in the simulated envelope, so there
are a priori no other length scales that we need to resolve. The numerical
challenge is that the simulations have to be run for $\sim 5000t_{0}$ and
$\sim 3\times 10^{6}$ timesteps to reach thermal equilibrium.

The grid configuration for Model F is as follows. The domain covers a spatial
volume of $(40\,r_{0})^{3}$, centered on $x=y=z=0$, and has a base resolution
of $128^{3}$ cells, which translates to a base cell size of $\Delta x=\Delta
y=\Delta z=0.3125$. We use static mesh refinement (SMR) to increase resolution
with decreasing $r$. For Model F, there are 5 levels of refinement above the
base resolution and the refinement transitions occur where $x$, $y$, and $z$
have the values $[\pm 0.625,\pm 1.25,\pm 3.75,\pm 7.5,\pm 10.0]$. At each of
these transitions, $\Delta x$, $\Delta y$, and $\Delta z$ each decrease by 2.
The highest refinement region, $-0.625<x,y,z<0.625$, has a cell size of
$\Delta x_{\rm min}=\Delta y_{\rm min}=\Delta z_{\rm
min}=0.3125/2^{5}=0.0098$, so that $a/\Delta x_{\rm min}=16.3$ for the
fiducial model.

Model R has the same domain configuration, but with a reduced base resolution
of $64^{3}$ cells instead of $128^{3}$. This translates to a cell size of
$\Delta x_{\rm min}=\Delta y_{\rm min}=\Delta z_{\rm min}=0.625/2^{5}=0.0195$
in the highest refinement region ($-0.625<x,y,z<0.625$). The Model A grid
structure is the same as Model F, except that we include a sixth level of
refinement where $-0.3125<x,y,z<0.3125$ so that $\Delta x_{\rm min}=\Delta
y_{\rm min}=\Delta z_{\rm min}=0.0049$.

Convection Model F provides the initial conditions for our collapse
simulations, which are summarized in Table 2. The collapse simulations vary
the sink size, $r_{s}$, as well as the start time, $t_{s}$, for the collapse.
The ‘Restart Time’ column indicates which snapshot from Model F was used to
start the collapse simulation. We activate the sink immediately, so $t_{s}$ is
also the time that collapse begins. The grid resolution of each collapse
simulation is inherited from Model F.

##  4 Convection Simulations

In this section, we present the results of our suite of convection
simulations, summarized in Table 2. We discuss our fiducial simulation in
detail in Section 4.1. We study the effects of changing the softening length,
$a$, and the grid resolution in Section 4.2. We compare our results to
analytical predictions in Section 4.3.

###  4.1 Fiducial Model

Our fiducial model (convection Model F) has $a=0.16r_{0}$ and the heating
region extends to $r_{\rm heat,out}=1.5a=0.24r_{0}$. The cooling radius is at
$r_{\rm cool}=5r_{0}$. From $r_{\rm heat,out}$ to $r_{\rm cool}$, our model
covers a factor of $\sim 20$ in radius in the convective zone. Details about
the computational domain were given in Section 3.7. We run the model to a
maximum time of $5000\,t_{0}$ which is $\sim$ 400 dynamical times at $r_{\rm
cool}$ and $\sim 3\times 10^{6}$ timesteps.

In Fig. 2, we plot the instantaneous, volume-integrated heating and cooling
rates, $L_{\rm heat}$ and $L_{\rm cool}\equiv\int Q_{\rm cool}dV$,
respectively. $L_{\rm heat}$ is a constant and integrates to eq. (16). The
cooling rate is time-dependent as it depends on the local temperature
difference in the cooling region (see eq. 18). At early times ($t\lesssim
2000t_{0}$), there is transient convective flow due to the temperature and
density drop-off at $\sim r_{\rm cool}$, resulting in a spike in $L_{\rm
cool}$. After this transient flow subsides, convection is driven by our
heating source and the flow is near thermal equilibrium with $L_{\rm
cool}\approx L_{\rm heat}$.

Figure 2: Volume-integrated heating rate (dashed line) and cooling rate (solid
line) for Model F. For $t\gtrsim 2000t_{0}$, the flow approaches thermal
equilibrium with $L_{\rm cool}\approx L_{\rm heat}$. The flow is qualitatively
similar for all $t\gtrsim 2000t_{0}$. For quantitative analysis, however, we
will restrict our attention to $4000<t/t_{0}<5000$, when the time average of
$L_{\rm cool}$ is very nearly equal to $L_{\rm heat}$ (see Fig. 3).

We consider the flow to be in thermal equilibrium if (1) the time- and
spherically-averaged convective luminosity, $\langle\overline{L}_{\rm
conv}\rangle$, is roughly independent of radius and (2) the time-averaged
cooling rate, $\langle L_{\rm cool}\rangle$, is equal to $L_{\rm
heat}$.222Here and in the remainder of this paper, we use $\overline{q}(r,t)$
or, simply, $\overline{q}$ to denote the radial profile of the 3D field
quantity $q(x,y,z,t)$ and use $\langle\,\rangle$ to denote time-averaged
quantities. All radial profiles are computed by finding the volume-weighted
mean of the quantity of interest over spherical shells.

We compute $\langle\overline{L}_{\rm conv}\rangle$ as follows. The radial
total energy flux is

$\mathcal{F}_{\rm tot}(x,y,z,t)=\bigg{(}\frac{1}{2}\rho
v^{2}+\frac{\gamma}{\gamma-1}P+\rho\phi\bigg{)}v_{r},$ (28)

with spatial mean $\overline{\mathcal{F}}_{\rm tot}(r,t)$. Following Parrish
et al. (2008), we decompose the density, temperature, and velocity fields into
a mean radial profile and a local deviation from the mean

$\displaystyle\rho$ $\displaystyle=\overline{\rho}(r,t)+\delta\rho$ (29)
$\displaystyle T$ $\displaystyle=\overline{T}(r,t)+\delta T$ (30)
$\displaystyle\mn@boldsymbol{v}$
$\displaystyle=\overline{\mn<EMAIL_ADDRESS>(31)

Inserting the decomposed fields into $\overline{\mathcal{F}}_{\rm tot}(r,t)$
yields advective terms proportional to $\overline{\rho}\cdot\overline{v_{r}}$
and $\overline{\delta\rho\delta v_{r}}$ (e.g., eq. (41) of Parrish et al.).
Subtracting these advective terms from $\overline{\mathcal{F}}_{\rm tot}(r,t)$
gives the radial convective flux profile $\overline{\mathcal{F}}_{\rm
conv}(r,t)$. The instantaneous convective luminosity profile is
$\overline{L}_{\rm conv}(r,t)=4\pi r^{2}\overline{\mathcal{F}}_{\rm
conv}(r,t)$ whose time average is denoted $\langle\overline{L}_{\rm
conv}\rangle$. In general, we find that the advective terms are small so
$\langle\overline{L}_{\rm tot}\rangle\approx\langle\overline{L}_{\rm
conv}\rangle$.

The purple, solid lines in Fig. 3 show $\langle\overline{L}_{\rm conv}\rangle$
for the time ranges listed in the lower legend. The constant heating rate,
$L_{\rm heat}$, is indicated with a dashed, black line. The cooling rate,
averaged over $4000\leq t/t_{0}\leq 5000$, is shown with the dotted, blue
line. After $t\sim 2000t_{0}$, $\langle\overline{L}_{\rm conv}\rangle$ is
roughly independent of radius and we will see in the next sections that the
flow is dynamically similar over these times. Our conditions for thermal
equllibrium are met after $t\sim 4000t_{0}$, when $\langle\overline{L}_{\rm
conv}\rangle\simeq\langle L_{\rm cool}\rangle\simeq L_{\rm heat}$. The black,
dotted lines in Fig. 1 show radial profiles of the simulation upon reaching
thermal equilibrium.

Figure 3: Time- and spherically-averaged convective luminosity (purple lines;
times used for each average are given in the legend), heating rate (black,
dashed line), and time-averaged cooling rate (blue, dotted line; time average
over $4000-5000$ $t_{0}$) for Model F. Grey shading is the same as in Fig. 1.
The convective luminosity is roughly independent of radius and the flow is
dynamically similar for $t>2000t_{0}$. The flow is in thermal equlibrium with
$\langle\overline{L}_{\rm conv}\rangle\simeq\langle L_{\rm cool}\rangle\simeq
L_{\rm heat}$ for $4000\lesssim t/t_{0}\lesssim 5000$. Figure 4: Velocity
field slices through the $z=0$ plane of our 3D simulation of convection at
three simulation times for Model F. The $z$-axis is arbitrary as our setup has
no preferred direction. Slices through other planes that include the origin
look similar. Top, middle, and bottom rows show the radial, polar, and
azimuthal components, respectively, of the fluid velocity field. With
$\theta=0$ along the $+z$ axis, $v_{\theta}>0$ indicates material flowing into
the page. The cooling radius, $r_{\rm cool}=5r_{0}$, begins the transition to
an exponential fall-off in gas density. In $v_{r}$, large plumes form and
subside throughout the simulation, bringing high entropy material out to
$r_{\rm cool}$. Cooled gas falls back towards the center in streams of varying
sizes that collide with material rising from the center. Flows perpendicular
to $r$ ($v_{\theta}$ and $v_{\phi}$) have a more obvious radial dependence
with smaller-scale structures at small $r$, where the scale height is
smallest.

Figure 4 shows snapshots of the total velocity field of the flow at $4000$,
$4500$, and $5000$ $t_{0}$ for fiducial Model F. Each panel is a slice through
the $z=0$ plane. The $z=0$ plane is an arbitrary choice; our setup has no
preferred axis and slices through other planes that include the origin are
similar. The top row shows radial velocity, $v_{r}(x,y,z=0)$, the second row
shows velocities in the $\theta$ direction, $v_{\theta}(x,y,z=0)$, and the
bottom row shows velocities in the $\phi$ direction, $v_{\phi}(x,y,z=0)$. Our
angles $\theta$ and $\phi$ are defined in the usual way, relative to the
(arbitrary) $z$-axis. Although our simulation domain extends to $\pm 20r_{0}$
in each Cartesian direction, we restrict the slices to $-6.5<x/r_{0}<6.5$ and
$-6.5<y/r_{0}<6.5$. The density falls off exponentially beyond $r\approx
6r_{0}$, so there is little momentum outside of this radius. The physically
relevant part of the domain is $r_{\rm heat,out}=0.24r_{0}\lesssim r\lesssim
r_{\rm cool}=5r_{0}$.

The panels of Figure 4 show turbulent flow at all scales. The radial flows
(top row) exhibit plumes that individually occupy a large fraction of the star
in both solid angle and radius. This morphology is consistent with other 3D
simulations of RSG envelopes (Chiavassa et al., 2009; Goldberg et al., 2021)
and with the large convective cells that have been inferred from observations
of Betelgeuse (Chiavassa et al., 2010; Dupree et al., 2020). The large, rising
plumes are often interrupted by narrower streams of cool material falling back
towards the origin from $r_{\rm cool}$. Though not apparent in this figure,
the infalling streams colliding at small $r$ often launch waves that move out
through the rising plumes. Features of the flow are qualitatively consistent
with analytical expectations. Given the density scale height in the envelope
of $H=r/b\approx 0.48r$, mixing-length-theory predicts large eddies at large
$r$. On average, the smallest-scale structures indeed exist at the smallest
$r$.

For a more quantitative analysis of our fiducial model, we compute radial
profiles over time to characterize the convective Mach number and angular
momentum content of the flow in the following subsections.

####  4.1.1 Mean Turbulent Mach Number

Figure 5: Spherically-averaged turbulent Mach number, eq. (32), as a function
of radius and time for Model F. The flow is dynamically similar for $t\gtrsim
2000t_{0}$ with Mach numbers of $\approx 0.05-0.25$. Comparing to Fig. 8,
times with larger-than-average $\mathcal{M}_{\rm turb}$ translate to higher-
than-average specific angular momentum.

We first consider profiles of the turbulent Mach number, $\mathcal{M}_{\rm
turb}$, which is the root-mean-square turbulent velocity profile divided by
the mean sound speed profile

$\mathcal{M}_{\rm turb}(r,t)=\frac{v_{\rm turb}(r,t)}{c_{s}(r,t)}.$ (32)

To find $v_{\rm turb}(r,t)$, we first compute the random velocity field

$\delta\mn@boldsymbol{v}(x,y,z,t)=\mn@boldsymbol{v}(x,y,z,t)-\overline{v_{r}}(r,t)\hat{r}$
(33)

where $\mn@boldsymbol{v}(x,y,z,t)$ is the total velocity field and
$\overline{v_{r}}(r,t)\hat{r}$ is the mean radial velocity profile. The
turbulent velocity profile, $v_{\rm
turb}(r,t)\equiv\overline{\|\delta\mn@boldsymbol{v}\|}(r,t)$, is the average
of the magnitude of the random velocity in each shell. We note that
subtracting $\overline{v_{r}}(r,t)\hat{r}$ is only necessary during the
initial transient phase (see. Fig. 2). After $\sim 1500t_{0}$,
$\overline{v_{r}}(r,t)$ is negligible.

Fig. 5 shows $\mathcal{M}_{\rm turb}(r,t)$ for our fiducial model for
$2000\leq t/t_{0}\leq 5000$, when $\langle\overline{L}_{\rm conv}\rangle$ is
independent of radius and when the kinetic energy in the box is constant to
within $\approx 10-15\%$. Each column in the figure plots the mean radial
profile $\mathcal{M}_{\rm turb}(r)$ at fixed time. We limit the $y$-axis to
$r\lesssim 5r_{0}$, which is the region of interest in our domain. The
$y$-axis is log-scaled to show features at small radii. The flow achieves
turbulent Mach numbers of $\approx 0.05-0.25$. The profiles are qualitatively
similar over the times shown, though there are periods of higher- and lower-
than-average Mach number.

Figure 6: Turbulent Mach number (top panel) and Brunt-Väisälä frequency
(bottom panel) profiles for our fiducial Athena++ model at the simulation
times listed in the legend (solid and dotted curves). The pink, dashed curves
show our MESA RSG model. The $x$-axes and grey shading are the same as in Fig.
1 and $N^{2}$ for the Athena++ model has been scaled to physical units as in
that figure. With our choice of $f_{\rm heat}=6\times 10^{-4}$, we capture the
convective Mach number profile of supergiant envelopes well.

The top panel of Fig. 6 compares $\mathcal{M}_{\rm turb}(r)$ for the
simulation times listed in the legend (solid, grey curves) to the convective
Mach number333The convective Mach number in MESA derives from the radial
convective flux from mixing length theory while our definition of $v_{\rm
turb}$ uses all three components of $\delta\mn@boldsymbol{v}$. profile of our
MESA RSG model (pink, dashed curve). The bottom panel of the figure compares
the Brunt-Väisälä frequency instead. The $r$ values of our simulation have
been scaled as in Fig. 1. Our simulation achieves envelope turbulent Mach
numbers that are similar to this MESA model. The $\mathcal{M}_{\rm turb}$
profile for the MESA model shown is representative of RSGs and YSGs and,
therefore, our simulation is broadly applicable to a wide range of supergiant
progenitors.

####  4.1.2 Total Angular Momentum

An important objective of our convection simulations is to quantify the (non-
zero) mean specific angular momentum at each radius that arises due to the
turbulent convective flows. Before turning to that question in Section 4.1.3,
it is important to first emphasize that the star is not rotating, so the net
angular momentum in the domain must be close to zero. We explore this in Fig.
7, where we plot the cumulative Cartesian components of the total angular
momentum vector, which are computed as follows.

Figure 7: Components of $\mn@boldsymbol{J}(r,t)$, the total angular momentum
enclosed at $r$, divided by the total gas mass in the domain. Although the
total angular momentum is very small, there is a non-zero angular momentum
enclosed at each $r$ $\lesssim 8r_{0}$ due to random convective flows.
Accretion of material during collapse leads to the accretion of finite angular
momentum as a function of time even though the star has roughly zero net
angular momentum.

Let $\ell_{i}=\rho(\mn@boldsymbol{r}\times\mn@boldsymbol{v})_{i}$ be the $i$th
component of the angular momentum density in each grid cell at time $t$. We
compute the mean radial profile of $\ell_{i}$ to give
$\overline{\ell}_{i}(r,t)$. Then we compute the $i$-th component of the total
angular momentum vector enclosed at a radius $r$ at fixed time,

$J_{i}(r,t)=\int_{0}^{r}4\pi r^{2}\overline{\ell}_{i}(r,t)dr.$ (34)

We also compute the total gas mass in the box at each time, $m(t)$. The top
three panels of Fig. 7 plot $J_{i}(r,t)/m(t)$ for $i=x,y,$ and $z$. The bottom
panel shows the magnitude $J/m=\sqrt{\sum_{i}J_{i}^{2}}/m$. The curves in each
panel correspond to the different simulation times listed in the legend. We
have divided by $m(t)$ in order to put the $y$ axis in units of specific
angular momentum ($j_{0}$), as in later figures.

Fig. 7 shows that although each component of the cumulative Cartesian total
angular momentum vector, $\mn@boldsymbol{J}(r,t)$, goes to zero444There is a
small non-zero total specific angular momentum in the box due to gas leaving
the domain. See Sec. 3.3.4. when integrated out to $r\approx 8.5r_{0}$, there
is a finite $\mn@boldsymbol{J}(r,t)$ at each $r\lesssim 8.5r_{0}$. Consider,
for example, the $t=5000t_{0}$ curve for $J_{z}(r,t)$ in the third panel. The
curve becomes large and positive out to $r\approx 3r_{0}$ as all of the inner
shells with positive $J_{z}$ are added to the sum. The curve then drops toward
zero as the material with negative $J_{z}$ is added to the integral.
Considering the $t=5000t_{0}$ curve in each panel, if only the material out to
$r\approx 3r_{0}$ were accreted, a net
$\boldsymbol{J}/m\approx(0.001\hat{x}+0.0005\hat{y}+0.004\hat{z})j_{0}$ with
magnitude $J/m=0.0045j_{0}$ would be available to feed a rotationally-
supported structure at small radii. Although infall of the entire envelope
implies a total specific angular momentum budget of
$\approx\mn@boldsymbol{0}$, this (perhaps surprisingly) does not mean that the
spin of the BH would be zero if the entire star collapsed. We return to this
in Section 6.2.1.

For this particular simulation, the magnitudes of each component of
$\mn@boldsymbol{J}$ are within a factor of a few of one another for the times
shown. This is just a coincidence. The turbulent velocity field is random, as
is the direction of $\mn@boldsymbol{J}$. At other times and in other
simulations, one or two of the Cartesian components is larger than the
other(s) by an order-of-magnitude. There is no preferred direction in our
setup.

####  4.1.3 Mean Specific Angular Momentum

We now turn to the central question of characterizing the angular momentum
profile of the convective envelope. We define the specific angular momentum
vector of the gas in each grid cell as $\mn@boldsymbol{j}_{\rm rand}$. The
angle between $\mn@boldsymbol{j}_{\rm rand}$ and the +$z$-axis is denoted
$\theta_{j}$ and the angle between $\mn@boldsymbol{j}_{\rm rand}$ and the
$+x$-axis (in the $x$-$y$ plane) is $\phi_{j}$. Finally, we define
$\hat{{\jmath}}_{\rm rand}$ to be the unit vector in the direction of
$\mn@boldsymbol{j}_{\rm rand}$.

Figure 8: Spherically-averaged specific angular momentum as a function of
radius and time for Model F. There is finite mean angular momentum at each
radius. When scaled to our MESA supergiant model, $\overline{\jmath}_{\rm
rand}(r,t)\approx(4-28)\times 10^{17}$ cm2 s-1 (see eq. (26) and accompanying
text) implying circularization radii of many hundreds times the ISCO radius
for a newly-formed BH (see Fig. 17). Figure 9: Polar (top panel) and azimuthal
(bottom panel) direction angle profiles of the spherically-averaged specific
angular momentum vector, $\mn@boldsymbol{\overline{\jmath}}_{\rm rand}$, over
time for Model F. In the top panel, white regions show where
$\mn@boldsymbol{\overline{\jmath}}_{\rm rand}$ lies near the $x$-$y$ plane
while red (blue) colors correspond to a positive (negative) $z$-component of
$\mn@boldsymbol{\overline{\jmath}}_{\rm rand}$. Two regions of similar
brightness but opposite color have anti-parallel $z$-components of
$\mn@boldsymbol{\overline{\jmath}}_{\rm rand}$. In the lower panel, colors
that differ by $\pi$ correspond to anti-parallel projections of
$\mn@boldsymbol{\overline{\jmath}}_{\rm rand}$ in the $x$-$y$ plane. Material
from different regions in $r$ and fixed $t$ show features with opposite
direction angles, as is required for $\mathbf{J}\approx\mathbf{0}$ when
integrated over the box. The direction of
$\mn@boldsymbol{\overline{\jmath}}_{\rm rand}$ at each radius is not fixed in
time but flips in sign over long timescales. The large-scale structures in
each panel show connected regions in space and time where the direction of
$\mn@boldsymbol{\overline{\jmath}}_{\rm rand}$ is changing slowly. This
coherence of $\mn@boldsymbol{\overline{\jmath}}_{\rm rand}$ over radius and
time suggests coherent accretion of $\mn@boldsymbol{\overline{\jmath}}_{\rm
rand}$ as material falls in towards the BH.

Fig. 8 shows $\overline{\jmath}_{\rm
rand}(r,t)\equiv\|\mn@boldsymbol{\overline{\jmath}}_{\rm rand}(r,t)\|$ for our
fiducial model. Although the net angular momentum in the envelope is nearly
zero, convection gives rise to fluid motions that are perpendicular to $r$,
resulting in a finite mean angular momentum at each radius. The distribution
of momentum as a function of $r$ varies in time. Each time (that is, each
column in the figure) is a specific realization of the stochastic angular
momentum distribution. The time that the star collapses is arbitrary with
respect to the flow, and the accretion rate of angular momentum over time will
depend on the state of the convective flow over the time that it takes for the
envelope to collapse. We will explore this idea in greater detail in the
context of our collapse simulations in Section 5.

The spherically-averaged specific angular momentum profiles of Fig. 8 vary
within $0.005\lesssim\overline{\jmath}_{\rm rand}/j_{0}\lesssim 0.032$, in our
code units with $j_{0}=\sqrt{GMr_{0}}$ (for comparison, the specific angular
momentum for rotation at break-up is $0.5-2.0$ $j_{0}$ between $0.24r_{0}$ and
$5r_{0}$; so the random flows have specific angular momenta that are $0.3-6$ %
of breakup). We will apply this result to a specific supergiant model in Fig.
17, but for now we make a rough comparison to eq. (3). Using eq. (26) with
$M_{\bullet}=10M_{\odot}$ and $r_{\rm ph}=1000R_{\odot}$, our simulated
envelope achieves specific angular momenta of $\overline{\jmath}_{\rm
rand}\approx(6-40)\times 10^{17}$ cm2 s-1 (note that $\overline{\jmath}_{\rm
rand}>j_{\rm ISCO}=1.5\times 10^{17}$ cm2 s-1 for a non-spinning $10M_{\odot}$
BH). The values of $j_{\rm rand}$ that we find are a factor of a few $-$ 10
times the estimate of eq. (3) that was derived by Quataert et al. (2019); we
return to this later in Fig. 11.

For comparison, the $15M_{\odot}$ model of Heger et al. (2005), which has a
birth rotational velocity of $200$ km s-1, has rotational specific angular
momentum of $5\times 10^{18}$ cm2 s-1 at its surface at birth. By the pre-SN
phase, the star has rotational specific angular momentum of $\sim 8\times
10^{17}$ cm2 s-1 at its surface, declining to $\sim 5\times 10^{16}$ cm2 s-1
near the base of the hydrogen envelope. The net rotational angular momentum is
small compared to the random angular momentum of the convection zone except at
the surface of the star where the two are comparable.

The direction of $\mn@boldsymbol{\overline{\jmath}}_{\rm rand}$ over time is
important for determining the dynamics at smaller radii when the envelope
eventually accretes onto the newly-formed BH. The upper and lower panels of
Fig. 9 plot $\overline{\theta}_{j}(r,t)$ and $\overline{\phi}_{j}(r,t)$,
respectively, in units of $\pi$. In the top panel, white regions show where
$\mn@boldsymbol{\overline{\jmath}}_{\rm rand}$ lies near the $x$-$y$ plane.
The large regions of dark red show where $\overline{\theta}_{j}\lesssim\pi/4$.
Similarly, the darkest blue regions of space and time are where the $-z$
component of $\mn@boldsymbol{\overline{\jmath}}_{\rm rand}$ is largest. In the
lower panel, the darkest colors correspond to $\overline{\phi}_{j}\approx\pi$
and lightest colors correspond to $\overline{\phi}_{j}\approx 2\pi$. To
interpret the implications of Fig. 9, we note that if $\hat{{\jmath}}_{\rm
rand}$ is rapidly varying on timescales of the dynamical time at the
circularization radius, then the gas could not all circularize in the same
plane and one might expect a more spherical accretion flow at scales of the
circularization radius, $r_{\rm circ}$. Additionally, rapidly varying
$\hat{{\jmath}}_{\rm rand}$ in the radial profile could mean greater
cancelation of $\mn@boldsymbol{j}_{\rm rand}$ as shells (and individual
parcels) each reach their respective $r_{\rm circ}$ and interact with one
another. Contrast this with the case of a rotating star in which
$\hat{{\jmath}}_{\rm rand}$ is constant and material can build up into a disk
as all of the rotationally-supported flow orbits the same axis.

For material with $j_{\rm rand}\sim 10^{-2}j_{0}$ (Fig. 8), the
circularization radius is $r_{\rm circ}\sim 10^{-4}r_{0}$ where the dynamical
time is $\sim 10^{-6}t_{0}$. Periodograms of $\overline{\jmath}_{\rm
rand}(r,t)$, $\overline{\theta}_{j}(r,t)$, and $\overline{\phi}_{j}(r,t)$ show
that the most power is on timescales larger than $80t_{0}$ $\gg 10^{-6}t_{0}$.
This implies that the direction of angular momentum is quite coherent on
timescales relevant to the circularization of the disk at small radii. In
Section 5.2, we will explore how much these mean profiles are modified as the
envelope collapses.

###  4.2 Dependence on $a$ and Resolution

To study the influence of the softening length, $a$, on the angular momentum
associated with convective flows, we run a simulation with $a=0.08r_{0}$
(Model A) for comparison to our fiducial run with $a=0.16r_{0}$ (Model F). The
grid is the same for the two models except that we add a sixth SMR level to
Model A (see Section 3.7). The two models are otherwise identical. To study
the influence of resolution, we run Model R, which is identical to Model F
except that it adopts a base resolution of $64^{3}$ cells instead of $128^{3}$
cells. Both simulations have 5 SMR levels above the base resolution with the
refinement transitions placed at the same $x$,$y$, and $z$ values. To compare
the different simulations, we take time averages of the instantaneous profiles
$\mathcal{M}_{\rm turb}(r,t)$ and $\overline{\jmath}_{\rm rand}(r,t)$ over the
simulation times that our conditions for thermal equilibrium are satisfied.
For Models F and A, these are $4000<t/t_{0}<5000$. For Model R, these are
$2500<t/t_{0}<3500$. The top panel of Fig. 10 shows the resultant time-
averaged convective Mach number profiles, $\langle\mathcal{M}_{\rm
turb}\rangle$, and the bottom panel shows the time-averaged specific angular
momentum profiles, $\langle\overline{\jmath}_{\rm rand}\rangle$. We restrict
each curve to $r_{\rm heat,out}<r<r_{\rm cool}$.

Comparing the solid curves for Model F to the dashed curves for Model A, both
simulations achieve nearly identical $\langle\mathcal{M}_{\rm turb}\rangle$
with a maximum difference of $5\%$ at $r\sim 2r_{0}$. This corresponds to
roughly a $20\%$ difference in $\langle\overline{\jmath}_{\rm rand}\rangle$ at
that radius. At small $r$, $\langle\overline{\jmath}_{\rm rand}\rangle$ is
identical between the two simulations.

For Model F and Model R (dotted curves), $\langle\mathcal{M}_{\rm
turb}\rangle$ differs from $5\%$ to $45\%$ as $r$ decreases from $\sim
2.5r_{0}$. Over this same region, $\langle\overline{\jmath}_{\rm rand}\rangle$
is lower in the lower-resolution simulation by $\approx 10-50\%$.

Our simulation with smaller $a$ (that is, simulating smaller radii closer to
the base of the convective envelope) shows very good convergence in
$\langle\overline{\jmath}_{\rm rand}\rangle$ for $r\lesssim 0.6r_{0}$ where
our grid is most well-resolved. Although we are not converged in
$\langle\overline{\jmath}_{\rm rand}\rangle$ at all radii, the more physical
simulation with higher resolution increases $\langle\overline{\jmath}_{\rm
rand}\rangle$, reenforcing our result that convective flows with
$\mathcal{M}_{\rm turb}$ similar to supergiant envelopes give rise to $j_{\rm
rand}>j_{\rm ISCO}$.

Figure 10: Time- and spherically-averaged convective Mach number profile (top)
and specific angular momentum profile (bottom) for our fiducial Model F (solid
lines), Model A with smaller Plummer smoothing length $a$ (dashed lines), and
Model R with $2\times$ lower base resolution (dotted lines). The time averages
are taken over $4000<t/t_{0}<5000$ (Models F and A) or $2500<t/t_{0}<3500$
(Model R), when the simulations are in thermal equilibrium. Comparing Model F
($a=0.16r_{0}$) to Model A ($a=0.08r_{0}$): $\langle\mathcal{M}_{\rm
turb}\rangle$ is very similar and where $\langle\mathcal{M}_{\rm turb}\rangle$
differs by $\approx 5\%$, $\langle\overline{\jmath}_{\rm rand}\rangle$ differs
by $\approx 20\%$. Comparing Model F to Model R (the half-resolution run): The
lower-resolution simulation has lower Mach numbers at small $r$, which results
in a $\approx 50\%$ lower specific angular momentum at small radii.

###  4.3 Comparison to the Analytical Estimate for $j_{\rm rand}$

To compare our measured values of specific angular momentum to the analytical
estimate of eq. (2), we use the instantaneous profiles $\overline{\jmath}_{\rm
rand}(r,t)$ and $v_{\rm turb}(r,t)$ to compute the time-averaged profiles,
$\langle\overline{\jmath}_{\rm rand}\rangle$ and $v_{c}\equiv\langle v_{\rm
turb}\rangle$, respectively, over $4000<t/t_{0}<5000$. Then, $v_{c}$ and
$H=r/b$ are used in eq. (2) to give the analytical estimate for $j_{\rm
rand}$.

Figure 11: Comparison of the time- and spherically-averaged specific angular
momentum profile, $\langle\overline{\jmath}_{\rm rand}\rangle$, of Model F
(solid line) and Model A (dashed line), to the local analytical estimate. The
dot-dashed curve uses the time-averaged turbulent velocity,
$v_{c}\equiv\langle v_{\rm turb}\rangle$, to estimate $j_{\rm rand}$ according
to eq. (2). We multiply this estimate by the factor of $1/2$ found by Quataert
et al. (2019). Instead of the $j_{\rm rand}\propto r$ scaling of eq. (2), we
find that $\langle\overline{\jmath}_{\rm rand}\rangle$ is flatter in radius,
suggesting that the convective flows are global and have roughly constant
angular momentum over a range of radii, rather than being set locally at each
radius (as is assumed in eq. (2)).

The solid line in Fig. 11 plots $\langle\overline{\jmath}_{\rm rand}\rangle$
versus radius over $r_{\rm heat,out}<r<r_{\rm cool}$ for Model F (solid line)
and Model A (dashed line). The dash-dotted line is the analytical estimate for
$\langle\overline{\jmath}_{\rm rand}\rangle$ using $v_{c}$ computed from Model
A. In the plot, we include the factor of $1/2$ that was found to describe the
local simulations of Quataert et al. (2019, see their equation 4). Eq. (2) is
a local estimate for $j_{\rm rand}$ in a given shell as it assumes the value
of $j_{\rm rand}$ just depends on the local values of $v_{c}$ and $H$. The
local estimate seems to capture our simulation results at large $r$ well.
However, the $H\propto r$ scaling in the power-law convective envelope gives
eq. (2) a strong radial dependence that is absent from our simulated
$\langle\overline{\jmath}_{\rm rand}\rangle$ over a large range in $r$.
Instead, $\langle\overline{\jmath}_{\rm rand}\rangle$ is essentially flat in
radius above $r\approx 0.6r_{0}$. At smaller $r$,
$\langle\overline{\jmath}_{\rm rand}\rangle$ falls off with decreasing $r$,
though with a shallower slope than the local estimate. By $r=0.12r_{0}$, there
is an order-of-magnitude difference between $\langle\overline{\jmath}_{\rm
rand}\rangle$ and the local estimate.

We interpret our simulation results to indicate that $j_{\rm rand}$ is being
set at large $r$ (where there is good agreement with the local estimate) and
$j_{\rm rand}$ is roughly conserved as material flows to smaller $r$. We noted
in the context of Fig. 4 that the convective cells in our simulation appear
coherent over a large fraction of the stellar radius. Our results suggest that
these large structures approximately conserve $j_{\rm rand}$, significantly
increasing the angular momentum content of the flows at smaller radii. The
recent radiation hydrodynamic simulations of Goldberg et al. (2021) find
similar results (see their sec. 4.2 for a detailed comparison).

As shown in Fig. 1, the base of the convective hydrogen envelope for our MESA
RSG model is located at $r=0.03r_{0}$ in our code units. If we extrapolate the
$\langle\overline{\jmath}_{\rm rand}\rangle$ curves of Fig. 11 to
$\log(r/r_{0})=-1.52$, we find $\langle\overline{\jmath}_{\rm
rand}\rangle=10^{-3}j_{0}$. Using eq. (26) with $r_{\rm ph}=840R_{\odot}$ and
$M_{\bullet}=6M_{\odot}$, then $\langle\overline{\jmath}_{\rm
rand}\rangle\approx 8.8\times 10^{16}$ cm2 s-1. For comparison, for a non-
spinning $6M_{\odot}$ BH, $j_{\rm ISCO}\approx 9\times 10^{16}$ cm2 s-1. This
extrapolation of our results shows that $j_{\rm rand}$ is important even near
the base of the hydrogen envelope. Material near the base of the convective
hydrogen zone is the most bound and the most likely to survive the weak
explosion in a FSN, though we reiterate that most of the convective envelope
we simulate will remain bound during a FSN (see panel (d) of Fig. 1).

The time-averaged profiles considered here and in the previous section are
useful for characterizing the properties of the flow, but the instantaneous
specific angular momentum in the envelope is what matters for the collapse of
the star. The collapse begins from a particular state of the flow and the
specific angular momentum accreting to small radii as a function of time
depends on the state of the envelope at the start of collapse, as we explore
with our collapse simulations in the next section.

##  5 Collapse Simulations

In this section, we present the results of our collapse simulations, which are
listed in Table 2. Each collapse run is initialized from a restart file that
was output during our fiducial convection run, Model F. The time associated
with the restart file, $t_{s}$, is the start time of the collapse run, at
which time the sink with radius $r_{s}$ is activated at the origin. The
physically relevant portion of our domain extends to $r=5r_{0}$ so we stop the
collapse simulation before shells that began at $r\gtrsim 5r_{0}$ begin to
accrete. The rarefaction wave that is launched when the sink is introduced
travels out at the local sound speed. Once the rarefaction wave reaches a
particular $r$, that shell begins to fall in towards the sink. The integrated
sound-crossing time from the sink radius to $r=5r_{0}$ is $10.8t_{0}$ (eq. 46
with $b=2.1$ and $\gamma=1.4762$). The infall time (from rest) from $5r_{0}$
back to $r_{s}$ is $14.1t_{0}$ for a total accretion time of $\approx 25t_{0}$
(see eq. 48). We thus run our collapse simulations until $t-t_{s}=25t_{0}$. To
be conservative, given that the material has variations in initial $v_{r}$, we
only use data to a time of $t-t_{s}=20t_{0}$.

Collapse Runs 1 and 2 are initialized at $t_{s}=4560$ and $4900$ $t_{0}$,
respectively, when the convective flow has somewhat higher specific angular
momentum with maximum $\overline{\jmath}_{\rm rand}(r,t_{s})$ of $\approx
0.028-0.032$ $j_{0}$ (see Fig. 8). Collapse Runs 3 and 4, with $t_{s}=4270$
and $4800$ $t_{0}$, respectively, begin from states of somewhat lower specific
angular momentum with maximum $\overline{\jmath}_{\rm rand}(r,t_{s})$ of
$\approx 0.015-0.02$ $j_{0}$. Section 5.1 presents the results of collapse
Runs 1-4. Collapse Run 1s is identical to collapse Run 1, except that we adopt
a smaller sink size and, in Section 5.3, we show that our results are
converged with respect to $r_{s}$. Section 5.2 considers the extent to which
material is restructured during the collapse by comparing the measured
accretion rates to rates predicted from assuming ballistic infall of the
convective material.

###  5.1 Flow and Instantaneous Accretion Rates

Figure 12: Slices of local Mach number, $\mathcal{M}=v/c_{s}$, through the
$z=0$ plane for a simulation with an absorbing sink at small radii that
initiates collapse of the convective envelope (collapse Run 1 in Table 2). The
snapshots show $5$, $10$, and $15t_{0}$ after the sink is activated at
$t_{s}=4560t_{0}$. The transonic rarefaction wave separates the subsonic
convective background flow from the supersonic inflow of material roughly in
free-fall towards the low-pressure sink. The somewhat non-spherical
rarefaction wave and infalling material carry the imprints of the initial
state of the convective flows.

Fig. 12 shows snapshots of Run 1 at $5$, $10$, and $15$ $t_{0}$ after the
collapse begins at $t_{s}=4560t_{0}$. Here we plot slices of the gas Mach
number, $\mathcal{M}=v/c_{s}$, through the $z=0$ plane. The outgoing
rarefaction wave (roughly the outer edge of the white region) moves out
through the subsonic convective flow, enclosing a region of supersonic infall
as material falls back towards the sink. The rarefaction wave is nearly but
not exactly spherical and instead reflects the differences in the local radial
velocity of the turbulent background. For example, the flow at positive $x$
shows the imprint of two streams of material that were already falling back
towards the origin before the introduction of the sink. As expected for
$\gamma<5/3$ the gas reaching the sink is supersonic.

Figure 13: Instantaneous accretion rates of mass and angular momentum for the
collapse of Model F starting at four different times. Panel (a): mass
accretion rates. Panels (b), (c), and (d): magnitude, polar angle, and
azimuthal angle, respectively, of $\boldsymbol{j}_{\rm
acc}\equiv\boldsymbol{{\dot{J}}}/\dot{m}$, the specific angular momentum
vector of the accreted material. The $x$-axis is the elapsed time since the
start of the collapse. Panels (a) and (b) can be scaled to physical units via
eqs. (27) and (26), respectively. The $\phi_{j}$ values for Runs 2 and 4 have
been given small phase shifts to avoid the discontinuity $\phi_{j}=0=2\pi$
where the values would otherwise cross the $x$-axis.

Fig. 13 shows the rates of mass and angular momentum accretion for Runs 1 - 4
as a function of elapsed time, $t-t_{s}$, after introduction of the sink. The
rates reported here are computed by measuring the fluxes through the sink
surface at $r_{\rm flux}=r_{s}=0.08r_{0}$. The mass accretion rate, panel (a),
falls off in time because the mean density profile falls off with radius as
$\propto r^{-2.1}$. With $j_{\rm rand}\ll j_{0}$ and $\mathcal{M}_{\rm
turb}<1$ neither centrifugal nor turbulent pressure modify the infall so there
is no difference in $\dot{m}$ between collapse runs started from different
convection states, even with a factor of 2 difference in maximum
$\overline{\jmath}_{\rm rand}(r,t_{s})$.

The remaining panels of Fig. 13 characterize the specific angular momentum
vector of the accreted material, $\mn@boldsymbol{j}_{\rm
acc}\equiv\dot{\boldsymbol{J}}/\dot{m}$, over time. Panels (b), (c), and (d)
show the magnitude, $j_{\rm acc}$, the polar angle, $\theta_{j}$, and the
azimuthal angle, $\phi_{j}$, respectively (the direction angles are defined as
in Section 4.1.3). Panel (b) shows that, critically, the infall of the
envelope has not restructured the material in such a way as to erase or
diminish the $j_{\rm rand}$ resulting from the convective motions in the
envelope. Indeed, $j_{\rm acc}$ over time for all runs is just as expected
from Fig. 8 with $0.005\lesssim j_{\rm acc}\lesssim 0.03$ for Runs 1 and 2,
and with $0.002\lesssim j_{\rm acc}\lesssim 0.02$ for Runs 3 and 4.

The direction of $\mn@boldsymbol{j}_{\rm acc}$ as a function of time is
important for determining the ultimate fate of material as it plunges to yet
smaller radii. Focusing on the thick, solid line for Run 1, panels (c) and (d)
show that the largest variations in $\theta_{j}$ and $\phi_{j}$ occur over
$\sim 5-10$ $t_{0}$ while smaller-amplitude changes take place on $\sim t_{0}$
timescales. Other curves show one or two large-amplitude changes on somewhat
shorter timescales. For example, Run 4 (dotted, pink lines) shows large swings
in $\theta_{j}$ and $\phi_{j}$ over the first $\sim 3t_{0}$ of the collapse
before transitioning to a phase where the angles are relatively constant.
Recall, however, that the dynamical time at $r_{\rm circ}$ is $10^{-6}t_{0}\ll
t_{0}$. Taken together, these curves show that the magnitude, direction, and
variability of $\mn@boldsymbol{j}_{\rm acc}$ depend on the state of the
envelope at the start of the collapse. Overall, however, the direction of
$\mn@boldsymbol{j}_{\rm acc}$ tends to be a slowly varying function over the
time it takes for a large fraction of the envelope to reach the sink.

###  5.2 Comparison to Semi-Analytical Predictions

An important goal of this work is to understand the extent to which the mean
angular momentum associated with convective flows is modified during infall.
We explore in this section whether the mean profiles from the convection
simulation, that is $\overline{\jmath}_{\rm rand}$ of Fig. 8 and
$\overline{\theta}_{j}$ and $\overline{\phi}_{j}$ of Fig. 9, can be used to
predict the accretion rates that are realized in the collapse simulation under
the assumption of ballistic infall of each shell. There is significant
dispersion in $v_{r}$, $j_{x}$, $j_{y}$ and $j_{z}$ in each radial shell about
the mean profiles that we have shown thus far (see Fig. 14, which shows
histograms of $j_{z}$ values of material about to enter the sink in our
collapse run). One can imagine that if material with the largest $v_{r}$
consistently has, e.g., $j_{z}$ that is far from the mean, then material
reaching the sink at a given time could be sampling the tail of the $j_{z}$
distribution of many different shells at the same time, and therefore would
not match our ballistic prediction using the mean profiles.

Figure 14: Histogram of $j_{z}$ for grid cells in a shell just outside the
sink ($r_{s}=0.08r_{0}$) during collapse simulation Run 1. Each histogram is
for a different simulation time (listed in the legend) and is computed using
all grid cells in the shell with $0.1<r/r_{0}<0.32$. All cells have the same
volume and the inflow is supersonic in this region. The times chosen are
spaced such that new material has entered the region (and the old material has
been accreted by the sink) when each histogram is computed. While the mean
values in each shell (vertical lines) are $-0.0024,0.0021,$ and $0.012$
$j_{0}$, the width of the distribution is a factor of 10 larger. The width of
the corresponding distribution in circularization radius is roughly $100$
times the mean circularization radius of the shell. The parcels of gas in each
shell will begin to circularize at a large range of radii, which will result
in a complicated flow as centrifugal support becomes important.

We now briefly explain how we predict accretion rates from the time-dependent
mean radial profiles of mass density, $\overline{\rho}(r,t)$, angular momentum
density, $\overline{\ell}(r,t)\equiv\overline{\rho}(r,t)\overline{\jmath}_{\rm
rand}(r,t)$, and direction angles $\overline{\theta}_{j}(r,t)$ and
$\overline{\phi}_{j}(r,t)$ (see Appendix A for a more detailed explanation and
a simpler test case in which the profiles are not changing in time). Upon
introduction of the sink, the outgoing rarefaction wave travels at the sound
speed and a given shell begins the infall once the rarefaction wave arrives.
The time it takes for the rarefaction wave to reach $r$ is the integrated
sound-crossing time, $t_{\rm wave}(r)$ given in eq. (46). The shell falls from
rest from $r$ over a time $1.14t_{\rm ff}(r)$ where $t_{\rm ff}(r)$ is the
integrated free-fall time, eq. (47). The factor of $1.14$ (supplied by Eric R.
Coughlin and based on the self-similar rarefaction wave solutions of Coughlin
et al. 2019), comes from the fact that infall begins from a profile in
hydrostatic equilibrium. There is thus a pressure gradient behind the
rarefaction wave, so the infall is not true zero-pressure ‘free-fall’ as
assumed by $t_{\rm ff}(r)$. The total time for the shell initially at $r$ to
accrete is the sum $t_{\rm acc}(r)=t_{\rm wave}(r)+{1.14}t_{\rm ff}(r)$.

In our ballistic prediction, we assume that the mass and angular momentum that
arrive at the sink at time $t_{\rm acc}(r_{i})$ are just the amount of mass
and angular momentum that were contained in the shell at $r_{i}$ when the
rarefaction wave arrived. These values are computed from the time-dependent
mean profiles from the convection simulation. We then assume the computed
value of the quantity of interest in that shell, e.g. the total mass, does not
change as that shell falls in over the subsequent time of ${1.14}t_{\rm
ff}(r_{i})$.

Figure 15: Comparison of the predicted accretion rates of mass and angular
momentum with those measured in the collapse simulation. Solid, black curves
are the instantaneous accretion results for Run 1 shown in Fig. 13. The
purple, dotted curves use the time-dependent, spherically-averaged profiles
from convection Model F to predict the accretion results by accounting for the
travel time of the rarefaction wave from the sink to a given shell and
assuming each shell then undergoes ballistic infall. The agreement between the
curves in each panel shows that there is little restructuring of
$\mn@boldsymbol{j}_{\rm rand}$ during this phase of infall in which $r>r_{\rm
circ}$.

The purple, dotted curves in Fig. 15 plot the predicted values of the mass
accretion rate, $\dot{m}$, as well as the predicted magnitude and direction
angles of the specific angular momentum vector of the accreted material for
Model F assuming collapse begins at $t_{s}=4560t_{0}$. For comparison, the
accretion results realized in the equivalent collapse simulation (Run 1) are
shown with solid, black lines in each panel.

Panel (a) shows that $\dot{m}$ is predicted extremely well and serves as
independent confirmation of the Coughlin et al. (2019) rarefaction wave
solutions. While there are some small differences between the two curves in
panels (b)-(d), the magnitude and direction angles of $\mn@boldsymbol{j}_{\rm
acc}$ are well-predicted using the mean profiles and the assumption of
ballistic infall from rest. This good agreement with the actual accretion
rates is non-trivial, especially given the large dispersion of $v_{r}$ and the
components of $\mn@boldsymbol{j}_{\rm rand}$ in each shell as well as the fact
that all of the material has non-zero $v_{r}$ rather than starting the infall
from rest, as we assumed.

We interpret the excellent agreement between measured and predicted
$\mn@boldsymbol{j}_{\rm acc}$ to mean that during the transition from
convection to infall and in the subsequent supersonic infall, very little
restructuring of the material takes place. This should remain the case until
another source of dynamical support becomes important, e.g. centrifugal
pressure. This will not occur until much smaller radii than simulated here
($r_{\rm circ}\ll r_{s}$).

Although we do not show the results here, we also computed predictions for
comparison to collapse Run 3 (that is, computing the prediction assuming
collapse of Model F from $t_{s}=4270t_{0}$). We found the same level of
agreement between predicted and actual accretion rates as in Fig. 15.

###  5.3 Dependence on Sink Size

In this section, we compare our measurement of the circularization radius of
the accreted material, $r_{\rm circ}(t)=j_{\rm acc}^{2}/GM$, for collapse Runs
1 and 1s which differ only in choice of sink radius ($r_{s}=0.08$ and $0.04$
$r_{0}$, respectively). For these simulations, we measure mass and angular
momentum accretion fluxes at both the sink radius, $r_{\rm flux}=r_{s}$, and
at a larger radius of $r_{\rm flux}=0.48r_{0}$.

The two solid lines in Fig. 16 plot $r_{\rm circ}$ of the accreted material as
measured at $r_{\rm flux}=0.48r_{0}$. The circularization radii measured at
$0.48r_{0}$ in the two simulations agree to better than $0.15\%$, indicating
that the flow has already converged with respect to sink size by
$r_{s}=0.08r_{0}$.

The dotted line in Fig. 16 shows the accretion rate for Run 1s measured at
$r_{\rm flux}=r_{s}=0.04r_{0}$. The curve is shifted by the finite travel time
from $0.48$ to $0.04$ $r_{0}$ ($\approx 0.18t_{0}$) but is otherwise the same
as the other two curves.

Figure 16: Comparison of the circularization radius of accreted material for
simulations with different sink sizes initiated from convection Model F at
$t_{s}=4560t_{0}$. The thick grey and thin black curves, which lie on top of
each other, adopt sink sizes of $r_{s}=0.04$ and $0.08$ $r_{0}$, respectively,
and show $r_{\rm circ}$ computed using the accretion rates measured at
$0.48r_{0}$ in both simulations. The dotted and thick grey lines show $r_{\rm
circ}$ for the simulation with $r_{s}=0.04r_{0}$ computed using the accretion
rates measured at $0.04$ and $0.48$ $r_{0}$, respectively. They are the same
except for a small shift in time due to the finite travel time from $0.48$ to
$0.04$ $r_{0}$.

##  6 Discussion

In this section we consider our results in the context of FSN, scaling our
convection results to a typical RSG and discussing the possible character of
accretion and outflows during the infall of the convective hydrogen envelope.

###  6.1 Application to Supergiants

In code units, the circularization radius of material with specific angular
momentum $j_{\rm rand}$ is

$r_{\rm circ}=({j_{\rm rand}}/{j_{0}})^{2}r_{0}.$ (35)

Making use of $r_{0}=r_{\rm ph}/6$ to write the previous expression in
physical units,

$\displaystyle r_{\rm circ}=1.16\times 10^{9}\bigg{(}\frac{j_{\rm
rand}}{10^{-2}\,j_{0}}\bigg{)}^{2}\bigg{(}\frac{r_{\rm
ph}}{1000R_{\odot}}\bigg{)}{\rm cm}.$ (36)

For a non-spinning BH, $r_{\rm ISCO}=6GM_{\bullet}/c^{2}$, so the ratio of
circularization radius to the radius of the ISCO is

$\frac{r_{{\rm circ}}}{r_{\rm ISCO}}=130\,\bigg{(}\frac{j_{\rm
rand}}{10^{-2}\,j_{0}}\bigg{)}^{2}\bigg{(}\frac{r_{\rm
ph}}{1000R_{\odot}}\bigg{)}\bigg{(}\frac{10M_{\odot}}{M_{\bullet}}\bigg{)}.$
(37)

Typical values for the fraction inside the first set of parenthesis are
roughly 1-3 (see Fig. 8 and panel (b) of Fig. 13).

Figure 17: Spherically-averaged circularization radius profile,
$\overline{r}_{\rm circ}(r,t)$, as a function of radius and time for our
fiducial convection model. The left- and right-hand colorbar labels report,
respectively, the values in code units (eq. 35) and as a multiple of the ISCO
radius (eq. 37). We take $r_{\rm ph}=840R_{\odot}$ and
$M_{\bullet}=6M_{\odot}$ in eq. (37), appropriate for the MESA model of Fig.
1. White patches in the plot indicate regions that fall below the colorbar
floor of $10$ $r_{\rm ISCO}$. For nearly all radii and at all times,
$\overline{r}_{\rm circ}(r,t)$ of each shell is $\sim 10-1500\times r_{\rm
ISCO}$.

Fig. 17 shows the mean profiles of circularization radius for the convective
material, $\overline{r}_{\rm circ}(r,t)$, which are computed using
$\overline{\jmath}_{\rm rand}(r,t)$ of Fig. 8 in eq. (35). The left side of
the colorbar is in code units. The right side converts those code values to
$\overline{r}_{\rm circ}/{r_{\rm ISCO}}$ using eq. (37) and adopting $r_{\rm
ph}\approx 840R_{\odot}$ and $M_{\bullet}=6M_{\odot}$, as in Fig. 1. The white
regions are where $\overline{r}_{\rm circ}$ is below the colorbar lower-limit
of $10r_{\rm ISCO}$ (if we set this lower bound to $r_{\rm ISCO}$ instead,
only a few white points would remain). The figure shows that turbulent Mach
numbers of $\sim 0.1$ (Fig. 5) correspond to circularization radii that are
many hundreds to a thousand times larger than $r_{\rm ISCO}$. In addition,
$\overline{r}_{\rm circ}\gtrsim 10{r_{\rm ISCO}}$ at nearly all radii and all
times. So although $\mn@boldsymbol{\overline{\jmath}}_{\rm rand}$ is a
function of both radius and time, the convective material begins the collapse
with $\overline{r}_{\rm circ}\gtrsim 10{r_{\rm ISCO}}$ almost independent of
when the collapse of the envelope begins.

We further showed in Section 5.2 that there is little restructuring of the
specific angular momentum as the matter infalls to smaller radii. Instead, we
found that both the magnitude and direction of $\mn@boldsymbol{j}_{\rm acc}$
are well-predicted by the mean profiles of the envelope prior to collapse
(Fig. 15). We therefore conclude that, although the times at which we started
the collapse in our collapse runs were arbitrary with respect to the flow,
Fig. 17 shows that $r_{\rm circ}$ of the convective material will be well
outside of $r_{\rm ISCO}$ no matter when the collapse begins.

###  6.2 Implications for Failed Supernovae

In Sections 4.3 and 6.1, we showed that the entire hydrogen convective zone
has $j_{\rm rand}\gtrsim j_{\rm ISCO}$. Based on our simulation results, we
argued that the local estimate of $\sim\frac{1}{2}$eq. (2) seems to set
$j_{\rm rand}$ at large radii in the convection zone and that $j_{\rm rand}$
is roughly conserved down to smaller radii. The helium layer is also
convective, though with smaller convective velocities, $v_{c}$. If we assume
that for the helium convective zone $j_{\rm rand}$ is related to $v_{c}$ by
$\frac{1}{2}$eq. (2), then $\overline{\jmath}_{\rm rand}\leq 5\times 10^{14}$
cm2 s-1 or $\overline{\jmath}_{\rm rand}\leq 10^{-2}j_{\rm ISCO}$ for our MESA
RSG model (see also fig. 3 of Quataert et al. 2019). This value is consistent
with simulations of convection in the helium layer performed by Gilkis & Soker
(2016) for which the mean profiles of the Cartesian components of
$\mn@boldsymbol{j}_{\rm rand}$ were $\lesssim 7\times 10^{14}$ cm2 s-1 in
magnitude or $\lesssim 0.03j_{\rm ISCO}$ for their assumed BH mass. The other
convective regions interior to the helium layer have even smaller
$\overline{\jmath}_{\rm rand}$. We conclude that, for non-rotating stars,
convective material interior to the hydrogen convection zone is likely to
accrete spherically onto the BH. We note that this conclusion still holds even
if the convective Mach number in the oxygen layer is a factor of $4-5$ larger
than captured in the MESA model, as found in the 3D simulations of Fields &
Couch (2020).

One subtlety related to the infall of the helium layer is the large dispersion
in $j_{\rm rand}$ about the mean $\overline{\jmath}_{\rm rand}$ in each radial
shell of the convective zone (see Fig. 14). Gilkis & Soker (2016) find that
some parcels in their helium layer have $j_{\rm rand}\sim 10^{16}$ cm2 s-1
$\sim j_{\rm ISCO}$; they argue that these individual parcels could generate
an outflow even though the majority of the infalling material has $j_{\rm
rand}\ll j_{\rm ISCO}$. It is unclear whether localized regions with $j_{\rm
rand}>j_{\rm ISCO}$ in the helium layer can reverse the inflow of the bound
envelope or are simply advected into the BH along with the bulk of the
material. This requires further study. In what follows, we focus on the infall
of the hydrogen envelope onto the newly-formed BH.

For our MESA RSG model, panel (d) of Fig.1 shows that the binding energy of
the hydrogen envelope is $E_{\rm bind}\approx 2\times 10^{48}$ erg. Some of
the hydrogen envelope may be ejected due to weak shocks induced by neutrino-
cooling during the proto-NS phase. The arrows in panel (d) of Fig. 1 show the
amount of mass that could be unbound based on the range of shock energies
found by Ivanov & Fernández (2021). For the lowest shock energies, only
$0.5M_{\odot}$ is likely to be unbound (open arrows) while for the highest
shock energy, $3.6M_{\odot}$ of the envelope could be ejected (filled arrows).
In either case, at least $6.9M_{\odot}$ of the convective hydrogen zone
remains bound to the BH, which represents the bulk of the convection zone in
our simulations.

For $r_{\rm circ}\approx 500r_{\rm ISCO}\approx 0.05R_{\odot}$, as is typical
in our simulations, $E_{\rm bind}$ for the entire hydrogen envelope could be
supplied by accretion of an amount of mass

$\delta m\sim\frac{r_{\rm circ}E_{\rm bind}}{GM_{\bullet}}\sim 3\times
10^{-3}M_{\odot}.$ (38)

(if inflow makes it to $r_{\rm ISCO}$, then $\delta m\approx
10^{-5}M_{\odot}$). For our MESA model, the accretion rates of Fig. 13
translate to $\dot{m}\sim 20M_{\odot}$ yr-1. At that rate, $\delta m\sim
10^{-3}M_{\odot}$ is accreted after roughly half an hour.

At these super-Eddington accretion rates of $\approx$ $10-100$ $M_{\odot}$
yr-1, the flow is optically thick and unable to cool by radiation. With
$\dot{m}\ll 10^{-3}M_{\odot}$ s-1 there is no neutrino cooling either
(Beloborodov, 2008), so accretion of the hydrogen envelope is an optically-
thick, advection-dominated accretion flow (optically-thick ADAF; Begelman &
Meier, 1982; Abramowicz et al., 1988; Narayan & Yi, 1994). As material falls
to small radii without the ability to cool, gravitational potential energy can
only be converted into kinetic and thermal energy and the material may just
return to large radii in an outflow. In addition to the collimated outflows
and disk winds that super-Eddington disks inevitably produce (Jiang et al.,
2019), accretion of stellar material with fixed rotation axis and roughly
uniform circularization radius sets up an accretion shock that sweeps through
the infalling material and can unbind the outer parts of the star (Lindner et
al., 2010; Murguia-Berthier et al., 2020).

If an accretion shock does not unbind all of the envelope and material can
circularize into a disk, then preferential outflow into a region along the
angular momentum axis would allow accretion to continue for longer periods
than could be possible for more spherical outflow. The relevant timescale in
this case is that over which the direction of $\mn@boldsymbol{j}_{\rm acc}$
varies. For example, if the coherence time for the direction of
$\mn@boldsymbol{j}_{\rm acc}$ is $\delta t\sim 2t_{0}\approx 25$ days (the
approximate time during which $\theta_{j}$ changes by $\lesssim\pi/2$ for all
curves in Fig. 13), then $\delta m\sim\dot{m}\delta t\sim 1.4M_{\odot}$ would
have time to fall in. If a fraction $f$ (possibly $\ll 1$) of this material
reaches $r_{\rm ISCO}$ for a $6M_{\odot}$ BH, the change in potential energy
would be $|\delta U|\sim 4\times 10^{53}f$ erg. Envelope gas within the funnel
region would easily be unbound as this energy drives outflows from small
radii. If the orientation of the disk changes on the same timescale as that of
$\mn@boldsymbol{j}_{\rm acc}$ ($\gtrsim 2t_{0}\sim 1$ month), then the
collimated outflow could sweep out a large fraction of $4\pi$ on this
timescale, ejecting most of the material and resulting in an energetic, long-
duration transient.

The timescale and energetics inferred in the previous paragraph are similar to
what is needed to power extremely long-duration gamma-ray transients such as
Swift 1644+57 (Bloom et al., 2011) and Swift J2058.4+0516 (Cenko et al., 2012)
by stellar core-collapse (e.g., Quataert & Kasen, 2012; Woosley & Heger,
2012). The biggest uncertainty in this application of our results is whether
the tenuously-bound hydrogen envelope is unbound before much of it can accrete
onto the newly-formed BH.

####  6.2.1 BH Spin

The efficiency of tapping accretion power to unbind the hydrogen envelope
determines not only how much material is returned to the star’s environment
and the nature of the transient produced, but it also determines the final
spin of the BH. We can make a rough estimate of the mass and final BH spin by
assuming the BH accretes mass at a rate of $f\dot{m}$, where $\dot{m}$ is the
rate that matter falls in from the envelope and $0<f\leq 1$, and by computing
the angular momentum accretion rate of the BH according to

$\dot{\boldsymbol{J}_{\bullet}}=\begin{cases}f\dot{m}\mn@boldsymbol{j}_{\rm
acc}&\text{if }j_{\rm acc}<j_{\rm ISCO}\\\ f\dot{m}j_{\rm
ISCO}\hat{\jmath}_{\rm acc}&\text{if }j_{\rm acc}\geq j_{\rm ISCO}\end{cases}$
(39)

(here $\hat{\jmath}_{\rm acc}$ is the unit vector in the direction of
$\mn@boldsymbol{j}_{\rm acc}$ and $\dot{m}(t)$ and $\mn@boldsymbol{j}_{\rm
acc}(t)$ are computed from the instantaneous profiles from our convection
simulations using the methods of Sec 5.2). By limiting $j_{\rm acc}$ to
$j_{\rm ISCO}$, this estimate accounts for the fact that an accretion disk or
outflow transports angular momentum to infinity to allow a fraction $f$ of the
infalling material to reach $r_{\rm ISCO}$.

To compute the spin of the BH, we assume that everything interior to the
convective hydrogen envelope is accreted and carries no angular momentum, so
we initialize the BH mass to $M_{\bullet}=6M_{\odot}$ and the BH angular
momentum to $\boldsymbol{J}_{\bullet}=\mathbf{0}$. Working outward from the
base of the hydrogen envelope, we update $M_{\bullet}$ given $\dot{m}(t)$ and
$\boldsymbol{J}_{\bullet}$ is increased using eq. (39). Throughout this
process, $j_{\rm ISCO}$ is evolved with the mass and spin of the BH using the
standard relations (evaluated at the equator) for a Kerr BH. We do not account
for the angle between the spin vector of the BH and the angular momentum
vector of the shell except to check whether the dot product is positive or
negative. If the dot product is positive or zero, we assume the material is on
a prograde orbit when computing $j_{\rm ISCO}$. For a negative dot product, we
assume a retrograde orbit. We also ensure that the spin of the BH never
exceeds the maximum value $J_{\rm max}=GM_{\bullet}^{2}/c$, though this ends
up not being necessary because $J_{\bullet}/J_{\rm max}$ is always $\lesssim
0.8$.

Figure 18: Dimensionless spin parameter of the BH versus the radius up to
which the convective envelope is incorporated into the BH. Here radius refers
to the location in the envelope before collapse begins and the spin at each
$r$ is the spin that the BH would have if everything exterior to $r$ were
unbound. Grey shading indicates where $r>r_{\rm cool}$. The mass and angular
momentum in each accreted shell are calculated using the instantaneous
profiles from the convection simulations as described in Sec. 5.2. The upper
panel assumes all of the mass in each shell is accreted ($f=1$ in eq. 39). The
legend gives the convection model and assumed start time for the collapse,
$t_{s}$, that were used to compute each curve. In all cases, the BH spin
approaches $\approx 80\%$ of maximum as more of the envelope is accreted
before declining (due to a sign change in the accreted
$\dot{\boldsymbol{J}}$). Although the total angular momentum in the envelope
is zero, the BH ends up with non-zero spin because only a small fraction of
the specific angular momentum in each shell can be accreted by the BH. For
partial accretion of the mass in each shell ($f<1$, lower panel) the BH
achieves smaller maximum spin at the peak of the curve. For $0.03<f<1$, the
maximum spins are $\sim 0.06-0.59$. Weak shocks like those studied by Ivanov &
Fernández (2021) imply BH spins near the peak of each curve with spins of
$\sim 0.5$ when $f\sim 1$ and spins of $\sim 0.1$ when $f\sim 0.1$.

In Fig. 18, we plot the dimensionless spin parameter of the BH,
$J_{\bullet}/J_{\rm max}$, as a function of the radius $r$ out to which the
envelope is incorporated into the BH.555The final density profile in the
Athena models differs a bit from the initial profile (see Fig. 1). In
evaluating the BH spin in Fig. 18, we renormalize the density by a factor of
0.75 to keep the total mass within $r=5r_{0}$ the same as the initial
condition. The spin calculated without this renormalization is very similar to
that shown in Fig. 18, so our conclusions are not sensitive to this choice.
The $r$ coordinate refers to the radius of the material before the collapse.
So at each $r$ in the figure, the value of the spin assumes material located
at radii $>r$ prior to collapse is unbound and a fraction $f$ of the material
that had radius $<r$ before collapse is accreted onto the BH. The curves in
the upper panel correspond to a different convection model and start time for
the collapse, as noted in the legend, for the case in which $f=1$. The spin
parameter increases slowly as material with $j\sim j_{\rm ISCO}$ is added to
the initially non-spinning BH. For all the models shown, the spin parameter
saturates near $\sim 0.8$ before slowly decreasing as material with opposite
spin direction is accreted. If material out to $r_{\rm cool}$ is accreted, the
final spin parameter is $0.3-0.7$. There is some angular momentum exterior to
$r_{\rm cool}$ (grey shaded region), so we extend the curves to $r=9r_{0}$,
where the cumulative profiles of Fig. 7 approach zero. If all of this material
is incorporated, the final spin parameter is $0.3-0.4$. The spin parameter
does not return to zero (as one might expect from Fig. 7) because only a small
fraction of the total available specific angular momentum in each shell can be
incorporated into the BH. As a result, there is uneven cancelation of the
angular momentum vectors across many different shells, unlike the integral
over the star at fixed time that is shown in Fig. 7.

The BH spin versus $r$ curves shown in the upper panel of Fig. 18 represent
upper limits on the spin for two reasons. First, we assumed that all mass
reaching small radii accretes. In reality, some of this mass will go into the
outflow and will not accrete onto the BH. Second, once an outflow from small
radii occurs, much of the hydrogen envelope will be blown away, thus modifying
the $\dot{m}$ that can be fed from large radii by the envelope. Batta &
Ramirez-Ruiz (2019) considered these two effects in more detail in the context
of rapidly rotating collapsars; here we fold both effects into the factor $f$
in eq. (39). To explore the influence of a reduced accretion rate on the final
BH mass and maximum spin, we repeated each of the calculations of the upper
panel of the figure with $f=0.03,0.1$, and $0.3$, instead of $f=1$. The lower
panel of Fig. 18 shows the resulting curves for the particular case of
collapse of Convect Model F from $t_{s}=4560t_{0}$. Across all models, the
lower values of $f$ give maximum spin parameters of $0.06-0.1$, $0.19-0.29$,
and $0.45-0.59$, respectively, and final BH masses (integrating out to $r_{\rm
cool}$) of $\sim 6.6$, $7.3$, and $9.4$ $M_{\odot}$, respectively. For
comparison, $f=1$ gives a final BH mass of $\sim 16.8M_{\odot}$ and maximum
spin parameters of $0.81-0.85$. If weak shocks, like those studied by Ivanov &
Fernández (2021), unbind $\sim$ few $M_{\odot}$ of material, the BH spin would
be $\sim 0.1-0.5$, near the peak of the curves in Fig. 18.

####  6.2.2 Timescale for Angular Momentum Redistribution

Our calculations neglect self-gravity; thus it is useful to check whether
internal gravitational torques can significantly reduce the mean specific
angular momentum of a shell before the shell reaches $r_{\rm circ}$. Using
snapshots from collapse Run 1, we computed the timescale for angular momentum
redistribution for individual shells of material outside of the sink due to
torques from the rest of the gas that is in supersonic infall. At
$t=4578t_{0}$, the material is supersonic out to $r\approx 4.5r_{0}$, so we
compute the vector torque on a shell with $0.1\leq r/r_{0}\leq 0.2$ due to the
gas with $0.08\leq r/r_{0}\leq 4$. We do similar calculations for different
shell radii, different regions of the background gas, and for different
snapshots ($t=4570$ and $4575$ $t_{0}$). In all cases, the timescale for
angular momentum redistribution is a few thousand times the free-fall time of
the shell. So gravitational torques cannot alter the angular momentum vector
of the infalling gas before centrifugal pressure becomes important at $\sim
r_{\rm circ}$.

In additional to centrifugal pressure, magnetic fields can play an important
role as the material continues to infall. RSGs may have $\sim 1$G magnetic
fields at their surface (Aurière et al., 2010; Petit et al., 2013; Tessore et
al., 2017). Assuming flux freezing holds, so that $B\propto r^{-2}$, then the
magnetic pressure is equal to the ram pressure of the infalling matter
(assumed to be in free-fall) at a radius of

$r_{B}\sim 0.09\bigg{(}\frac{B_{0}}{1{\rm
G}}\bigg{)}^{\frac{4}{3}}\bigg{(}\frac{r_{\rm
ph}}{1000R_{\odot}}\bigg{)}^{\frac{8}{3}}\bigg{(}\frac{10M_{\odot}/{\rm
yr}}{\dot{m}}\bigg{)}^{\frac{2}{3}}\bigg{(}\frac{10M_{\odot}}{M_{\bullet}}\bigg{)}^{\frac{1}{3}}R_{\odot}$
(40)

where $B_{0}$ is the field strength at $r_{\rm ph}$. For our $18M_{\odot}$ RSG
and for the range of accretion rates in Fig. 13, $r_{B}\sim 0.01-0.07$
$R_{\odot}$ which is roughly $500r_{\rm ISCO}\sim 0.05R_{\odot}$, comparable
to $r_{\rm circ}$ realized in our models. Magnetic fields may thus become
dynamically important during the infall. We note that the field in the
interior of the star may be much larger than 1$G$. For example, if $\rho
v_{\rm turb}^{2}\sim B^{2}/8\pi$, then in the outer $500R_{\odot}$ of our MESA
RSG model, $B\sim 200-800$ $G$, meaning that magnetic fields could be yet more
important during collapse.

##  7 Summary and Conclusions

A fraction of RSGs and YSGs may end their lives in failed supernovae (FSNe),
in which core collapse does not lead to a successful, neutrino-powered SN
explosion. Even if weak shocks are launched by radiation of neutrino energy
prior to the NS collapsing to a BH, much of the hydrogen envelope remains
bound and will fall into the newly-formed BH (Ivanov & Fernández, 2021); see
Fig. 1 panel (d) and associated discussion. Previous work by Gilkis & Soker
(2014, 2016) and Quataert et al. (2019) has shown that, even when the star has
zero net angular momentum, the random velocity field in the convective
hydrogen envelope gives rise to finite specific angular momentum (at each
radius) that is larger than that associated with the ISCO of the BH. This
suggests that accretion disks generically form during the infall of the
hydrogen envelope in FSN.

We perform two sets of 3D hydrodynamical simulations to study the random
angular momentum associated with convective flows in the context of FSN. We
first simulate convection in polytropic models that are applicable to the
convective envelopes of RSGs and YSGs. Our convection simulations provide the
initial conditions for a set of collapse simulations, which follow the infall
of the convective material after we introduce a low-pressure sink at the
origin that mimics core collapse. Our convection simulations extend the work
of Quataert et al. (2019) by simulating a significant fraction of the hydrogen
envelope and in the spherical geometry of a star. We confirm their finding
that the specific angular momentum associated with convective flows is larger
than the specific angular momentum of the ISCO of a $10M_{\odot}$ BH, although
we find a different scaling with radius than inferred from their results,
described in more detail below. We further show with our collapse simulations
that the specific angular momentum of the convective flows is largely
conserved during the infall of the material at radii larger than the
circularization radius of the gas. The direction of the specific angular
momentum vector is slowly varying over the timescales relevant to the region
where the material will circularize. All of this implies that, during the
collapse of a supergiant in a FSN, the random angular momentum of the
convective hydrogen envelope is likely to lead to the formation of
centrifugally-supported gas at small radii that could drive outflows and
generate an observable transient. The following summarizes our main results:

  * •

Convective flows in supergiant envelopes give rise to finite angular momentum
at each radius even when the total angular momentum of the envelope is zero
(Fig. 7).

  * •

For polytropic models with convective Mach numbers of $0.05<\mathcal{M}_{\rm
turb}<0.25$, consistent with the hydrogen envelopes of RSGs and YSGs (Fig. 6),
the convective flows give rise to specific angular momenta of $j_{\rm
rand}\approx(4-28)\times 10^{17}$ cm2 s-1, where we have used eq. (26) to
scale our results to a RSG with photosphere radius $840R_{\odot}$ and an
assumed BH mass of $6M_{\odot}$. These values correspond to circularization
radii relative to the BH ISCO of 10 $\lesssim r_{\rm circ}/r_{\rm
ISCO}\lesssim$ 1500 (Fig. 17).

  * •

At the largest radii in the convection zone in our simulation, $j_{\rm rand}$
agrees with the local simulations and analytical scaling (given in our eq. 2)
of Quataert et al. (2019). However, we find a different scaling with $r$ than
suggested by this local estimate. In our simulations, the convective flows
roughly conserve $j_{\rm rand}$ down to small radii. This results in our
simulated $j_{\rm rand}$ being larger than the estimate of eq. (2) over a
large fraction of the envelope (Fig. 11). This is an important result because
the outer radii of the convective zone are more likely to be unbound by weak
shocks generated during collapse. Our simulations imply larger $j_{\rm rand}$
and a larger likelihood of circularization at precisely the radii that are
more likely to remain bound and accrete onto the BH.

  * •

The specific angular momentum of the convective envelope is not significantly
altered during collapse. Instead, mean profiles of the specific angular
momentum vector $\mn@boldsymbol{j}_{\rm rand}$ (in magnitude and direction)
from the convection simulations can be used to predict the specific angular
momentum vector of the accreted material measured in the collapse simulations
($\mn@boldsymbol{j}_{\rm acc}\equiv\dot{\boldsymbol{J}}/\dot{m}$) by assuming
ballistic infall and accounting for the finite sound-travel and infall time
for each shell (Fig. 15).

  * •

The specific angular momentum of the accreted material depends on the state of
the envelope at the start of collapse. However, even when collapse begins at
one of the lowest angular momentum states (Run 4 with collapse start time
$t_{s}=4800t_{0}$), $r_{\rm circ}>r_{\rm ISCO}$. Indeed, because the mean
profiles from the convection simulation can be used to accurately predict the
accretion rates of mass and angular momentum flowing to small radii in our
collapse calculations, we can conclude from Fig. 17 that $r_{\rm circ}$ of the
accreted hydrogen envelope is always greater than $\approx 10$ $r_{\rm ISCO}$
for a $6M_{\odot}$ BH.

  * •

The direction of $\mn@boldsymbol{j}_{\rm acc}$ is slowly varying (Fig. 13) on
the timescales that we are sensitive to. These timescales, which are set by
convective flow times in the envelope, are very long relative to the dynamical
time at the circularization radius of the material. Because the direction of
the mean specific angular momenta of the convective envelope are not
significantly modified during the infall, it is likely that the direction of
$\mn@boldsymbol{j}_{\rm acc}$ remains coherent down to $r_{\rm circ}$ scales,
allowing for the presence of coherent rotation at small radii. This supports
the idea that accretion at small radii will generate outflows that can carry
energy to large radii. The resulting accretion dynamics at small radii are
likely to be very different from standard accretion simulations because there
is a very large dispersion in accreted $\mn@boldsymbol{j}_{\rm acc}$ in each
spherical shell (Fig. 14). As material reaches $r_{\rm circ}$, the
distribution of $r_{\rm circ}$ within each shell will drive interactions
between individual parcels, which will generate mixing, torques, and shock-
heating as individual parcels are deflected from spherical infall. The outcome
of these flows requires further study.

  * •

Accretion of the entire envelope by the BH leads to finite BH spin, even
though the total angular momentum of the envelope is zero. The BH spin is
$\sim 0.5$ if most of the envelope accretes, but less if outflows at small
radii remove significant mass (Fig. 18). The finite BH spin occurs because
$j_{\rm acc}>j_{\rm ISCO}$, the specific angular momentum of the ISCO, in
nearly all of the infalling shells. The BH cannot accrete more than $j_{\rm
ISCO}$ from each shell, so infalling shells with identical but oppositely-
oriented specific angular momenta can never cancel.

Our results imply that much of the hydrogen envelope that remains bound in FSN
is likely to become rotationally-supported during collapse rather than fall
directly into the BH. This should apply both to the initial collapse and to
any subsequent fallback. Sufficient mass is available to accrete onto the
newly-formed BH that it could power a month $-$ year long transient with
energies up to $\sim 10^{52}$ ergs (Sec. 6.2). This could be the origin of
unusually long-duration gamma-ray transients such as Swift 1644+57 (Bloom et
al., 2011) and Swift J2058.4+0516 (Cenko et al., 2012) and unusually energetic
Type II supernovae (e.g., OGLE-2014-SN-073 Moriya et al., 2018). The key
uncertainty is that the binding energy of the hydrogen envelope is only $\sim
10^{48}$ ergs and so a small amount of accretion may unbind the envelope
before sufficient mass accretes to power an energetic transient. The FSN
candidate discovered by Gerke et al. (2015) showed a modest brightening up to
$\sim 10^{6}L_{\odot}$ for a year before fading to a luminosity significantly
below that of the progenitor star (Adams et al., 2017). The initial
brightening prior to the star ‘disappearing’ is consistent with a $\sim$ few
$\times 10^{48}$ erg shock unbinding the hydrogen envelope (e.g., Lovegrove &
Woosley, 2013) but not with a more energetic transient. Given the highly
super-Eddington accretion rates and significant rotational support implied by
our results, it is also unclear how rapidly the star is expected to fade away
after core collapse, or whether it might remain a near-Eddington source for an
extended period of time powered by super-Eddington fallback accretion.
Understanding this better is important for understanding the FSN candidate
discovered by Gerke et al. (2015) and for identifying and interpreting
additional candidates in ongoing and future transient surveys. We expect a
large diversity in the observational manifestation of RSG collapse to a BH
given the stochastic nature of the angular momentum in the convective
envelope.

There are several key aspects of supergiant convection and FSN that are not
captured in our study. The random convective angular momentum in simulations
is a more subtle effect than what is usually studied in simulations of
convection. Increasing resolution somewhat increases the time-averaged $j_{\rm
rand}$ but it is not known what resolution is required to fully converge.
Additionally, we have neglected the phase of the FSN when a low-energy shock
wave sweeps through the convective hydrogen envelope. The outgoing shock could
restructure the convective material and will unbind some fraction of the
hydrogen envelope. Simulations that follow the propagation of the weak shock
through the convective envelope, the subsequent expansion of the still-bound
material to apocenter, and the ultimate infall and accretion onto the BH are
needed to self-consistently model FSN. Finally, future work will extend our
simulations by following the inflow to radii where rotational support becomes
dynamically important. Such studies are critical for understanding the extent
to which the random angular momentum associated with convection leads to flows
at small radii that can drive energy and momentum to large radii and may be
capable of producing a transient.

## Acknowledgements

We thank the anonymous referee, Lars Bildsten, Eric R. Coughlin, Paul C.
Duffell, Drummond Fielding, Jared A. Goldberg, Yan-Fei Jiang, Wenbin Lu,
Morgan MacLeod, Philipp Moesta, Enrico Ramirez-Ruiz, Sean Ressler, and Stephen
Ro for useful conversations and suggestions. A.A. gratefully acknowledges
support from the University of California, Berkeley Fellowship, the Cranor
Fellowship at U.C. Berkeley, and the National Science Foundation Graduate
Research Fellowship under Grant No. DGE 1752814. E.Q. was supported in part by
a Simons Investigator award from the Simons Foundation. This work benefited
from workshops supported by the Gordon and Betty Moore Foundation through
Grant GBMF5076.

The simulations presented in this article were performed on computational
resources managed and supported by Princeton Research Computing, a consortium
of groups including the Princeton Institute for Computational Science and
Engineering (PICSciE) and the Office of Information Technology’s High
Performance Computing Center and Visualization Laboratory at Princeton
University. Preliminary studies used the Savio computational cluster resource
provided by the Berkeley Research Computing program at the University of
California, Berkeley (supported by the UC Berkeley Chancellor, Vice Chancellor
for Research, and Chief Information Officer). This project was made possible
by the following publicly available software: Astropy666http://www.astropy.org
(Astropy Collaboration et al., 2013, 2018), Athena++ (Stone et al., 2020),
matplotlib (Hunter, 2007), MESA (Paxton et al., 2011; Paxton et al., 2013,
2015, 2018, 2019), MESA SDK (Townsend, 2020), NumPy (Harris et al., 2020),
yt777https://github.com/yt-project/yt (Turk et al., 2011).

## Data Availability

Data related to the results of simulations in this article will be shared on
reasonable request to the corresponding author via e-mail.

## References

  * Abramowicz et al. (1988) Abramowicz M. A., Czerny B., Lasota J. P., Szuszkiewicz E., 1988, ApJ, 332, 646
  * Adams et al. (2017) Adams S. M., Kochanek C. S., Gerke J. R., Stanek K. Z., Dai X., 2017, MNRAS, 468, 4968
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, aj, 156, 123
  * Aurière et al. (2010) Aurière M., Donati J. F., Konstantinova-Antova R., Perrin G., Petit P., Roudier T., 2010, A&A, 516, L2
  * Basinger et al. (2020) Basinger C. M., Kochanek C. S., Adams S. M., Dai X., Stanek K. Z., 2020, arXiv e-prints, p. arXiv:2007.15658
  * Batta & Ramirez-Ruiz (2019) Batta A., Ramirez-Ruiz E., 2019, arXiv e-prints, p. arXiv:1904.04835
  * Begelman & Meier (1982) Begelman M. C., Meier D. L., 1982, ApJ, 253, 873
  * Beloborodov (2008) Beloborodov A. M., 2008, in Axelsson M., ed., American Institute of Physics Conference Series Vol. 1054, Cool Discs, Hot Flows: The Varying Faces of Accreting Compact Objects. pp 51–70 (arXiv:0810.2690), doi:10.1063/1.3002509
  * Bloom et al. (2011) Bloom J. S., Butler N. R., Cenko S. B., Perley D. A., 2011, GRB Coordinates Network, 11847, 1
  * Burrows et al. (2020) Burrows A., Radice D., Vartanyan D., Nagakura H., Skinner M. A., Dolence J. C., 2020, MNRAS, 491, 2715
  * Cenko et al. (2012) Cenko S. B., et al., 2012, ApJ, 753, 77
  * Chiavassa et al. (2009) Chiavassa A., Plez B., Josselin E., Freytag B., 2009, A&A, 506, 1351
  * Chiavassa et al. (2010) Chiavassa A., Haubois X., Young J. S., Plez B., Josselin E., Perrin G., Freytag B., 2010, A&A, 515, A12
  * Couch et al. (2015) Couch S. M., Chatzopoulos E., Arnett W. D., Timmes F. X., 2015, ApJ, 808, L21
  * Coughlin et al. (2018a) Coughlin E. R., Quataert E., Fernández R., Kasen D., 2018a, MNRAS, 477, 1225
  * Coughlin et al. (2018b) Coughlin E. R., Quataert E., Ro S., 2018b, ApJ, 863, 158
  * Coughlin et al. (2019) Coughlin E. R., Ro S., Quataert E., 2019, ApJ, 874, 58
  * Dupree et al. (2020) Dupree A. K., et al., 2020, ApJ, 899, 68
  * Ertl et al. (2016) Ertl T., Janka H. T., Woosley S. E., Sukhbold T., Ugliano M., 2016, ApJ, 818, 124
  * Fernández et al. (2018) Fernández R., Quataert E., Kashiyama K., Coughlin E. R., 2018, MNRAS, 476, 2366
  * Fields & Couch (2020) Fields C. E., Couch S. M., 2020, ApJ, 901, 33
  * Gerke et al. (2015) Gerke J. R., Kochanek C. S., Stanek K. Z., 2015, MNRAS, 450, 3289
  * Gilkis & Soker (2014) Gilkis A., Soker N., 2014, MNRAS, 439, 4011
  * Gilkis & Soker (2016) Gilkis A., Soker N., 2016, ApJ, 827, 40
  * Goldberg et al. (2021) Goldberg J. A., Jiang Y.-F., Bildsten L., 2021, arXiv e-prints, p. arXiv:2110.03261
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Heger et al. (2005) Heger A., Woosley S. E., Spruit H. C., 2005, ApJ, 626, 350
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Ivanov & Fernández (2021) Ivanov M., Fernández R., 2021, ApJ, 911, 6
  * Jiang et al. (2019) Jiang Y.-F., Stone J. M., Davis S. W., 2019, ApJ, 880, 67
  * Lindner et al. (2010) Lindner C. C., Milosavljević M., Couch S. M., Kumar P., 2010, ApJ, 713, 800
  * Linial et al. (2021) Linial I., Fuller J., Sari R., 2021, MNRAS, 501, 4266
  * Lovegrove & Woosley (2013) Lovegrove E., Woosley S. E., 2013, ApJ, 769, 109
  * Lovegrove et al. (2017) Lovegrove E., Woosley S. E., Zhang W., 2017, ApJ, 845, 103
  * MacFadyen & Woosley (1999) MacFadyen A. I., Woosley S. E., 1999, ApJ, 524, 262
  * Moriya et al. (2018) Moriya T. J., Terreran G., Blinnikov S. I., 2018, MNRAS, 475, L11
  * Murguia-Berthier et al. (2020) Murguia-Berthier A., Batta A., Janiuk A., Ramirez-Ruiz E., Mandel I., Noble S. C., Everson R. W., 2020, ApJ, 901, L24
  * Nadyozhin (1980) Nadyozhin D. K., 1980, Ap&SS, 69, 115
  * Narayan & Yi (1994) Narayan R., Yi I., 1994, ApJ, 428, L13
  * Neustadt et al. (2021) Neustadt J. M. M., Kochanek C. S., Stanek K. Z., Basinger C. M., Jayasinghe T., Garling C. T., Adams S. M., Gerke J., 2021, arXiv e-prints, p. arXiv:2104.03318
  * O’Connor & Ott (2011) O’Connor E., Ott C. D., 2011, ApJ, 730, 70
  * Papish & Soker (2011) Papish O., Soker N., 2011, MNRAS, 416, 1697
  * Parrish et al. (2008) Parrish I. J., Stone J. M., Lemaster N., 2008, ApJ, 688, 905
  * Paxton et al. (2011) Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton B., et al., 2013, ApJS, 208, 4
  * Paxton et al. (2015) Paxton B., et al., 2015, ApJS, 220, 15
  * Paxton et al. (2018) Paxton B., et al., 2018, ApJS, 234, 34
  * Paxton et al. (2019) Paxton B., et al., 2019, ApJS, 243, 10
  * Perna et al. (2018) Perna R., Lazzati D., Cantiello M., 2018, ApJ, 859, 48
  * Petit et al. (2013) Petit P., Aurière M., Konstantinova-Antova R., Morgenthaler A., Perrin G., Roudier T., Donati J. F., 2013, Magnetic Fields and Convection in the Cool Supergiant Betelgeuse. p. 231, doi:10.1007/978-3-642-30648-8_9
  * Piro (2013) Piro A. L., 2013, ApJ, 768, L14
  * Powell & Müller (2020) Powell J., Müller B., 2020, MNRAS, 494, 4665
  * Quataert & Kasen (2012) Quataert E., Kasen D., 2012, MNRAS, 419, L1
  * Quataert et al. (2019) Quataert E., Lecoanet D., Coughlin E. R., 2019, MNRAS, 485, L83
  * Ro et al. (2019) Ro S., Coughlin E. R., Quataert E., 2019, ApJ, 878, 150
  * Soker (2010) Soker N., 2010, MNRAS, 401, 2793
  * Stone et al. (2008) Stone J. M., Gardiner T. A., Teuben P., Hawley J. F., Simon J. B., 2008, ApJS, 178, 137
  * Stone et al. (2020) Stone J. M., Tomida K., White C. J., Felker K. G., 2020, ApJS, 249, 4
  * Sukhbold & Adams (2020) Sukhbold T., Adams S., 2020, MNRAS, 492, 2578
  * Sukhbold et al. (2016) Sukhbold T., Ertl T., Woosley S. E., Brown J. M., Janka H. T., 2016, ApJ, 821, 38
  * Sukhbold et al. (2018) Sukhbold T., Woosley S. E., Heger A., 2018, ApJ, 860, 93
  * Tessore et al. (2017) Tessore B., Lèbre A., Morin J., Mathias P., Josselin E., Aurière M., 2017, A&A, 603, A129
  * Townsend (2020) Townsend R., 2020, MESA SDK for Linux, doi:10.5281/zenodo.3706650, https://doi.org/10.5281/zenodo.3706650
  * Turk et al. (2011) Turk M. J., Smith B. D., Oishi J. S., Skory S., Skillman S. W., Abel T., Norman M. L., 2011, The Astrophysical Journal Supplement Series, 192, 9
  * Ugliano et al. (2012) Ugliano M., Janka H.-T., Marek A., Arcones A., 2012, ApJ, 757, 69
  * Woosley (1993) Woosley S. E., 1993, ApJ, 405, 273
  * Woosley & Heger (2012) Woosley S. E., Heger A., 2012, ApJ, 752, 32
  * Yadav et al. (2020) Yadav N., Müller B., Janka H. T., Melson T., Heger A., 2020, ApJ, 890, 94
  * Zhang et al. (2008) Zhang W., Woosley S. E., Heger A., 2008, ApJ, 679, 639

##  Appendix A Method for Predicting Accretion Rates

In this section we discuss how mass and angular momentum accretion rates can
be predicted from the instantaneous snapshot of the convective flow at the
start of collapse. Our interest here is in cases where the gas is roughly in
hydrostatic equilibrium at the start of collapse (i.e this is not free-fall-
from-rest collapse; the gas has pressure support). Here we use a simple test
problem to illustrate the procedure.

For the test problem, we use a power-law density profile in hydrostatic
equilibrium with a point mass potential. Functionally, the potential, the
density profile, and the adiabatic sound speed profile are, respectively,

$\displaystyle\Phi(r)$ $\displaystyle=-GM/r$ (41) $\displaystyle\rho(r)$
$\displaystyle=\rho_{0}(r_{0}/r)^{b}$ (42) $\displaystyle c_{s}(r)$
$\displaystyle=\sqrt{\frac{\gamma}{(b+1)}\frac{GM}{r}},$ (43)

where $r_{0}$, $\rho_{0}$, and $b$ are constants and $\gamma$ is the adiabatic
index of the gas. Code units are the same as defined in Section 3.6. We use
$b=2.5$ and $\gamma=4/3$ for our test problem. At initialization, we give all
of the gas in the domain the same angular velocity, $\Omega$, about the
z-axis. This gives a $z$-angular momentum density profile of
$\frac{2}{3}\Omega r^{2}\rho(r)$. The mass contained in the shell at radius
$r$ and thickness $\Delta r$ is

$\Delta m(r)=4\pi r^{2}\rho(r)\Delta r$ (44)

and the z-angular momentum in the shell is

$\Delta J_{z}(r)=4\pi r^{2}\bigg{[}\frac{2}{3}\Omega
r^{2}\rho(r)\bigg{]}\Delta r=\frac{8\pi}{3}\Omega r^{4}\rho(r)\Delta r.$ (45)

To predict the accretion rate, we compute the time it takes the shell at $r$
to reach the sink (with radius $r_{s}$) after the sink is activated. Call this
function $t_{\rm acc}(r)$. We also need to know the time over which that shell
accretes, call that $\Delta t_{\rm acc}(r)$. When we activate the sink, a
rarefaction wave begins at $r_{s}$ and propagates outward in radius. The time
it takes the rarefaction wave to reach $r$ is the integrated sound-crossing
time from $r_{s}$ to $r$,

$t_{\rm
wave}(r)=\int_{r_{s}}^{r}\frac{dr^{\prime}}{c_{s}(r^{\prime})}=\frac{2}{3}\bigg{(}\frac{b+1}{\gamma
GM}\bigg{)}^{1/2}\bigg{[}r^{3/2}-r_{s}^{3/2}\bigg{]}.$ (46)

If the pressure behind the rarefaction wave were zero, then the shell would
fall from rest from its initial position at $r$ to the sink radius $r_{s}$
over the integrated free-fall time,

$t_{\rm
ff}(r)=\bigg{(}\frac{r^{3}}{2GM}\bigg{)}^{1/2}\Bigg{\\{}\frac{\pi}{2}-\sin^{-1}\bigg{(}\sqrt{\frac{r_{s}}{r}}\bigg{)}+\sqrt{\frac{r_{s}}{r}\Big{(}1-\frac{r_{s}}{r}\Big{)}}\Bigg{\\}}.$
(47)

There is, in fact, a small pressure gradient behind the rarefaction wave
(Coughlin et al., 2019), so the time for the shell to fall back after the
arrival of the rarefaction wave is actually $1.14\,t_{\rm ff}(r)$ for our
simulated values of $b$ and $\gamma$ (E. R. Coughlin, private communication).
The total time required for the shell at $r$ to reach the sink is thus

$t_{\rm acc}(r)=1.14\,t_{\rm ff}(r)+t_{\rm wave}(r).$ (48)

and

$\Delta t_{\rm acc}(r)=t_{\rm acc}(r+\Delta r/2)-t_{\rm acc}(r-\Delta r/2).$
(49)

We want accretion rates as a function of time, so we invert $t_{\rm acc}(r)$
to obtain $r_{\rm acc}(t)$. The analytical accretion rate predictions are thus

$\displaystyle\dot{m}(t)$ $\displaystyle={\Delta m(r_{\rm acc})}/{\Delta
t_{\rm acc}(r_{\rm acc})}$ (50) $\displaystyle\dot{J}_{z}(t)$
$\displaystyle={\Delta J_{z}(r_{\rm acc})}/{\Delta t_{\rm acc}(r_{\rm acc})}.$
(51)

Figure 19 compares eqs. (50) and (51) to the accretion rates measured in the
simulation. The predictions reproduce the simulated rates very well. The
bottom panel shows the $\%$ error between the predicted and simulated rates.