indicated by the Soltan ( 1982 ) argument. During BH accretion, the infalling matter forms an accretion disc around the BH and loses angular momentum and energy via viscous dissipation. A fraction of this energy is thermalized via collisional processes resulting in the production of substantial radiation. During high accretion this will result in very high luminosity (up to ∼10 48 erg s -1 ) radiation emanating from a very compact region of milliparsec size in the centre of the galaxy. These objects, called Active Galactic Nuclei (AGNs), are some of the most powerful sources of radiation in the Universe. Injection of energy and momentum feedback from AGN to the surrounding interstellar medium known as 'AGN feedback' is thought to be responsible for regulating star formation in massive galaxies (Benson et al. 2003 ;Di Matteo, Springel & Hernquist 2005 ;Dubois et al. 2010 ;Bourne & Sijacki 2017 ). It is also postulated to be the reason behind the observed correlations between the central SMBH mass and host galaxy properties (Sijacki et al. 2007 ;Di Matteo et al. 2008 ). There is strong observational evidence for AGN feedback in the form of galactic outflows with velocities of the order of ∼1000 km s -1 and outflow rates (of the order of ∼1000 M yr -1 ) much larger than the star formation rates (Rupke & Veilleux 2011 ;Sturm et al. 2011 ). Such high velocities and outflow rates are difficult to achieve through stellar feedback alone (Fabian 2012 ). Furthermore, the strongest outflow velocities are found in the central regions of the galaxies (Rupke & Veilleux 2011 ). These clearly indicate the presence of AGN activity.

In the past few decades there has been much effort to model galaxy formation and evolution with numerical simulations (see Vogelsberger et al. 2020 , for a technical re vie w). Numerous recent cosmological hydrodynamics simulations have been able to produce observationally consistent results (e.g. Genel et al. 2014 ;Vogelsberger et al. 2014a , b ;Schaye et al. 2015 ;McAlpine et al. 2017 ;Nelson et al. 2018 ;Pillepich et al. 2018b ;Dav é et al. 2019 ) and also identify AGN feedback as the main mechanism for regulating star formation in massive galaxies (Teyssier et al. 2011 ;Dubois et al. 2013 ;Sijacki et al. 2015 ;Weinberger et al. 2018 ). One of the major challenges in modelling BH accretion and feedback is the large range of scales involved. Fueling SMBHs requires large inflows of gas from galactic scales to much smaller length scales such as the Bondi radius (of the order of ∼ pc). Due to resolution limitations, galactic scale simulations cannot directly follow the formation of accretion discs around BHs. Instead, BH accretion and feedback are modelled using subgrid (or subresolution) prescriptions, which are designed to capture large-scale effects of unresolved dynamical processes (e.g. Sijacki et al. 2015 ). BHs are usually treated as sink particles that swallow nearby gas according to accretion rates calculated based on gas properties at tens of pc-kpc scales. These scales are larger than the radius of influence of BHs (distance out to which the gravitational potential of the black hole dominates the gravitational potential of the host galaxy, which is typically around a few pc to tens of pc) and hence the accretion rates are being estimated without resolving the important dynamics that ultimately delivers gas to the black hole, potentially leading to inaccuracies.

The subgrid models of the interstellar medium (ISM) can also affect the accuracy of BH accretion rates. Many hydrodynamic simulations of galaxy evolution use an effective equation of state (eEOS) approach to describe the ISM (Springel & Hernquist 2003 ;Vogelsberger et al. 2013 ; see Section 2.1 for a brief description). In the model we adopt, gas cooling, star formation, and stellar feedback occurring on unresolved scales are approximated via a two-phase medium in pressure equilibrium. Because these models generally do not consider gas cooling and collapse beyond a certain limit (typically ∼10 4 K), they produce an overly smooth gas distribution. Apart from this, the geometrical structure of the ISM is thought to play an important role in determining how the AGN feedback couples to the ISM (Wagner, Bicknell & Umemura 2012 ;Wagner, Umemura & Bicknell 2013 ;Bieri et al. 2017 ;Torrey et al. 2020 ). A dense and uniform ISM will allow the feedback to do more work on the environment compared to a porous ISM (Faucher-Gigu ère & Quataert 2012 ; Torrey et al. 2020 ). The geometric structure of the ISM is mostly shaped by stellar feedback and hence depends on the subgrid models of stellar evolution. Nuclear scale simulations with a multiphase ISM that resolve sub parsec scale injection of AGN feedback also show that AGNs can launch po werful gas outflo ws which strongly suppress nuclear star formation and BH growth (Hopkins et al. 2016 ). Thus, it is important to study BH accretion and feedback using an ISM model with a well-resolved multiphase structure and local injection of stellar feedback.

There exist several models in the literature that treat the ISM explicitly. Some of the recent ones are the Feedback In Realistic Environments (FIRE; Hopkins et al. 2018 ) in the mesh-free magnetohydrodynamics code GIZMO (Hopkins 2015 ), the Agertz et al. ( 2013 ) treatment of the ISM in the Eulerian adaptive mesh refinement code RAMSES (Teyssier 2002 ), and Stars and MUltiphase Gas in GaLaxiEs (SMUGGLE; Marinacci et al. 2019 ) in the moving-mesh magnetohydrodynamics code AREPO (Springel 2010 ). Cosmological simulations using FIRE have been able to reproduce many observed correlations of galaxy properties (e.g. Hopkins et al. 2014 ;Ma et al. 2016 ;Wetzel et al. 2016 ;Muratov et al. 2017 ;Sparre et al. 2017 ;Orr et al. 2018 ). Idealized galaxy simulations using the SMUGGLE model has been able to produce well-resolved multiphase ISM with observationally consistent star formation rates and galactic outflows (Marinacci et al. 2019 ), realistic star cluster properties (Li et al. 2020(Li et al. , 2022 ) ), and Hydrogen emission brightness profiles (Smith et al. 2022 ;Tacchella et al. 2022 ) and constant density cores in dwarf galaxies (Jahn et al. 2021 ).

Recently the FIRE model has also been used to study SMBH physics. Angl és-Alc ázar et al. ( 2017b ) used a gravitational torque based accretion prescription (Hopkins & Quataert 2011 ) to model BH growth in cosmological zoom simulations and showed that bursty stellar feedback limits the fueling of high redshift AGNs. Angl és-Alc ázar et al. ( 2021) implemented a hyper-Lagrangian refinement scheme in these simulations with which they were able to resolve the sub-parsec scale gas inflows on to the SMBH for short ∼10 Myr timescales. These and related studies ( C ¸atmabacak et al. 2022 ;Tillman et al. 2022 ) of SMBHs using the FIRE model have also shown that stellar feedback and galaxy mergers play an important role in shaping the BH-galaxy scaling relations. The FIRE model was also used to simulate AGN feedback to study its impact on host galaxies and the intracluster medium (Torrey et al. 2020 ;Su et al. 2021 ;Wellons et al. 2022 ).

Observational and theoretical studies have improved our understanding of black hole accretion and AGN physics. Nevertheless there are still many open questions in this field. These are some of the questions we are trying to address through this work: How does black hole fueling depend on the surrounding ISM conditions such as the availability of gas and strength of stellar feedback? What are the time-scales and variability of AGN activity? To what extent do galaxy mergers trigger AGN activity?

In this paper, we study BH fueling by simulating idealized galaxies with AREPO , using SMUGGLE to model an explicit, multiphase ISM. We investigate the nature of BH fueling in a diverse range of galactic environments and o v er long (few Gyr) time-scales. We combine our high resolution simulations with a super-Lagrangian refinement scheme which impro v es the gas mass resolution in the immediate neighbourhood of the BHs. This feature allows us to resolve gas dynamics near the BH on scales closer to its radius of influence. We highlight the key differences between BH accretion in the explicit ISM versus an ef fecti ve equation of state ISM. After establishing the resolution convergence of our model we look at BH growth in idealized galaxy mergers. Galaxy mergers are thought to be important triggers of some AGN activity, because gravitational torques during merger can funnel gas to the central regions of the galaxies (e.g. Sanders et al. 1988 ;Hernquist 1989 ;Barnes & Hernquist 1991, 1996 ;Mihos & Hernquist 1996 ;Springel, Di Matteo & Hernquist 2005 ;Hopkins et al. 2009 ;Cotini et al. 2013 ;Weston et al. 2016 ;Blumenthal & Barnes 2018 ;Ellison et al. 2019 ). The role of mergers in producing AGNs is still a debated topic, ho we ver. Some observ ations sho w that there is no strong correlation between galaxy mergers and AGN activity

Downloaded from https://academic.oup.com/mnras/article/517/4/4752/6747143 by Oregon State University user on 23 March 2023 MNRAS 517, 4752-4767 (2022) (Grogin et al. 2005 ;Pierce et al. 2007 ;Kocevski et al. 2011 ;Villforth et al. 2019 ). Thus, it is important to better understand mechanisms for BH fueling in both isolated and merging environments.

This paper is organized as follows. In Section 2 we describe the numerical methods including the subgrid models used. In Section 3 we present the results of our simulations and analyse them. In Section 4 we discuss the implications of our results.

Our simulations use AREPO , which is a moving-mesh magnetohydrodynamics code (Springel 2010 ;Pakmor et al. 2016 ). AREPO uses a finite volume method to solve hydrodynamics equations on an unstructured mesh defined by a Voronoi tessellation of discrete points. The mesh points are free to mo v e with the fluid flo w allo wing a continuous and automatic adjustment of spatial resolution (for details see Springel 2010 ). Simulations of galaxy formation and evolution involve both dark matter and baryonic matter. In AREPO , dark matter and star particles are modelled using the collisionless Boltzmann equation coupled to Poisson's equation, while the gas is modelled as an ideal gas following Euler's equations (Vogelsberger et al. 2020 ). The baryonic physics of galaxy evolution involves a large number of dynamical processes such as star formation, stellar feedback, black hole accretion, gas cooling and heating, among others, spanning a wide range of spatial and temporal scales. Because these scales are beyond the resolution limits of current simulations, these processes are implemented as subgrid models.

In this paper, we will focus on simulations using the SMUGGLE model. For comparison, we will also run some of the simulations with the eEOS model. In the following sections we briefly discuss these two models.

In ef fecti ve equation of state models of the ISM (Springel & Hernquist 2003 ), the small-scale processes go v erning the dynamics of the ISM are assumed to be in an equilibrium state where the temperature of the gas can be approximated as a function of only the gas density. These models do not explicitly resolve the small-scale multiphase structure of the ISM gas but instead model it as a two phase medium composed of cold gas clouds in pressure equilibrium with a hot gas. A relation of the form T ∝ ρ γ between the local gas temperature and density is imposed (i.e. an equation of state). Stars are allowed to form stochastically from gas cells which are abo v e (below) a given density (temperature) threshold. Due to resolution limitations, individual stars cannot be resolved and each star particle spawned from gas cells represents a stellar population. An initial mass function (IMF) is used to model stellar evolution within each star particle (for instance, the galaxy formation model of Vogelsberger et al. 2013 , used in the Illustris simulations, uses the Chabrier 2003 IMF). In the local version of stellar feedback model, the energy and mass released by the stars, calculated using the IMF, is injected into a prescribed number of neighbouring gas cells of the star particle in a kernelweighted fashion. Ho we ver in the vast majority of cosmological simulations using the eEOS model, hydrodynamically decoupled wind particles are used to deposit the energy outside the dense ISM instead of injecting it locally. These winds are also modelled to capture the galactic outflows in cosmological simulations. We do not use the non-local winds in our simulations. This eEOS model has been used by Illustris (Vogelsberger et al. 2014b ), IllustrisTNG (Weinberger et al. 2017 ;Pillepich et al. 2018a ) and Auriga (Grand et al. 2017 ) simulations.

In order to resolve the ISM structure on scales below a few hundred parsecs, which is the scale of molecular clouds, the eEOS has to be replaced with a more detailed modelling of the cold dense gas.

Current state-of-the-art galactic scale simulations have sufficient resolution to resolve gas cooling, local star formation, and feedback on much smaller scales; such simulations do not require an eEOS approach. SMUGGLE is a recently developed explicit ISM and stellar feedback model for AREPO , which is suitable for simulations with resolution of ∼10 4 M or better (Marinacci et al. 2019 ). The SMUGGLE model includes gas heating and cooling mechanisms such as low temperature atomic and molecular cooling, cosmic rays, and photoelectric heating, which can create temperatures from ∼10 K to ∼10 8 K and allow a natural development of the multiphase ISM. SMUGGLE also models star formation and stellar feedback in a localized fashion. Star particles are allowed to form only from very cold and dense gas that is also gravitationally bound. As in Marinacci et al. ( 2019 ), we adopt a density threshold of 100 cm -3 which is in the range of densities of giant molecular clouds, compared to a threshold of 0 . 13 cm -3 typically adopted in the eEOS model of Springel & Hernquist ( 2003 ). Stellar feedback inputs are calculated assuming the Chabrier IMF and include supernovae (SN) energy and momentum injection, radiative feedback from young, massive stars, and energy and momentum injection from AGB and OB winds. The number of SN explosions and the mass released are calculated by sampling from a Poisson distribution in order to mimic the discrete nature of these events. SN implementation also takes into account the momentum boost due to unresolved Sedov-Taylor expansion phases. All stellar feedback injections happen locally into the nearest neighbouring gas cells of the star particle. A cubic spline kernel is used to calculate weights and the number of neighbours is set to 32 with a tolerance of ±1. Galaxy properties simulated using SMUGGLE was shown to be resolution convergent at ∼10 4 M resolution. Details of the numerical methods and implementation can be found in Marinacci et al. ( 2019 ). Some of the key parameters used in our simulations using the eEOS and SMUGGLE models are listed in Table 1 .

Similar to the ISM physics, it is impossible to resolve the flow of gas into the central SMBH because of the extremely small spatial scales. In our simulations the subgrid model of BH accretion is based on Eddington-limited Bondi-Hoyle prescription. In this model the accretion rate is given by,

where Ṁ Bondi is the Bondi-Hoyle accretion rate,

and Ṁ Edd is the Eddington accretion rate,

Here ρ is the local gas density and c s is the sound speed near the BH. The factor G is Newton's constant, c is the speed of light in vacuum, r = 0.2 is the radiative efficiency of accretion, m p is the mass of a proton, σ T is the Thompson cross-section, and α is a dimension-less scaling parameter. Equation ( 2) is obtained by assuming spherically symmetric and steady accretion of ideal gas on to the BH. The Eddington accretion rate represents the maximum accretion rate beyond which the radiation pressure will o v ercome the gravitational pull on the gas. The gas density ρ and sound speed c s in equation ( 2) are obtained by averaging over 64 (with a tolerance of ±1) gas cells nearest to the black hole in a kernel-weighted fashion.

A cubic spline kernel is used to calculate the weight of a cell as a function of its distance from the BH. The BH accretion rates do not depend strongly on the number of neighbours used in the calculation as long as it is not too large (i.e. ࣠ 100). We will often parametrize the accretion rate in terms of the Eddington ratio, χ Edd ≡ Ṁ Bondi / Ṁ Edd .

The simulations in this study do not include AGN feedback. In the absence of feedback from BHs, the central gas densities can become many orders of magnitude higher than they would be in the presence of feedback. This issue is exacerbated by the SMUGGLE ISM model, which produces gas clouds with much higher densities than can be achieved in eEOS ISM models. In practice, we find that this results in long periods of Eddington-limited accretion that depletes the central gas reservoir in a few hundred megayear period.

In order to study BH fueling o v er longer time-scales without introducing a complicated dependence on the highly complex interplay between BH fueling and feedback, we adopt a non-standard approach: we use the parameter α to scale down the Bondi accretion rate by a factor of 10 -5 . This value was chosen based on the Eddington ratios of the unscaled Bondi accretion rate, which we found to have typical values of ∼10 5 . With this scaling of the accretion rate, we find that the Eddington ratios in our simulations vary from 10 -12 to 1, with mean values (av eraged o v er the entire simulation) of 0.05-0.3 for different simulation setups. This scaling is used only in the simulations with the SMUGGLE model, where the higher gas densities create these unphysical accretion rates. As described below, we also carry out a set of eEOS simulations for comparison. In all the eEOS simulations, α is set to 1, which we find yields similar average accretion rates as in the SMUGGLE model.

Historically, this α parameter has been used to scale up the Bondi accretion rate by an arbitrary factor (often chosen to be ∼100; e.g. Springel et al. 2005 ;Booth & Schaye 2009 ;Hayward et al. 2014 ). This was justified based on the assumption that gas densities would continue to increase on subresolution spatial scales, such that the Bondi rate would be underestimated in limited-resolution simulations. Essentially, this is indeed what we find in our high-resolution SMUGGLE simulations that allow high-density gas clouds to form on small scales. In reality, AGN feedback would regulate the gas density in the vicinity of the BH, such that a large reduction in the Bondi rate would not be necessary to achieve reasonable accretion rates. We note that this issue is exacerbated by the BH repositioning scheme in our runs where the BH is pinned to the local gravitational potential minimum. This issue is discussed in the next section. While the use of a scaling parameter of α = 10 -5 in place of AGN feedback is admittedly artificial, it enables us to achieve the primary goal of this study: conducting the first analysis of nuclear gas inflows around BHs within the SMUGGLE ISM model. The use of the α parameter is discussed further in Sections 3 and 4 .foot_1

In AREPO , gas cells are routinely refined and de-refined when their mass deviates from a fixed target mass by more than a factor of two. This ensures that all gas cells in the simulation have roughly a constant mass and hence a uniform mass resolution. In addition to this, one can also modify the refinement criterion to lower the target mass in specific regions of interest, in order to achieve a higher resolution locally with minimal increase in the o v erall CPU cost. Such a refinement scheme has already been used in AREPO to study cosmological gas accretion from the circumgalactic medium (Suresh et al. 2019 ). We implemented an analogous scheme in AREPO to study BH accretion. We lower the target gas cell mass by a factor F inside a sphere of radius r min centred at the black hole. Beyond this radius the target mass is linearly interpolated to its original value at a radius r max as follows,

where m 0 is the uniform resolution target mass and r is the distance of the gas cell from the black hole. The parameters F , r min , and r max can be varied. We set r min = 2 . 86 kpc, r max = 14 . 3 kpc and vary F between 3 and 30 for different runs. The value of r min has been chosen based on the largest values of the kernel radius, which is defined as the radius of a sphere centred around the BH such that the weighted number of gas cells inside the sphere is equal to the specified number of neighbours. Note that we do not switch off star formation and stellar feedback in the refinement region. The slow transition between the two target masses o v er a radial separation of 11.44 kpc minimizes any numerical errors due to the interaction between particles/cells of significantly different masses. In our refinement runs only a very small fraction ( < 0 . 5 per cent ) of the star particles in the refinement re gion hav e masses lar ger than ∼5 times the local tar get mass. We note that a different refinement prescription that modifies the gas cell radius rather than the gas cell target mass was implemented in AREPO by Curtis & Sijacki ( 2015 ). In that scheme, the radius of the gas cells inside the refinement region is a linear function of the distance from the BH. Curtis & Sijacki ( 2016b ) used this approach to formulate a modified Bondi prescription incorporating the angular momentum of the gas and showed that the angular momentum can limit BH growth and delay the quasar phase following a galaxy merger. It has also been used in several AGN feedback studies (Curtis & Sijacki 2016a ;Bourne & Sijacki 2017 ;Koudmani et al. 2019 ). The impro v ed resolution offered by the refinement scheme was shown to make significant changes in the coupling efficiency of AGN feedback relative to identical uniform resolution simulations. (Beckmann, Devriendt & Slyz 2019 ) to study BH accretion. Resolution of gas cells near the BH can also have an impact on the dynamics of BHs. When the mass of gas cells and other particles are comparable to the mass of the BH due to resolution limitations, twobody interactions can cause the BH to wander artificially about the galactic nucleus. In order to mitigate this issue, the BH is often re-positioned to the local potential minimum at every time-step. Refining gas cells near the black hole will minimize such numerical errors and better resolve dynamical friction due to gas, making the black hole dynamics significantly more stable relative to the uniform resolution simulations. Ho we ver, we find that in some of our merger simulations the position of the BHs have a few kpc scale fluctuations without the repositioning scheme even at high resolutions. Because of this, we use the repositioning scheme in all our simulations. We note that in the SMUGGLE model, especially in the absence of AGN feedback, using the repositioning scheme can lead to some enhancement of BH accretion rates due to frequent repositioning of the BH to the high density gas clouds. We find that in simulations where the BH position fluctuates by a similar amount with and without the repositioning scheme, BH masses do not differ by more than a factor of a few, indicating that our results are not dominated by the effects of the repositioning scheme. We discuss this issue in more detail in Section 4 .

Our isolated galaxy simulations involve three separate initial disc g alaxies: A Milky-Way type g alaxy (MW), a Small Magellanic Cloud lik e dw arf galaxy (SMC), and a LIRG-like galaxy (Sbc). Each galaxy initial conditions (ICs) contain a dark matter halo, central SMBH, stellar bulge, stellar and gaseous discs. The dark matter halo and stellar bulge are modelled using the Hernquist ( 1990 ) profile and the stellar and gaseous discs follow an exponential profile radially. The stellar disc has a sech 2 ( z) vertical profile and the gaseous disc has a profile set by hydrostatic equilibrium. The stellar bulge and disc particles are tracked as distinct particle types. The IC parameters are given in Table 2 . Our ICs are similar to the ones used by Hopkins, Quataert & Murray ( 2012 ) and Hayward et al. ( 2014 ). We also carry out equal-mass galaxy merger simulations, using the same initial conditions as progenitor galaxies (MW-MW, Sbc-Sbc, & SMC-SMC). In each case, the two galaxies are placed on a parabolic orbit. The initial separation is chosen such that the pericentric passage happens at 0.57 Gyr. The initial angular momenta of the two galaxies are at angles (30 • , 60 • ) and ( -30 • , 45 • ) with respect to the orbital angular momentum. This orbit is chosen as a well-studied example of a strongly interacting major merger. We run simulations at four different mass resolutions (resx0. 1, resx0.3, resx1, resx3). The lowest resolution level (resx0.1) has a target gas mass of 2.5 × 10 4 M , dark matter particle mass of 2.64 × 10 5 M , stellar bulge, and disc particles with mass of 2.59 × 10 4 M and 2.78 × 10 4 M . In the higher resolution runs, these masses are decreased by factors of 3, 10, and 30, respectively. These parameters are summarized in Table 3 . In some of the runs we combine low resolution with the refinement scheme described abo v e to obtain a higher resolution in the central re gions. F or e xample, the resx0.1 run with a refinement factor of 30 has the same gas mass resolution in the central regions as the resx3 run. In such cases only the target gas mass in the refinement region changes, whereas the masses of DM, stellar disc, and bulge particles remain unchanged.

We finally note that, owing to the smooth gas density profile in the initial conditions, an initial relaxation period is needed in our simulations before star formation and BH accretion reach steady state. We therefore do not analyse any data from the first 0.36 Gyr of each simulation. In addition, we do not allow BH accretion during the initial relaxation period. In the refinement runs the super-Lagrangian refinement scheme is used throughout the simulation including the initial relaxation period.

In this section, we look at the main differences between isolated galaxy simulations using the SMUGGLE and eEOS model of Springel & Hernquist ( 2003 ). Fig. 1 shows the face-on gas density projections of our fiducial-resolution MW disc galaxy simulated using the SMUGGLE (on the left) and eEOS model (on the right). One of the main differences in the structure of the ISM is the smooth gas distributions in the eEOS ISM compared to the complex and turbulent structure of the SMUGGLE ISM. The local injection of stellar feedback in SMUGGLE leads to large ( ∼kpc scale) low density cavities in the ISM. SMUGGLE is also able to produce gas clouds with much higher densities than the peak densities produced by the eEOS. This is due to the explicit modelling of the cold and dense phase of gas. Fig. 2 compares the time evolution of the total star formation rate, the gas density near the BH, and the BH Eddington ratio ( χ Edd = Ṁ Bondi / Ṁ Edd ) in the two simulations. Due to the turbulent and stochastic nature of the SMUGGLE ISM, all three quantities are much more stochastically variable than in the eEOS ISM. Unsurprisingly, the environment near the BH is much more sensitive to local fluctuations in the gas density than is the global SFR. The SFR fluctuates by a factor of a few o v er ∼Myr time-scales, while the BH Eddington ratio fluctuates by more than three orders of magnitude. In the eEOS ISM, the time evolution of these quantities is much smoother. The average SFR in the eEOS run is slightly higher than that of the SMUGGLE run because of the absence of stellar winds in the eEOS run. This does not make any qualitative changes to our results. The average gas densities near the BH in the SMUGGLE run are more than an order of magnitude larger than that of the BH in the eEOS run. The Eddington ratios are comparable between the two runs, ho we ver, which is a direct result of the α = 10 -5 scaling Table 3. Target gas mass m g , softening length of baryonic particles and dark matter DM , dark matter particle mass m DM , stellar disc particle mass m d and stellar bulge particle mass m b at different resolution levels of our simulations. The last column shows the refinement factors ( F ) used (in addition to the uniform resolution runs) at different resolution levels. In all refinement runs r min and r max are set to 2.86 and 14.3 kpc, respectively. These parameters are the same for all three initial conditions given in

Table 2 . Name m g (M ) (pc) DM (pc) m DM (M ) m d (M ) m b (M ) Refinement factor(s) resx0.1 2.5 × 10 4 45.7 98.6 2.64 × 10 5 2.78 × 10 4 2.59 × 10 4 3, 10, 30 resx0.3 7.7 × 10 3 30.6 65.9 8.82 × 10 4 9.26 × 10 3 8.62 × 10 3 10 resx1 2.5 × 10 3 21.4 45.7 2.64 × 10 4 2.78 × 10 3 2.59 × 10 3 -resx3 7.7 × 10 2 14.3 30.7 8.82 × 10 3 9.26 × 10 2 8.62 × 10 2 - factor used in the SMUGGLE runs, versus α = 1 in the eEOS runs.

The actual Eddington ratio of the BH (before applying the scaling factor) in the SMUGGLE run is several orders of magnitude higher than that of the eEOS run. As a result, Fig. 2 is useful primarily as a qualitative comparison between the two simulation setups.

Next we look at the resolution convergence of BH mass growth in the isolated MW galaxy simulations with the SMUGGLE ISM, including the simulations with super-Lagrangian refinement around the central BH. In Fig. 3 we plot the mass of BHs in the MW isolated galaxies at four different resolution levels, including four uniform-resolution and four refinement runs. The dashed lines correspond to runs with refinement, and the solid lines correspond to uniform-resolution runs.

In our colour scheme, runs with the same gas mass resolution in the central region are given the same colour.

For the first several hundred Myr after accretion is turned on, the BHs in all simulations grow at about the same rate. After t = 1.0 Gyr, the BH masses begin to diverge somewhat. Owing to the M 2 BH scaling of the Bondi accretion rate (equation 2 ), we expect to see any differences in BH growth amplified o v er time. By comparing the uniform resolution runs, we see that the BH growth rates decrease somewhat with increasing resolution, but BH masses converge within a factor of two for the fiducial (resx1) and higher resolutions. Further increasing the resolution does not make any significant changes. The resx3 run has a slightly higher BH mass than the resx1 run, but this difference is within the amount of fluctuations expected in these runs due to the stochasticity of the ISM. The intermediate resolution run (resx0.3) has a factor of a few higher BH growth and the lowest resolution run (resx0.1) has almost an order of magnitude higher final BH mass compared to the fiducial resolution run. This trend is consistent with the resolution convergence of SFR (not shown here) which decreases slightly with decreasing gas mass resolution. Note that resx0.3 simulations are consistent with the lowest resolution at which SMUGGLE has previously been shown to produce reasonably well converged results for stellar and ISM evolution (Marinacci et al. 2019 ).

By comparing the solid and dashed lines in Fig. 3 , we see that the low-resolution runs with refinement produce BH growth very close to the corresponding higher, uniform-resolution run. The refinement runs are also convergent at the fiducial gas mass resolution in the central re gion. We hav e v erified that increasing the refinement factor beyond 30 (we tested up to F = 100) does not make any significant changes in the BH growth. Since at such high resolutions ( m gas ∼ 250 M ) we are reaching the limits of the stellar evolution model,foot_3 we restrict the refinement factor to a maximum of 30 (equi v alent to the resx3 level). We have also included a run with a refinement factor of 10 in the resx0.3 run to show that changing both the background resolution and refinement will produce the same results, as long as the gas in the central region has the same resolution.

In Fig. 4 we compare the face on view of the gas density distribution of the resx0.1 (left-hand panel), resx0.1 F10 (middle), and the resx1 (right-hand panel) runs. In the top plots, we see that the resx0.1 and resx0.1 F10 runs look very similar but less resolved compared to the resx1 which has finer structures. The lower plots zoom into the central 4 kpc region of the top plots. We see that the gas distribution in the central region of the lowest resolution run Figure 5. Gas properties near the BH in MW isolated galaxy runs at resx0.1 (uniform low), resx1 (uniform high), and resx0.1 F10 (low refined) resolution levels. Top and middle panels show the kernel-weighted gas density and sound speed near the BH, and the bottom panel shows the average gas density within 1 kpc radius of the BH. The kernel-weighted gas density and sound speed have been averaged over 200 time-steps, which is roughly 10 Myr, and the average gas density within 1 kpc is plotted at every 7.1 Myr. The kernel-weighted gas density and sound speed of all three runs agree well. Ho we v er, the av erage gas density in the central region of the resx0.1 run is slightly higher than the resx1 run, whereas the resx0.1F10 run is in very good agreement with resx1 run.

is very smooth on kpc scales whereas the central region of the refinement run has more structure and looks closer to that of the higher resolution run. Note that the smooth central region ( r < 0 . 5 kpc) in the resx0.1F10 run is just a transient feature which lasts for around ∼200 Myr. Overall, the gas in the central region of the resx0.1F10 and the resx1 runs hav e v ery similar density profiles and morphology. The gas densities, star formation rates, and stellar mass in this region of the two simulations also agree well. Thus, by using the localized, BH-based refinement scheme we can accurately reproduce the central gas dynamics of the corresponding higher uniform-resolution runs.

In Fig. 5 we compare the kernel-weighted gas density and sound speed near the BH and the average gas density within a 1 kpc radius of the BH in the isolated MW simulations with resolutions of resx0.1, resx0.1 F10, and resx1. Despite the factor of ∼10 difference in the final BH masses between the low resolution and the higher resolution runs, the local gas densities and sound speeds near the BH of the three runs are converged within a factor of ∼2-3. Ho we ver, by comparing the average gas densities in the central region of the galaxies we see that the refinement run has a much better agreement with the higher resolution run compared to the low resolution run. In the resx0.1 F10 run we see some fluctuations in the central gas density M BH ( t) ) 2 , where Ṁ BH and M BH are the Bondi accretion rates and BH mass of the runs shown in Fig. 3 . The same line styles are used to differentiate each simulation, but the y -axis scale is much smaller here. All the runs have very similar mass growth, indicating that M 2 BH dependence of Bondi accretion rate amplifies small differences o v er time resulting in the large differences seen in Fig. 3 . and sound speed at ∼1 . 1 and ∼2 . 1 Gyr. The density spike at 1.1 Gyr leads to a spike in the central SFR. The SN explosions in the star particles created during the spike leads to the formation of a large cavity in the central region. This leads to a large drop in the central gas density and a corresponding increase in the sound speed which lasts for ∼200 Myr. We note that fluctuations of this type are not anomalous in the SMUGGLE model since the ISM is intermittently burst and turbulent (Torrey et al. 2017 ). Moreover, these quantities are calculated by averaging over a very small region of 10-100 pc size in the central region making them v ery sensitiv e to small changes in the ISM.

As noted abo v e, the Bondi accretion rate depends on the square of the BH mass. This means that even if the average ISM properties are convergent, small differences in the instantaneous BH accretion rate can, o v er time, lead to fairly large differences in the BH mass. To disentangle this non-linearity and to isolate the effects of gas properties, we re-scale the accretion rate by keeping the BH mass fixed and equal to the initial mass. Thus the modified accretion rate is Ṁ mod ( t) = Ṁ BH ( t)( M BH (0)

M BH ( t) ) 2 where Ṁ BH and M BH are the Bondi accretion rates and BH mass of the runs shown in Fig. 3 . We then integrate this modified accretion rate with time. This gives the gas that would have been accreted by a BH of fixed mass if it were embedded in the same environment as the actual BH in the simulation. In Fig. 6 we compare this quantity at different resolution and refinement levels. Note that this is a post-processing analysis of the same simulations shown in Fig. 3 , in a similar format. As expected we see a much better convergence here, with only a negligible difference in the total amount of gas accreted at different resolution levels (including the lowest resolution level). There is also much less gas accretion o v erall in this modified model-note the very different scales between Figs 3 and 6 . This indicates that the differences in ISM properties determining the BH accretion rates are very small in these runs, but the non-linear dependence of the Bondi rate on the BH mass exacerbates these differences. We would therefore expect to see even better resolution convergence in accretion models with a weaker dependence on BH mass (e.g.

DeBuhr et al. 2010 ; Hobbs et al. 2011 ; Hopkins & Quataert 2011 ). Downloaded from https://academic.oup.com/mnras/article/517/4/4752/6747143 by Oregon State University user on 23 March 2023 (blue) and SMC (green) isolated galaxies at the fiducial resolution. Bottom row: Eddington ratio of BHs in MW (left-hand panel), Sbc (middle), and SMC (right-hand panel) isolated galaxies at the fiducial resolution. The thick translucent lines show the Eddington ratios averaged over 150 Myr window. The higher disc gas fraction of SMC and Sbc galaxies leads to more bursty star formation which results in large fluctuations in the central gas distribution. This leads to much larger fluctuations in BH accretion rates compared to that of MW.

The main moti v ation for using the super-Lagrangian refinement scheme is to attain a higher resolution near the BH without a significant increase in the run time of the simulations. The refinement runs are much faster compared to the corresponding uniform resolution runs with the differences increasing with increasing resolution. F or e xample, the resx0.1 F10 run is roughly 10 times faster than the fiducial resolution run (resx1) and the highest refinement run (resx0.1 F30) is approximately 25 times faster than the highest uniform resolution run (resx3) and five times faster than the fiducial resolution run. Thus, in studies focusing mainly on SMBHs, using the refinement scheme can save significant computational resources while still being able to produce accurate results. We note that global properties like the SFR do not change significantly in runs with refinement.

We now compare results from the isolated MW disc simulation to results from the other two isolated simulations: the Sbc and SMC initial conditions. In Fig. 7 we compare the SFR, specific SFR (sSFR), BH mass and Eddington ratios of the Sbc and SMC galaxies to the MW. The top left-hand panel shows the total SFR and the sSFR of the three galaxies. The most rele v ant dif ference between these three simulated galaxies is that the Sbc and SMC have much higher disc gas fractions compared to MW, which leads to ele v ated and bursty star formation with fluctuations that at times exceed an order of magnitude. As expected based on the higher gas content, the sSFRs of the Sbc and SMC simulations are very similar to each other, and both are greater than that of the MW by a factor of ∼30. Due to the higher disc gas fraction, the Sbc and SMC galaxies also have much higher peak gas densities near the BH compared to MW. The bursty star formation in the Sbc and SMC runs leads to large fluctuations in the central gas densities, which is reflected in the Eddington ratio evolution plots. The Eddington ratio of the BH in the MW simulation fluctuates by two orders of magnitude with an average value of 0.16, whereas the corresponding fluctuations in Sbc and SMC span ∼9 orders of magnitude. These extreme fluctuations essentially mean that the BH is actively accreting for only short periods of time ( ∼ a few Myr) when the BH resides in a high-density gaseous region, and in between these bursts of growth it is essentially quiescent. None the less, the average Eddington ratios of the BHs in all three galaxies are similar (0.15 for MW, 0.13 for Sbc, and 0.11 for SMC). The median Eddington ratios are 0.06 for MW, 1.8 × 10 -7 for Sbc, and 2.2 × 10 -5 for SMC, owing to the much lower minimum accretion rates in those simulations. The similarity of the average Eddington ratios in Sbc and SMC indicates that the BHs spend similar amounts of time in the on state. The low values of the median Eddington ratio in Sbc and SMC are due to the large fluctuations. After BH accretion has been allowed to proceed in the simulation for a total of 2.5 Gyr, the BH masses in MW, Sbc, and SMC grow by factors of 22, 18, and 10, respectively.

In this section, we show the results of our idealized merger simulations, focusing first on mergers between the MW progenitor galaxies. As with the isolated runs, we have verified that the super-Lagrangian refinement merger simulations are in reasonably good agreement with the corresponding uniform high resolution runs. Accordingly, we do not present further detailed comparison between the refinement and uniform-resolution merger simulations. Fig. 8 shows both faceon and edge-on gas column density projection of the fiducialresolution SMUGGLE MW merger at four different epochs: before the first pericentric passage, during the first pericentric passage, between the first and second pericentric passages, and during the final coalescence. In the top images the line of sight is aligned with the angular momentum of the galaxy being plotted (and the BH is at the centre of the box). In the bottom plots the line of sight is rotated by 90 degrees. By comparing the first and third snapshots, we see that the close interaction between the galaxies strips some of the gas from the outer regions of the galaxies. None the less, the galaxies still retain the disc structure and higher gas densities in the central region, which is clear in both the edge-on and face-on views. The fourth snapshot shows that the final coalescence completely deforms the shape of the galaxies. A large fraction of the gas ( ∼35 per cent ) in the galaxies is now clustered in a very small ( ∼5 kpc) region at the centre with densities decreasing very rapidly with distance from the centre. After this snapshot strong feedback from SN explosions expel most of the gas from the central region of the galaxies.

In Fig. 9 we compare the MW-MW merger simulations at the fiducial resolution (resx1) with the SMUGGLE (on the left) and eEOS (on the right) models. In the bottom panel, we plot the distance between the BHs in the two galaxies. In both models, the galaxies go through two pericentric passages before the final coalescence 3 at t ∼ 1 . 7 Gyr. In the first and second panels we show the total SFR and the gas content in the central 3 kpc of the galaxies in these runs. The green and grey lines correspond to the SMUGGLE and eEOS 3 The BHs of the two galaxies are merged when they are within the neighbour search radius. isolated galaxies at the fiducial resolution. The SFRs of the isolated galaxies are multiplied by two for easier comparison.

In the SMUGGLE run, after the pericentric passages and the coalescence, there is a spike in the central gas content of the merging galaxies relative to the isolated galaxies. This indicates that the interaction between the merging galaxies drives an inflow of gas towards the nuclear region. The increase in nuclear gas content is also correlated with a spike in the star formation rates of these galaxies due to increased gas densities. After the first pericentric passage the SFR of the merging galaxies goes up by a factor of ∼10 followed by another increase by a factor of ∼10 after the coalescence. After ∼1 . 7 Gyr, the SFR fluctuates by 2-3 orders of magnitude and decreases rapidly as most of the central gas is expelled by stellar feedback. In the eEOS merger, we see a similar increase in the nuclear gas content and SFR after the pericentric passages and coalescence, but the SFR varies in a smooth fashion, whereas the SMUGGLE galaxies are significantly more bursty. After the final coalescence in the eEOS runs there is no secular decrease in the SFR and central gas content.

In rows 3-6 of Fig. 9 we compare the time evolution of the BH accretion rates, gas density, and sound speed near the BHs, and the BH masses in these runs. We focus first on the SMUGGLE results. The gas densities near the BHs in the SMUGGLE merger run stay roughly the same with small fluctuations until the first pericentric passage. After the passage, the peak densities increase by 2-3 orders of magnitude, followed by a further increase by a factor of few after the final coalescence. Ho we ver, the strong stellar feedback due to the ele v ated and bursty star formation after the pericentric passage T op panel: T otal SFR of the two MW galaxies in the merger run and twice the SFR of the isolated MW galaxy with SMUGGLE (left-hand panel) and eEOS (right-hand panel) models at the fiducial resolution. Panel 2: gas mass within the central 3 kpc region of the two galaxies in the merger run (red) and isolated galaxy (gre y). P anel 3: Accretion rate of the BHs in the merging and isolated galaxies. Panels 4 and 5: Kernel-weighted gas densities and sound speeds near the BHs in the merging and isolated runs. Panel 6: Masses of the BHs in the merging and isolated MW galaxies. Bottom panel: Distance between the two BHs in the merging galaxies. Sound speeds and gas densities near the BHs are plotted at every time-step of the simulation. All other quantities are plotted at 7.1 Myr intervals. The interaction of the galaxies trigger gas inflow towards the central region resulting in elevated SFR and BH accretion rates. Due to the large fluctuations in the central gas distribution in the SMUGGLE run, the merger induced enhancement of BH mass growth is much smaller than that of the eEOS runs.

leads to the formation of large low density cavities in the central regions of the galaxies. This results in large fluctuations (4-5 orders of magnitude) in the central gas density and limits the average gas density as the BHs end up spending significant amounts of time in these cavities. The maximum central sound speeds also increase by 2-3 orders of magnitude, as the same feedback that dro v e the formation of the cavities also heated the gas. Note that the minima of the sound speed are correlated with the peaks of density . Additionally , the minimum central sound speeds in the merger run are lower by a factor of a few relative to the equivalent isolated galaxy.

The higher densities and lower sound speeds lead to large spikes in the Bondi accretion rates, which are suppressed by the α scaling factor and by the Eddington limit. During the density peaks, the BHs accrete near the Eddington limit for short durations, in between which the accretion rates are low. Between the first pericentric passage and the final coalescence there is a factor of few increase in the average accretion rates. During and after the coalescence, the accretion rates increase further, to a factor of ∼10 abo v e their initial average values. Thus, although local stellar feedback modulates the BH accretion rate on short time-scales as in the isolated galaxy, here the interaction between MW-type galaxies drives cold gas to the galactic nuclei, thereby enhancing the total BH growth. In the fiducial-resolution SMUGGLE runs the total BH mass growth o v er a period of ∼1 . 9 Gyr in the merger run is enhanced by a factor of ∼4 relative to that of the isolated run.

In the eEOS merger, there is a steady increase in the gas density near the BH after the first pericentric passage. As in the SMUGGLE merger, this results from the nuclear inflow of gas driven by gravitational perturbations. But in this case, the smooth gas distribution in the central region of the galaxies (as seen in Fig. 1 ) yields steady, Eddington-limited accretion for ∼500 Myr, after which the gas densities and accretion rates decrease as the galaxies run out of gas in the central region. During the second pericentric passage and the final coalescence, more gas is funnelled into the central region, and as a result the central gas density and accretion rates spike again. The central sound speed increases slightly after the pericentric passages and the coalescence. Although the time-averaged central densities are comparable in both the eEOS and SMUGGLE mergers, the factor of ∼10 difference in the minimum sound speeds make a significant difference in the accretion rates in the two models due to the c -3 s dependence of the Bondi accretion rate. The Eddington ratio of the BHs in the eEOS merger has an average value of 0.65 and a median of 0.59 during 250 Myr period immediately after the first passage. During the coalescence the average is 0.63 and median is 0.8. In the SMUGGLE merger the average Eddington ratio is 0.23 during both phases, but the median is lower during these periods ( ∼2-3 × 10 -6 ), again reflecting the large fluctuations between active accretion episodes and inactive periods. After ∼2 Gyr the total BH growth in the eEOS merger run is enhanced by a factor of ∼40 relative to the eEOS isolated run, which is an order of magnitude larger enhancement than that in the SMUGGLE run.

Note that α = 1 in the eEOS runs, meaning that no scaling factor is applied to the Bondi rate. As indicated by the density and sound speed plots (Fig. 9 ), the BH growth in the SMUGGLE runs would be much larger than that in the eEOS runs if no scale factor were applied to the Bondi rate. We also note that since most of the BH growth in the eEOS merger run happens during the steady, Eddington-limited accretion, the final BH mass is not expected to change significantly for α in the range ∼0.001-1. It is difficult to predict exactly how these comparisons would change in SMUGGLE runs with α = 1 that also included some form of BH feedback. Depending on the nature and strength of the feedback, and the galactic environment, we might see either a larger or a smaller difference in BH growth in merging galaxies relative to the isolated galaxies.

Our results from the eEOS run are in general consistent with past studies of BH growth during mergers (e.g. Blecha et al. 2011 ;DeBuhr, Quataert & Ma 2011 ;Thomas et al., in preparation). In the SMUGGLE model the large scale (few kpc) gas inflow driven by the interaction of the galaxies is very similar to that in the eEOS model. But the fueling of BHs that happens at much smaller scales ( ∼10 pc) is dramatically different due to the stochastic stellar feedback and the presence of cold and dense gas clouds. Even though the BH fueling in the SMUGGLE run is enhanced during the merger, the level of enhancement is much more modest compared to that in the eEOS model.

In Fig. 10 we show the SFR, BH accretion rate, BH mass, and separation between the BHs in the Sbc-Sbc (left-hand panel) and SMC-SMC (right-hand panel) merger simulations with the SMUGGLE model. In the SMC run the SFR and nuclear gas content spikes for a short duration at the time of the first pericentric passage and decreases afterwards. This is different from the case of MW merger, where the SFR and central gas content increases almost steadily following the first passage. As seen in Fig. 7 , star formation is significantly more bursty in the SMC run compared to the MW. The close interaction of the galaxies during the first passage leads to a further increase of SFR and hence stronger stellar feedback. Strong winds driven by the feedback expel a large amount of gas from the central region of the galaxies, thereby limiting the gas inflow and lowering the SFR. This is seen as the sudden drop of M gas at ∼0 . 8 Gyr. This results in lowered accretion rates until the second pericentric passage. The final coalescence brings more gas to the central region of the galaxies. After the coalescence, the central gas content increases by a factor of ∼5 and the global SFR increases by roughly two orders of magnitude with large fluctuations o v er ∼100 Myr time-scales. This increases the local gas density near the BHs and consequently the time-averaged accretion rate increases by an order of magnitude. The final BH mass in the SMC merger at t = 2.4 Gyr is ∼4.3 times larger than twice the final BH mass in the isolated galaxy. The enhancements in the SFR and the BH accretion rates after the final coalescence in the SMC merger are similar to that of the MW merger.

In the Sbc merger the global SFR and nuclear gas content behave very similarly to SMC until the second pericentric passage due to the similar bursty star formation of Sbc. The spike in SFR during the first passage is smaller compared to that of SMC. The expulsion of gas by stellar feedback from the galactic centres after the first passage slightly lowers the accretion rates. After the second pericentric passage and the final coalescence the SFR and nuclear gas content increase by less that an order of magnitude and a factor of ∼3, respectiv ely. The time-av eraged BH accretion rates increases only by a factor of few after the coalescence, and the total BH mass growth is enhanced only by a factor of 1.25 relative to the isolated run. These enhancements are much smaller than the corresponding enhancements of SFR and BH accretion rates in the MW and SMC merger.

Thus, the interaction of the galaxies enhance the BH accretion rates in all three merger simulations. The effect of the first pericentric passage on the SFR and BH accretion rates is more pronounced in the MW merger, whereas the SMC and Sbc mergers have virtually no nuclear gas inflows at this stage. After the coalescence all three systems become significantly more bursty, with the MW merger Downloaded from https://academic.oup.com/mnras/article/517/4/4752/6747143 by Oregon State University user on 23 March 2023 Figure 10. T op panel: T otal SFR of the merging galaxies and twice the SFR of the isolated galaxy with Sbc on the left and SMC on the right. Panel 2: Gas mass within 3kpc of the two merging galaxies and the isolated galaxy. Panel 3: Kernel-weighted gas density near the merging and isolated BHs at each time-step of the simulations. Panels 4 and 5: Accretion rates and masses of the BHs in the merging and isolated galaxies. Panel 6: Distance between the two BHs in the merging galaxies. Gas densities near the BHs are plotted at every time-step of the simulation. All other quantities are plotted at 7.1 Myr intervals. In both SMC and Sbc mergers, strong stellar feedback after the 1st pericentric passage expels gas from the central region resulting in lower BH accretion until the 2nd passage. The final coalescence leads to gas inflows resulting in ele v ated accretion rates. The o v erall enhancement of BH growth in SMC is similar to that of MW, whereas Sbc has only a very small enhancement.

showing the largest fluctuations in SFR. BH accretion rates also increase, in correlation with the SFR. The accretion rates in Sbc and SMC mergers are lower that that of the isolated runs until the coalescence and increases to a higher value after the coalescence resulting in an o v erall enhancement of BH mass growth. Thus, our results indicate that galaxy mergers can play an important role in triggering AGN activity, in agreement with previous work, but we also find that the amount and time-scale of merger-induced BH fueling episodes depend strongly on the nature of the surrounding ISM.

In this paper, we have investigated the nature of gas inflows and BH fueling in AREPO hydrodynamics simulations using the multiphase ISM and explicit stellar evolution model SMUGGLE. Our simulation suite included initial conditions for three different progenitor g alaxies: a g as-poor MW type g alaxy and two g as rich g alaxies (Sbc and SMC). These were each evolved in isolation, and also in a set of idealized galaxy merger simulations, o v er long time-scales of ∼2-3 Gyr. Most of our simulations resolve gas dynamics at ∼10 -100 pc scales.

We find that the clumpy and turbulent ISM of the SMUGGLE model results in stochastic fueling of BHs with orders of magnitude fluctuations in the accretion rates o v er ∼ Myr time-scales. Star formation is also stochastic in the SMUGGLE model, but we find that BH accretion rates are much more sensitive to the stochastic variations in the local ISM density. Gas-rich galaxies with bursty star formation have larger fluctuations in the BH accretion rates due to the formation of large low density cavities in the ISM. The stochastic nature of BH accretion shown in our work is similar to the behaviour seen in FIRE simulations (Angl és-Alc ázar et al. 2017b ) with a gravitational torque based accretion prescription (Angl és-Alc ázar et al. 2017a ), although the type of galaxies in their study were different from ours. They also found bursty BH accretion modulated by stellar feedback, especially at high redshifts.

By varying the resolution levels, we have demonstrated that the Bondi-Hoyle accretion model coupled to SMUGGLE is resolution convergent for gas cell masses ࣠ 3 × 10 3 M . The BH masses after 2.5 Gyr agree to within a factor of ∼4 for all except the lowest resolution run. When the M 2 BH dependence is scaled out of the Bondi accretion rate to compare the amount of gas accreted by BH of fixed mass, all the runs agree within a factor of ∼1.5. Additionally, we implemented a super-Lagrangian refinement scheme which increases the gas mass resolution in the immediate neighbourhood (within 2.86 kpc) of the BHs. Using this scheme in the lowest resolution simulations ( m gas = 2.5 × 10foot_4 M ) produced central gas dynamics and BH growth in good agreement with that of the high-resolution runs. This allows high resolution studies of BH dynamics with only a minimal increase in CPU cost. Isolated galaxy runs with refinement factors of 10 and 30 were approximately 10 and 25 times faster than the corresponding high uniform resolution runs. Often, the highest resolution achieved by cosmological and zoom in simulations is comparable to the lowest resolution considered here. This implies that using the refinement scheme in these simulations can accurately resolve BH accretion where the BH growth will otherwise be o v erestimated due to the low resolution. However, further testing is required to ensure the validity of the refinement scheme for refinement factors abo v e 30.

In our isolated disc simulations, we find that the higher gas fraction of Sbc and SMC galaxies leads to much more stochastic star formation and BH fueling relative to that of the MW galaxy. The BH in the MW-like simulation has an average Eddington ratio of 0.16 with roughly two orders of magnitude fluctuations o v er Myr timescales whereas the BHs in Sbc and SMC like galaxies have average Eddington ratios of ∼ 0.1 with nine orders of magnitude fluctuations. Thus, in the latter cases the BH is essentially switching back and forth between an active fueling state and a quiescent state, dictated by whether it is instantaneously located in a high-density gas cloud or a low-density region evacuated by local stellar feedback. In reality, the accretion rate should vary more gradually as the accretion disc is consumed by the BH o v er the viscous time-scale, and AGN feedback would further modulate the fueling environment. None the less, these results demonstrate the dramatic impact that the multiphase ISM can hav e o v er short time-scales on the cold gas supply available to the central BH.

We also studied BH fueling during galaxy mergers. The impact of different stages of the merger on BH accretion depends strongly on the type of galaxies. In the MW merger, the first pericentric passage leads to large inflow of gas to the central region of the galaxies. This increases the peak gas densities near the BHs by several orders of magnitude while also leading to large fluctuations due to the bursty stellar feedback. The net effect is an increase in average accretion rates. In the Sbc and SMC mergers, owing to their high sSFR, strong stellar feedback following the first pericentric passage expels a large amount of gas from the central region of the galaxies resulting in lower accretion rates. Thus the first pericentric passage produces opposite effects on BH fueling in Sbc and SMC runs compared to the MW run. Ho we ver, the second pericentric passage and the final coalescence leads to inflows in all three systems resulting in enhanced BH accretion. The inflows in MW and SMC mergers are much stronger compared to the inflows in the Sbc merger. The total BH mass growth in the MW and SMC merger runs is enhanced by a factor of ∼4 relative to that of the corresponding isolated runs whereas in the Sbc merger run there is only a factor of ∼1.25 enhancement. Thus the o v erall enhancement of BH growth depends strongly on the initial conditions. We also compared our SMUGGLE runs with simulations of the same initial conditions using the eEOS model, the ISM and stellar evolution model used in Illustris, IllustrisTNG, and Auriga simulations. We find that the merger induced enhancement of BH fueling in the SMUGGLE model is much smaller compared to that in the eEOS model. This is because of the bursty stellar feedback in the SMUGGLE runs which results in high variability at short time-scales of the central gas distribution counteracting the inflows. In the MW merger with the eEOS model, the smooth ISM and the inflow of gas leads to a steady increase in BH accretion rate with a factor of ∼40 enhancement in the total BH mass growth relative to the isolated runs.

These results imply that major mergers of the type of galaxies considered here can trigger AGN acti vity. Ho we ver, simulations that do not explicitly model stellar feedback and resolve the multiphase structure of the ISM can significantly o v erestimate merger induced enhancements in AGN activity.

Although the initial conditions of our simulations represent diverse types of galactic environments, to reach more general conclusions regarding the nature of BH accretion flows in a multiphase ISM, more e xhaustiv e study of the parameter space would be needed. Such type of studies will be easier to perform using simulations of a cosmological volume, an avenue that we aim to pursue in future work.

The subgrid prescriptions for BH accretion, feedback, and dynamics are also important areas for future work. The simulations presented here do not include any AGN feedback, in order to focus on nuclear gas inflows separately from their highly non-trivial coupling to AGN outflows. Because the explicit ISM model allows gas to condense to very high peak densities, particularly in galactic nuclei, the Bondi accretion model produces unrealistic BH mass growth o v er time in the absence of AGN feedback. As discussed in Section 2.3 , in order to simulate long periods (several Gyr) of galactic evolution in isolated and merging systems, we used the α parameter in the Bondi accretion formula to artificially scale down the accretion rates by a large factor (10 -5 ). With this scaling, the average and median values of the uncapped Eddington ratios in the isolated discs are (0.31, 0.06) in MW, (0.75, 2.2 × 10 -5 ) in SMC, and (110, 1.8 × 10 -7 ) in Sbc runs. In the merger runs the average and median Eddington ratios are (260, 2.5 × 10 -6 ) in MW, (23, 1.1 × 10 -5 ) for SMC, and (150, 1.7 × 10 -7 ) in Sbc. Given the stochastic, highly variable nature of BH accretion within the SMUGGLE ISM model, these appear to be reasonable values for the range of AGN simulated here. 4 In reality, BH growth is expected to be self-regulated through AGN feedback.

In future simulations, we plan to study this self-regulated growth directly.

BH dynamics, including dynamical friction, cannot be fully resolv ed ev en in the highest resolution re gimes of our simulation suite. To prevent spurious wandering of the BH due to numerical noise, we therefore employ the common approach of repositioning the BH on the potential minimum at every time-step. As mentioned in Section 2.4 , we find that this BH repositioning scheme can have some significant effects on BH dynamics in the SMUGGLE ISM. Because this model directly traces the formation of high-density gas clouds and local stellar feedback injection within a clumpy, multiphase ISM, the gas distribution varies dramatically o v er short time and spatial scales. As a result, the repositioning scheme has a tendency to localize the BH within high density gas clouds, which can lead to some artificial enhancement of the accretion rate. Ho we ver, we MNRAS 517, 4752-4767 (2022) find that not using the repositioning scheme in merger simulations can artificially suppress the BH accretion rates, particularly in the post-merger stage, as the BH is subject to spurious dynamical kicks from comparable-mass particles. In future work, it will be useful to consider a more realistic dynamical friction prescription in place of this repositioning scheme. Ho we ver, we note that this issue will also be mitigated by the inclusion of AGN feedback, as the feedback will quickly disperse the high density gas clouds and decrease the gas density in the central region of the galaxies.

The interplay between BH fueling, star formation, and stellar feedback is a crucial but complex aspect of BH and galaxy coevolution. We have undertaken a novel numerical study of BH fueling in an explicit, multiphase ISM. We implemented a novel super-Lagrangian refinement scheme and studied gas inflows on to BHs o v er long time-scales and in diverse merging and isolated galactic environments. Our results clearly demonstrate the profound impact that the local ISM can have on BH fueling rates and duty cycles, while confirming previous findings that global perturbations (such as galaxy interactions) also modulate nuclear gas inflows and trigger AGN activity. Our work lays the groundwork for future studies of BH/galaxy co-evolution, including studies that also consider the role of AGN feedback within an explicit, multiphase ISM.