# In LIGO’s Sight? Vigorous Coherent Gravitational Waves from Cooled Collapsar
Disks

Ore Gottlieb<EMAIL_ADDRESS>Center for Computational
Astrophysics, Flatiron Institute, 162 5th Avenue, 6th floor, New York, NY
10010  Physics Department and Columbia Astrophysics Laboratory, Columbia
University, 538 West 120th Street, New York, NY 10027  Amir Levinson The
Raymond and Beverly Sackler School of Physics and Astronomy, Tel Aviv
University, Tel Aviv 69978, Israel  Yuri Levin Physics Department and Columbia
Astrophysics Laboratory, Columbia University, 538 West 120th Street, New York,
NY 10027  Center for Computational Astrophysics, Flatiron Institute, 162 5th
Avenue, 6th floor, New York, NY 10010  Department of Physics and Astronomy,
Monash University, Clayton, VIC 3800, Australia

###### Abstract

We present the first numerical study of gravitational waves (GWs) from
collapsar disks, using state-of-the-art 3D general-relativistic
magnetohydrodynamic simulations of collapsing stars. These simulations
incorporate a fixed Kerr metric for the central black hole (BH) and employ
simplified prescriptions for disk cooling. We find that cooled disks with an
expected scale height ratio of $H/R\gtrsim 0.1$ at $\sim 10$ gravitational
radii induce Rossby instability in compact, high-density rings. The trapped
Rossby vortices generate vigorous coherent emission regardless of disk
magnetization and BH spin. For BH mass of $\sim 10\,M_{\odot}$, the GW
spectrum peaks at $\sim 100\,{\rm Hz}$ with some breadth due to various non-
axisymmetric modes. The spectrum shifts toward lower frequencies as the disk
viscously spreads and the circularization radius of the infalling gas
increases. Weaker-cooled disks with $H/R\gtrsim 0.3$ form a low-density
extended structure of spiral arms, resulting in a broader, lower-amplitude
spectrum. Assuming an optimistic detection threshold with a matched-filter SNR
of 20 and a rate similar to supernovae Ib/c, LIGO-Virgo-KAGRA (LVK) could
detect $\lesssim 1$ event annually, suggesting that GW events may already be
hidden in observed data. Third-generation GW detectors could detect dozens to
hundreds of collapsar disks annually, depending on the cooling strength and
the disk formation rate. The GW amplitudes from collapsar disks are $\gtrsim
100$ times higher with a substantially greater event rate than those from
core-collapse supernovae, making them potentially the most promising burst-
type GW class for LVK and Cosmic Explorer. This highlights the importance of
further exploration and modeling of disk-powered GWs, promising insights into
collapsing star physics.

##  1 Introduction

The past decade has marked the beginning of the gravitational wave (GW) era,
during which the LIGO-Virgo-KAGRA (LVK; Accadia et al., 2012; Acernese et al.,
2015; LIGO Scientific Collaboration et al., 2015; Abbott et al., 2018; Akutsu
et al., 2021) collaboration has detected cosmic ripples arising from
inspiraling sources (Acernese et al., 2019; Tse et al., 2019; Buikema et al.,
2020). As we anticipate the upcoming phase with advanced LVK runs, followed by
third-generation GW detectors such as Cosmic Explorer (CE; Abbott et al.,
2017; Reitze et al., 2019; Evans et al., 2023), Einstein Telescope (ET;
Punturo et al., 2010a, b; Hild et al., 2011), and Laser Interferometer Space
Antenna (LISA; Babak et al., 2021), there is a keen interest in discovering
new GW sources, particularly those distinct from inspirals. The exploration of
new types of GWs holds the promise of unveiling new insights into the physics
of their progenitors and facilitating a deeper understanding of the diverse GW
events occurring throughout the Universe.

Similar to inspiralling sources, burst-type sources are astrophysical events
that generate transient GWs. However, their GW emission is less coherent,
posing challenges for creating deterministic waveform models and,
consequently, for their detection. A promising site for such burst-type
sources is powerful explosions such as core-collapse supernovae (CCSNe).
However, numerical simulations have shown that due to their quasi-spherical
structure, only a minute fraction of the CCSN energy is released in GWs, where
the dominant mechanism is stochastic excitation of the protoneutron star (PNS)
f- and g-modes, which are strongly coupled to gravitational radiation (e.g.,
Nakamura et al., 2016; Andresen et al., 2017, 2019, 2021; Yakunin et al.,
2017; Morozova et al., 2018; Radice et al., 2019; Abdikamalov et al., 2020;
Mezzacappa et al., 2020, 2023; Vartanyan et al., 2023; Mezzacappa & Zanolin,
2024; Powell & Müller, 2024; Richardson et al., 2024). The predicted small
excitation amplitudes imply that even the third-generation GW detectors can
only detect GWs powered by PNS accretion in the Milky Way and satellite
galaxies at design sensitivity (Srivastava et al., 2019; Evans et al., 2023;
Corsi et al., 2024).

If the progenitor star has sufficient angular momentum, its collapse will form
an accretion disk around a newly formed black hole (BH, “collapsar”).
Accretion disks have been considered promising GW emitters due to various
instabilities (Kobayashi & Mészáros, 2003; van Putten & Levinson, 2003;
Levinson et al., 2015; van Putten et al., 2019). If a global instability
develops, non-axisymmetric modes could give rise to density waves, introducing
a quadrupole moment that emits GWs. One criterion to examine the stability of
the disk is by comparing its pressure and gravity. The pressure gradient can
counteract the gravitational forces trying to pull material inward, thereby
preventing excessive collapse or disruption. If the disk becomes
gravitationally unstable (Chen & Beloborodov, 2007), it may fragment into
high-density clumps (Gammie, 2001), which might give rise to strong GW
emission (Piro & Pfahl, 2007). However, it remains to be seen whether such
accretion disks exist in astrophysical systems and the extent of their
asymmetry, highlighting the need for general-relativistic magnetohydrodynamic
(GRMHD) simulations.

Previous numerical attempts to assess the GW emission from accretion disks
employed a preset self-gravitating torus setup (Kiuchi et al., 2011; Korobkin
et al., 2011; Wessel et al., 2021). These studies found that non-axisymmetric
modes emerge in hydrodynamic simulations due to Papaloizou-Pringle
instabilities (PPI), which are quenched by magnetorotational instabilities
(MRI) when magnetic fields are present (Bugli et al., 2018; Wessel et al.,
2023). However, the relevance of torus studies to astrophysical systems such
as collapsars is questionable for several reasons: (1) Ad-hoc prescribed tori
evolve differently than accretion disks that form self-consistently during the
collapse; (2) Torus simulations lack the constant gas supply onto the disk and
the disk-envelope interaction; (3) Previous torus studies have not considered
cooling effects, which are expected in the high-density collapsar disks due to
neutrino emission (e.g., Narayan et al., 2001; van Putten & Levinson, 2003;
Kohri et al., 2005; Chen & Beloborodov, 2007; Just et al., 2022). Cooling
alters the disk thermodynamics and increases its density, potentially making
it the primary driving mechanism for the disk to become gravitationally
unstable (Batta & Lee, 2014), leading to vigorous GW emission.

Another potentially important source of density perturbations is the Rossby
wave instabilities (RWIs Lovelace et al., 1999; Tagger & Pellat, 1999; Li et
al., 2000, 2001). RWIs emerge due to sharp gradients in the disk, which act as
a potential well that can trap and amplify coherent structures like Rossby
vortices. Therefore, RWIs are expected to be particularly strong in
transiently forming disks. This insight motivates detailed first-principles
simulations of self-consistently forming disks without considering self-
gravity. Indeed, we find that the overdense Rossby vortices appear robustly in
our simulations, creating a substantial mass quadrupole, responsible for
generating GWs that may be detectable in LVK. The vigorous RWI-powered GWs
were overlooked in previous studies because their setups lacked cooling
mechanisms, necessary to increase the disk mass density and facilitate the
formation of RWIs.

It is worth mentioning that there may be other mechanisms by which collapsars
may produce strong GWs. For example, if a favorable magnetic field
configuration is present, the central BH will power a relativistic jet that
will generate a long gamma-ray burst (lGRB; Levinson & Eichler, 1993; Woosley,
1993), and a hot turbulent cocoon, powered by the jet-star interaction. The
aspherical nature and enormous energy of these outflows render the jet (Sago
et al., 2004; Birnholtz & Piran, 2013; Leiderschneider & Piran, 2021) and the
cocoon (Gottlieb et al., 2023b) interesting GW sources for LISA (Amaro-Seoane
et al., 2017) and CE/ET (Gottlieb et al., 2023b), respectively. While lGRB
jets and cocoons are more promising GW emitters compared to quasi-spherical
CCSNe, their scarcity reduces their GW detection prospects. Therefore, the
detectable GW signals from accretion disks may be appreciably more abundant
than those from relativistic outflows.

In this _letter_ , first in a series of papers, we begin to build the case for
serious and intense searches for GW bursts from collapsar disks by current and
future GW interferometer teams. The structure of the letter is as follows. We
outline the numerical setup and the GW calculation methodology in §2. In §3,
we describe the disk evolution and the emergence of instabilities that give
rise to non-axisymmetric modes. Subsequently, in §4, we present the resultant
vigorous GW signal, followed by a discussion on the expected event rates in
LVK and third-generation detectors in §5. We deliberate on the implications of
our results and conclude in §6.

##  2 Numerical simulations

###  2.1 Numerical setup

Model | Setup | $H/R$ | $\beta_{p}$ | ${\rm max}(\sigma_{0})$ | $\,{a_{\rm BH}}$ | $\,{M_{\rm BH}}\,[\,{M_{\odot}}]$ | $T_{s}\,[\rm s]$ | $R_{\rm max}\,[r_{g}]$ | $\,{E_{\rm GW}}\,[\,{{\rm erg}}]$ |  $\varrho$ (LVK) |  $\varrho$ (CE) | LVK rate [${\rm yr^{-1}}$] | $\,{f_{\rm GW}}\,[{\rm Hz}]$
---|---|---|---|---|---|---|---|---|---|---|---|---|---
$B$ | BNS merger | 0.1 | $10^{4}$ | - | 0.68 | 2.67 | $0.3$ | $10^{3}$ | $2\times 10^{46}$ | 0.1; 0.2 | 0.4; 1.4 | 0 | 500-2000
$C$ | Collapsar | 0.1 | - | $10^{-3}$ | 0.8 | 10 | $24$ | $10^{5}$ | $7\times 10^{50}$ | 25; 46 | 390; 750 | $\lesssim 1$ | 30-300
$Ca$ | Collapsar | 0.1 | - | $10^{-3}$ | 0.1 | 10 | 6.2 | $10^{5}$ | $1.5\times 10^{49}$ | 16; 33 | 180; 360 | $\lesssim 1$ | 200-300
$Cb$ | Collapsar | 0.1 | - | $10^{-4}$ | 0.8 | 10 | 7.4 | $10^{5}$ | $8\times 10^{49}$ | 27; 54 | 350; 690 | $\lesssim 1$ | 100-200
$Cc$ | Collapsar | 0.3 | - | $0$ | 0.8 | 10 | 5.4 | $10^{5}$ | $1.5\times 10^{49}$ | 5; 9 | 74; 130 | $\lesssim 10^{-2}$ | 100-200
Table 1:  A summary of the models’ parameters. The model names stand for the
system type: BNS merger ($B$) and collapsar ($C$): $C$ for our canonical
setup; $Ca$ for low $\,{a_{\rm BH}}$; $Cb$ for lower $b$; and $Cc$ for weaker
cooling. $H/R$ is the imposed disk scale height ratio which is effectively
higher in models $Ca$ and $Cb$, $\beta_{p}$ is the characteristic plasma beta
in the disk, ${\rm max}(\sigma_{0})$ is the maximum initial magnetization in
the star, $\,{a_{\rm BH}}$ is the BH dimensionless spin, $\,{M_{\rm BH}}$ is
the BH mass, $T_{s}$ is the simulation time, and $R_{\rm max}$ is the outer
radius of the simulation grid in units of gravitational radius, $r_{g}$.
Quantities measured from simulations: $\,{E_{\rm GW}}$ is the GW energy, and
$\varrho$ is the matched-filter SNR for edge-on and face-on observers,
respectively, shown for LVK A+ and CE at $D=10\,{\rm Mpc}$, and corrected for
accretion disk lifetime of $t_{\rm acc}\approx 100\,\,{{\rm s}}$ (see §5). The
LVK rate is calculated assuming the model represents all collapsar disks and
that their emergence rate is similar to SNe Ib/c rate. $\,{f_{\rm GW}}$ is the
characteristic GW frequency range.

A robust GW signal may emerge from compact, dense accretion disks. Such
conditions are prevalent in accretion disks around BHs formed due to the
gravitational collapse of a rotating stellar core. Similar conditions also
exist in disks formed following a merger of two compact objects, at least one
of which is an NS. Both of these systems exhibit high mass accretion rates of
$\dot{M}\gtrsim 10^{-2}\,\,{M_{\odot}}$ s-1, making them susceptible to
neutrino cooling. Analytic models (Chen & Beloborodov, 2007) and numerical
simulations incorporating neutrino transport (Siegel et al., 2019; Foucart,
2023) suggest that neutrino-cooled disks cool down to a characteristic ratio
of disk height to radius of $0.1\lesssim H/R\lesssim 0.3$, depending on the
effective viscosity and BH spin. $H/R$ tends to be smaller, $H/R\approx 0.1$,
in the innermost disk regions, where, as we shall see, most of the GW
radiation is emitted. Neutrino transport schemes are computationally expensive
and thus do not allow following the disk evolution over long timescales.
Therefore, our simulations do not incorporate full neutrino transport, but a
cooling prescription that artificially maintains a fixed $H/R$ (Noble et al.,
2009). This prescription mimics radiative cooling by reducing the internal
energy of the disk over the local Keplerian timescale until it matches a
target disk scale height.

The fiducial collapsar simulations (model $C$ in table 1) adhere to the setup
outlined in Gottlieb et al. (2022b), employing a stellar mass of
$M_{\star}=20\,\,{M_{\odot}}$, and maximum initial magnetization
($\sigma_{0}\equiv b^{2}/4\pi\rho c^{2}$, where $b$ is the comoving magnetic
field strength and $\rho$ is the comoving mass density) within the stellar
core capped at ${\rm max}(\sigma_{0})=10^{-3}$. The angular momentum profile
provides a fixed circularization radius of $r_{\rm circ}=25\,r_{g}$, similar
to what is suggested by stellar evolution models during the first several
seconds (Gottlieb et al. 2024, in prep.). The disk scale height is imposed to
maintain $H/R=0.1$ at all times. The BH mass and dimensionless spin parameter
are $\,{M_{\rm BH}}=10\,\,{M_{\odot}}$, and $\,{a_{\rm BH}}=0.8$,
respectively. In addition to our fiducial model $C$, we also consider a
slowly-spinning BH with $\,{a_{\rm BH}}=0.1$ (model $Ca$) and a weaker
$B$-field (model $Cb$). Both exhibit effectively thicker disks which increase
in height from $H/R=0.1$ close to the BH to $H/R\gtrsim 0.2$ at $R\gtrsim
20\,r_{g}$, consistent with neutrino-transport collapsar simulations (Issa et
al. 2024, in prep.). These disks are thicker because less gravitational energy
(due to increased ISCO in $Ca$) or magnetic dissipation ($Cb$) is available to
heat, making the cooling process less efficient. We also consider a weakly-
cooled disk with $H/R=0.3$ in model $Cc$ which increases to $H/R\approx 0.4$
at $r=20\,r_{g}$, similar to a collapsar disk without cooling in Gottlieb et
al. (2022a). We find that GWs from weakly-cooled magnetized disks with
$H/R=0.3$ are outshone by the cocoon-powered GWs. Therefore, to investigate
the GW signature of the weakly-cooled disk, we impose $b=0$ in model $Cc$ such
that no electromagnetically-driven jet and cocoon form. We outline the list of
simulations in Table 1.

In Appendix A, we present a post-merger simulation, which mirrors that
conducted by Gottlieb et al. (2023c). In short, it begins with a merger of NSs
of masses $1.06\,\,{M_{\odot}}$ and $1.78\,\,{M_{\odot}}$, described by the
LS220 equation of state for NSs (Lattimer & Swesty, 1991). The merger product
is a BH with a mass of $\,{M_{\rm BH}}=2.67\,\,{M_{\odot}}$ and $\,{a_{\rm
BH}}=0.68$. Surrounding the BH is a substantial accretion disk weighing
$M_{d}=0.096\,\,{M_{\odot}}$. The vector potential is defined by:

$A=A_{\phi}\propto\begin{cases}\rho&\rho/\rho_{\rm max}\geq 0.2\\\ &\\\
0&\rho/\rho_{\rm max}<0.2\\\ \end{cases}~{},$ (1)

and is imposed to be divergence-free. The potential normalization is chosen so
the gas-to-magnetic pressure ratio is $\beta_{p}\approx 10^{4}$.

We conduct the simulations using the 3D GPU-accelerated code h-amr (Liska et
al., 2022), employing an ideal equation of state with an adiabatic index of
$4/3$, suitable for radiation-dominated gas. The simulations adopt spherical
polar coordinates, $r$, $\theta$, $\phi$ with a logarithmic cell distribution
along the $r$-direction, spanning from just inside the event horizon to
$R_{\rm max}$, which for collapsars is located outside of the star radius (see
Tab. 1). Along $\theta$\- and $\phi$-directions, the cells are distributed
uniformly. In total, the base grid comprises $N_{r}\times N_{\theta}\times
N_{\phi}$ cells in the $r$-, $\theta$\- and $\phi$-directions, respectively.
For the collapsar and BNS merger simulations, $N_{r}=320,192$;
$N_{\theta}=96,96$; and $N_{\phi}=192,96$, respectively. We incorporate local
adaptive time-stepping and static mesh refinement. One (two) level of
refinement is applied in the innermost $100\,(50)\,r_{g}$ in all collapsar
(merger) simulations. We ensure that the wavelength of the fastest-growing MRI
mode in the magnetohydrodynamic simulations is adequately resolved when the
disk reaches a steady state evolution. In particular, we find that the MRI
quality factor $Q_{\rm MRI}$, denoting the number of cells per the MRI
wavelength, satisfies $Q_{\rm MRI}\gtrsim 30$ in both the $\hat{\theta}$ and
$\hat{\phi}$ directions, surpassing the minimum requirement of $Q_{\rm
MRI}\sim 10$ (e.g., Hawley et al., 2011).

###  2.2 GW approximations

The GW fields can be computed using the Green’s function method:

$h^{\mu\nu}(r,t)\approx\frac{4G}{c^{4}}\int{\frac{T_{eff}^{\mu\nu}\left(\boldsymbol{r}^{\prime},t-\lvert\boldsymbol{r}-\boldsymbol{r}^{\prime}\rvert/c\right)}{\lvert\boldsymbol{r}-\boldsymbol{r}^{\prime}\rvert}d^{3}r^{\prime}}\,,$
(2)

where $T_{eff}^{\mu\nu}=\sqrt{-g}T^{\mu\nu}+t^{\mu\nu}$, $T^{\mu\nu}$ is the
contravariant stress-energy tensor in Einstein equations, and $t^{\mu\nu}$ is
a pseudo-tensor composed of nonlinear combinations of $h_{\mu\nu}$ and their
derivatives (see, e.g., Sec. V.A in Thorne 1980 for details). Eq. (2)
constitutes an integral equation for $h^{\mu\nu}$ and is therefore not
practically useful for numerical computations. In the linear approximation,
where internal stresses of the source can be neglected, Eq. (2) can be
rendered applicable by omitting the term $t^{\mu\nu}$. However, this approach
is invalid for self-gravitating systems, such as the accretion disks under
consideration. For sufficiently slow systems, for which the reduced wavelength
largely exceeds the GW source’s characteristic size, a multipole emission
expansion can be sought. In the weak gravity limit, wherein $t^{0\mu}\ll
T^{0\mu}$, useful formulae for the mass quadrupole, mass octupole, and current
quadrupole moments, that account for gravitational stresses, can be derived
(see St. Germaine-Fuller et al., 2013). The mass quadrupole formula they
obtained can be reduced, in this regime, to the standard form (see also Misner
et al. 1973 for a similar derivation) that we use to compute the low-frequency
emission [see Eq. (6) below]. For the high-frequency emission, the multipole
expansion fails because the reduced wavelength $\lambda/(2\pi)$ becomes
comparable to the source size. We therefore have to use Eq. (2) directly to
compute the time-dependent perturbations to the metric. The right-hand side of
this equation contains the pseudo-tensor representing gravitational stresses
in the system; their direct computation goes beyond the scope of this paper.
Instead, we use Eq. (2) in the linear approximation and neglect $t^{\mu\nu}$.
We argue that this might be justified for unbound material, however, how
accurate is the high-frequency spectrum thereby obtained remains to be
assessed. It is reassuring however that the Green’s function and the
quadrupole approximation give consistent results in the intermediate range of
frequencies, as we show in Appendix B.

As stated above, the multipole expansion approximation is only applicable when
the reduced wavelength is larger than the size of the source along the line of
sight, e.g., for edge-on observers:

$2\pi\,{f_{\rm GW}}R_{d}\lesssim c\,,$ (3)

where $R_{d}$ is the density wave orbital radius. For the angular frequency
$\Omega$ of the emitting gas in the quasi-Keplerian disk:

$\Omega=\sqrt{\frac{G\,{M_{\rm BH}}}{R_{d}^{3}}}\,,$ (4)

the characteristic GW peak frequency is:

$\,{f_{\rm GW}}=\frac{\tilde{m}\Omega}{2\pi}\approx
180\tilde{m}\left(\frac{\,{M_{\rm
BH}}}{10\,\,{M_{\odot}}}\right)^{1/2}\left(\frac{100\,\,{{\rm
km}}}{R_{d}}\right)^{3/2}\,{\rm Hz}\,.$ (5)

where $\tilde{m}$ is the non-axisymmetric mode. Eq. (5) implies that the
multipole expansion remains valid when $R_{d}\gtrsim\tilde{m}^{2}r_{g}$. For
the quadrupole moment ($\tilde{m}=2$), this corresponds to $\,{f_{\rm
GW}}\lesssim 800\,{\rm Hz}$.

####  2.2.1 Comparison to CCSN studies

Previous CCSN studies (e.g., Nakamura et al., 2016; Andresen et al., 2017,
2019, 2021; Yakunin et al., 2017; Radice et al., 2019; Mezzacappa et al.,
2020, 2023; Vartanyan et al., 2023; Powell & Müller, 2024; Richardson et al.,
2024) have employed advanced numerical simulations to track the SN shock
expansion. By applying the quadrupole approximation in post-processing, they
have estimated the resulting GW signal. These studies have identified distinct
GW signals at both low and high frequencies.

The high-frequency GW emission ($\,{f_{\rm GW}}\gtrsim{\rm kHz}$) originates
from accretion onto the PNS and convection within the PNS. The compactness of
the PNS ensures that the quadrupole approximation remains valid even at high
frequencies [Eq. (3)]. This was recently confirmed by Longo Micchi et al.
(2023), who compared the high-frequency emission obtained by the quadrupole
approximation with the Neuman-Penrose formalism and found that the two methods
are consistent around the kilohertz peak.

The low-frequency emission ($\,{f_{\rm GW}}\lesssim 0.2\,{\rm kHz}$) arises
from the neutrino heating region below the SN shock wave. In these
simulations, the shock wave expands to distances of $\sim 10^{3}\,\,{{\rm
km}}$, for which Eq. (3) dictates that the quadrupole approximation is only
valid at $f\lesssim 50\,{\rm kHz}$. The invalidity of the quadrupole
approximation for the low-frequency peak calculation underscores the need for
more reliable calculations to verify the low-frequency emission.

####  2.2.2 Post-process calculation

We solve the GWs from the simulations by combining two methods: The quadrupole
approximation at frequencies below $\sim 0.2\,{\rm kHz}$, and the flat
spacetime Green’s function at high frequencies (see Appendix B). For the
quadrupole calculation, we follow Vartanyan et al. (2023) to calculate the
first time-derivative of the quadrupole111To validate the quadrupole
approximation, we verify that the octupole component in St. Germaine-Fuller et
al. (2013) can be neglected at $f\lesssim 800\,{\rm Hz}$.:

$\dot{Q}_{m,ij}(r,t)=\int{\rho^{\rm ret}(v_{i}^{\rm ret}x_{j}+v_{j}^{\rm
ret}x_{i}-\frac{2}{3}\delta_{ij}v_{r}^{\rm
ret}r)dx_{i}^{\prime}dx_{j}^{\prime}dx_{k}^{\prime}}\,,$ (6)

where all retarded quantities are defined as $q^{\rm
ret}=q\left(r^{\prime},t-\lvert r-r^{\prime}\rvert/c\right)$. Then, the strain
tensor components are:

$h_{\mu\nu}(r,t)\approx\frac{2G}{c^{4}D}\frac{d\dot{Q}_{m,\mu\nu}(r,t)}{dt}\,,$
(7)

where $D$ is the distance to the source.

For the Green’s function calculation, we follow Gottlieb et al. (2023b) who
calculated the emission from unbound ejecta, neglecting gravity. The strain
tensor components are determined using the far field limit of Eq. (2) in flat
spacetime, specifically:

$h_{\mu\nu}(r,t)\approx\frac{4G}{c^{4}D}\int{T_{\mu\nu}\left(\boldsymbol{r^{\prime}},t-\frac{\lvert\boldsymbol{r}-\boldsymbol{r^{\prime}}\rvert}{c}\right)dx_{i}^{\prime}dx_{j}^{\prime}dx_{k}^{\prime}}\,,$
(8)

where $T_{\mu\nu}$ is the covariant stress-energy tensor, and the integration
spans the volume of the disk.

Without the loss of generality, for an observer along the $\hat{x}$-axis, the
two independent solutions constituting the plane wave polarizations in the
transverse-traceless (TT) gauge metric are given by:

$\begin{split}h_{+}=h_{yy}-h_{zz}\,;\\\ h_{\times}=2h_{yz}\,.\end{split}$ (9)

Analogous solutions can be derived for observers situated along the
$\hat{y}$\- and $\hat{z}$-axes.

The characteristic strain is defined as (Moore et al., 2015):

$h_{c}(f)=2f\lvert\tilde{h}(f)\rvert\,,$ (10)

where

$\tilde{h}(f)=\sqrt{\tilde{h}_{+}^{2}(f)+\tilde{h}_{\times}^{2}(f)}\,,$ (11)

and $\tilde{h}_{\times}(f)$ and $\tilde{h}_{+}(f)$ represent the strain
polarizations in the frequency domain. To evaluate the signal detectability
and event rate, we estimate the matched-filter signal-to-noise ratio (SNR) as
(Flanagan & Hughes, 1998; Moore et al., 2015):

$\varrho=\sqrt{\int_{0}^{\infty}\frac{4|\tilde{h}_{A}(f)|^{2}}{S(f)}df}\,,$
(12)

where $S(f)$ denotes the noise power spectral density of the GW detector, for
which we use the PyCBC package (Usman et al., 2016; Nitz et al., 2021) for the
noise power-spectral densities. $\tilde{h}_{A}(f)$ is the dimensionless GW
strain, defined by:

$\tilde{h}_{A}(f)=F_{+}\tilde{h}_{+}(f)+F_{\times}\tilde{h}_{\times}(f)\,,$
(13)

where $F_{+},F_{\times}$ are functions of the antenna pattern of the detector,
which depend on the sensitivity to each polarization and the sky location of
the observer. For the SNR calculations, we adopt $F_{+}=1$ and $F_{\times}=0$.

##  3 Hydrodynamic evolution

Figure 1: 3D renderings illustrate the logarithmic mass density (in .c.g.s.)
of collapsar accretion disks. The length scale can be measured with respect to
the BH event horizon radius, represented by a black circle. (a) The early
formation stages of the disk ($\sim 1\,\,{{\rm s}}$) feature low-density
spiral density waves. (b) With the accumulation of infalling gas, a high-
density ring forms at $R\approx 15\,r_{g}$. RWIs in the disk give rise to non-
axisymmetric modes. Over time, the disk viscously spreads and its density
progressively drops. (c) In the absence of strong cooling, the disk exhibits
an extended structure [note that the horizon size in panels (a) and (b) is
considerably larger than that in panel (c)], and attains lower densities
compared to strongly-cooled disks. Consequently, weakly-cooled disks exhibit
spiral arms at all times, giving rise to a weaker GW signal. Full movies
showcasing the evolution of all collapsar accretion disks are available at
http://www.oregottlieb.com/disk_gw.html.

The first stellar shell with a circularization radius exceeding the innermost
stable circular orbit (ISCO) radius, will form an accretion disk. We simulate
the collapse of a spherically symmetric mass distribution star, incorporating
a solid body rotation profile that facilitates the rapid formation of an
accretion disk. Figure 1 presents 3D renderings depicting the logarithmic mass
density at different phases of strongly-cooled disk formation (panels a,b) and
for weakly-cooled disks (c). The strongly-cooled accretion disks (models
$C,Ca,Cb$) exhibit similar evolution. We identify two stages in the formation
of strongly-cooled disks in collapsars. As illustrated in Figure 1(a), during
the initial stages of disk formation ($\lesssim 1\,\,{{\rm s}}$), the disk is
composed of low-density plasma. In low-density environments, pressure and
density waves propagate freely, enabling differential rotation to amplify
small perturbations into large, coherent spiral arms. Consequently, at the
onset of disk formation, the radial expansion of the spiral waves results in a
roughly flat and attenuated GW spectrum.

Fig. 1(b) indicates that over time, the accumulation of gas in the disk
increases the disk density, leading to sharp pressure, density, and magnetic
field gradients that render the disk Rossby unstable (e.g., Tagger & Pellat,
1999) – The centrifugal force acts against the gradient to trap a pronounced
high-density Rossby vortex within the disk. Thus, even without self-gravity in
our simulations, the RWIs generate high-density “banana”-shaped vortices that
contain a substantial fraction of the disk’s mass, as shown in Fig. 1(b). The
stability of the flow in the disk against the development of turbulence due to
shear and density gradients can be assessed through the Richardson number:

$R_{i}=g\left(\frac{1}{\gamma}\frac{d{\rm ln}p}{dr}-\frac{d{\rm
ln}\rho}{dr}\right)\left(r\frac{d\Omega}{dr}\right)^{-2}\,,$ (14)

where $\gamma$ is the adiabatic index, and $g$ denotes the effective gravity:

$g=\frac{GM}{r^{2}}-r\Omega^{2}\,.$ (15)

If $R_{i}\lesssim 0.25$, the flow becomes unstable. We find that in our
simulations $R_{i}<0$ in the high-density regions in the innermost $\sim
20r_{g}$, supporting the conditions for the development of local turbulence.
The turbulent mixing and differential rotation may facilitate the formation of
sharp gradients, which are essential for the emergence of Rossby vortices. We
leave a detailed analysis of the RWI criterion in collapsar disks for future
work.

Figure 2: (a) Velocity streamlines in the co-rotating frame of the high-
density wave at $(x,y)\approx(-150,100)\,\,{{\rm km}}$ portray a trapped
Rossby vortex. (b) The mass accretion rate on the BH horizon remains roughly
constant in time and across models. (c) The jet efficiency on the BH horizon
is substantial for rapidly spinning BHs with strong magnetic fields.

Figure 2(a), which presents the velocity streamlines in the disk, offers
visual evidence of these high-density clumps behaving as Rossby vortices
within their circular orbit. Namely, in the corotating frame, the vortex flow
consists of banana-shaped closed streamlines and it is trapped between
prograde and retrograde circular flows on either side of the vortex. These
modes emerge stochastically, with the $\tilde{m}$-mode number constantly
changing throughout the simulation (see movies at
http://www.oregottlieb.com/disk_gw.html).

Stellar evolution models suggest that the circularization radius remains
roughly constant during the first several seconds at $R\sim 20\,r_{g}$
(Gottlieb et al. 2024, in prep.), as employed in our simulations.
Nevertheless, the presence of strong magnetic fields facilitates MRI, which
transports angular momentum outwards, thereby driving disk expansion and
reducing disk density. This implies that later falling gas will merge with the
outer parts of the disk before reaching its circularization radius. After
several seconds, stellar evolution models (e.g., Woosley & Heger, 2006)
suggest that the circulation radius of the gas increases as a function of the
radius at the time of collapse, and thus the disk expansion will become
faster. We leave a detailed exploration of the impact of the stellar
rotational profile on the disk structure for future work.

Fig. 1(c) depicts the extended accretion disk in model $Cc$, characterized by
weak cooling with $H/R=0.3$. The increased disk size results in lower mass
density in the disk, similar to the hydrodynamic conditions at early times of
strongly-cooled disks. The low-density conditions in the disk set the stage
for developing the spiral arm structure with mode $\tilde{m}=2$.

Fig. 2(b) indicates that all models exhibit a comparable constant mass
accretion rate on the BH, as expected from free-falling gas with an initial
radial power-law index of $-1.5$ (Gottlieb et al., 2023a). The subtle
differences between the models stem from the outflow feedback on accretion,
with stronger outflows impeding accretion more. Fig. 2(c) depicts the jet
efficiency on the BH horizon, $\eta_{j}$, defined as:

$\eta_{j}=\dot{M}^{-1}\int_{r_{h},\sigma>1}\sqrt{-g}(-T^{r}_{t}-\rho
u^{r})\mathrm{d}\theta\mathrm{d}\varphi\,,$ (16)

where the integration is calculated on the BH horizon $r_{h}$, and only for
fluid elements with magnetization $\sigma>1$, $g$ is the metric determinant,
$T^{r}_{t}$ denotes the radial energy flux density component of the mixed
stress-energy tensor $T$, and $u^{\mu}$ is the four-velocity such that $\rho
u^{r}$ represents the radial mass-energy flux density. Model $C$ exhibits high
jet efficiency through the Blandford-Znajek mechanism (BZ; Blandford & Znajek,
1977). In contrast, slowly rotating and weakly-magnetized BHs can eject merely
$\eta_{j}\lesssim 1\%$ of the accreted power due to a smaller angular momentum
reservoir or weaker magnetic flux on the horizon. We note that the emergence
of collapsar jets in weakly-magnetized disks is observed only when cooling is
strong (see Gottlieb et al., 2022a, for the absence of jets in weakly-
magnetized thick disks).

##  4 GW emission

###  4.1 Analytic estimates

We begin by providing an analytic approximation for the characteristics of the
GW signal. Given that $2\pi\,{f_{\rm GW}}R_{d}/c=m\sqrt{r_{g}/R_{d}}\lesssim
1$ at the peak frequency, we can employ the gravitational quadrupole moment
($\tilde{m}=2$) to estimate the GW strain magnitude as:

$h=\frac{2G}{Dc^{4}}\frac{d^{2}Q_{m}}{dt^{2}}\approx\frac{2G}{Dc^{4}}{E_{d}\epsilon}\,,$
(17)

where $\epsilon$ signifies the degree of asymmetry introduced by the azimuthal
modes and $E_{d}$ is the disk energy, estimated as:

$E_{d}\approx\frac{1}{2}\frac{G\,{M_{\rm BH}}M_{d}}{R_{d}}\approx
10^{52}\frac{\,{M_{\rm
BH}}}{10\,\,{M_{\odot}}}\frac{M_{d}}{0.1\,\,{M_{\odot}}}\frac{100\,\,{{\rm
km}}}{R_{d}}\,\,{{\rm erg}}\,,$ (18)

where $M_{d}$ is the mass of the disk. We find the disk mass and radius in
model $C$ to grow linearly with time such that $E_{d}\approx{\rm const}\approx
7\times 10^{51}\,\,{{\rm erg}}$. Plugging Eq. (18) into Eq. (17), the
characteristic strain of the disk can be estimated as:

$h_{c}\approx 7\times 10^{-23}\epsilon\frac{10\,{\rm Mpc}}{D}\frac{\,{M_{\rm
BH}}}{10\,\,{M_{\odot}}}\frac{M_{d}}{0.1\,\,{M_{\odot}}}\frac{100\,\,{{\rm
km}}}{R_{d}}\,.$ (19)

Since $\,{M_{\rm BH}}\sim R_{d}$, $h_{c}$ is independent of the BH mass and
scales as $M_{d}$, which depends on the disk thickness, the circularization
radius of the gas and the envelope mass.

Comparing Eq. (19) with our numerical calculations, we find that
$\epsilon\approx 0.3$ in model $C$ and consistently $\epsilon\gtrsim 0.2$ in
all strongly-cooled disks, but $\epsilon\approx 0.1$ when cooling is weak in
model $Cc$. Suppose the high-density wave that dominates the GW emission is in
a quasi-Keplerian orbit. In that case, the characteristic GW peak frequency is
governed by the angular frequency $\Omega$ [Eq. (4)] and the $\tilde{m}$-mode,
given by Eq. (5). Higher subdominant modes (multipoles) will broaden the
spectrum above the peak frequency. As larger orbits emerge with the disk
expansion, the emission will shift to lower frequencies. If the high-density
wave exhibits radial dependence, such as in spiral waves, it will span varying
distances from the BH. This implies that the emission originates from
different orbits, resulting in a broader spectrum and a weaker GW signal (see
also Levinson et al., 2015; van Putten et al., 2019).

###  4.2 Simulation results

Figure 3:  Cross (a,c) and plus (b,d) polarizations for face-on observers
(a,b) and edge-on observers (c,d) feature two characteristic timescales. Short
timescales, characterized by high variability, are dictated by the orbit of
non-axisymmetric density waves close to the BH. Modulations on longer
timescales ($\sim 0.1\,\,{{\rm s}}$) in model $Cc$ imply that modes at larger
orbits contribute to the GW spectrum in extended disks.

Figure 3 depicts the two GW polarizations observed by face-on (a, b) and edge-
on (c, d) observers, as described in Eq. (9). All polarizations feature high
variability stemming from the orbital timescale of the high-density waves. As
discussed earlier, the amplitude of the short timescale variability is lower
in model $Cc$, where $H/R=0.3$, compared to the strongly-cooled disks with
$H/R\gtrsim 0.1$. Additionally, in model $Cc$, all observers would detect a
lower frequency peak over $\sim 0.1\,\,{{\rm s}}$ seen by the sinusoidal waves
(purple lines). This phenomenon arises from larger orbits in the extended
weakly-cooled disk.

We observe a phase shift of $\pi/4$ between cross and plus polarizations for
face-on observers (i.e., circular polarization), and an anti-phase between
cross and plus polarizations for edge-on observers (linear polarization). As
we shall see next, the GW spectra are rather similar for different observers.
This suggests that detecting the elliptical polarization could serve as a
probe for the inclination angle of the disk relative to the Earth’s line of
sight. We plan to explore this further in future work.

Figure 4:  Top: edge-on (red) and face-on (blue) spectra for models $C$ (first
column from left), $Ca$ (second column), $Cb$ (third column), and $Cc$ (fourth
column) are presented, with different rows depicting the spectra at various
times. The black lines represent the sensitivity of various GW detectors.
During the disk formation stage (first row), all models exhibit a quasi-flat
spectrum. When cooling is strong (first three columns), Rossby vortices
develop, leading to pronounced GW peak frequencies as the disk evolves. A
movie showcasing the evolution of the GW spectrum is available at
http://www.oregottlieb.com/disk_gw.html. Bottom: The spectrograms for face-on
observers over a moving time-window of $0.4\,\,{{\rm s}}$ illustrate the
signal amplification after $\lesssim 1\,\,{{\rm s}}$. The spectrogram of model
$Ca$ resembles that of $C$.

Figure 4 depicts the integrated GW spectra at various time windows for the
different models. During the disk formation stage within the first $\sim
1\,\,{{\rm s}}$ [Fig. 1(a)], low-density spiral waves predominantly contribute
to the GW emission. These spirals, spanning various radii, translate to a
roughly flat GW spectrum, characterized by a relatively weak signal (first row
in Fig. 4). Over time, accretion proceeds, the disk mass grows and the GW
emission intensifies. For strongly-cooled disks, high-density rings form at
$R_{d}\approx 15\,r_{g}$ [Fig. 1(b)]. Instabilities within these rings give
rise to azimuthal modes induced by high-density waves [Fig. 2(a)], which yield
coherent GW emission. Eq. (5) suggests that during this phase, the GW signal
peaks at $\sim 100\,\tilde{m}\,{\rm Hz}$. As $\tilde{m}$ changes
stochastically over time, higher multipoles contribute to the spectra some
breadth at $f\gtrsim 100\,{\rm Hz}$.

Comparing the GW spectra of strongly cooled disks, we find that all such disks
evolve similarly, regardless of their magnetic field strength and BH spin. For
example, models $C$ and $Ca$ exhibit similar evolution during the first $\sim
6\,\,{{\rm s}}$, as shown by their spectra at $5\,\,{{\rm s}}<t<6\,\,{{\rm
s}}$. Over time, the disk viscously expands and the spectrum broadens to lower
frequencies, as seen in the spectrum of model $C$ at $23\,\,{{\rm
s}}<t<24\,\,{{\rm s}}$, resulting in a broader spectrum in model $C$ due to
its longer $T_{s}$.

We estimate the GW energy using:

$\,{E_{\rm
GW}}=\sum_{\mu\nu}\frac{c^{3}}{20G}\int_{0}^{t}\left(\frac{dh_{\mu\nu}D}{dt}\right)^{2}dt~{}.$
(20)

Fig. 5 illustrates that all models exhibit a similar GW luminosity evolution
in time, indicating that the differences in $\,{E_{\rm GW}}$ in Tab. 1 and in
the integrated spectra over $T_{s}$ between the models are attributed to the
differences in $T_{s}$. The increase in the GW luminosity can be explained by
the accumulation of additional mass in the disk, which both increases the GW
amplitude [Eq. (19)] and facilitates the emergence of non-axisymmetric modes.
The asymptotic value of $L_{\rm GW}\gtrsim 10^{49}\,\,{{\rm erg}}\,\,{{\rm
s}}^{-1}$ suggests that a few percent of the total disk energy is radiated
away as GWs. In nature, the lifetime of collapsar disks might be $\sim
100\,\,{{\rm s}}$, on the order of lGRB timescales, implying the integrated GW
energy may reach $\,{E_{\rm GW}}\sim 10^{51}\,\,{{\rm erg}}$, as indicated by
model $C$ (see Tab. 1).

Figure 5: The GW luminosity as a function of time, $L_{\rm GW}=dE_{\rm
GW}/dt$, illustrates that the GW power grows similarly in all models.

In model $Cc$, characterized by weak cooling, some non-axisymmetric modes
emerge, similar to other models. However, this model does not feature a
distinct high-density ring but spiral arms (§3), resulting in a weaker
integrated signal, also increasing linearly with time. In model $Cc$, the disk
extends to $\sim 10^{3}\,\,{{\rm km}}$. Therefore, low-frequency modes are
also present due to the extended disk orbital frequency. The weakly-cooled
disks have amplitudes consistent with what has been found by Wessel et al.
(2023) for non-cooled disks.

The bottom panels of Fig. 4 depict the spectrograms in models $C,Cb,Cc$. In
strongly-cooled disks, the azimuthal density wave at $R_{d}\approx
200\,\,{{\rm km}}$ emerges at $t\gtrsim 2\,\,{{\rm s}}$, giving rise to peak
GW frequencies at $\,{f_{\rm GW}}\approx 100-300\,{\rm Hz}$. The stochastic
emergence of various non-axisymmetric modes and the disk spreading result in a
broad spectrum. In the weakly-cooled disk, no high-density rings are present,
and the disk is more extended. This leads to a broader frequency range with
lower characteristic frequencies, which are not fully captured in the
spectrogram moving time windows of $0.4\,\,{{\rm s}}$.

With the strong GW emission in cooled disks, one may expect that the GW losses
could alter the system dynamics. We make a crude estimate of this effect by
treating the orbiting vortex as a member of a binary, with the central BH
being the other member, and by computing the radiation-reaction timescale for
such a binary. The latter is given by (e.g., Shapiro & Teukolsky, 1983):

$t_{rr}=\frac{5}{256}\frac{c^{5}}{G^{3}}\frac{a^{4}}{M^{2}\mu}\,,$ (21)

where $a$ is the semi-major axis, well-approximated by the orbital radius of
the vortex, $M\approx M_{\rm BH}$ is the total mass of the binary that is
well-approximated by the BH mass, and $\mu$ is the reduced mass of the binary
that is well-approximated by the vortex mass $m$. We find:

$t_{rr}=\frac{5}{256}\frac{c^{5}}{G^{3}}\frac{R^{4}}{\,{M_{\rm
BH}}^{2}m}\approx 10\left(\frac{10\,\,{M_{\odot}}}{\,{M_{\rm
BH}}}\right)^{2}\frac{10^{-2}\,\,{M_{\odot}}}{m}\left(\frac{R}{150\,\,{{\rm
km}}}\right)^{4}\,\,{{\rm s}}\,.$ (22)

Gravitational radiation reaction is important if it acts on a shorter
timescale than the viscous timescale, which is given by:

$t_{v}=\frac{R^{2}}{\alpha c_{s}H}\approx 0.75\frac{R_{d}}{150\,\,{{\rm
km}}}\left(\frac{\alpha}{0.1}\frac{c_{s}}{10^{9}\,{\rm cm\,\,{{\rm
s}}^{-1}}}\frac{H/R}{0.2}\right)^{-1}\,\,{{\rm s}}\,,$ (23)

where $\alpha$ is the effective viscosity coefficient and $c_{s}$ is the sound
speed.

The normalization values in Eqs. (22), (23) are based on typical values found
in the simulations. The order-of-magnitude ratio suggests that the
gravitational radiation reaction does not dramatically alter the disk dynamics
and the Rossby vortices. The $\alpha$ viscosity parameter is not calculated in
the simulation and its value can vary by more than an order of magnitude.
However, when the magnetic fields are strong, as considered in collapsars, its
value is typically $\alpha\gtrsim 0.1$ (Kotko & Lasota, 2012; Ju et al.,
2017), thereby unlikely to affect this result.

##  5 Detectability

Reliable computations of the SNR, a customary measure of detectability in GW
astronomy, encounter significant uncertainties stemming from various factors:
(1) The accuracy of our findings is contingent upon our simplistic
prescription of disk cooling; (2) The absence of capturing global
instabilities driven by self-gravity in our simulations presents a challenge
to the robustness of the observed disk hydrodynamics; (3) Although our
simulations are the longest 3D GRMHD collapsar simulations ever performed,
they are constrained to only a fraction of the expected disk’s lifetime.
Consequently, late-time disk characteristics may diverge from early-time
conditions, e.g., due to non-linear feedback from outflows; (4) The disk’s
lifetime and spatial dimensions are intricately linked to the progenitor’s
rotational and density profiles, which are poorly constrained observationally
and fraught with uncertainties in theoretical modeling, particularly regarding
angular momentum transport within rotating stars. Similarly, the masses of the
stellar envelope and the BH affect $M_{d}$ and $R_{d}$, and thus the GW
amplitude and frequency. Unlike the progenitor profiles, these quantities are
inferred from observations but span a wide range of values. (5) The fraction
of stellar core collapses that lead to the formation of accretion disks
remains poorly constrained observationally.

Advanced simulations should be designed to assess the uncertainty arising from
(1) and (2). We incorporate the remaining uncertainties using the following
approximations:

(3,4) The lifespan of an accretion disk could be deduced from the rotational
profile in stellar evolution models of rapidly rotating stars. Rapid rotation
leads to the stripping of the outermost stellar envelope, resulting in a star
radius of $R\approx 10^{11}\,{\rm cm}$ (Woosley & Heger, 2006), corresponding
to an accretion timescale of hundreds of seconds. However, the launch of
outflows impedes accretion, shortening the accretion disk timescale. Our best
estimate is that the simulated accretion disk lifetime is approximately an
order of magnitude shorter than what simple infall time calculations suggest.

The expected disk lifetime is longer than our current simulation runs, and it
would be appealing to introduce a simple scaling of SNR with $\,{t_{\rm
acc}}$. However, the increase in SNR with $\,{t_{\rm acc}}$ is not
straightforward to infer. As shown in Fig. 5, the GW luminosity grows with
mass accumulation in the disk over time. On the other hand, as the disk
expands, the GW emission is anticipated to weaken and shift towards lower
frequencies, where the sensitivity of LVK is diminished and that of CE
increases [see Eq. (19)]. The frequency is also contingent upon $\,{M_{\rm
BH}}$ and $R_{d}$ as per Eq. (5). A BH could potentially form with a mass as
low as $\,{M_{\rm BH}}\approx 3\,\,{M_{\odot}}$, resulting in a smaller disk
radius and higher frequencies – $\,{f_{\rm GW}}\sim\,{M_{\rm BH}}^{-1}$ [see
Eq. (5)], deviating from the peak sensitivity range of both LVK and CE.
Nevertheless, this might be offset by smaller BHs exhibiting an extended
accretion timescale, which scales as $\,{t_{\rm acc}}\sim\,{M_{\rm
BH}}^{-0.5}$. Ultimately, the BH mass growth with $\dot{M}\approx
0.1\,\,{M_{\odot}}\,\,{{\rm s}}^{-1}$ [Fig. 2(b)], and the increasing orbital
radius due to the accretion of shells with larger circularization radius and
the disk viscous spreading will augment $R_{d}$, likely culminating in a
signal peaking at $\,{f_{\rm GW}}\approx 100\,{\rm Hz}$, unless BH-driven
outflows impede accretion in later stages. We consider a conservative
potential enhancement of the numerically calculated SNR values [Eq. (12)] at
$D=10\,{\rm Mpc}$ by a factor of $100\,\,{{\rm s}}/T_{s}$, given that
$\,{t_{\rm acc}}$ might be hundreds of seconds while the SNR grows sub-
linearly. Overall, collapsar disks power a GW signal with SNR in LVK (CE) of
dozens (hundreds) at $D=10\,{\rm Mpc}$.

(5) We constrain the formation rates of collapsar accretion disks by
establishing lower and upper limits. For the lower limit, we argue that as the
formation of lGRB jets necessitates the presence of both an accretion disk and
a robust magnetic field, collapsar disks are at least as common as lGRBs. The
observed rate of the latter is ${\cal{R}_{\rm lGRB}}\sim 1\,{\rm
Gpc^{-3}}~{}{\rm yr}^{-1}$ (Wanderman & Piran, 2010), and thus for a typical
jet opening angle of $\theta_{j}\approx 0.1$, the collapsar rate is
${\cal{R}_{\rm Collapsar}}\approx 100~{}{\rm Gpc^{-3}}~{}{\rm yr}^{-1}$. It is
worth noting that if the jet exhibits wobbling behavior, the inferred lower
limit on the rate could be lower by up to an order of magnitude (Gottlieb et
al., 2023b). For the upper limit, it has been proposed that most jets do not
generate typical lGRBs (Bromberg et al., 2012), but rather power a separate
class of low-luminosity GRBs (llGRBs) driven by a mildly-relativistic shock
breakout (Tan et al., 2001; Waxman et al., 2007; Nakar, 2015). Such mildly-
relativistic shocks necessitate the formation of jets, which are likely
associated with collapsar disks. Therefore, we adopt the high llGRB rate
${\cal{R}_{\rm llGRB}}\approx 10^{4}~{}{\rm Gpc^{-3}}~{}{\rm yr}^{-1}$
(Soderberg et al., 2006), which is also comparable with the rate of stripped-
envelope SNe Ib/c of ${\cal{R}_{\rm Ib/c}}\approx 2.6\times 10^{4}~{}{\rm
Gpc^{-3}}~{}{\rm yr}^{-1}$, as an upper limit on the rate of collapsar disks.

Following the above approximations and assuming that the simulated collapsar
configurations represent the broader collapsing star population, we are ready
to discuss the detection prospects of GWs from collapsar disks. Setting an
optimistic matched-filter SNR threshold for detection at
$\varrho_{c}=20\varrho_{20}$ (estimating $20\lesssim\varrho_{c}\lesssim 50$;
Macquet, private communication), we project that LVK O4 holds the potential to
detect GWs emanating from accretion disks up to distances of a few dozen Mpc.
Conversely, CE is anticipated to discern such signals at distances extending
to hundreds of Mpc. In the CE scenario, even weakly-cooled disks hold
detectability prospects spanning dozens of Mpc.

If collapsar disk rates are similar to those of lGRBs, then the expected LVK
event rate is $\sim 10^{-2}\,\varrho_{20}^{-3}\,{\rm yr}^{-1}$. However, if
disks emerge at rates similar to SNe Ib/Ic, a few GW events might already have
been detectable by LVK, as indicated by more than a dozen SNe Ib/c events
within $D\lesssim 40\,{\rm Mpc}$ in the Zwicky Transient Facility (ZTF; Graham
et al., 2019; Masci et al., 2019), including events as close as $D\approx
15\,{\rm Mpc}$ (e.g., Jacobson-Galán et al., 2020; Kilpatrick et al., 2021;
Rho et al., 2021; Ho et al., 2023; Kuncarayakti et al., 2023). In CE, several
events per year might be detectable even if the cooling efficiency is low and
collapsar disks are as common as lGRBs. If a substantial fraction of SNe Ib/c
gives rise to accretion disks, hundreds of events could become detectable with
third-generation GW detectors, depending on the disk cooling efficiency. In
either scenario, both the strength and the rate of detectable GW signals from
accretion disks are orders of magnitude more promising than those expected
from CCSNe.

Electromagnetic counterparts have the potential to significantly enhance
temporal and spatial localization, thereby aiding the search for GWs emanating
from collapsar disks. Relativistic outflows such as lGRBs can provide valuable
constraints on time localization, up to the propagation of the jet in the
stellar envelope of $\lesssim 100\,\,{{\rm s}}$. However, the detectability of
such jets is often low, estimated at only $\lesssim 1\%$ due to their beamed
emission (Frail et al., 2001). Non-relativistic electromagnetic counterparts
are anticipated in most SN Ib/c events but offer weaker time localization,
spanning from hour-long shock breakout to week-long radioactive decay emission
post-collapse. Consequently, collapsars may not benefit from precise time
localization, presenting a challenge in conducting targeted searches for GW
emission from collapsar disks. Nevertheless, regardless of whether the
electromagnetic signal aids in the detection of the GW signal, all GW events
are expected to be accompanied by a detectable electromagnetic counterpart, at
least those detectable by present-day GW detectors, rendering these events
multi-messenger phenomena.

##  6 Discussion

In this _letter_ , we underscore the potential of collapsar accretion disks as
promising isotropic GW sources, even for present-day GW detectors. To
investigate this, we conducted the first numerical simulations tracking the
evolution of collapsar disks. We explored various disk magnetizations, BH
spins, and cooling strengths to assess their respective impacts on the emitted
GWs. We found a consistent picture across simulations with cooling exerting
the most significant influence on the resulting GW emission. While post-merger
disks are also subject to strong cooling, their relatively short lifetimes and
low masses render their detection prospects unlikely (see Appendix A).

We implemented artificial cooling that operates instantaneously. Specifically,
the cooling timescale is nearly zero, implying that cooling directly induces
disk instability rather than through an increase in surface density (Gammie,
2001). When cooling is strong, the dominant instability observed is the
Rossby-Wave Instability, which generates a Rossby vortex within a high-density
ring in the innermost part of the disk at $R_{d}\approx 15\,r_{g}$ (see e.g.,
Tagger & Pellat, 1999). The number of high-density clumps dictates the
azimuthal mode number, consequently governing the GW frequency. In our
simulations, the characteristic GW frequency begins with $\,{f_{\rm
GW}}\gtrsim 100\,{\rm Hz}$ with some breadth due to contribution from various
azimuthal modes. As time progresses, the disk expands via viscous spreading
and the larger circularization radius of the stellar outer shells, gradually
shifting the spectrum towards lower frequencies. When the disk cooling is
weaker, the large radial extent of the disk reduces the mass density in the
disk. This supports the development of spiral arms, giving rise to a broader
and attenuated GW spectrum. In §5, we discussed potential uncertainties in the
progenitor rotational profiles, but our understanding is that the qualitative
behavior of the Rossby-Wave Instability should persist regardless of the
specific model, provided that the cooling is sufficiently strong.

Our findings of a robust signal motivate the exploration of GWs from collapsar
accretion disks in observed LVK data. Abbott et al. (2021) conducted an
extensive all-sky search for long-duration ($t\gtrsim 1\,\,{{\rm s}}$) GW
transients in LVK O3 run data. They established upper limits on the GW strain
for various transient sources, including accretion disks. Using our simulated
spectra, we estimate these limits correspond to weak constraints on the non-
detection of accretion disks at $D\lesssim 1\,{\rm Mpc}$. However, previous
searches did not utilize templates but focused on detecting excess power in
spectrograms (Thrane et al., 2011). Consequently, accurately assessing the
appropriate matched-filter SNR for detection remains challenging. Our
conservative yet optimistic estimate of a matched-filter SNR of
$\varrho_{c}=20$ for detection corresponds to a maximum detection distance of
a few dozen Mpc for GWs from collapsar disks.

In addition to uncertainties in $\varrho_{c}$, the unknown rate with which
collapsing stars form accretion disks presents a challenge in estimating the
expected detection rate of these objects. Assuming that a significant fraction
of stripped-envelope SNe is accompanied by disk formation and
$\varrho_{c}=20$, a rough estimate of the detection rate suggests that a few
events could already be detectable using present-day GW detectors,
highlighting the need for follow-up modeled and targeted searches. The number
of GW detections from collapsar disks is anticipated to increase by orders of
magnitude in future detectors such as Cosmic Explorer or Einstein Telescope.
Even if the GW emission from collapsar disks is lower than our estimates, e.g.
due to a lower cooling efficiency, dozens of events are expected to be
detectable. To explicitly state our worst-case scenario expectations: if
accretion disks form at a relatively low rate, closer to that observed for
lGRBs, or a large fraction of them is not strongly concentrated and cooled, we
may need to await future advancements in GW detectors to detect these GW
events.

Collapsar accretion disks have emerged as highly promising sources of GW
signals, boasting significantly greater strength by $\gtrsim 2$ orders of
magnitude compared to extensively studied CCSN events. While exhibiting a
comparable strength to cocoon-powered GWs, accretion disks are considerably
more abundant, increasing their detection prospects. Furthermore, among all
suggested GW sources from collapsars, accretion disks are anticipated to
exhibit the most coherent pattern, potentially offering an opportunity to
develop phenomenological templates for their GW bursts.

Detecting GWs emanating from accretion disks presents a novel opportunity to
study the inner workings of collapsing stars, as the disk time evolution
reflects the radial structure of the collapsing star. The duration of the GW
signal holds valuable insights into the accretion timescale dictated by the
balance between the collapse and the outflows. Furthermore, the GW spectrum
can help determine important disk properties such as mass, radial extent, and
energy. The GW elliptical polarization may reveal the orientation of the disk
with respect to the line of sight. As statistical data accumulates, this
approach aims to elucidate the intricate relationship between accretion disks,
CCSNe, jet launching, and lGRBs.

In summary, considering the relatively high expected event rate, we argue that
collapsar accretion disks might be the most promising burst-type GW sources to
date. We therefore urge both the theoretical and observational communities to
embark on exploration and modeling of GWs originating from collapsar accretion
disks. On the theoretical side, enhancing the numerical experiments with
neutrino-cooling schemes and self-gravity, as well as utilizing progenitors
from stellar evolution models, will place theoretical predictions on firmer
ground. On the observational side, it seems important to comprehensively
characterize GW signals from various collapsar models and develop detailed
templates for these signals. The GWs might be imminently detectable in
current-day detectors and could be targeted by their electromagnetic
counterparts. Such efforts will facilitate a more efficient GW search and
enrich the depth of information attainable from a multi-messenger signal
concerning the physical properties of the source.

We thank Will Farr, Anthony Mezzacappa, David Radice, Katerina Chatziioannou,
Eric Thrane, Neil Cornish, Adrian Macquet, Florent Robinet, Adam Burrows, and
Sharan Banagiri for useful discussions. OG is supported by the Flatiron
Research Fellowship. This work was supported by a grant from the Simons
Foundation (MPS-EECS-00001470-04) to AL and YL. YL’s work on this subject is
also supported by Simons Investigator Grant 827103.

## References

  * Abbott et al. (2017) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, Classical and Quantum Gravity, 34, 044001, doi: 10.1088/1361-6382/aa51f4
  * Abbott et al. (2018) —. 2018, Living Reviews in Relativity, 21, 3, doi: 10.1007/s41114-018-0012-9
  * Abbott et al. (2021) Abbott, R., Abbott, T. D., Acernese, F., et al. 2021, Phys. Rev. D, 104, 102001, doi: 10.1103/PhysRevD.104.102001
  * Abbott et al. (2023) —. 2023, Physical Review X, 13, 011048, doi: 10.1103/PhysRevX.13.011048
  * Abdikamalov et al. (2020) Abdikamalov, E., Pagliaroli, G., & Radice, D. 2020, arXiv e-prints, arXiv:2010.04356, doi: 10.48550/arXiv.2010.04356
  * Accadia et al. (2012) Accadia, T., Acernese, F., Alshourbagy, M., et al. 2012, Journal of Instrumentation, 7, 3012, doi: 10.1088/1748-0221/7/03/P03012
  * Acernese et al. (2015) Acernese, F., Agathos, M., Agatsuma, K., et al. 2015, Classical and Quantum Gravity, 32, 024001, doi: 10.1088/0264-9381/32/2/024001
  * Acernese et al. (2019) Acernese, F., Agathos, M., Aiello, L., et al. 2019, Phys. Rev. Lett., 123, 231108, doi: 10.1103/PhysRevLett.123.231108
  * Akutsu et al. (2021) Akutsu, T., Ando, M., Arai, K., et al. 2021, Progress of Theoretical and Experimental Physics, 2021, 05A101, doi: 10.1093/ptep/ptaa125
  * Amaro-Seoane et al. (2017) Amaro-Seoane, P., Audley, H., Babak, S., et al. 2017, arXiv e-prints, arXiv:1702.00786, doi: 10.48550/arXiv.1702.00786
  * Andresen et al. (2021) Andresen, H., Glas, R., & Janka, H. T. 2021, MNRAS, 503, 3552, doi: 10.1093/mnras/stab675
  * Andresen et al. (2017) Andresen, H., Müller, B., Müller, E., & Janka, H. T. 2017, MNRAS, 468, 2032, doi: 10.1093/mnras/stx618
  * Andresen et al. (2019) Andresen, H., Müller, E., Janka, H. T., et al. 2019, MNRAS, 486, 2238, doi: 10.1093/mnras/stz990
  * Babak et al. (2021) Babak, S., Hewitson, M., & Petiteau, A. 2021, arXiv e-prints, arXiv:2108.01167, doi: 10.48550/arXiv.2108.01167
  * Baiotti et al. (2008) Baiotti, L., Giacomazzo, B., & Rezzolla, L. 2008, Phys. Rev. D, 78, 084033, doi: 10.1103/PhysRevD.78.084033
  * Batta & Lee (2014) Batta, A., & Lee, W. H. 2014, MNRAS, 437, 2412, doi: 10.1093/mnras/stt2061
  * Birnholtz & Piran (2013) Birnholtz, O., & Piran, T. 2013, Phys. Rev. D, 87, 123007, doi: 10.1103/PhysRevD.87.123007
  * Blandford & Znajek (1977) Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433, doi: 10.1093/mnras/179.3.433
  * Bromberg et al. (2012) Bromberg, O., Nakar, E., Piran, T., & Sari, R. 2012, ApJ, 749, 110, doi: 10.1088/0004-637X/749/2/110
  * Bugli et al. (2018) Bugli, M., Guilet, J., Müller, E., et al. 2018, MNRAS, 475, 108, doi: 10.1093/mnras/stx3158
  * Buikema et al. (2020) Buikema, A., Cahillane, C., Mansell, G. L., et al. 2020, Phys. Rev. D, 102, 062003, doi: 10.1103/PhysRevD.102.062003
  * Chen & Beloborodov (2007) Chen, W.-X., & Beloborodov, A. M. 2007, ApJ, 657, 383, doi: 10.1086/508923
  * Corsi et al. (2024) Corsi, A., Barsotti, L., Berti, E., et al. 2024, arXiv e-prints, arXiv:2402.13445.  https://arxiv.org/abs/2402.13445
  * Ecker et al. (2024) Ecker, C., Gorda, T., Kurkela, A., & Rezzolla, L. 2024, arXiv e-prints, arXiv:2403.03246, doi: 10.48550/arXiv.2403.03246
  * Evans et al. (2023) Evans, M., Corsi, A., Afle, C., et al. 2023, arXiv e-prints, arXiv:2306.13745, doi: 10.48550/arXiv.2306.13745
  * Fernández et al. (2020) Fernández, R., Foucart, F., & Lippuner, J. 2020, MNRAS, 497, 3221, doi: 10.1093/mnras/staa2209
  * Flanagan & Hughes (1998) Flanagan, É. É., & Hughes, S. A. 1998, Phys. Rev. D, 57, 4535, doi: 10.1103/PhysRevD.57.4535
  * Foucart (2023) Foucart, F. 2023, Living Reviews in Computational Astrophysics, 9, 1, doi: 10.1007/s41115-023-00016-y
  * Frail et al. (2001) Frail, D. A., Kulkarni, S. R., Sari, R., et al. 2001, ApJ, 562, L55, doi: 10.1086/338119
  * Gammie (2001) Gammie, C. F. 2001, ApJ, 553, 174, doi: 10.1086/320631
  * Gottlieb et al. (2023a) Gottlieb, O., Jacquemin-Ide, J., Lowell, B., Tchekhovskoy, A., & Ramirez-Ruiz, E. 2023a, ApJ, 952, L32, doi: 10.3847/2041-8213/ace779
  * Gottlieb et al. (2022a) Gottlieb, O., Lalakos, A., Bromberg, O., Liska, M., & Tchekhovskoy, A. 2022a, MNRAS, 510, 4962, doi: 10.1093/mnras/stab3784
  * Gottlieb et al. (2022b) Gottlieb, O., Liska, M., Tchekhovskoy, A., et al. 2022b, ApJ, 933, L9, doi: 10.3847/2041-8213/ac7530
  * Gottlieb et al. (2023b) Gottlieb, O., Nagakura, H., Tchekhovskoy, A., et al. 2023b, ApJ, 951, L30, doi: 10.3847/2041-8213/ace03a
  * Gottlieb et al. (2023c) Gottlieb, O., Metzger, B. D., Quataert, E., et al. 2023c, ApJ, 958, L33, doi: 10.3847/2041-8213/ad096e
  * Graham et al. (2019) Graham, M. J., Kulkarni, S. R., Bellm, E. C., et al. 2019, PASP, 131, 078001, doi: 10.1088/1538-3873/ab006c
  * Hawley et al. (2011) Hawley, J. F., Guan, X., & Krolik, J. H. 2011, ApJ, 738, 84, doi: 10.1088/0004-637X/738/1/84
  * Hild et al. (2011) Hild, S., Abernathy, M., Acernese, F., et al. 2011, Classical and Quantum Gravity, 28, 094013, doi: 10.1088/0264-9381/28/9/094013
  * Ho et al. (2023) Ho, A. Y. Q., Sollerman, J., & Perley, D. 2023, Transient Name Server AstroNote, 57, 1
  * Jacobson-Galán et al. (2020) Jacobson-Galán, W. V., Margutti, R., Kilpatrick, C. D., et al. 2020, ApJ, 898, 166, doi: 10.3847/1538-4357/ab9e66
  * Ju et al. (2017) Ju, W., Stone, J. M., & Zhu, Z. 2017, ApJ, 841, 29, doi: 10.3847/1538-4357/aa705d
  * Just et al. (2022) Just, O., Goriely, S., Janka, H. T., Nagataki, S., & Bauswein, A. 2022, MNRAS, 509, 1377, doi: 10.1093/mnras/stab2861
  * Kilpatrick et al. (2021) Kilpatrick, C. D., Drout, M. R., Auchettl, K., et al. 2021, MNRAS, 504, 2073, doi: 10.1093/mnras/stab838
  * Kiuchi et al. (2011) Kiuchi, K., Shibata, M., Montero, P. J., & Font, J. A. 2011, Phys. Rev. Lett., 106, 251102, doi: 10.1103/PhysRevLett.106.251102
  * Kobayashi & Mészáros (2003) Kobayashi, S., & Mészáros, P. 2003, ApJ, 589, 861, doi: 10.1086/374733
  * Kohri et al. (2005) Kohri, K., Narayan, R., & Piran, T. 2005, ApJ, 629, 341, doi: 10.1086/431354
  * Korobkin et al. (2011) Korobkin, O., Abdikamalov, E. B., Schnetter, E., Stergioulas, N., & Zink, B. 2011, Phys. Rev. D, 83, 043007, doi: 10.1103/PhysRevD.83.043007
  * Kotko & Lasota (2012) Kotko, I., & Lasota, J. P. 2012, A&A, 545, A115, doi: 10.1051/0004-6361/201219618
  * Kuncarayakti et al. (2023) Kuncarayakti, H., Sollerman, J., Izzo, L., et al. 2023, A&A, 678, A209, doi: 10.1051/0004-6361/202346526
  * Lattimer & Swesty (1991) Lattimer, J. M., & Swesty, F. D. 1991, Nucl. Phys., A535, 331, doi: 10.1016/0375-9474(91)90452-C
  * Lehner et al. (2016) Lehner, L., Liebling, S. L., Palenzuela, C., & Motl, P. M. 2016, Phys. Rev. D, 94, 043003, doi: 10.1103/PhysRevD.94.043003
  * Leiderschneider & Piran (2021) Leiderschneider, E., & Piran, T. 2021, Phys. Rev. D, 104, 104002, doi: 10.1103/PhysRevD.104.104002
  * Levinson & Eichler (1993) Levinson, A., & Eichler, D. 1993, ApJ, 418, 386, doi: 10.1086/173397
  * Levinson et al. (2015) Levinson, A., van Putten, M. H. P. M., & Pick, G. 2015, ApJ, 812, 124, doi: 10.1088/0004-637X/812/2/124
  * Li et al. (2001) Li, H., Colgate, S. A., Wendroff, B., & Liska, R. 2001, ApJ, 551, 874, doi: 10.1086/320241
  * Li et al. (2000) Li, H., Finn, J. M., Lovelace, R. V. E., & Colgate, S. A. 2000, ApJ, 533, 1023, doi: 10.1086/308693
  * LIGO Scientific Collaboration et al. (2015) LIGO Scientific Collaboration, Aasi, J., Abbott, B. P., et al. 2015, Classical and Quantum Gravity, 32, 074001, doi: 10.1088/0264-9381/32/7/074001
  * Liska et al. (2022) Liska, M. T. P., Chatterjee, K., Issa, D., et al. 2022, ApJS, 263, 26, doi: 10.3847/1538-4365/ac9966
  * Longo Micchi et al. (2023) Longo Micchi, L. F., Radice, D., & Chirenti, C. 2023, MNRAS, 525, 6359, doi: 10.1093/mnras/stad2420
  * Lovelace et al. (1999) Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805, doi: 10.1086/306900
  * Masci et al. (2019) Masci, F. J., Laher, R. R., Rusholme, B., et al. 2019, PASP, 131, 018003, doi: 10.1088/1538-3873/aae8ac
  * Mezzacappa & Zanolin (2024) Mezzacappa, A., & Zanolin, M. 2024, arXiv e-prints, arXiv:2401.11635, doi: 10.48550/arXiv.2401.11635
  * Mezzacappa et al. (2020) Mezzacappa, A., Marronetti, P., Landfield, R. E., et al. 2020, Phys. Rev. D, 102, 023027, doi: 10.1103/PhysRevD.102.023027
  * Mezzacappa et al. (2023) —. 2023, Phys. Rev. D, 107, 043008, doi: 10.1103/PhysRevD.107.043008
  * Misner et al. (1973) Misner, C. W., Thorne, K. S., & Wheeler, J. A. 1973, Gravitation
  * Moore et al. (2015) Moore, C. J., Cole, R. H., & Berry, C. P. L. 2015, Classical and Quantum Gravity, 32, 015014, doi: 10.1088/0264-9381/32/1/015014
  * Morozova et al. (2018) Morozova, V., Radice, D., Burrows, A., & Vartanyan, D. 2018, ApJ, 861, 10, doi: 10.3847/1538-4357/aac5f1
  * Nakamura et al. (2016) Nakamura, K., Horiuchi, S., Tanaka, M., et al. 2016, MNRAS, 461, 3296, doi: 10.1093/mnras/stw1453
  * Nakar (2015) Nakar, E. 2015, ApJ, 807, 172, doi: 10.1088/0004-637X/807/2/172
  * Narayan et al. (2001) Narayan, R., Piran, T., & Kumar, P. 2001, ApJ, 557, 949, doi: 10.1086/322267
  * Nitz et al. (2021) Nitz, A., Harry, I., Brown, D., et al. 2021, gwastro/pycbc: 1.18.0 release of PyCBC, v1.18.0, Zenodo, Zenodo, doi: 10.5281/zenodo.4556907
  * Noble et al. (2009) Noble, S. C., Krolik, J. H., & Hawley, J. F. 2009, ApJ, 692, 411, doi: 10.1088/0004-637X/692/1/411
  * Piro & Pfahl (2007) Piro, A. L., & Pfahl, E. 2007, ApJ, 658, 1173, doi: 10.1086/511672
  * Powell & Müller (2024) Powell, J., & Müller, B. 2024, arXiv e-prints, arXiv:2406.09691.  https://arxiv.org/abs/2406.09691
  * Punturo et al. (2010a) Punturo, M., Abernathy, M., Acernese, F., et al. 2010a, Classical and Quantum Gravity, 27, 084007, doi: 10.1088/0264-9381/27/8/084007
  * Punturo et al. (2010b) —. 2010b, Classical and Quantum Gravity, 27, 194002, doi: 10.1088/0264-9381/27/19/194002
  * Radice et al. (2019) Radice, D., Morozova, V., Burrows, A., Vartanyan, D., & Nagakura, H. 2019, ApJ, 876, L9, doi: 10.3847/2041-8213/ab191a
  * Reitze et al. (2019) Reitze, D., Adhikari, R. X., Ballmer, S., et al. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 35, doi: 10.48550/arXiv.1907.04833
  * Rho et al. (2021) Rho, J., Evans, A., Geballe, T. R., et al. 2021, ApJ, 908, 232, doi: 10.3847/1538-4357/abd850
  * Richardson et al. (2024) Richardson, C. J., Andresen, H., Mezzacappa, A., et al. 2024, arXiv e-prints, arXiv:2404.02131, doi: 10.48550/arXiv.2404.02131
  * Sago et al. (2004) Sago, N., Ioka, K., Nakamura, T., & Yamazaki, R. 2004, Phys. Rev. D, 70, 104012, doi: 10.1103/PhysRevD.70.104012
  * Shapiro & Teukolsky (1983) Shapiro, S. L., & Teukolsky, S. A. 1983, Black holes, white dwarfs and neutron stars. The physics of compact objects, doi: 10.1002/9783527617661
  * Shibata & Taniguchi (2006) Shibata, M., & Taniguchi, K. 2006, Phys. Rev. D, 73, 064027, doi: 10.1103/PhysRevD.73.064027
  * Siegel et al. (2019) Siegel, D. M., Barnes, J., & Metzger, B. D. 2019, Nature, 569, 241, doi: 10.1038/s41586-019-1136-0
  * Soderberg et al. (2006) Soderberg, A. M., Kulkarni, S. R., Nakar, E., et al. 2006, Nature, 442, 1014, doi: 10.1038/nature05087
  * Srivastava et al. (2019) Srivastava, V., Ballmer, S., Brown, D. A., et al. 2019, Phys. Rev. D, 100, 043026, doi: 10.1103/PhysRevD.100.043026
  * St. Germaine-Fuller et al. (2013) St. Germaine-Fuller, J., Drasco1, S., Ott, C., & Mösta, P. 2013, https://dcc.ligo.org/LIGO-T1400235/public
  * Tagger & Pellat (1999) Tagger, M., & Pellat, R. 1999, A&A, 349, 1003, doi: 10.48550/arXiv.astro-ph/9907267
  * Tan et al. (2001) Tan, J. C., Matzner, C. D., & McKee, C. F. 2001, ApJ, 551, 946, doi: 10.1086/320245
  * Thorne (1980) Thorne, K. S. 1980, Reviews of Modern Physics, 52, 299.  https://api.semanticscholar.org/CorpusID:122954742
  * Thrane et al. (2011) Thrane, E., Kandhasamy, S., Ott, C. D., et al. 2011, Phys. Rev. D, 83, 083004, doi: 10.1103/PhysRevD.83.083004
  * Tse et al. (2019) Tse, M., Yu, H., Kijbunchoo, N., et al. 2019, Phys. Rev. Lett., 123, 231107, doi: 10.1103/PhysRevLett.123.231107
  * Usman et al. (2016) Usman, S. A., Nitz, A. H., Harry, I. W., et al. 2016, Classical and Quantum Gravity, 33, 215004, doi: 10.1088/0264-9381/33/21/215004
  * van Putten & Levinson (2003) van Putten, M. H. P. M., & Levinson, A. 2003, ApJ, 584, 937, doi: 10.1086/345900
  * van Putten et al. (2019) van Putten, M. H. P. M., Levinson, A., Frontera, F., et al. 2019, European Physical Journal Plus, 134, 537, doi: 10.1140/epjp/i2019-12932-3
  * Vartanyan et al. (2023) Vartanyan, D., Burrows, A., Wang, T., Coleman, M. S. B., & White, C. J. 2023, Phys. Rev. D, 107, 103015, doi: 10.1103/PhysRevD.107.103015
  * Wanderman & Piran (2010) Wanderman, D., & Piran, T. 2010, MNRAS, 406, 1944, doi: 10.1111/j.1365-2966.2010.16787.x
  * Waxman et al. (2007) Waxman, E., Mészáros, P., & Campana, S. 2007, ApJ, 667, 351, doi: 10.1086/520715
  * Wessel et al. (2021) Wessel, E., Paschalidis, V., Tsokaros, A., Ruiz, M., & Shapiro, S. L. 2021, Phys. Rev. D, 103, 043013, doi: 10.1103/PhysRevD.103.043013
  * Wessel et al. (2023) —. 2023, Phys. Rev. D, 107, 123031, doi: 10.1103/PhysRevD.107.123031
  * Woosley (1993) Woosley, S. E. 1993, ApJ, 405, 273, doi: 10.1086/172359
  * Woosley & Heger (2006) Woosley, S. E., & Heger, A. 2006, ApJ, 637, 914, doi: 10.1086/498500
  * Yakunin et al. (2017) Yakunin, K. N., Mezzacappa, A., Marronetti, P., et al. 2017, arXiv e-prints, arXiv:1701.07325, doi: 10.48550/arXiv.1701.07325

##  Appendix A Compact binary mergers

During the merger of a BNS or a BH-NS system, some of the NS crust is
disrupted and circularized around the merger remnant. The tidal disruption of
the NS naturally gives rise to non-axisymmetric modes. This phenomenon is
illustrated in Figure 6(a), which presents a 3D rendering of the logarithmic
mass density around the post-merger BH $10\,{\rm ms}$ after a BNS merger. The
spiral structure persists for a few dozen milliseconds before transitioning
into an axisymmetric configuration within a more uniformly distributed post-
merger disk [Fig. 6(b)]. The compact nature of these post-merger disks
(Fernández et al., 2020) indicates that their characteristic GW frequency is
notably higher compared to those observed in collapsars.

Fig. 6(c-f) present the integrated GW spectra at different times after the
merger. During the first $\sim 30\,{\rm ms}$, the spiral arm generates strong
GW emission due to a significant deviation from axisymmetry. The contribution
from different radii results in a roughly flat spectrum. As the disk starts to
be accreted and become more axisymmetric, its GW emission decays with time. At
$t\gtrsim 100\,\,{{\rm ms}}$, the spiral arm vanishes, and density waves at
$R_{d}\approx 25\,\,{{\rm km}}$ govern the GW emission. The small radius of
post-merger accretion disks results in a characteristic GW frequency of
$\,{f_{\rm GW}}\approx 1.5\,{\rm kHz}$ [right panels, see Eq. (5)].

The SNR indicated in Tab. 1 is likely too low for detection, even if
considering the upper limit of BNS merger rate of222Post-merger massive disks
with $M_{d}\approx 0.1\,\,{M_{\odot}}$ power lGRBs (Gottlieb et al., 2023c).
If standard short GRBs are more abundant in the Universe than lGRBs from
mergers, then our simulated setup is uncommon among mergers, rendering this
rate over-optimistic. ${\cal{R}_{\rm BNS}}\sim 10-1700\,{\rm Gpc\,yr^{-3}}$
(Abbott et al., 2023) by future-generation GW detectors. Additionally, the GWs
from the disk might be overshadowed by the GWs from the BH ringdown during the
first several ms after the merger. This is particularly true considering
unequal mass binaries tend to feature more massive disks for promising GW
detectability, but also longer ringdown timescales. If the post-merger remnant
is a meta-stable NS that peaks at $\,{f_{\rm GW}}\gtrsim 1\tilde{m}\,{\rm
kHz}$ (Shibata & Taniguchi, 2006; Baiotti et al., 2008; Lehner et al., 2016),
then the ringdown can even last dozens of ms (Ecker et al., 2024).

Figure 6: 3D renderings illustrate the logarithmic mass density (in .c.g.s.)
of post-merger accretion disks. The length scale can be measured with respect
to the size of the BH in the center. (a) The early formation stages of the
disk ($\sim 10\,{\rm ms}$) feature high-density spiral arms resulting from the
disrupted NS. (b) As the matter axisymmetrized, the compact accretion disk
becomes homogenous. The full movie showcasing the evolution of the post-merger
accretion disk is available at http://www.oregottlieb.com/disk_gw.html. (c-f)
edge-on (red) and face-on (blue) spectra for model $B$, with different rows
depicting the spectra at various times. The dotted black lines represent the
CE sensitivity curve.

##  Appendix B Comparison of Green’s function

Figure 7 shows a comparison between the Green’s function calculation and the
quadrupole approximation using the spectrum of model $C$. Around the peak
frequency, $\,{f_{\rm GW}}\approx 100\,{\rm Hz}$, the two methods coincide.
However, at low frequencies, the Green’s function overestimates the spectrum,
whereas, at high frequencies, the quadrupole approximation becomes invalid,
resulting in a second peak at $f>{\rm kHz}$. For the spectra shown in Fig. 4,
we use the minimum of the two methods.

Figure 7: The spectra in model $C$ as calculated using the quadrupole
approximation (blue) and the Green’s function (red) for face-on observers.