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

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 ≳ 0 . 1 at ∼ 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 ∼ 10 M ⊙ , the GW spectrum peaks at ∼ 100Hz with some breadth due to various nonaxisymmetric 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 ≳ 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 signal-to-noise ratio of 20 and a rate similar to Type Ib/c supernovae, LIGO–Virgo–KAGRA (LVK) could detect ≲ 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


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 ogottlieb@flatironinstitute.orgfacilitating a deeper understanding of the diverse GW events occurring throughout the Universe.
Similar to inspiraling 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 (CC-SNe).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 proto-neutron star (PNS) f -and g-modes, which are strongly coupled to gravitational radiation (e.g., Nakamura et al. 2016;Andresen et al. 2017Andresen et al. , 2019Andresen et al. , 2021;;Yakunin et al. 2017;Morozova et al. 2018;Radice et al. 2019;Abdikamalov et al. 2020;Mezzacappa et al. 2020Mezzacappa et al. , 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, nonaxisymmetric 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;Siegel et al. 2022).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 nonaxisymmetric modes emerge in hydrodynamic simulations owing to Papaloizou-Pringle instabilities, 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 owing 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. 2000Li et al. , 2001)).RWIs emerge as a result of 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 selfgravity.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, the 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 nonaxisymmetric 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.

Numerical setup
A robust GW signal may emerge from compact, dense accretion disks.Such conditions are prevalent in accretion disks around BHs formed owing 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 a neutron star (NS).Both of these systems exhibit high mass accretion rates of Ṁ ≳ 10 −2 M ⊙ 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 ≲ H/Rx ≲ 0.3, depending on the effective viscosity and BH spin.
B BNS merger 0.  1) adhere to the setup outlined in Gottlieb et al. (2022b), employing a stellar mass of M ⋆ = 20 M ⊙ , and maximum initial magnetization (σ 0 ≡ b 2 /4πρc 2 , where b is the comoving magnetic field strength and ρ is the comoving mass density) within the stellar core capped at max(σ 0 ) = 10 −3 .The angular momentum profile provides a fixed circularization radius of r circ = 25 r g , similar to what is suggested by stellar evolution models during the first several seconds (Gottlieb et al. 2024).The disk scale height is imposed to maintain H/R = 0.1 at all times.The BH mass and dimensionless spin parameter are M BH = 10 M ⊙ , and a BH = 0.8, respectively.In addition to our fiducial model C, we also consider a slowly-spinning BH with a 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 ≳ 0.2 at R ≳ 20 r g , consistent with neutrinotransport collapsar simulations (Issa et al. 2024, in prep.).These disks are thicker because less gravitational energy (due to increased innermost stable circular orbit (ISCO) in Ca) or magnetic dissipation (Cb) is available to heat, making the cooling process less efficient.We also consider a weaklycooled disk with H/R = 0.3 in model Cc which increases to H/R ≈ 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 out-shone 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 ⊙ and 1.78 M ⊙ , described by the LS220 equation of state for NSs (Lattimer & Swesty 1991).The merger product is a BH with a mass of M BH = 2.67 M ⊙ and a BH = 0.68.Surrounding the BH is a substantial accretion disk weighing M d = 0.096 M ⊙ .The vector potential is defined by: and is imposed to be divergence-free.The potential normalization is chosen so that the gas-to-magnetic pressure ratio is β p ≈ 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 radiationdominated gas.The simulations adopt spherical polar coordinates, r, θ, ϕ with a logarithmic cell distribution along the r-direction, spanning from just inside the event horizon to R max , which for collapsars is located outside of the star radius (see Tab. 1).Along the θand ϕ-directions, the cells are distributed uniformly.In total, the base grid comprises N r × N θ × N ϕ cells in the r-, θand ϕ-directions, respectively.For the collapsar and binary NS (BNS) merger simulations, N r = 320, 192; N θ = 96, 96; and N ϕ = 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 MRI , denoting the number of cells per the MRI wavelength, satisfies Q MRI ≳ 30 in both the θ and φ directions, surpassing the minimum requirement of Q MRI ∼ 10 (e.g., Hawley et al. 2011).

GW approximations
The GW fields can be computed using the Green's function method: where T µν e f f = √ −gT µν + t µν , T µν is the contravariant stressenergy tensor in Einstein equations, and t µν is a pseudotensor composed of nonlinear combinations of h µν and their derivatives (see, e.g., Sec.V.A in Thorne 1980 for details).Eq. ( 2) constitutes an integral equation for h µν 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 µν .However, this approach is invalid for selfgravitating 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µ ≪ T 0µ , 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 highfrequency emission, the multipole expansion fails because the reduced wavelength λ/(2π) 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 µν .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π where R d is the density wave orbital radius.For the angular frequency Ω of the emitting gas in the quasi-Keplerian disk: the characteristic GW peak frequency is: where m is the nonaxisymmetric mode.Eq. ( 5) implies that the multipole expansion remains valid when R d ≳ m2 r g .For the quadrupole moment ( m = 2), this corresponds to f GW ≲ 800 Hz.

Comparison to CCSN studies
Previous CCSN studies (e.g., Nakamura et al. 2016;Andresen et al. 2017Andresen et al. , 2019Andresen et al. , 2021;;Yakunin et al. 2017;Radice et al. 2019;Mezzacappa et al. 2020Mezzacappa et al. , 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 GW ≳ 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 GW ≲ 0.2 kHz) arises from the neutrino heating region below the SN shock wave.In these simulations, the shock wave expands to distances of ∼ 10 3 km, for which Eq. ( 3) dictates that the quadrupole approximation is only valid at f ≲ 50 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.

Post-process calculation
We solve the GWs from the simulations by combining two methods: The quadrupole approximation at frequencies below ∼ 0.2 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 quadrupole1 : where all retarded quantities are defined as Then, the strain tensor components are: 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: where T µν 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 x-axis, the two independent solutions constituting the plane wave polarizations in the transverse-traceless (TT) gauge metric are given by: Analogous solutions can be derived for observers situated along the ŷand ẑ-axes.
The characteristic strain is defined as (Moore et al. 2015): with h× ( f ) and h+ ( f ) representing the strain polarizations in the frequency domain.To evaluate the signal detectability and event rate, we estimate the matched-filter signal-tonoise ratio (SNR) as (Flanagan & Hughes 1998;Moore et al. 2015): 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.Parameter hA ( f ) is the dimensionless GW strain, defined by: where F + , F × 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 × = 0.

HYDRODYNAMIC EVOLUTION
The first stellar shell with a circularization radius exceeding the 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 (≲ 1 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 highdensity "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: where γ is the adiabatic index, and g denotes the effective gravity: If R i ≲ 0.25, the flow becomes unstable.We find that in our simulations R i < 0 in the high-density regions in the innermost ∼ 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), 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 bananashaped closed streamlines and it is trapped between prograde and retrograde circular flows on either side of the vortex.These modes emerge stochastically, with the 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 ∼ 20 r g (Gottlieb et al. 2024), as employed in our simulations.Nevertheless, strong magnetic fields facilitate 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.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, η j , defined as: where the integration is calculated on the BH horizon r h , and only for fluid elements with magnetization σ > 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 µ is the four-velocity such that ρ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 η j ≲ 1% of the accreted power owing 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).

Analytic estimates
We begin by providing an analytic approximation for the characteristics of the GW signal.Given that 2π f GW R d /c = m r g /R d ≲ 1 at the peak frequency, we can employ the gravitational quadrupole moment ( m = 2) to estimate the GW strain magnitude as: where ϵ signifies the degree of asymmetry introduced by the azimuthal modes and E d is the disk energy, estimated as: 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 ≈ const ≈ 7 × 10 51 erg.Plugging Eq. ( 18) into Eq.( 17), the characteristic strain of the disk can be estimated as: Since M BH ∼ 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 ϵ ≈ 0.3 in model C and consistently ϵ ≳ 0.2 in all strongly-cooled disks, but ϵ ≈ 0.1 when cooling is weak in model Cc.Suppose the high-density wave dominating the GW emission is in a quasi-Keplerian orbit.In that case, the characteristic GW peak frequency is governed by the angular frequency Ω [Eq.( 4)] and the 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).

Simulation results
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 vari- ability is lower in model Cc, where H/R = 0.3, compared to the strongly-cooled disks with H/R ≳ 0.1.Additionally, in model Cc, all observers would detect a lower frequency peak over ∼ 0.1 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 π/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 Earth's line of sight.We plan to explore this further in future work.
Figure 4 depicts the integrated GW spectra at various time windows for the different models.During the disk formation stage within the first ∼ 1 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 ≈ 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 ∼ 100 m Hz.As m changes stochastically over time, higher multipoles contribute to the spectra some breadth at f ≳ 100 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 ∼ 6 s, as shown by their spectra at 5 s < t < 6 s.Over time, the disk viscously expands and the spectrum broadens to lower frequencies, as seen in the spectrum of model C at 23 s < t < 24 s, resulting in a broader spectrum in model C owing to its longer T s .
We estimate the GW energy using: Fig. 5 illustrates that all models exhibit a similar GW luminosity evolution in time, indicating that the differences in E 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 nonaxisymmetric modes.The asymptotic value of L GW ≳ 10 49 erg 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 ∼ 100 s, on the order of lGRB timescales, implying the integrated GW energy may reach E GW ∼ 10 51 erg, as indicated by model C (see Tab. 1).
In model Cc, characterized by weak cooling, some nonaxisymmetric 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 ∼ 10 3 km.Therefore, low-frequency modes are also present owing 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 ≈ 200 km emerges at t ≳ 2 s, giving rise to peak GW frequencies at f GW ≈ 100 − 300 Hz.The stochastic emergence of various nonaxisymmetric 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 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 radiationreaction timescale for such a binary.The latter is given by (e.g., Shapiro & Teukolsky 1983): where a is the semi-major axis, well-approximated by the orbital radius of the vortex, M ≈ M BH is the total mass of the binary, is well-approximated by the BH mass, and µ is the reduced mass of the binary, well-approximated by the vortex mass m.We find: Gravitational radiation reaction is important if it acts on a shorter timescale than the viscous timescale, which is given by: where α 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-ofmagnitude ratio suggests that the gravitational radiation reaction does not dramatically alter the disk dynamics and the Rossby vortices.The α 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 α ≳ 0.1 (Kotko & Lasota 2012;Ju et al. 2017), thereby unlikely to affect this result.

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, latetime 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 factors 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 ≈ 10 11 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 acc .However, the increase in SNR with t 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 on M BH and R d as per Eq. ( 5).A BH could potentially form with a mass as low as M BH ≈ 3 M ⊙ , resulting in a smaller disk radius and higher frequenciesf GW ∼ M 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 acc ∼ M BH −0.5 .Ultimately, the BH mass growth with Ṁ ≈ 0.1 M ⊙ 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 GW ≈ 100 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 Mpc by a factor of 100 s/T s , given that t 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 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 R lGRB ∼ 1 Gpc −3 yr −1 (Wanderman & Piran 2010), and thus for a typical jet opening angle of θ j ≈ 0.1, the collapsar rate is R Collapsar ≈ 100 Gpc −3 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 R llGRB ≈ 10 4 Gpc −3 yr −1 (Soderberg et al. 2006), which is also comparable to the rate of stripped-envelope Type Ib/c SNe (SNe Ib/c) of R Ib/c ≈ 2.6 × 10 4 Gpc −3 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 ϱ c = 20ϱ 20 (estimating 20 ≲ ϱ c ≲ 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 sce-nario, 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 ∼ 10 −2 ϱ −3 20 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 SN Ib/c events within D ≲ 40 Mpc in the Zwicky Transient Facility (ZTF; Graham et al. 2019;Masci et al. 2019), including events as close as D ≈ 15 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 ≲ 100 s.However, the detectability of such jets is often low, estimated at only ≲ 1% owing to their beamed emission (Frail et al. 2001).Nonrelativistic 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 multimessenger phenomena.

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 life-times 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 ≈ 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 GW ≳ 100 Hz with some breadth due to contributions 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 ≳ 1 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 ≲ 1 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 ϱ 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 ϱ 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 strippedenvelope SNe is accompanied by disk formation and ϱ 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 ≳ 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 accumulate, 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 neutrinocooling 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.Quadrupole Green Figure 7.The spectra in model C as calculated using the quadrupole approximation (blue) and the Green's function (red) for face-on observers.

Figure 1 .
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 (≈ 1 s) feature low-density spiral density waves.(b) At t ≳ 1 s, with the accumulation of infalling gas, a high-density ring forms at R ≈ 15 r g .RWIs in the disk give rise to nonaxisymmetric modes (shown at t ≈ 5 s).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 shown at t ≈ 2 s), 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.

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 m = 2. Fig. 2(b)

Figure 2 .
Figure 2. (a) Logarithmic mass density map of model C at t = 3.4 s with velocity streamlines in the corotating frame of the high-density wave at (x, y) ≈ (−150, 100) 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 3 .
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 nonaxisymmetric density waves close to the BH.Modulations on longer timescales (∼ 0.1 s) in model Cc imply that modes at larger orbits contribute to the GW spectrum in extended disks.

Figure 4 .
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 s illustrate the signal amplification after ≲ 1 s.The spectrogram of model Ca resembles that of C.

Figure 5 .
Figure5.The GW luminosity as a function of time, L GW = dE GW /dt, illustrates that the GW power grows similarly in all models.

Figure 6 .
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 (∼ 10 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.

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 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, β p is the characteristic plasma beta in the disk, max(σ 0 ) is the maximum initial magnetization in the star, a BH is the BH dimensionless spin, M BH is the BH mass, T s is the simulation time, and R max is the outer radius of the simulation grid in units of gravitational radius, r g .Quantities measured from simulations: E GW is the GW energy, and ϱ is the matched-filter SNR for edge-on and face-on observers, respectively, shown for LVK A+ and CE at D = 10 Mpc, and corrected for accretion disk lifetime of t acc ≈ 100 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 GW is the characteristic GW frequency range.