Context. The detection of supermassive black holes SMBHs in highredshift luminous quasars may require a phase of rapid accretion and as a precondition substantial gas influx toward seed black holes BHs from kiloparsec or parsec scales. Our previous research demonstrated the plausibility of such gas supply for BH seeds within starforming giant molecular clouds GMCs with high surface density 10SUP4SUP MSUBSUB pcSUP2SUP facilitating hyperEddington accretion via efficient feeding by dense clumps which are driven by turbulence and stellar feedback. Aims. This article presents an investigation of the impacts of feedback from accreting BHs on this process including radiation mechanical jets and highly relativistic cosmic rays. Methods. We ran a suite of numerical simulations to explore diverse parameter spaces of BH feedback including the subgrid accretion model feedback energy efficiency mass loading factor and initial metallicity. Results. Using radiative feedback models inferred from the slim disk we find that hyperEddington accretion is still achievable yielding BH bolometric luminosities of as high as 10SUP41SUP  10SUP44SUP ergs depending on the GMC properties and specific feedback model assumed. We find that the maximum possible mass growth of seed BHs MSUPmaxSUP SUBBHSUB is regulated by the momentumdeposition rate from BH feedback SUBfeedbackSUBSUBBHSUBc which leads to an analytic scaling that agrees well with simulations. This scenario predicts the rapid formation of 10SUP4SUP MSUBSUB intermediatemassive BHs IMBHs from stellarmass BHs within 1 Myr. Furthermore we examine the impacts of subgrid accretion models and how BH feedback may influence star formation within these cloud complexes.
['stars: black holes', 'stars: formation', 'quasars: general', 'quasars: supermassive black holes', 'Astrophysics - Astrophysics of Galaxies', 'Astrophysics - High Energy Astrophysical Phenomena']
Feedbackregulated seed black hole growth in starforming molecular clouds and galactic nuclei
{'Feedback-regulated Seed Black Hole Growth in Star-Forming Molecular Clouds and Galactic Nuclei': 'Yanlong Shi 1 , 2 ⋆ , Kyle Kremer 1 , and Philip F. Hopkins 1  \n- 1 TAPIR, MC 350-17, California Institute of Technology, Pasadena, CA 91125, USA\n- 2 Canadian Institute for Theoretical Astrophysics, University of Toronto, Toronto, ON M5S 3H8, Canada  \nReceived ; accepted', 'ABSTRACT': "Context. The detection of supermassive black holes (SMBHs) in high-redshift luminous quasars may require a phase of rapid accretion, and as a precondition, substantial gas influx toward seed black holes (BHs) from kilo-parsec or parsec scales. Our previous research demonstrated the plausibility of such gas supply for BH seeds within star-forming giant molecular clouds (GMCs) with high surface density ( ∼ 10 4 M ⊙ pc -2 ), facilitating 'hyper-Eddington' accretion via e ffi cient feeding by dense clumps which are driven by turbulence and stellar feedback.  \nAims. This article investigates the impacts of feedback from accreting BHs on this process, including radiation, mechanical jets, and highly relativistic cosmic rays.  \nMethods. We run a suite of numerical simulations to explore diverse parameter spaces of BH feedback, including the sub-grid accretion model, feedback energy e ffi ciency, mass loading factor, and initial metallicity.  \nResults. Utilizing radiative feedback models inferred from the slim disk, we find that hyper-Eddington accretion is still achievable, yielding BH bolometric luminosities as high as 10 41 - 10 44 erg / s, depending on the GMC properties and specific feedback model assumed. We find the maximum possible mass growth of seed BHs ( ∆ M max BH ) is regulated by the momentum deposition rate from BH feedback, ˙ p feedback / ( ˙ M BH c ), which leads to an analytic scaling that agrees well with simulations. This scenario predicts the rapid formation of ∼ 10 4 M ⊙ intermediate-massive BHs (IMBHs) from stellar-mass BHs within ∼ Myr. Furthermore, we examine the impacts of sub-grid accretion models and how BH feedback may influence star formation within these cloud complexes.  \nKey words. (galaxies:) quasars: supermassive black holes - (galaxies:) quasars: general - stars: black holes - stars: formation", '1. Introduction': "Formation of supermassive black holes (SMBHs), especially those weighing ∼ 10 9 M ⊙ at very high redshift like z ≳ 7 (e.g., recent observations by Bañados et al. 2018; Yang et al. 2020; Wang et al. 2021), has long been an intriguing astrophysical problem (e.g., recent reviews by Inayoshi et al. 2020; Volonteri et al. 2021),. Observations have shown that these SMBHs grow from lighter 'seed' BHs (Yu & Tremaine 2002). The seeds, typically in the mass range of intermediate-mass black holes (IMBHs; ∼ 100 - 10 6 M ⊙ ), are proposed to form in multiple astrophysical scenarios like the direct collapse of the pristine gas (e.g., Bromm & Loeb 2003; Latif et al. 2022), Population III star remnants (e.g. Madau & Rees 2001; Ryu et al. 2016), stellar-mass BHs that undergo hyper-Eddington accretion (e.g., Lupi et al. 2016; Pezzulli et al. 2016; Regan et al. 2019; Shi et al. 2023; Lupi et al. 2024), and runaway stellar mergers in dense star clusters (e.g., Portegies Zwart et al. 2004; Kremer et al. 2020; Shi et al. 2021). If we use a reference radiative e ffi ciency of 0.1, the e -folding time scale at the Eddington accretion rate is ∼ 45 Myr; on the other hand, due to the limited time allowed for SMBH formation (especially for z ≳ 7 quasars that formed only ≲ 10 9 yr after the Big Bang), a sustained phase of fast accretion is inevitable, possibly at super- or hyper-Eddington ( ≳ 500 ˙ M Edd) accretion rates (Inayoshi et al. 2020). 1  \nA number of theoretical works on small-scale BH accretion physics have shown that super-Eddington accretion is achievable and sustainable in principle within the BH accretion disk (Begelman 1979; Quataert & Gruzinov 2000; Blandford & Begelman 2004; Inayoshi et al. 2016). This is also supported by simulations with RHD (Ohsuga et al. 2005), RMHD (Jiang et al. 2014, 2019) and GRMHD (S˛adowski et al. 2015). However, these simulations typically embed the BH inside a gas reservoir with a su ffi -ciently high mass supply rate from large radii ( ≳ 10 3 ˙ M Edd). This level of gas feeding is not necessarily possible in more realistic star-forming (much 'larger-scale' compared to BH disks) astrophysical environments, since the gas medium may fragment due to self-gravity and result in star formation; and feedback from newly-formed stars, including radiative pressure, photoionization / heating, and winds, will develop turbulence or bulk motion in the medium (e.g., see star formation simulations Grudi'c et al. 2018, 2021a). Strong stellar feedback near an accreting BH may deplete the available fuel and challenge the feasibility of superEddington accretion (Dubois et al. 2015; Habouzit et al. 2017; Bower et al. 2017). Light seeds, like IMBHs or stellar-mass BHs, are especially prone to the accretion challenges inherent to these environments since regimes influenced / dominated by their gravity (i.e., the Bondi-Hoyle radius, Hoyle & Lyttleton 1939; Bondi 1952) are small compared with galactic scales.  \nTo reconcile the problem, 'larger-scale' (kpc or pc scale) simulations are not only required to model star formation and feedback in su ffi cient details, but also to resolve the some gas dynamics below the gravitational capture scale (Bondi-Hoyle ra-  \ndius) of at least some rapid-accreting BHs. In Shi et al. (2023) (Paper I hereafter), we made some first steps towards addressing this problem, using simulations of star formation and feedback in giant molecular clouds (GMCs) similar to that in Grudi'c et al. (2018), and studying BH seed growth in the absence of BH feedback. Although the stellar feedback makes the GMC uneven and turbulent, it also generates dense clumps and shocks (Klessen 2000), which may have low internal velocity dispersion (Mac Low & Klessen 2004; McKee & Ostriker 2007) - these factors make the simulation possible to resolve the gravitational capture of gas for some 'lucky' BHs encountering dense gas clumps with small relative velocity, since the Bondi-Hyle radii are large under these circumstances. Similar to the star formation e ffi -ciency explored in the literature (Grudi'c et al. 2018; Kim et al. 2018; He et al. 2019; Fukushima et al. 2020; Grudi'c et al. 2021a; Kim et al. 2021; Chevance et al. 2023), BH accretion is also regulated by the surface density ( Σ 0 ≡ M cl /π R 2 cl ) of the cloud: higher Σ 0 means higher self-gravity ( ∼ GM cl / R 2 cl ∼ π G Σ 0) which keeps the cloud bound rather than disrupted by stellar feedback. Paper I found that for clouds with high surface density ( Σ 0 ≈ 10 4 M ⊙ pc -2 ), the gas supply reaching the BH's accretion disk could be su ffi cient to support hyper-Eddington accretion with f Edd ≡ ˙ M BH / ˙ M Edd ∼ 1000 for at least a small but non-negligible amount of seeds.  \nMore realistically, BH feedback mechanisms, like electromagnetic radiation (Fabian 2012; Hopkins et al. 2016), winds / jets (Silk & Rees 1998; Murray et al. 2005; Di Matteo et al. 2005; Ostriker et al. 2010; S˛adowski et al. 2016; Torrey et al. 2020; Hu et al. 2022; Massonneau et al. 2023), and cosmic rays (CRs; Sijacki et al. 2008; Guo & Mathews 2012; Zweibel 2017; Su et al. 2021; Ishibashi & Fabian 2023) will deposit energy and mechanical momentum back to the gas reservoir. The energy and momentum heat and deplete the gas fuel, or even prevent its infall to the BH (Dubois et al. 2013), potentially challenging the hyper-Eddington accretion scenario. Still, close to the BH ( ≲ 1000 r Sch), at su ffi ciently high accretion rates and opacities, photons are e ff ectively trapped with the inflow, which means the radiated energy (Watarai et al. 2000; Madau et al. 2014, and other references above) will only grow sub-linearly with the accretion rate (measured as ˙ M BH / ˙ M Edd). This is also true for the mass outflow powered by radiation (Hu et al. 2022). As a result, the feedback strength at high accretion rates is suppressed and hyper-Eddington accretion could remain possible.  \nThough BH feedback is not included in the simulations of Paper I, there is an estimate in terms of radiative feedback based on the density criteria for hyperEddington accretion predicted by Inayoshi et al. (2016): n ∞ ≳ 10 5 ( M BH / 10 4 M ⊙ ) -1 ( T ∞ / 10 4 K) 3 / 2 cm -3 , where n ∞ and T ∞ are the density and temperature of gas 'near' the BH (at sub-pc scales). For BHs in simulations that grow rapidly, especially those inside GMCs with Σ 0 ≈ 10 4 M ⊙ pc -2 , the ambient gas density is above this critical value, which suggests hyper-Eddington accretion may be achievable even if there is radiative feedback. Still, there is a necessity to include BH radiative feedback in simulations like Paper I for completeness since BH feedback may still a ff ect GMC behaviors at larger scales. The influence of other forms of BH feedback, namely, mechanical feedback and cosmic rays on accretion is also to be determined.  \nIn this article, we extend the discussion of Paper I and explore the impact of multiple BH feedback mechanisms, including radiation, mechanical outflow, and CRs. In particular, we are interested in how the BH feedback mechanisms determine the maximum mass a BH can reach from accretion, and how they scale in di ff erent GMCs. To reach this goal, we run a suite of  \nFig. 1. Mass and energy flow of BH accretion and feedback. This study extends the framework of star formation in Grudi'c et al. (2018) and BH accretion without feedback in Paper I. To implement BH feedback, we use a model of mass transfer from the disk to the BH. With only a fraction of the accretion flow going to the BH, other portions are deposited back to the GMC in the form of radiation, jets / winds, and cosmic rays. We note that since the feedback energy e ffi ciencies are typically small ( ≲ 0 . 1), ˙ M BH GLYPH<17> ˙ M BH , 0.  \n<!-- image -->  \nnumerical simulations by adopting the basic setups and physics involved in Paper I where the authors showed hyper-Eddington accretion was possible without BH feedback, but with a new sub-grid BH accretion and feedback model, as well as related physics.  \nThe article is organized as follows. In § 2, we introduce the background and implementations of BH accretion and feedback, as well as the simulation initial conditions. In § 3 , we present the result with the fiducial feedback models and explore the effects of di ff erent accretion / feedback physics. In § 4, we discuss the result and important caveats of our simulations. Finally, we conclude in § 5.", '2. Background and Method': 'This study follows the basic framework of Paper I: each simulation is initialized with a number of seed BHs inside a precollapse pre-star-formation parsec-scale GMC, and then evolves self-consistently with MHD + gravity and additional physics including star formation / feedback, BH accretion / feedback, etc.', '2.1. Overview': "Our simulation follows the same numerical framework for star formation and feedback as that in previous star-formation simulations (Grudi'c et al. 2018, 2021b) and the previous study of this series (Paper I). As a short summary, we use the meshless, Lagrangian, Godunov MHD code GIZMO (Hopkins 2015, 2016; Hopkins & Raives 2016) in its Meshless Finite Mass (MFM) mode, including the physics of self-gravity, radiative cooling, star formation, and feedback. In particular, star formation and feedback are based on the FIRE-3 implementation of the Feedback In Realistic Environments 2 (FIRE) framework (Hopkins et al. 2018b, 2023). Each 'star' particle in the simulation represents an IMF-averaged ensemble of stars, which evolves and deposits mass outflow, metals, and radiation back to the ISM as its feedback (Hopkins et al. 2018a). This treatment has been shown to reasonably reproduce many of the properties of star formation GMCs, like the star formation e ffi ciency, star cluster dynamics, and cluster mass distribution (Grudi'c et al. 2018, 2021b).  \nTable 1. Initial conditions involved in this study. We chose some gas complexes in the 'high surface-density' ( Σ 0 = 13000 M ⊙ pc -2 ) group of Paper I, with the initial radius ( R cl) of 5 pc or 50 pc. Each kind of gas complex is simulated at high or low resolution. Starting from the third column, we show the initial free-fall time ( t ff ), gas mass resolution ( m gas), softening radius for stars ( r star soft ) and BHs ( r BH soft ), characteristic jet energy e ffi ciency ( ϵ w , c) and velocity ( v w , c) as described in §2, number of BH seeds ( N BH), and initial BH mass ( M ini BH ) for di ff erent kinds of simulations.  \nThe simulation also follows the same scheme for BH accretion from gravitational capture as that in Paper I, based on that stated in Grudi'c et al. (2021a): gas cells are captured by a BH if (i) they are within a preset 'sink radius' near the BH, (ii) they are gravitationally bound to the BH, and (iii) their individual Keplerian orbit is within the sink radius (Bate et al. 1995). Moreover, we choose the softening radius of BHs to be the same as its sink radius and set the value to make sure that the Bondi-Hoyle accretion radius is resolved (i.e., greater than the sink radius) for fast-accreting BHs, which is extensively discussed in Paper I.  \nFor this study, we extend the BH accretion scheme with a sub-grid model of BH accretion to account for BH feedback. To help explain this, we present a schematic view of the simulation in Fig. 1. Accretion onto a BH can be decomposed into two successive phases: (i) the gravitational capture of gas from the surrounding ISM (as a reservoir), or from ∼ pc to sub-pc scales, as studied in Paper I and resolved in our simulations here; (ii) the mass transfer from the BH accretion disk (as a secondary reservoir) to the event horizon, or from sub-pc to ∼ km (or r Sch) scales, which we do not resolve in our simulations. Along with the accretion, gravitational energy is transformed into other forms of energy, like radiation, winds / jets, and cosmic rays.  \nAt the second phase (smaller scales), not all of the mass / energy flows reach the BH: (i) some fraction of the accretion flow will be powered by kinetic / radiative pressure and form jets / winds; (ii) some fraction of the energy 'leaks' from the disk and is deposited into the GMC, in the form of radiation and CRs. As a result, intuitively BH feedback mechanisms typically serve as negative factors that suppress the BH accretion in both two ways: (i) reducing the final mass transfer rate onto the BH due to energy and mass loss; (ii) preventing additional gas in GMC from accreting.  \nMoreover, there is thus a huge mismatch among the scales and physics involved, and the second process is beyond the resolution limit ( ∼ 0 . 1 pc) of our MHD simulation. So we treat this with a sub-grid model: the BH sink particle accretes gas from the parent simulation on resolved scales, stores the mass inside a gas reservoir (or disk), then transfers the mass to the BH at a rate ˙ M dep, which is also used to calculate mass and energy flow rates for di ff erent feedback mechanisms.  \nDetails of the sub-grid model are expanded below, following the schematic 'order' of mass flow from the disk to the BH shown in Fig. 1. However, we note that the 'order' is only for defining quantities like the mass loading factor and energy e ffi ciency, not in the chronological sense. Moreover, since the feedback energy e ffi ciency (defined as the ratio between the outflow energy and ˙ M BH c 2 ) is typically very small (at most 0 . 1) in the simulation, we treat ˙ M BH , 0 and ˙ M BH interchangeably in the simulation and this article, without impacting our results at least in the order-of-magnitude sense. For example, the radia-  \nn power is defined as L bol ≡ ϵ r M BH , 0 c 2 in Fig. 1, the relation L bol GLYPH<17> ϵ r M BH c 2 still holds as a very good approximation.", '2.2.1. Disk mass depletion': "As shown in the flowchart Fig. 1, gas is captured by sink particles, via resolved gravitational accretion. It is then added to a 'reservoir' of mass M d, which we will call the 'disk', surrounding the BH of mass M BH. We require a sub-grid model for the accretion rate ˙ M d from the disk to the BH, since this determines the feedback properties, so we define ˙ M d = M d / t dep in terms fo the depletion time t dep.  \nThe classical way to describe the BH accretion disk is the thin disk model (Shakura & Sunyaev 1973), where the accretion rate follows  \n˙ M ∼ 2 πα Σ d c 2 s / Ω K . (1)  \nHere Σ d, c s, and Ω K are respectively the surface density, sound speed, and Keplerian frequency at a certain radius of the disk; α is a dimensionless constant that characterizes the e ff ective viscosity. Motivated heuristically by Eq. (1), we set  \nt dep = t dyn , sink / (2 πα ) . (2)  \nNote that we redefined α by absorbing the factor of ( Ω K r sink / c s) 2 and other constants. Here t dyn , sink ≡ [ r 3 sink / G ( M BH + M d)] 1 / 2 is the dynamical time scale of at the sink radius, which is the fastest possible mass-depletion time. By varying α we may bracket different conditions including slow and fast mass depletion rates.  \nWe caution that for realistic BH accretion, the mass inflow rate can be dependent on the radius since there is mass outflow in the form of winds (Blandford & Begelman 1999; Hu et al. 2022). Although in Eq. (1) we actually assume a scaleindependent mass inflow rate, the e ff ect of BH wind feedback is indeed considered (the next subsection), where we assume only a fraction of the mass can reach the BH while the remaining are deposited back to the GMC in the form of wind and jet.  \nIn the simulation, we also add two constraints to the BH mass transfer. We first set a maximum mass of the accretion disk to represent the e ff ects like fragmentation limiting accretion if mass 'stalls' in the disk, which is quantified by the upper limit of M d / M BH: once M d reaches this limit, no gas cells will be absorbed by the BH in the next time step. Secondly, there is also an upper limit on the Eddington ratio which is defined as f Edd ≡ ˙ M BH / ˙ M Edd. We vary the upper limit M d / M BH. Due to our particular interest in hyper-Eddington accretion which is predicted to be extreme (Inayoshi et al. 2020), we fix the upper limit of f Edd to be 1000, but obtained qualitatively similar results so long as this limit is ≫ 1, as we discuss below.  \n≡  \nTable 2. Free parameters in our BH feedback sub-grid model. Along with the fiducial setup, we also study the impact of di ff erent physical processes by varying the corresponding quantities.  \nFig. 2. Sub-grid models for radiative feedback (top panel) and mechanical feedback (bottom two panels). Top : energy and momentum ejection rates as a function of accretion rates, including 'log-form' models inferred from the slim disk ( solid ), and simple ones with constant ϵ r ≡ L / ( ˙ M BH c 2 ) ( dashed ). For realistic models the feedback strength Bottom : The momentum ejection rates for mechanical feedback with constant ϵ w ≡ ˙ E kin / ( ˙ M BH c 2 ) ( dashed ), and the model in Hu et al. (2022) - though the released kinetic energy is almost constant ( E kin ≈ 4 × 10 -4 ˙ M Edd c 2 ), the momentum injection rate grows linearly ( ˙ p ≈ 0 . 03 f Edd ˙ M Edd c ) due to the variable mass loading factor.  \n<!-- image -->", '2.2.2. Mechanical feedback': "In this simulation, we treat the mass and energy flow as shown in Fig. 1. For a given ˙ M dep, we assume some fraction of this will be ejected in the format of winds / outflows / jets, which we combine here as 'mechanical' outflow and will simply refer to as 'winds' throughout. We define the jet mass loading factor as η w ≡ ˙ M w / ˙ M BH , 0, where ˙ M w is the jet mass outflow rate, so ˙ M w = ˙ M dep · η w / (1 + η w). The kinetic luminosity of the winds can be parameterized with a coe ffi cient ϵ w, such that L w = ϵ w ˙ M BH , 0 c 2 , or with e ff ective outflow velocity v w, then L w = M w v 2 w / 2.  \nFor fiducial simulations, we set η w = 1 and a characteristic ϵ w such that the mechanical feedback energy is capable of disrupting the whole cloud, i.e., the accumulated jet energy fills the potential well of the GMC: R ϵ w ˙ M BH , 0 c 2 d t ∼ GM 2 cl / R cl. Treating ϵ w as a constant, we integrate over BH accretion history and find  \nthe LHS turns into ϵ w ∆ M BH , 0 c 2 , where ∆ M BH , 0 is the total mass of gas reaching near the event horizon. We choose a characteristic value that ∆ M BH , c = 10 4 M ⊙ . Then the characteristic ϵ w is ϵ w , c = M cl / ∆ M BH , c( v cl / c ) 2 . Here v 2 cl = GM cl / R cl is the characteristic circular velocity of the cloud. Following the energy argument we assumed above, if ϵ w ≪ ϵ w , c then the mechanical feedback is insignificant; while if ϵ w ≫ ϵ w , c then the strong mass outflow may disrupt the GMC quickly. It also implies the critical jet velocity at the critical jet luminosity: v w , c = v cl p M cl / ∆ M BH , c.  \nMotivated by simulations of hyper-Eddington accretion (S˛adowski et al. 2016), we assume fiducial η w = 1 but vary v w widely. Moreover, some small-scale simulations fitted sub-grid models which vary η w as a function of the accretion rate. In this work, we implemented a version presented in Hu et al. (2022), which roughly sets η w = f Edd -1 (where f Edd ≡ ˙ M BH / ˙ M Edd) while keeping a nearly constant kinetic energy output E kin ≈ 4 × 10 -4 ˙ M Edd c 2 at the super-Eddington regime of f Edd ≳ 3.  \nNumerically, the winds / jets in the simulation are implemented as high-resolution gas cells that are ejected in a bipolar fashion from the BH (Torrey et al. 2020). In the rest frame of the BH, these gas cells are ejected along the spin of the BH, which is further determined by the accumulated angular momentum from the gas. A similar implementation was also used and described in Su et al. (2021); Grudi'c et al. (2021a).", '2.2.3. Radiative feedback': 'Given some remaining accretion rate ˙ M BH , 0 (after the loss of winds from ˙ M dep), we can calculate the bolometric accretion disk luminosity L bol in terms of some radiative e ffi ciency: L bol = ϵ r ˙ M Edd c 2 . Given L bol, we explicitly model the outgoing BH spectrum based on the template in Shen et al. (2020) and its e ff ects in our multi-band (X-ray, EUV, FUV, NUV, OIR, far-IR) radiation treatment (as that in Hopkins et al. 2020), including its e ff ects on photo-ionization and photo-heating, Compton heating, and radiation pressure on gas and dust, as described in Hopkins et al. (2023).  \nAnalytical studies and numerical simulations have shown that due to the photon-trapping e ff ect on small scales near BHs, the bolometric luminosity grows slowly with increasing accretion rate (rescaled to f Edd ≡ ˙ M BH , 0 / ˙ M Edd), typically in a logarithm form. Due to this reason, super-Eddington accretion is not limited by the radiative energy loss (Madau et al. 2014).  \nIn our simulation, we therefore compare various constantϵ r choices with two models for a variable ϵ r (as shown in Fig. 2). The first one by Watarai et al. (2000) is analytically derived from the thin disk model, which features a linear form at low accretion rate (with contact ϵ r = 0 . 1) and a logarithm form at high  \nFig. 3. Visualization of simulations with and without BH feedback physics, each from a snapshot when a BH is undergoing rapid growth. Left two zoom-in plots : Gas morphology for the 10 6 M ⊙ GMCs without ( top ) or with ( bottom ) BH feedback. From left to right, we zoom in at a BH undergoing hyper-Eddington accretion. As the smallest scale, each panel shows the velocity field near the BH, with the circular velocity ( √ GM cl / R cl) presented as a gauge. Bottom two zoom-in plots : the same as the left plots, but for the 10 8 M ⊙ GMC.  \n<!-- image -->  \naccretion rate. Another form is fitted by Madau et al. (2014) from original simulations of relativistic slim disks in S˛adowski (2009), which also features dependence on the BH spin parameter a . We set a = 0 for this form in our simulation since Madau et al. (2014) found that the super-Eddington accretion e -folding time is nearly independent of a .', '2.2.4. Cosmic rays': "Previous simulations have argued for the importance of cosmic rays (ultra-relativistic particles) accelerated around BHs or jets, e.g., Sijacki et al. (2008); Guo & Mathews (2012); Su et al. (2021). In addition, CRs also lose energy due to Coulomb interactions, ionization, and catastrophic losses, which also heat the ambient gas (e.g., Guo & Oh 2008).  \nIn this study, we simulate CRs based on implementations described in Hopkins et al. (2022), representing them as a 'single-bin' (relativistic ∼ 1 -10 GeV protons) which obey a two-moment streaming + di ff usion + acceleration equation, with Alfvenic streaming and a fiducial scattering rate set to ν ∼ 10 9 s -1 calibrated to solar system and Milky Way cosmic ray observations (from Voyager, Fermi, and AMS-02, see Hopkins et al. 2022). All relevant cosmic ray loss / cooling and gas coupling / scattering / ionization / heating details are included as described therein.  \nIn our simulations, CR feedback from BHs is coupled alongside the mechanical feedback (i.e., the newly spawned jet cells). We vary the strength of the BH CR feedback through its energyloading coe ffi cient ϵ CR, defined by the injection rate L CR = ϵ CR ˙ M BH , 0 c 2 . Additionally, for each simulation with the BH CR feedback, we set a constant CR energy ejection rate of 10% for supernovae to account for the CR feedback from stars. This particular e ff ect is not considered in other sets of simulations without CRs.", '2.3. Simulation initial conditions': "We follow the same initial condition (IC) for GMCs as Paper I, but only choose GMCs with the highest initial surface density ( Σ 0 = 13 , 000 M ⊙ / pc 2 ) with di ff erent radii (5 or 50 pc) and therefore, di ff erent masses (10 6 or 10 8 M ⊙ ). These are chosen as the clouds in Paper I where BHs accrete significantly without BH  \nFig. 4. Mass evolution and accretion rate for BH particles in di ff erent GMCs. In the first and third panels (from top to bottom), we select 5 BHs with significant mass growth from high-resolution simulations of the two complexes and plot the evolution of the mass bound to the BH (BH + disk, solid ) and the BH-only mass ( dashed ). In the second and fourth panels, we show the accretion rate for the same BHs, each rescaled to its Eddington accretion rate ( f Edd = ˙ M BH / ˙ M Edd).  \n<!-- image -->  \nff  \nfeedback, while other cases would be unlikely to see significant mass growth with the feedback physics. The clouds are initially with the solar metallicity ( Z ⊙ ) in the fiducial case. The ICs are listed in Tab. 1. The time limit for each simulation is 2 . 5 t ff ,  \n0  \nwhere t ff = [ R 3 cl / (8 GM cl)] 1 / 2 is the initial free-fall time of the cloud. For each mass group, there are low (initially with 64 3 equal-mass gas cells) and high (128 3 ) resolution runs. Force softening is set as described in Paper I.  \nEach simulation is initialized with a number of BH seeds with random mass, position, and velocity, which are sampled following the same method as Paper I: BHs are randomly distributed following uniform spatial distribution within the volume; the initial velocity magnitude is confined to be below the local circular velocity (for discussions of these choices, see Shi et al. 2023). In particular, the BH mass follows log-uniform distribution within 10 (100) - 1000 M ⊙ (10 4 M ⊙ ) for the 10 6 M ⊙ (10 8 M ⊙ ) cloud, which covers the range of stellar-mass BHs and IMBHs.  \nThe independent parameters we vary to quantify the impact of BH feedback are listed in Tab. 2. For each parameter, we set a fiducial value as the 'baseline' and vary each parameter in turn.", '3.1. Fiducial results: radiative-inefficient models': 'Compared with Paper I, we include a \'log-form\' radiativeine ffi cient model (Madau et al. 2014) to account for the radiative feedback from slim disks and a constant wind energy e ffi ciency. This naturally causes a suppression in the mass accretion rate of seed BHs (as will be described in §3.2). Here we focus on the two simulations of 10 6 M ⊙ and 10 8 M ⊙ cloud complexes with \'high resolution\' (initially with 128 3 gas cells, see Tab. 1).  \nAs a first impression, Fig. 3 compares the gas morphology of simulations with or without BH feedback at di ff erent scales. From left to right, we zoom in on a BH experiencing significant mass accretion we select. The particular BH which accreted the most mass throughout the simulation time range, and the snapshot shown is also chosen when that BH is in its fastestaccreting phase. Each panel shows the line-of-sight averaged density, which is defined as R ρ ( s )d s / R d s . At the smallest scales (panels in the right column), we evaluate v gas -v BH and show the velocity field within a thin layer of 1 / 8 of the field of view (centered at the BH) on top of the density field.  \nAs was also suggested in Paper I, the fast-growing BH seed is often near a dense gas clump ( ρ ≳ 10 5 M ⊙ / pc 3 , typically more than 100 times denser than the mean density). With a small relative velocity | v gas -v BH | < | v BH | which ensures e ffi cient gravitational capture of gas. Despite feedback from BHs, these features are also observed in corresponding simulations of both 10 6 M ⊙ and 10 8 M ⊙ complexes, which also fit the expectations of the Bondi-Hoyle accretion.  \nIn the absence of BH feedback, the gas in the simulations is relatively dense. Additionally, near the accreting BH without feedback (see the first and third row of Fig. 3), there are features like disks and spiral arms at ∼ pc scales at late times, which means that there is coherent gas inflow due to the potential well of the GMC. At this time, there is also star formation in the 10 6 M ⊙ GMC near the BH, creating a cavity in the gas distribution (see the first row of Fig. 3), though not observed in the 10 8 M ⊙ GMCthrough the simulation time. In contrast, when BH feedback is turned on, additional bubbles, outflows, sheets, and clumps appear and the inflow is relatively incoherent.  \nFig. 4 shows the BH mass growth history in both fiducial simulations, by selecting 5 seed BHs with the highest final-toinitial mass ratios. For each BH, both the BH-only mass and the sink-particle (BH + disk) mass are shown. As in Paper I, single  \nFig. 5 shows the cumulative distribution function (CDF) of ∆ M BH (relative mass growth at the end of the simulation of each BH) measured from simulations. As expected in Paper I, only a modest fraction ( ≲ 10%) of BH seeds have appreciable mass accretion as the resolution allows. Comparing the two GMCs, we find that BHs in the 10 8 M ⊙ GMC have much more significant mass growth, which is again expected from Paper I. Moreover, we compare the low- and high-resolution simulations with the same BH feedback recipes in the same figure and found a qualitatively similar distribution.  \n<!-- image -->  \n/circledot  \nFig. 5. Cumulative distribution function (CDF) of the accreted mass, Eddington ratio, and bolometric luminosity for BH populations in simulations with the fiducial BH feedback setups. In the right two panels, we show the distribution of the \'peak\' value throughout the evolution. Both simulations of high ( solid ) and low ( dashed ) resolutions are displayed and show qualitatively similar trends. The artificial flat \'plateau\' at the left of some lines reflects the mass resolution limit of that simulation.  \nBHs grow rapidly when they encounter dense clumps, while the accretion rate is low most of the time.  \nWe see little di ff erence in the BH-only mass and the sink particle mass at the early stage of the evolution when the BH accretion rate is low. If the BH feedback is significantly strong, the overall mass accretion is suppressed, and the BH accretion is determined by the material that is already bound to its gravity. Due to the timescale for mass transfer from the disk to the BH, there is a \'phase lag\' between M BH and M sink, but the e ff ect is not significant in our fiducial model.  \nFor single BHs, the BH mass accretion rates ˙ M BH (the actual accretion rate arriving at the BH) are also shown in Fig. 4. We scale the rate in the unit of the Eddington rate for each BH. Here we show the same 5 BHs selected in Fig. 4. These samples typically reach hyper-Eddington accretion abruptly at some time in their evolution and maintain the status for a short period of time ( ∼ 0 . 1 t ff ), during which time ˙ M BH is capped by the upper limit set in the code. Once the fast accretion phase is terminated, the mass transfer from the disk to the BH declines with the characteristic scale of t dep as defined in § 2.  \nThus, from Fig. 4, there is evidence that hyper-Eddington accretion is achievable for BHs in the dense complexes we simulated, even with some BH feedback model. However, the BH feedback can have a negative impact on accretion, which we will quantify in § 3.2.  \nThroughout the simulation, we find these BHs with resolved accretion activity typically reach super-Eddington accretion at some stage of the evolution, and a small fraction ( ∼ 1%) of the total BHs can even reach the preset cap ( f Edd ≥ 1000). We also check the bolometric luminosity, defined as L bol = ϵ r( f Edd) ˙ M BH c 2 . For the 10 6 M ⊙ complex, BHs can emit radiation at ∼ 10 42 erg / s at some stage of their evolution, while the value can be ∼ 10 42 - 10 43 erg / s for BHs in the 10 8 M ⊙ complex.  \nv ini BH [ √ GM cl /R cl ] j ini BH [ √ GM cl R cl ] Fig. 6. Dependence of BH accretion on initial conditions of BHs. From our high-resolution fiducial simulations, we mark BHs with resolved accretion in a phase space of the accreted mass ( ∆ M BH) versus The initial mass ( M ini BH ), distance to the cloud\'s center ( r ini BH ), velocity ( v ini BH ), and specific angular momentum ( j ini BH ). We find only shallow dependence on these initial conditions. We may also clearly see the mass resolution limit from the lower \'cuto ff \' in ∆ M BH for the 10 8 M ⊙ complex.  \n<!-- image -->  \nFor BHs in these two complexes, the CDFs of ∆ M BH and L bol are very similar despite di ff erent mass resolutions and number of sampling points (i.e., the number of seed BHs in each simulation). We find that the distribution function of ∆ M BH can be well-fitted with a log-normal distribution: f ( ∆ M BH) ∝ exp h -log( ∆ M BH /µ ) 2 / (2 s 2 ) i . For both fiducial runs of the 10 6 M ⊙ and 10 8 M ⊙ complexes, the best-fit parameters are approximately µ = 10 -2 . 5 M ⊙ , s = 2 . 25. Explicitly, the CDF of ∆ M BH [M ⊙ ] is  \nP ( ∆ M BH > x ) = 1 2 " 1 -erf log x + 2 . 5 2 . 25 √ 2 !# . (3)  \nIn Paper I, we demonstrated that there is weak dependence of the results on the initial properties of BH seeds, like their initial mass, position, and velocity magnitude, as long as they are gravitationally bound to the cloud. Similarly, we explore such dependence for our new simulations by checking ∆ M BH for each BH at the end of the simulation. In each panel of Fig. 6, we list BH\'s initial mass, radii from the GMC\'s center, velocity magnitude, and the specific angular momentum. Again, for both the 10 6 M ⊙ and 10 8 M ⊙ cloud, there is no (or little) correlation with ∆ M sink. The result means that BH accretion in a turbulent environment is stochastic and much information from the initial condition is smeared in the process. This is why we simulated N BH ≫ 1 seeds in each cloud - there is only a ∼ 1 -10% probability of a seed being \'lucky\' and encountering a dense clump for accretion.  \nFig. 7. Parameter survey of the accretion model, which utilizes an effective α -parameter to quantify the accretion rate. Here we vary α and the upper limit of disk mass. In the fiducial case ( dashed ), α = 0 . 1 and M d / M BH ≤ 10.  \n<!-- image -->  \n/circledot  \nFrom the same figure, BH accretion has no strong dependence on the initial position or velocity. However, from the last panel, there is stronger evidence that BHs with significant mass growth tend to have low initial specific angular momentum ( ≲ 0 . 2 √ GM cl R cl). This is expected from the evolution of the cloud, since the gravitational collapse brings dense gas inward, and only BH seeds with low specific angular momentum can reach the inner part of the complex.', '3.2. Parameter survey: effects of different physics': 'To check the e ff ects of di ff erent physics, we perform a parameter survey by varying parameters in our sub-grid model, as listed in Tab. 2. For considerations in computational costs, we perform these simulations at low resolution (see Tab. 1).', '3.2.1. Sub-grid accretion model': "In Fig. 7 we show the e ff ect of free parameters in the disk mass depletion: the upper limit of M d / M BH, and the e ff ective α parameter for viscosity. In this plot, we show the CDF of the change in BH + disk and the BH-only mass, as well as the CDF of the Eddington ratio for BHs through evolution. Comparing the fiducial case, enforcing a smaller upper limit M d / M BH ≤ 1 leads to accreted mass lower by a factor of ∼ 10, and fewer BHs can reach hyper-Eddington accretion. This is a natural consequence of the smaller disk mass, which means that less mass can reach the BH + disk system. Because the mass comes in rapidly, it must either be able to accrete extremely rapidly through the disk ( α > 1), or be able to avoid fragmentation and star formation ( M d / M BH > 1) to reach the 'full' potential growth. But we stress that we still clearly see hyper-Eddington accretion for M d / M BH ≤ 1.  \nThen we test higher (1) or lower (0.01) values of the α parameter, in contrast with the fiducial one (0.1). As shown in Eq. (1), smaller α means slower mass transfer from the disk to  \nFig. 8. Parameter survey of the energy e ffi ciency in di ff erent feedback mechanisms based on experiments with the 10 6 M ⊙ complex. From left to right, the three columns show the CDF of accreted mass by BHs, as well as the highest possible f Edd and L bol throughout the evolution. Each row represents a set of simulations with varying models of a specific feedback mechanism (labeled in each left panel), with the fiducial model marked in black dashed lines.  \n<!-- image -->  \n/circledot  \nthe BH, this is true in our tests: by fixing M d / M BH ≤ 10, there is a larger deviation between the BH + disk curve and the BH-only curve in the mass CDF once α is smaller, and the Eddington ratio is typically lower for smaller α . The trend is also true for tests with fixed M d / M BH ≤ 1.  \nThe parameter α also a ff ects the mass accretion onto the BH + disk system in the simulation. If α is low, the feedback is also weaker due to less mass reaching the BH; if α is high, mass depletion from the disk is e ffi cient and the disk mass can be supplemented quickly in a very gas-rich environment. Both e ff ects are positive factors for accumulating more mass in the BH + disk system. This argument is also reflected in Fig. 7 when comparing tests with di ff erent α , where α = 1 and α = 0 . 01 tests sometimes have more mass accretion than the fiducial α = 0 . 1 case.  \nStill, we will show below that none of these e ff ects in the sub-grid accretion model is as dominate as the largest variations in the strength of feedback.", '3.2.2. Radiative feedback and metallicity': "We study the e ff ect of radiative feedback by varying its energy e ffi ciency ϵ r, including two 'log-form', accretion-ratedependent models from Watarai et al. (2000) and Madau et al. (2014), as well as fixed ϵ r at 10 -9 - 0 . 1. Similar to the layout of Fig. 5, the top panels of Fig. 8 and Fig. 9 show the CDF of accreted mass, the Eddington ratio, and the bolometric luminosity of the BH population.  \nComparing fixed-value energy e ffi ciencies, there is a clear trend that once ϵ r ≳ 10 -3 (10 -1 ) for the 10 6 M ⊙ (10 8 M ⊙ ) GMC, the accretion onto BHs is suppressed. The same behavior appears in the Eddington ratio. Despite relatively strong feedback, there are 1% of BHs reaching the preset cap in the Eddington ratio in many simulations, while the fraction of these BHs drops when the feedback is strong. Finally, the luminosity of BHs typically increases with radiative e ffi ciency as expected.  \nAs also shown in Fig. 2, the radiative e ffi ciency at high accretion rates ( f Edd = 1000) can be as low as 3 × 10 -2 for Madau et al. (2014) and 10 -2 for Watarai et al. (2000). Comparing simulations with fixed-value e ffi ciencies, we find that these radiativeine ffi cient models behave closely like those with fixed ϵ r = 10 -3 and ϵ r = 10 -2 respectively. The feedback strength is thus mainly dominated by the hyper-Eddington regime for these models. As the result, the Madau et al. (2014) model is less feedbackdominated.  \nWe further study the impact of the GMC's initial metallicity. In the no-BH-feedback limited studied in Paper I, experiments with di ff erent initial metallicity found no (or little) correlation between BH accretion and the initial metallicity. However, higher metallicity means tighter coupling between the radiation and wind due to more dust grains in the ISM, which may imply stronger radiative pressure (e.g., Larson & Starrfield 1971), which is one of the radiative feedback mechanisms (alongside with heating and photoionization) implemented in the simulation. In this study, we vary the initial metallicity for two kinds  \nFig. 9. Parameter survey of the energy e ffi ciency in di ff erent feedback mechanisms based on experiments with the 10 8 M ⊙ cloud complex. Conventions are the same as that in Fig. 8.  \n<!-- image -->  \n/circledot  \nof radiative models: the fiducial Madau et al. (2014) model and the fixed-value ϵ r = 0 . 01 model (which is in the strong feedback limit). In Fig. 10, we plot the maximum accretion in the BH population as a function of the initial metallicity. For most possible combinations of the initial GMC masses and feedback models, there is still no strong dependence on the metallicity. One exception is the 10 6 M ⊙ GMC with fixed ϵ r = 0 . 01, there is a clear trend that BH accretion drops as the metallicity increases.  \nThis can be explained with the aid of Fig. 8 and Fig. 9. For the 10 8 M ⊙ GMC, both radiative-feedback models above are in the weak-feedback limit even at the solar metallicity (since the CDF of BH mass accretion is close to that with ϵ r = 10 -6 ). However, for the 10 6 M ⊙ GMC, the ϵ r = 0 . 01 model is beginning to more into a 'strong' feedback limit as its CDF of mass accretion significantly deviates from the fiducial case. As a result, the GMCmetallicity becomes important for that particular case, but in most cases (where we are not very close to this boundary) its e ff ects are weak.", '3.2.3. Mechanical feedback (jets and winds)': "Mechanical feedback is also presented in the second rows of Figs. 8 and 9. For models with fixed mass loading factoring η w = 1, we again see the trend that the accretion onto black holes is suppressed once the energy e ffi ciency is su ffi ciently high, this happens at ϵ w ≳ 10 -5 for both complexes. These 'transition points' happen at lower energy than in the radiative feedback simulations. In terms of f Edd, we find the BHs can not reach the maximum value once ϵ w ≈ 1. All these simulations assume the  \nradiative feedback e ffi ciency following Madau et al. (2014), we find the bolometric luminosity decreases when the mechanical feedback is stronger, as expected.  \nThe model from Hu et al. (2022) has variable η w and ϵ w. Although its energy e ffi ciency is 10 -3 - 10 -2 depending on the accretion rate, the accretion rate is highly suppressed: the maximum accretion ( ∆ M max BH ) is even lower than the strong-feedback case with ϵ w = 0 . 1. Despite that the energy e ffi ciency is relatively low, the model has a strong momentum outflow that scales like ˙ P ≥ 0 . 01 f Edd ˙ M Edd c (see Fig. 2), resulting very strong BH feedback. Moreover, the fraction of mass flow reaching the BH is low (1 / f Edd in the hyper-Eddington phase), making the BH hard to grow. Instead, a much more significant fraction of the material is ejected as mechanical outflow.  \nAnother important quantity in the mechanical feedback model is the jet loading factor η w, which determines the fraction of mass flow that goes into the BH ( f acc = 1 / (1 + η w)). We present the result of this experiment in Fig. 11, where we vary η w among 1, 10, and 100, corresponding to f acc equals to 0.5, 0.1, and 0.01, while for each variation we keep ϵ w fixed at the fiducial value of each cloud complex (i.e., change the wind velocity v w). For both cloud complexes, we find that the fraction of BHs with mass accretion (with the mass resolution of gas) is almost unchanged despite di ff erent η w, while the maximum accreted mass ( ∆ M max BH ) di ff ers. As a rough estimate from constantϵ w models in Fig. 11, ∆ M max BH ∝ 1 / √ η w ≈ p f acc. We find that this is in agreement with an analytic scaling that will be covered in §4.  \n/circledot  \nFig. 10. The impact of the initial metallicity on BH accretion, which is made from a series of simulations varying the initial metallicity (10 -6 Z ⊙ -Z ⊙ ) and BH radiative-feedback models (fiducial or ϵ r = 0 . 01).  \n<!-- image -->", '3.2.4. Cosmic rays': "Fig. 12 shows gas morphology and CR energy density for selected simulations with high CR energy e ffi ciency (a 10 6 M ⊙ GMC with ϵ CR = 10 -3 and a 10 8 M ⊙ GMC with ϵ CR = 10 -2 ), and zooms in toward a rapidly-accreting BH. As also suggested in Fig. 3, the BH is located at a dense clump in the GMC. Due to strong CR feedback, there is also strong outflow near the BH. Due to CR feedback, there is a high-energy CR 'bubble' near the BH, whose energy density is comparable with CRs generated through stellar feedback. However, these bubbles are of the size of ∼ 0 . 25 pc to ∼ 2 . 5 pc, significantly smaller than the scale of stellar CRs (typically ∼ R cl). As a result, there is a higher energy (pressure) gradient from the high-energy CR bubbles near the BH, which may then expel gaseous material from the BH or even disrupt the whole GMC.  \nFrom Figs. 8 and 9 we find that BH CR feedback can be important for BH accretion in these dense GMCs we simulated if the energy e ffi ciency is high. For example, BH accretion is significantly suppressed when ϵ CR ≳ 10 -3 for the 10 6 M ⊙ GMC, or ϵ CR ≳ 10 -1 for the 10 8 M ⊙ GMC. We also see a similar trend in f Edd and L bol. Despite that CR feedback from SNe is also included in this set of simulations, we find no clear di ff erence when comparing the ϵ CR = 10 -6 run and the run with fiducial setups (without SNe CR feedback), which means this particular stellar feedback mechanism does not bring new e ff ects to BH accretion / feedback.", '4.1. BH Feedback and BH accretion': 'For all the feedback mechanisms simulated, we see (as expected) that su ffi ciently strong BH feedback suppresses BH accretion. To further quantify the e ff ect of feedback, we show the relationship between the feedback strength (indicated by energy e ffi ciencies ϵ r, ϵ w, and ϵ CR) and the maximum of the accreted mass of the BH population in Fig. 13. Besides the fixed-value models, we also show some specific sub-grid models for reference, in which energy e ffi ciency is estimated at the hyper-Eddington regime (at f Edd = 1000).  \nFor all mechanisms, there is a plateau at the weak-feedback (low energy-e ffi ciency) end, which simply represents the regime where that particular feedback mechanism does not regulate BH growth significantly. However, when feedback is stronger,  \nFig. 11. The e ff ect of jet / wind mass loading factor ( η w) variations, which is made from a series of simulations with constant η w = 1 , 10 , 100 (but with ϵ w fixed at the fiducial value), and the variable model in Hu et al. (2022).  \n<!-- image -->  \n/circledot  \nthere is a drop in ∆ M max BH . We find that this feedback-regulated limit can be approximated with simple analytic arguments at the order-of-magnitude level.', '4.1.1. The feedback-regulated limit: scalings': "We consider a simple self-regulated model akin to e.g., Silk & Rees (1998) or Murray et al. (2005); Ostriker et al. (2010). First, consider just radiative and mechanical feedback, since the fiducial model neglects CRs. If accretion is momentumlimited, then BH growth will cease when the 'feedback force' F feedback = ˙ p feedback (momentum injection rate due to feedback) equals the gravitational force on gas in the parent complex, F feedback ∼ F grav ∼ GM 2 cl / R 2 cl .  \nFor radiative feedback, the radiation pressure is dominated by single-scattering (with at most O (1) correction from IR multiple scattering) even in our dense clouds, so ˙ p rad ∼ L / c ∼ ϵ r ˙ M BH c . Similarly, for mechanical feedback, ˙ p mech ∼ ˙ Mwv w ∼ p 2 η w ϵ w ˙ M BH c .  \nFor CRs, the momentum injection is somewhat less straightforward because it arises via partial confinement / scattering of CRs generating a CR pressure gradient which accelerates the gas. From the simple analytic spherically-symmetric CRpressure-driven wind solutions in Ji et al. (2020); Hopkins et al. (2022), we can approximate this around a dominant point source as ˙ p CR ∼ R d 3 x ˆ r · ∇ P CR ∼ ℓ 2 P CR ∼ ˙ E cr / ¯ v stream , e ff , where ˙ E CR ∼ ϵ CR ˙ M BH c 2 and ¯ v stream , e ff is some e ff ective bulk CR streaming speed 3 ( ∼ κ/ 3 ℓ in terms of some characteristic CR gradient scale-length ℓ , in the di ff usion-dominated regime). So ˙ p CR ∼ ( ϵ CR α c / ¯ v stream , e ff ) ˙ M BH c where α collects the order-unity uncertainties in the scalings above. Inserting our assumed κ and assuming ℓ ∼ R cl, or checking the CR properties in the simulations directly, we see the CRs are not strongly confined, so the prefactor α c / ¯ v stream , e ff ∼ 10, roughly.  \nSo taking ˙ p feedback = ˙ p rad + ˙ p mech + ˙ p CR, we find ˙ p feedback ∼ ( ϵ r + p 2 η w ϵ w + 10 ϵ CR) ˙ M BH c . Equating this feedback force to the self-gravity of the cloud complex will thus determine the maximum ˙ M BH possible. From Fig. 4, we see the BH growth is dominated by phases with a high accretion rate, which typically reaches f Edd ∼ 1000. The peak accretion rate can then be pa-  \n5  \nFig. 12. Visualization of the GMCs at strong cosmic-ray feedback, each from a snapshot with rapid BH growth. Left two zoom-in plots : the lineaveraged mass density ( top ) and CR energy density ( bottom ) near a BH undergoing runaway accretion, based on a 10 6 M ⊙ GMC. We zoom in near the BH, from the left to the right panel, similar to that in Fig. 3. The line-averaged density is evaluated within the box of each panel. We see a high-energy CR bubble generated by the BH. Right two zoom-in plots : the same as the left plots, but embedded in the 10 8 M ⊙ GMC.  \n<!-- image -->  \nrameterized as ˙ M BH ∼ f Edd M BH / t Sal. Note that ∆ M max BH ∼ M max BH is reached for some BHs with the most significant mass growth, we have  \n∆ M max BH ∼ M max BH ∼ GM 2 cl R 2 cl t Sal f Edd c 1 ϵ r + p 2 η w ϵ w + 10 ϵ CR (4) ∼ 3 × 10 4 M ⊙ M cl 10 6 M ⊙ ! 2 R cl 5 pc ! -2 1 f Edd( ϵ r + p 2 η w ϵ w + 10 ϵ CR) .  \nWith Eq. (4), we plot the scaling relations for each feedback mechanism variation in Fig. 13. First, we check models varying radiative and mechanical feedback e ffi ciency, i.e., ϵ CR = 0. As described previously, for the fiducial mechanical feedback model, we input ϵ w = 10 -6 for the 10 6 M ⊙ complex and ϵ w = 10 -3 for the 10 8 M ⊙ one; for the fiducial radiative feedback, we input ϵ r = 3 × 10 -3 (c.f. Fig. 2) since f Edd = 1000 is reached. Varying the feedback e ffi ciency, we find both scaling relations for the radiative and mechanical feedback are in good agreement with the simulations.  \nConsidering CRs, we also find the scaling relation agrees quite well with the simulations, though from our analytic scaling, we see that this is more sensitive to other uncertainties besides ϵ CR, namely CR transport physics. It remains deeply theoretically uncertain how well-coupled GeV CRs are in dense, neutralphase ISM gas (see e.g. Zweibel 2017; Krumholz et al. 2020; Hopkins et al. 2021), and the CR transport model we adopted is purely phenomenological, so the details of this scaling should be taken with more caution. Nonetheless, it provides some order-ofmagnitude guidance.", '4.2. BH feedback and star formation': "Fig. 14 shows the star formation e ffi ciency (SFE, here defined as the ratio between the stellar mass at the end of the simulation and the initial total gas mass) for di ff erent kinds of feedback mechanisms and various energy e ffi ciencies. When the energy e ffi ciency is su ffi ciently low, we find that in both the 10 6 M ⊙ and 10 8 M ⊙ GMCs, the SFE is ∼ 0 . 3. This is expected since the two GMCs have the same initial mean surface density (Grudi'c et al. 2018, and more citations in § 1). Then once the energy e ffi ciency increase beyond ∼ 10 -3 , the SFE for the 10 6 M ⊙ GMCdrops below 0 . 1. However, the impact on the 10 8 M ⊙ is not significant.  \nFig. 13. Dependence of BH accretion on the energy e ffi ciency of specific feedback mechanisms for the 10 6 M ⊙ ( top panel ) and 10 8 M ⊙ ( bottom panel ) GMCs. We summarize experiments with fixed-value energy e ffi ciencies ( solid ) and show the scaling relation ( dashed ) based on Eq. (4). Several more sophisticated sub-grid models are also displayed for reference, for each of them the energy e ffi ciency is calculated at the hyper-Eddington regime of f Edd = 1000.  \n<!-- image -->  \nA very special case is the ϵ CR = 0 . 1 run, in which the SFE drops to a very low value ( < 10 -3 ). With further inspection of the run, we find that the cloud was disrupted by strong CR winds at ∼ 0 . 3 t ff , which is much shorter than the characteristic time scale of star formation at the low-BH-feedback limit (from simulations it is t SF 1 / 2 ∼ 1 . 25 t ff , defined as the time when the SFE reaches half of the final value). As a result, the SFE almost freezes at the (extremely low ) star-formation level at the time of GMC's disruption.  \nWe may roughly explain this quantitatively. From our simulations, a typical BH growth history is dominated by a few 'bursty' super-Eddington accretions which happen in a short period of time (typically less than 0 . 05 t ff ). At a time scale much longer than the bursty growth, the BH mass is then anticipated to be M BH ∼ M ini BH exp( ⟨ f Edd ⟩ t / t Sal). Applying this to the BH with the most significant mass growth, we find the disruption time of the GMC due to strong BH feedback to be  \nt d ∼ t Sal ⟨ f Edd ⟩ ln       1 + ∆ M max BH M ini BH       . (5)  \nFor each simulation, M max BH and ⟨ f Edd ⟩ are presented in the middle columns of Figs. 5, 8, and 9; and M max BH (or ∆ M max BH ) is also analytically evaluated from Eq. (4). In general, we find that ⟨ f Edd ⟩ is approximately 10 - 100, and there is M max BH ≫ M ini BH so t d is generally much larger than the free-fall time of the two GMCs we simulated. However, when there is very strong BH feedback, ∆ M max BH is low. Thus there are chances that ∆ M max BH ≪ M ini BH and t d ∼ t Sal ∆ M max BH / ( ⟨ f Edd ⟩ M ini BH ) is smaller than (or comparable with) the free-fall time. For the ϵ CR = 0 . 1 run of the 10 6 M ⊙ complex, we measure and find M ini BH = 702 M ⊙ and ∆ M max BH = 740 M ⊙ (thus ∆ M max BH / M ini BH ≈ 0 . 05), resulting t d ∼ 2 / ⟨ f Edd ⟩ Myr. The time scale is comparable to (or substantially smaller than if ⟨ f Edd ⟩ is large) the free-fall time of the GMC (0 . 19 Myr).  \nWe also note that the argument is dependent on the initial mass of the BH, thus t d is essentially short only when the 'lucky' seed BH is initially massive. As a result, the extreme suppression of SFE does not always happen even if the energy e ffi ciency is high, which is true for most experiments in Fig. 14 with high energy e ffi ciency.  \nIn summary, we find that BH feedback will not a ff ect the star formation e ffi ciency unless it is in a highly feedbackdominated / regulated regime and relatively lower-mass cloud complex with high mass 'active' seeds ( M BH ≳ 10 3 M ⊙ with ϵ ≳ 10 -3 for the 10 6 M ⊙ GMC). The significance of the impact is largely decided by how soon the winds induced by BH feedback 'interrupt' star formation (which happens at the free-fall time scale) - in some particular models (e.g., very strong feedback from a massive BH seed), even if BH accretion is limited, the strong BH feedback is still able to disrupt the whole GMC at the very early stage of the star-formation history and result in very low star formation e ffi ciency.", '4.3. Caveats and outlook': "Because of the huge mismatch in scales of the BH accretion / feedback problem (from kpc / pc scales to the Schwarzschild radius), it is di ffi cult to simulate every aspect in great detail. In this study, we concentrate on larger scales, where the gas inflow towards BHs' gravitational capture radii is resolved but dynamics of small scales (gas cells that are considered 'gravitationally bounded to the BH' including, of course, the accretion disks themselves) are described by sub-grid models with a few free parameters (as listed in Tab. 2). However, the setup is su ffi -cient to answer the question we were asking: whether there could be enough fuel from large scales (kpc / pc) reaching small scales (sub-pc) to fuel rapid BH growth in the presence of BH feedback. The setup is thus also able to study the impact of di ff erent sub-grid parameters on BH accretion, as an attempt to bridge the dynamics of small and large scales.  \nStill, the sub-grid models and resolution of our simulations do not allow us to connect the large and small scales in great  \nFig. 14. The star formation e ffi ciency (SFE) in our experiments with various fixed-value BH-feedback energy e ffi ciencies. For the 10 6 M ⊙ GMC, there is (sometimes extreme) suppression of the SFE at high energy e ffi ciency.  \n<!-- image -->  \nself-consistency. For example, the energy e ffi ciencies of di ff erent feedback mechanisms are treated as arbitrary inputs which are not dependent on BH ambient boundary conditions like magnetic field and angular momentum flow. There is thus plenty of space for improvement.  \nRecent zoom-in simulations of BH accretion or feedback may shed light on this problem (e.g., Talbot et al. 2021; Guo et al. 2023; Hopkins et al. 2024b,a). With the technique of superLagrangian refinement, small-scale structures (e.g., BH accretion disks or jets) and dynamics (e.g., the magnetic field in disks, jet launching mechanisms, BH-disk interactions) may also be simulated with higher level of self-consistency.  \nFinally, as also described in Paper I, next-generation simulations which resolve individual stars (Grudi'c et al. 2021a; Guszejnov et al. 2021) may enable other important processes like BH seed formation, stellar merging, and BH accretion / feedback in the same simulation.", '5. Conclusions': "This study is a parameter-space survey of seed BH accretion and feedback in star-forming GMCs, with key parameters listed in Tab. 2. We focus on high surface-density clouds which are ideal environments for runaway BH accretion in the weak BH feedback limit (Paper I), and find that BH feedback self-regualtes BH growth. Our major conclusions are listed below and illustrated in Fig. 15.  \n- 1. Even with feedback, the BH accretion scenario is the same as that presented in Paper I. Rapid BH accretion happens when BHs semi-randomly encounter dense gas clumps with low relative velocities. This happens only when the dynamics of both BH and gas are dominated by GMC's self-gravity, i.e., when the cloud has high enough surface density that stellar feedback is ine ffi cient, and the BH has low enough mass that BH feedback is (initially) weak. Other GMC properties like metallicity have weak e ff ects (unless in some models with strong radiative feedback). Also, like the condition without BH feedback (Paper I), initial information about black holes like initial mass, position, and velocity does not significantly correlate with BH accretion (as long as BHs are gravitationally bound to the GMC initially), i.e., there is significant stochasticity in the process.  \nFig. 15. Overview of seed BH accretion and feedback in star-forming GMCs. (a). Initial condition, where the seed BH is inside a pre-starformation, turbulent GMC. (b). If the self-gravity of the GMC (characterized by its mean surface density), star formation and feedback dominate and disrupt the GMC, suppressing BH accretion. (c). If the GMC's self-gravity is strong, gas is kept bound, creating a proper environment for BH accretion and feedback. (d). If the BH feedback is weak, there is steady accretion flow towards the BH if it is at the center of the potential well. (e). (f). (g). Strong feedback in di ff erent forms may suppress both BH accretion and star formation, this process is dominated by momentum outflow from BHs and can be approximated with the analytic argument as in Eq. (4).  \n<!-- image -->  \n- 2. With our fiducial feedback model ('log-form' radiative efficiency and 'critical' mechanical feedback e ffi ciency, we find there is suppression of BH accretion in both simulations with 10 6 M ⊙ and 10 8 M ⊙ GMCs. Significant BH accretion is more likely in the 10 8 M ⊙ GMC as expected, where BHs grow up to ∼ 5 × 10 4 M ⊙ , even though we assume relatively e ffi cient radiative ⟨ ϵ r ⟩ ∼ 0 . 003 and ⟨ ϵ w ⟩ ∼ 0 . 001 (so ⟨ ϵ w ⟩ / ⟨ ϵ r ⟩ ∼ 1 / 3). There are short phases of hyper-Eddington accretion ( f Edd ∼ 1000) for both simulations.\n- 3. For feedback physics involved in this study, we find that the maximum possible accreted mass for BHs ∆ M max BH declines as the feedback e ffi ciency increases (Fig. 13), in a manner consistent with momentum-regulated BH growth. The simulations agree well with a simple analytic scaling summarized in Eq. (4). Depending on the cloud's properties and the feedback e ffi ciency, the scenario can produce 10 4 M ⊙ IMBHs quickly in a few free-fall times (typically ∼ Myr). These IMBHs can be the 'massive seeds' of subsequent galactic SMBHs.\n- 4. For a given feedback e ffi ciency model, the amount BHs can grow increases both with the surface density / acceleration scale of the parent complex, and the gas mass available. So dense, massive environments such as high-redshift galaxy nuclei will (unsurprisingly) be most favorable to rapid IMBH and even SMBH-level growth. Owing to the di ff erent physical mechanisms by which energy is converted into momentum, the same 'energetic' e ffi ciency ϵ for di ff erent physical  \nmechanisms (radiation, mechanical outflows, cosmic rays) will not produce the same level of self-regulation (Eq. 4). Physically-motivated models of a variable radiative e ffi -ciency ϵ r (e.g. Watarai et al. 2000; Madau et al. 2014) tend to predict relatively low ϵ r at high Eddington ratios, so it seems unlikely that this is actually the dominant self-regulation mechanism on scales like what we resolve here. In contrast, some recent models for broad-opening angle, non-relativistic winds from super-Eddington disks (Hu et al. 2022) have argued for this producing e ffi ciently momentum-loaded winds with large ( η w ≫ 1) mass-loading, which would expel most of the gravitationally captured mass before it can reach the BH, and so in our simulations here much more e ffi ciently suppress BH growth.  \n- 5. Even when BH accretion overall is suppressed by BH feedback with high e ffi ciency, there are still short phases of very high accretion rate ( ≳ 100 ˙ M Edd) for a subset of BHs in the simulation ( ∼ 10% for the 10 6 M ⊙ GMC and ∼ 1% for the 10 8 M ⊙ one).\n- 6. When varying the BH feedback e ffi ciency, we find that the massive complex's global star formation e ffi ciency is generally not suppressed by BH feedback, unless the BH feedback e ffi ciency is very high ( ≳ 0 . 01) such that the could is disrupted within the first free-fall time.\n- 7. Other parameters are related to the sub-grid (unresolved) mass flow to the BH, like the e ff ective α , the upper limit of BH accretion rate, and the jet mass loading factor η w are largely degenerate, but do a ff ect the BH mass growth (as expected by definition in our models).  \nSome of these findings are similar to other recent works. For example, Lupi et al. (2016) simulated the super-critical accretion onto stellar-mass BH seeds in gaseous circum-nuclear disks (CNDs), with constant ϵ r = 0 . 1 or the slim disk model (Madau et al. 2014) for radiative e ffi ciency. Pezzulli et al. (2016) studied the growth of the central BH (initially with ∼ 100 M ⊙ seed mass) of high redshift quasars with a semi-analytic framework including the slim-disk model and AGN wind feedback e ffi ciency of ∼ 10 -3 (but not CRs). Our results are similar to theirs for simulations with similar feedback parameters.  \nExtension of this work may include more sophisticated subgrid models of BH accretion, which connect the BH accretion with the properties of the ambient gas. Another direction is to develop zoom-in simulations that attain higher resolution near the BH seed, which is more self-consistent but computationally more expensive. We leave these possibilities for future work.  \nAcknowledgements. Support for the authors was provided by NSF Research Grants 1911233, 20009234, 2108318, NSF CAREER grant 1455342, NASA grants 80NSSC18K0562, HST-AR15800. Numerical calculations were run on the Caltech compute cluster 'Wheeler,' allocations AST21010 and AST20016 supported by the NSF and TACC, and NASA HEC SMD-16-7592. A public version of the GIZMO code is available at http://www.tapir.caltech.edu/ ~phopkins/Site/GIZMO.html . YS acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC), [funding reference number 568580].", 'References': "Bañados, E., Venemans, B. P., Mazzucchelli, C., et al. 2018, Nature, 553, 473 Bate, M. R., Bonnell, I. A., & Price, N. M. 1995, MNRAS, 277, 362 Begelman, M. C. 1979, MNRAS, 187, 237 Blandford, R. D. & Begelman, M. C. 1999, MNRAS, 303, L1 Blandford, R. D. & Begelman, M. C. 2004, MNRAS, 349, 68 Bondi, H. 1952, MNRAS, 112, 195 Bower, R. G., Schaye, J., Frenk, C. S., et al. 2017, MNRAS, 465, 32  \nBromm, V. & Loeb, A. 2003, ApJ, 596, 34  \nChevance, M., Krumholz, M. R., McLeod, A. F., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII,  \ned. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 1 Di Matteo, T., Springel, V., & Hernquist, L. 2005, Nature, 433, 604 Dubois, Y., Pichon, C., Devriendt, J., et al. 2013, MNRAS, 428, 2885 Dubois, Y., Volonteri, M., Silk, J., et al. 2015, MNRAS, 452, 1502 Fabian, A. C. 2012, ARA&A, 50, 455  \nFukushima, H., Yajima, H., Sugimura, K., et al. 2020, MNRAS, 497, 3830 Grudi'c, M. Y., Guszejnov, D., Hopkins, P. F., et al. 2018, MNRAS, 481, 688 Grudi'c, M. Y., Guszejnov, D., Hopkins, P. F., O ff ner, S. S. R., & FaucherGiguère, C.-A. 2021a, MNRAS, 506, 2199  \nGrudi'c, M. Y., Kruijssen, J. M. D., Faucher-Giguère, C.-A., et al. 2021b, MNRAS, 506, 3239  \nGuo, F. & Mathews, W. G. 2012, ApJ, 756, 181  \nGuo, F. & Oh, S. P. 2008, MNRAS, 384, 251  \nGuo, M., Stone, J. M., Kim, C.-G., & Quataert, E. 2023, ApJ, 946, 26  \nGuszejnov, D., Grudi'c, M. Y., Hopkins, P. F., O  \nff  \nner, S. S. R., & Faucher-  \nGiguère, C.-A. 2021, MNRAS, 502, 3646  \nHabouzit, M., Volonteri, M., & Dubois, Y. 2017, MNRAS, 468, 3935  \nHe, C.-C., Ricotti, M., & Geen, S. 2019, MNRAS, 489, 1880  \nHopkins, P. F. 2015, MNRAS, 450, 53 Hopkins, P. F. 2016, MNRAS, 462, 576  \n- Hopkins, P. F., Chan, T. K., Ji, S., et al. 2021, MNRAS, 501, 3640  \nHopkins, P. F., Grudi'c, M. Y., Wetzel, A., et al. 2020, MNRAS, 491, 3702  \n- Hopkins, P. F. & Raives, M. J. 2016, MNRAS, 455, 51\n- Hopkins, P. F., Squire, J., & Butsky, I. S. 2022, MNRAS, 509, 3779  \nHopkins, P. F., Squire, J., Quataert, E., et al. 2024a, The Open Journal of Astrophysics, 7, 20  \nHopkins, P. F., Squire, J., Su, K.-Y., et al. 2024b, The Open Journal of Astrophysics, 7, 19  \nHopkins, P. F., Torrey, P., Faucher-Giguère, C.-A., Quataert, E., & Murray, N. 2016, MNRAS, 458, 816  \nHopkins, P. F., Wetzel, A., Kereš, D., et al. 2018a, MNRAS, 477, 1578  \nHopkins, P. F., Wetzel, A., Kereš, D., et al. 2018b, MNRAS, 480, 800  \nHopkins, P. F., Wetzel, A., Wheeler, C., et al. 2023, MNRAS, 519, 3154  \nHoyle, F. & Lyttleton, R. A. 1939, Proceedings of the Cambridge Philosophical  \nSociety, 35, 405 Hu, H., Inayoshi, K., Haiman, Z., Quataert, E., & Kuiper, R. 2022, ApJ, 934, 132 Inayoshi, K., Haiman, Z., & Ostriker, J. P. 2016, MNRAS, 459, 3738 Inayoshi, K., Visbal, E., & Haiman, Z. 2020, ARA&A, 58, 27 Ishibashi, W. & Fabian, A. C. 2023, MNRAS, 519, 1931 Ji, S., Chan, T. K., Hummels, C. B., et al. 2020, MNRAS, 496, 4221 Jiang, Y.-F., Stone, J. M., & Davis, S. W. 2014, ApJ, 796, 106 Jiang, Y.-F., Stone, J. M., & Davis, S. W. 2019, ApJ, 880, 67 Kim, J.-G., Kim, W.-T., & Ostriker, E. C. 2018, ApJ, 859, 68 Kim, J.-G., Ostriker, E. C., & Filippova, N. 2021, ApJ, 911, 128 Klessen, R. S. 2000, ApJ, 535, 869 Kremer, K., Spera, M., Becker, D., et al. 2020, ApJ, 903, 45 Krumholz, M. R., Crocker, R. M., Xu, S., et al. 2020, MNRAS, 493, 2817 Larson, R. B. & Starrfield, S. 1971, A&A, 13, 190  \nLatif, M. A., Whalen, D. J., Khochfar, S., Herrington, N. P., & Woods, T. E. 2022, Nature, 607, 48  \nLupi, A., Haardt, F., Dotti, M., et al. 2016, MNRAS, 456, 2993  \nLupi, A., Quadri, G., Volonteri, M., Colpi, M., & Regan, J. A. 2024, A&A, 686,  \nA256 Mac Low, M.-M. & Klessen, R. S. 2004, Reviews of Modern Physics, 76, 125 Madau, P., Haardt, F., & Dotti, M. 2014, ApJ, 784, L38 Madau, P. & Rees, M. J. 2001, ApJ, 551, L27  \nMassonneau, W., Volonteri, M., Dubois, Y., & Beckmann, R. S. 2023, A&A,  \n670, A180 Ostriker, J. P., Choi, E., Ciotti, L., Novak, G. S., & Proga, D. 2010, ApJ, 722,  \nMcKee, C. F. & Ostriker, E. C. 2007, ARA&A, 45, 565 Murray, N., Quataert, E., & Thompson, T. A. 2005, ApJ, 618, 569 Ohsuga, K., Mori, M., Nakamoto, T., & Mineshige, S. 2005, ApJ, 628, 368 642  \nPezzulli, E., Valiante, R., & Schneider, R. 2016, MNRAS, 458, 3047 Portegies Zwart, S. F., Baumgardt, H., Hut, P., Makino, J., & McMillan, S. L. W.  \n2004, Nature, 428, 724 Quataert, E. & Gruzinov, A. 2000, ApJ, 539, 809 Regan, J. A., Downes, T. P., Volonteri, M., et al. 2019, MNRAS, 486, 3892 Ryu, T., Tanaka, T. L., Perna, R., & Haiman, Z. 2016, MNRAS, 460, 4122 Shakura, N. I. & Sunyaev, R. A. 1973, A&A, 24, 337 Shen, X., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2020, MNRAS, 495, 3252  \nShi, Y., Grudi'c, M. Y., & Hopkins, P. F. 2021, MNRAS, 505, 2753 Shi, Y., Kremer, K., Grudi'c, M. Y., Gerling-Dunsmore, H. J., & Hopkins, P. F. 2023, MNRAS, 518, 3606  \nSijacki, D., Pfrommer, C., Springel, V., & Enßlin, T. A. 2008, MNRAS, 387, 1403  \nSilk, J. & Rees, M. J. 1998, A&A, 331, L1  \nS˛adowski, A. 2009, ApJS, 183, 171  \nS˛adowski, A., Lasota, J.-P., Abramowicz, M. A., & Narayan, R. 2016, MNRAS, 456, 3915  \nS˛adowski, A., Narayan, R., Tchekhovskoy, A., et al. 2015, MNRAS, 447, 49 Su, K.-Y., Hopkins, P. F., Bryan, G. L., et al. 2021, MNRAS, 507, 175 Talbot, R. Y., Bourne, M. A., & Sijacki, D. 2021, MNRAS, 504, 3619  \nTorrey, P., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2020, MNRAS, 497, 5292  \nVolonteri, M., Habouzit, M., & Colpi, M. 2021, Nature Reviews Physics, 3, 732 Wang, F., Yang, J., Fan, X., et al. 2021, ApJ, 907, L1 Watarai, K.-y., Fukue, J., Takeuchi, M., & Mineshige, S. 2000, PASJ, 52, 133 Yang, J., Wang, F., Fan, X., et al. 2020, ApJ, 897, L14 Yu, Q. & Tremaine, S. 2002, MNRAS, 335, 965  \nZweibel, E. G. 2017, Physics of Plasmas, 24, 055402"}