title: Cosmological simulations of quasar fueling to sub-parsec scales using Lagrangian hyper-refinement
authors: Daniel Anglés-Alcázar Eliot Quataert Philip Hopkins Rachel Somerville Christopher Hayward Claude-André Faucher-Giguère Greg Bryan Dušan Kereš Lars Hernquist James Stone
text: 
  
  
  
  
  
  Daniel Anglés-Alcázar
  angles-alcazar@uconn.edu
  
  
  Department of Physics, University of Connecticut, 196 Auditorium Road, U-3046, Storrs, CT 06269-3046, USA
  Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA
  
  
  Department of Astronomy and Theoretical Astrophysics Center, University of California Berkeley, Berkeley, CA 94720, USA
  Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA 
  
  
  TAPIR, Mailcode 350-17, California Institute of Technology, Pasadena, CA 91125, USA
  
  
  Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA
  Department of Physics and Astronomy, Rutgers University, 136 Frelinghuysen Rd, Piscataway, NJ 08854, USA
  
  
  Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA
  
  
  CIERA and Department of Physics and Astronomy, Northwestern University, 2145 Sheridan Road, Evanston, IL 60208, USA
  
  
  Department of Astronomy, Columbia University, 550 West 120th Street, New York, NY 10027, USA
  Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA
  
  
  Department of Physics, CASS, University of California at San Diego, 9500 Gilman Drive, La Jolla, CA 92093, USA
  
  
  Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA
  
  
  School of Natural Sciences, Institute for Advanced Study, 1 Einstein Drive, Princeton, NJ 08540
  
  
  
  
  
  
  We present cosmological hydrodynamic simulations of a quasar-mass halo (M_ halo≈ 10^12.5  M_⊙  at z=2) that for the first time resolve gas transport down to the inner 0.1 pc surrounding the central massive black hole.  We model a multi-phase interstellar medium including stellar feedback by supernovae, stellar winds, and radiation, and a hyper-Lagrangian refinement technique increasing the resolution dynamically approaching the black hole.  We do not include black hole feedback.
  We show that the sub-pc inflow rate (1) can reach ∼6  roughly in steady state during the epoch of peak nuclear gas density (z∼ 2), sufficient to power a luminous quasar, (2) is highly time variable in the pre-quasar phase, spanning 0.001–10  on Myr timescales, and (3) is limited to short (∼2 Myr) active phases (0.01–0.1 ) followed by longer periods of inactivity at lower nuclear gas density and late times (z∼1), owing to the formation of a hot central cavity.  
  Inflowing gas is primarily cool, rotational support dominates over turbulence and thermal pressure, and star formation can consume as much gas as provided by inflows across 1 pc–10 kpc.  Gravitational torques from multi-scale stellar non-axisymmetries dominate angular momentum transport over gas self-torquing and pressure gradients, with accretion weakly dependent on black hole mass.  
  Sub-pc inflow rates correlate with nuclear (but decouple from global) star formation and can exceed the Eddington rate by ×10.  
  The black hole can move ∼10 pc from the galaxy center on ∼0.1 Myr. 
  Accreting gas forms pc-scale, rotationally supported, obscuring structures often misaligned with the galaxy-scale disk.  
  These simulations open a new avenue to investigate black hole–galaxy co-evolution.
   
  
  
  
  
  
  
  
  
  
  
  § INTRODUCTION
  
  
  
  The inflow of gas from large scales down to galactic nuclei plays a key role in galaxy formation, driving the growth of central massive black holes and a variety of related phenomena: from bright quasars (QSOs) that outshine their host galaxies <cit.> to active galactic nuclei (AGN) “feedback” in the form of fast nuclear outflows <cit.>, galaxy-scale winds <cit.>
  and radio-emitting jets <cit.> that may have a significant impact on galaxy evolution <cit.>.    
  The scaling relations between central black hole mass and properties of their host galaxies <cit.> and the similarity between the global cosmic histories of star formation and black hole accretion <cit.> further suggest some form of black hole–galaxy co-evolution over cosmological timescales.  
  It is thus crucial to understand the mechanisms driving gas inflows from galactic scales down to the black hole accretion disk in a full cosmological context, which remains a major challenge.
  
  
  Modern large volume cosmological hydrodynamic simulations such as Magneticum <cit.>, Horizon-AGN <cit.>, Eagle <cit.>, Illustris <cit.>, BlueTides <cit.>, Romulus <cit.>, IllustrisTNG <cit.>, and SIMBA <cit.> have implemented sub-grid models of black hole growth and feedback with increasing success at reproducing global galaxy properties and black hole observables <cit.>.  
  However, despite much recent progress, the typical ∼kpc scale resolution available 
  (orders of magnitude from resolving accretion disk feeding on sub-pc scales)
  makes sub-grid black hole modeling very schematic by necessity and a major source of uncertainty.  With strong degeneracies between accretion and feedback, the successes of most current models rely on carefully tuning sub-grid parameters that significantly limit their predictive power <cit.>.
  For example, simulations implementing Bondi accretion invariably predict that AGN feedback self-regulation drives the black hole–galaxy scaling relations <cit.>, while simulations implementing gravitational torque-driven accretion emphasize the role of a common gas supply for star formation and black hole growth <cit.>. 
  The inferred impact of AGN feedback depends not only on the assumed efficiency and accretion model, but also on the sub-grid treatment of the interstellar medium (ISM), star formation, and stellar feedback.
  
  
  Cosmological “zoom-in” simulations provide a higher resolution framework to study massive black holes in galaxies <cit.>.  The Feedback In Realistic Environments (FIRE) project[<http://fire.northwestern.edu>.] <cit.> implements stellar feedback on the scale of star-forming regions, injecting energy, momentum, mass, and metals from supernovae (SNe), stellar winds, and radiation directly based on the predictions of stellar population synthesis models.  
  This yields a realistic ISM with multi-phase structure and self-consistent generation of galactic winds <cit.> while reproducing a variety of galaxy properties <cit.> and circumgalactic medium (CGM) observables <cit.>.  
  The increased resolution and detailed stellar feedback have important implications for black hole scaling relations, AGN demographics, and galaxy quenching:  stellar feedback can significantly suppress early black hole growth by efficiently evacuating gas from galactic nuclei <cit.>.  
  Qualitatively similar effects are seen in simulations implementing Bondi accretion <cit.> but the details crucially depend on resolution, ISM physics, and black hole parameterization.
  
  
  Idealized (non-cosmological) simulations offer the possibility to study gas transport in galaxies in a more controlled setup and at much higher resolution than typically available in cosmological simulations.  Early models showed that large-scale tidal torques induced by galaxy mergers or bar/spiral wave instabilities in self-gravitating disks can lead to angular momentum transport and rapid inflow of gas to the central sub-kpc of galaxies <cit.>, confirmed by more recent simulations <cit.>.  Gas inflow to <100 pc scales requires additional processes, as large-scale torques become less efficient <cit.>.  Subsequent secondary instabilities occurring at smaller scales <cit.> are a promising mechanism to transport gas down to the black hole accretion disk, where turbulent MHD processes and/or MHD winds are expected to dominate <cit.>. Idealized simulations with sub-pc resolution have shown that such cascades of instabilities do occur under certain conditions <cit.>.  
  
  
  These idealized studies have greatly informed our understanding of black hole fueling, providing the basis for more physically motivated accretion prescriptions that can be implemented in cosmological simulations <cit.>.
  However, idealized simulations are limited by artificial initial conditions that do not reflect the complexity seen in cosmological simulations, where the more realistic larger scale environment includes cold filamentary accretion <cit.>, pressure-supported gas in hot halos <cit.>, fountain flows <cit.>, and global mass transport in disks <cit.>. 
  Idealized models are also more likely to experience artificial transients (in e.g. morphological and star formation properties), cannot reach steady state owing to gas consumption/evacuation and the lack of larger scale inflow, and cannot self-consistently determine the origin of gravitational torques (and if e.g. instabilities are initiated by mergers or internal processes) owing to their limited dynamic range of scales.
  In addition, many previous idealized simulations have been limited by the use of simplified sub-grid ISM models inherited from simulations with ∼kpc resolution, missing the small-scale multi-phase structure of the ISM in the nuclear region.
  
  
  In this work, we aim to model black hole growth in a full cosmological galaxy formation context while reaching sub-pc resolution, bridging the crucial gap between idealized models and cosmological hydrodynamic simulations.
  Our goal is to model explicitly the inflow of gas down to the scales ∼0.1 pc below which the accretion disk is expected to form around the central black hole in a massive galaxy <cit.>, 
  including multi-scale gravitational torques between the gas and stellar components, gas consumption by star formation, turbulence and gas ejection driven by stellar feedback, and non trivial black hole dynamics in the nuclear potential.  We achieve this by implementing multi-phase ISM physics from the FIRE project, a novel hyper-Lagrangian refinement technique that increases the resolution dynamically closer to the black hole, and a sub-pc scale treatment of black hole growth that reduces many of the uncertainties introduced by parameterized sub-grid models on 0.1–1 kpc scales.
  
  
  As a first application of this methodology, we focus on understanding the mechanisms responsible for gas transport across spatial scales and black hole growth during qualitatively distinct phases of a massive galaxy (M_⋆∼ 6–20 × 10^10 ) before, during, and after its peak of nuclear gas density at z∼ 2.5→1.  
  This allows us to investigate the physical conditions conducive to a broad range of black hole accretion properties, from large QSO-like inflow rates down to extended periods of inactivity, which represents a unique challenge for current models.  Given the complex multi-physics involved, we choose to neglect the effects of black hole feedback in an attempt to reduce model uncertainties and avoid making prior assumptions about the efficiency of AGN feedback.
  The results presented here provide the baseline for future simulations including AGN feedback.
   
  
  We describe our methodology in <ref>, focusing on the galaxy formation model (<ref>) and the use of pre-existing FIRE cosmological zoom-in simulations as the initial conditions (<ref>) to implement our new hyper-Lagrangian refinement technique (<ref>). 
  We present an overview of our simulations in <ref>, illustrating the radical increase in resolution relative to previous models (<ref>) and their unique dynamic range (<ref>).  
  We then quantify the galaxy radial structure and thermodynamic state of gas over five orders of magnitude in spatial scales (<ref>), which provides significant insight into the physics driving black hole growth and the connection between accretion and galaxy properties on larger scales.
  Our main results are presented in <ref>, where we measure for the first time resolved accretion rates at 0.1 pc under different conditions at the center of a massive galaxy, reaching peak inflow rates sufficient to power a luminous quasar (<ref>).  We study the connection between star formation and gas inflow rate across the full range of spatial scales 0.1 pc–10 kpc (<ref>), with important implications for black hole fueling and the star formation-AGN connection.
  We quantify the angular momentum properties of gas across scales (<ref>) and evaluate the relative contributions of gravitational torques and pressure gradients to angular momentum transport (<ref>), a crucial step toward identifying the dominant mechanisms driving gas inflows.
  We further study the intrinsic dependence of accretion rate on black hole mass (<ref>), with important implications for sub-grid black hole accretion prescriptions.
  We discuss our results and conclude in <ref>.
  Additional details about the robustness of our methodology are presented in <ref>, including systematic tests of resolution convergence (<ref>) and other numerics (<ref>).
  
  
  This work has been developed with the support of the Simulating Multiscale Astrophysics to Understand Galaxies (SMAUG) consortium[<https://www.simonsfoundation.org/flatiron/center-for-computational-astrophysics/galaxy-formation/smaug/>], which attempts to reduce the uncertainties associated with sub-grid parameterizations of key baryonic processes in galaxy formation to improve the accuracy of fundamental cosmological measurements. 
  Other work in SMAUG has centered on star formation and stellar feedback <cit.>, the physics of the CGM <cit.>, and the connection between semi-analytic models and cosmological hydrodynamic simulations <cit.>.  Our work complements these by providing a key step towards developing better sub-grid models for black hole accretion in galaxy formation, which play a key role in the evolution of massive galaxies and large-scale structure.
  
  
  
  
  
  § METHODS
  
  
  
  
  
   §.§ Cosmological galaxy formation model
  
  
  We use the FIRE-2 galaxy formation model presented in <cit.>, which is an updated numerical implementation of the original FIRE simulations <cit.> with more accurate hydrodynamics, gravitational softening, and supernova coupling algorithms.  FIRE-2 employs the GIZMO[<http://www.tapir.caltech.edu/ phopkins/Site/GIZMO.html>] N-body+hydrodynamics code <cit.> in its “meshless finite mass" (MFM) mode.  MFM is a mesh-free, Lagrangian finite-volume Godunov method that combines the advantages of particle-based and grid-based methods: MFM provides automatic adaptivity in resolution while minimizing advection and angular momentum conservation errors relative to AMR codes and exhibits superior performance in shock-capturing and fluid-mixing problems while avoiding the low-order errors inherent to smooth particle hydrodynamics (SPH).
  As shown in <cit.>, MFM performs well with different particle masses unlike many SPH methods where the errors increase when having unequal-mass particles interacting, which is crucial for our hyper-refinement simulations.
  
  
  We adopt a “standard" flat ΛCDM cosmology with parameters H_0 = 69.7 km s^-1Mpc^-1, Ω_ M = 1-Ω_Λ = 0.2821, Ω_ b = 0.0461, σ_8 = 0.817, and n_s = 0.9646 <cit.>. 
  Gravitational forces are computed using an improved version of the tree-particle-mesh gravity solver of the GADGET-3 code <cit.>, including adaptive gravitational softenings for the gas component.  Following <cit.>, we assume that the mass distribution corresponding to each element has the same functional form as the interaction kernel in MFM.  The time integration is fully adaptive with a power-of-two hierarchy for assigning individual time-steps <cit.>, including a time-step limiter to handle strong feedback events <cit.>.
  
  
  
  
  
    §.§.§ Cooling, star formation, and stellar feedback
  
  
  We incorporate radiative cooling and heating processes from T = 10–10^10 K as detailed in <cit.>, including free-free, photoionization/recombination, Compton, photo-electric, metal-line, molecular, fine-structure, dust collisional, and cosmic ray[We include approximate cosmic ray heating in the dense ISM but do not include full cosmic ray transport and feedback as in e.g. <cit.> and <cit.>.] processes. 
  Star formation occurs only in gas which is (1) locally self-gravitating following the criterion of <cit.>, (2) molecular with self-shielding fraction given by the local Sobolev approximation <cit.>, (3) Jeans unstable with gas element mass larger than the thermal Jeans mass, and (4) above a minimum hydrogen number density n_ H≥ 1000 cm^-3. 
  Eligible gas elements are converted into collisionless star particles with a probability set by their SFR density ρ̇_⋆ integrated over a timestep, where 
  ρ̇_⋆ = ρ_ mol  / t_ ff and t_ ff is the local free-fall time.  This is consistent with the rate at which small, locally self-gravitating clumps fragment in higher-resolution simulations of turbulent clouds <cit.> and analytic models based on turbulent fragmentation <cit.>, which may only be a small fraction of the dense gas mass.  The global efficiency of star formation is self-regulated by stellar feedback at ∼1–10% per free-fall time and is insensitive to the details of the star formation model <cit.>.
  
  
  
  Each star particle represents a single stellar population with known mass, age, and metallicity.  
  Stellar feedback is implemented at the scale of star-forming regions following explicitly the time dependence of several different mechanisms, including energy, momentum, mass, and metal injection from (1) Type Ia and Type II Supernovae (SNe), (2) continuous stellar mass-loss from OB and AGB winds, (3) photoionization and photoelectric heating, and (4) local and long-range momentum flux from radiation pressure. 
  All feedback quantities and their time dependence are tabulated directly from the stellar population synthesis model starburst99 <cit.> assuming a <cit.> initial mass function (IMF) and non-rotating, non-binary stellar populations.
  
  The radiative feedback implementation is described in detail in <cit.>.
  The age and metallicity-dependent IMF-averaged spectrum of each star particle is locally extinguished by the surrounding gas using a Sobolev approximation with frequency and metallicity-dependent opacities.  The emerging luminosities (including re-radiated dust emission) are propagated through a tree structure in the optically thin limit, yielding a long-range incident flux which is then corrected for local extinction at the location of the gas element.  We also include a uniform, redshift-dependent photo-ionizing background <cit.> with self-shielding.  The resulting incident flux is used to compute the ionization state of the gas, radiative heating/cooling rates, and radiation pressure from photon absorption including UV/optical single-scattering and re-radiated infrared photons. 
  
  The SNe feedback algorithm, fully described in <cit.>, is constructed to ensure conservation of mass, energy, and momentum as well as isotropic injection in the rest frame of the star particle regardless of the geometry and dynamics of the surrounding gas distribution (non-trivial in Lagrangian codes).  
  Each individual SN is treated independently, with IMF-averaged ejecta mass of 1.4  and 10.5  per SNe Ia and SNe II explosion, respectively, each with 10^51 erg of ejecta energy.
  The correct relative amount of momentum and thermal energy is injected depending on the resolved coupling distance relative to the SN cooling radius.  This algorithm reproduces converged solutions in both energy and momentum independent of resolution, which is crucial for the multi-mass resolution simulations presented here.  The same feedback implementation has been successfully used in cosmological zoom-in simulations with mass resolution ranging from m_ b∼ 20  in dwarfs <cit.> to m_ b∼ 3×10^4  in massive galaxies <cit.>.   
  At the highest resolution reached here (m_ b∼ 15 ), individual star particles are still massive enough to accommodate a full SN II mass ejection, but the same algorithm performs well down to significantly lower particle masses, conserving mass in a time-average sense <cit.>.
  Stellar winds follow a similar implementation but with continuous rather than impulsive injection of mass, momentum, energy, and metals into surrounding gas particles.
  The discrete nature of individual SNe is crucial in generating a multi-phase ISM and driving galactic winds, but IMF sampling on radiative feedback and stellar winds has weak effects on galaxy-scale properties <cit.>.
  
  
  
  
  
    §.§.§ Black hole physics
  
  
  We follow the growth of a massive black hole located at the center of the main halo in the high resolution region of a pre-existing cosmological zoom-in simulation from <cit.>, which implemented a sub-grid accretion model based on gravitational torques <cit.> evaluated on ∼100 pc scales.  The radical increase in resolution achieved by the hyper-Lagrangian refinement technique introduced here allows us to study black hole growth without the need for a sub-grid parameterization for the inflow rate from larger scales down to the accretion disk. 
  
    
  
  The black hole is represented by a collisionless particle with mass M_ BH which is 
  many orders of magnitude larger than the mass of high resolution gas and star particles that dominate the potential in the central ∼100 pc.
  This allows us to track the black hole dynamics explicitly owing to resolved gravitational forces without the need for a sub-grid dynamic friction prescription <cit.> or the artificial drag forces and re-positioning algorithms required in large volume cosmological simulations <cit.>.   
  The initial black hole mass in the hyper-refinement simulations is chosen to be M_ BH = 10^8  (see <ref> for simulations with different initial M_ BH), which is thus generally different from the evolving black hole in the pre-existing cosmological zoom-in simulation (<ref>).
  
  
  
  We model accretion by direct gravitational capture of gas by the central black hole. 
  Gas particles located within the accretion radius R_ acc are captured if (1) their velocity relative to the black hole is lower than the escape velocity and (2) the apocentric radius of the particle relative to the black hole is also within R_ acc.  We choose the accretion radius R_ acc≡ 0.1 pc so that it roughly corresponds to the scale below which the accretion disk is expected to form, R_ acc∼ 10^4 R_ s <cit.>, where R_ s = 2 G M_ BH/c^2 ≈ 10^-5 pc is the Schwarzschild radius for M_ BH = 10^8 .  We do not allow the black hole to accrete star particles.
  A similar explicit gravitational capture approach was used in <cit.> and <cit.> to study black hole growth in idealized nuclear scale simulations.  
  For simplicity, we refer to the measured gas inflow rate through R_ acc as “black hole accretion rate" or “Ṁ_ BH” but it should be interpreted as the instantaneous feeding rate onto the black hole accretion disk.  
  This represents an upper limit to the actual black hole growth owing to e.g. mass loss in accretion disk winds <cit.>.
  We neglect the effects of black hole feedback in an attempt to understand the mechanisms responsible for mass transport without making any assumptions about the efficiency of AGN feedback.
  
  
  
  
  
   §.§ Zoom-in simulations as initial conditions
  
  
  We use the FIRE-2 zoom-in simulations of massive halos (M_ halo∼ 10^12.5-13  at z=1) presented in <cit.> to identify interesting redshift snapshots for re-simulation at ultra-high resolution.  The initial conditions correspond to a subset of the A series of MassiveFIRE halos that were originally simulated with the FIRE-1 model <cit.>.  For this work, we focus on halo A4 (studied also in  and ), which was evolved from early times down to z=1 including an on-the-fly treatment of black hole growth  but not black hole feedback.  The mass resolution employed was m_ b = 3.3 × 10^4  and m_ DM = 1.7 × 10^5  for the baryonic and dark matter components.  The force softenings were ϵ_ gas^ min = 0.7 pc, ϵ_⋆ = ϵ_ BH = 7 pc, and ϵ_ DM = 57 pc, where ϵ_ gas^ min is the minimum adaptive force softening for gas (with ϵ_ gas identical to the kernel smoothing scale) and the softenings for the stellar (ϵ_⋆), black hole (ϵ_ BH), and dark matter (ϵ_ DM) components are fixed in physical units at z<9. 
  
  
  
  
  
  
  
  Figure <ref> shows the global stellar mass (M_⋆; top) and SFR (middle) as a function of redshift for the central galaxy of halo A4.  Throughout this paper, global galaxy properties refer to integrated quantities within 0.1 R_ vir of the host halo.  We identify halos using the Amiga Halo Finder <cit.>, modified such that the virial radius R_ vir corresponds to the evolving overdensity definition of <cit.>.  
  We also show the nuclear gas surface density (Σ_ gas; bottom) measured within the black hole kernel, which is defined to contain 256 gas resolution elements and has typical values in the range R_0 ∼ 10–100 pc.
  Global galaxy properties are computed from the snapshots available at ∼20–25 Myr intervals, but we take advantage of black hole-specific data outputs that recorded the physical conditions within R_0 (including Σ_ gas) for every timestep of the simulation, reaching time resolution dt∼ 100 yr.  
  
  
  The evolution of the central galaxy in halo A4 is qualitatively similar to that of halo A2 described in detail in <cit.>. 
  During the first ∼ 3 Gyr, the stellar mass grows to M_⋆∼ 5×10^10  while undergoing strong bursts of star formation reaching SFR ∼ 10–300 .  Correlated SNe drive large scale winds that can evacuate a large fraction of ISM gas and temporarily shut down star formation, while recycling of wind material provides fuel for subsequent bursts of star formation <cit.>.  At later times, star formation becomes less bursty and galactic winds less efficient <cit.>, with the galaxy growing more steadily to M_⋆≳ 2×10^11  down to z=1.
  This transition from bursty to steady star formation has a significant impact on the nuclear gas reservoir: Σ_ gas fluctuates by orders of magnitude at early times but there is a more steady nuclear gas reservoir at late times.  <cit.> showed that this has important implications for black hole growth, which can be significantly suppressed at early times relative to the host galaxy growth (see also ; ; Byrne et al. in prep.).
  
  We focus on three qualitatively distinct host galaxy conditions to study black hole feeding at ultra-high resolution, indicated in Figure <ref> by vertical lines of different colors: 
  
  
    * We identify z = 2.25 (orange) as the time at which galaxy A4 experienced the highest nuclear gas surface density, with Σ_ gas(<100  pc) > 10^11  and global SFR ∼200 .  These extreme conditions are the most likely to trigger gas inflow rates sufficient to power a luminous quasar, which we denote as “” phase.
  
  
  
    * We investigate the conditions leading to the  phase by going ∼40 Myr back in time to z = 2.28 (green), with lower nuclear gas density Σ_ gas(<100  pc) ∼ 10^10  but highly variable conditions reaching higher global SFR ∼ 500  (the peak SFR is missed in the parent cosmological simulation owing to limited output time resolution).  We denote these as the “" initial conditions. 
  
  
  
    * We compare the z>2 conditions near the peak of star formation with the late time conditions at z = 1.10 (blue) with significantly lower activity, where SFR ∼10  and Σ_ gas(<100  pc) < 10^9 .  These conditions are presumably representative of more common, low-luminosity AGN, which we denote as “".
  
  As we will show throughout the paper, the  phase reaches a quasi-steady state over the whole hyper-refinement simulation time while accretion in the  and  phases is highly time variable.
  
  
   
  lccccccccc
  1
  Simulation parameters (units are physical): 
  (1) Name: simulation designation. 
  (2) z: initial redshift of hyper-refinement.
  (3) M_ vir: halo virial mass.
  (4) M_⋆: galaxy stellar mass.
  (5) M_ gas: galaxy gas mass.
  (6) SFR: galaxy star formation rate.
  (7) Δ t: total duration (physical run-time) of the simulation at the full hyper-refinement level.
  (8) dt_ min: minimum timestep of black hole particle.
  (9) m_ g: minimum mass of gas resolution element.
  (10) h_ sml^ min: minimum gas smoothing length.
  0pt
  
  1lName    
  1cz    
  1cM_ vir []    
  1cM_⋆ []    
  1cM_ gas []    
  1cSFR []    
  1cΔ t [Myr]    
  1cdt_ min [yr]    
  1cm_ g []    
  1ch_ sml^ min [pc]
  
       2.28      2.31e12       5.90e10       2.15e10       523        3         0.15      15      0.08 
  
       2.25       2.37e12       8.21e10       1.87e10       222       4       0.42       30      0.08 
  
       1.10       4.44e12       2.31e11       9.63e9         13       20       2.55        15      0.15 
  
  
  
   
  
  
  
   §.§ Hyper-Lagrangian refinement
  
  
  
  
  
  
  Lagrangian hydrodynamic methods are naturally adaptive, populating high density regions with a larger number of gas resolution elements.  In order to achieve super-Lagrangian refinement in MFM, we split gas resolution elements progressively as they approach the central black hole.
  The distance R from each gas particle to the black hole is computed while traversing the gravity tree 
  <cit.>, which allows us to define a mass refinement factor χ_ ref(R) as an analytic function of R.
  Gas resolution elements split into two new elements if m_ g≥χ_ ref(R)× 2  m_ b for their current mass m_ g and distance R, where m_ b is the baryonic mass resolution in the original zoom-in simulation.  The mass of each new split element is 0.5 m_ g and they are located at a distance dr =  min(0.25 h_ sml, 0.35 d_ ngb) from the position of the parent resolution element, in opposite sides along a random direction.  Here, h_ sml is the smoothing length, d_ ngb is the distance to the nearest neighbor, and the numerical factors are chosen to minimize perturbations to hydrodynamic quantities and avoid overlapping of fluid elements.    
  The velocity of split elements is identical to that of the parent gas element, conserving total energy, momentum, and angular momentum in the splitting process.  We have also experimented with a first-order correction according to the existing velocity gradient, providing a more accurate reconstruction of the velocity field while introducing second order errors in the energy and angular momentum; we find no significant differences.  All other fluid quantities such as density and internal energy per unit mass remain unchanged, but are recomputed immediately after particle splitting following the standard formulation in MFM hydrodynamics for the new configuration of gas elements <cit.>.
  Split elements are allowed to merge onto higher mass elements when they are significantly under-massive relative to the resolution requirement, if m_ g < χ_ ref(R)× m_ b / 200, which is infrequent for the simulations presented here.
  Our fiducial simulations only refine the gas component but we have also implemented splitting/merging of collisionless particles (see Appendix <ref> for details). 
  In either case, newly formed stars have the same mass as the progenitor gas element regardless of the level of resolution.
  
  
  Figure <ref> shows the mass of gas resolution elements as a function of R obtained for separate simulations using different levels of refinement for our fiducial χ_ ref(R) functional form:
  
      χ_ ref = 
      χ_ ref^ min× (1 - R/R_ ref) + (R/R_ ref)^2     if R ≤ R_ ref
       
       1    if R > R_ ref
  
  where χ_ ref^ min is the desired refinement factor near the black hole, which we vary from χ_ ref^ min = 1/2 (yellow) to χ_ ref^ min = 1/2048 (purple), and R_ ref = 1 kpc is the distance at which super-Lagrangian refinement begins to operate.  
  The original baryonic mass resolution m_ b = 3.3 × 10^4  is retained at R>1 kpc, decreasing rapidly down to m_ g≈χ_ ref× m_ b closer to the black hole.  In our highest resolution simulations we reach m_ g≈ 15 , which allows us to follow the inflowing gas well within the black hole radius of influence and down to the accretion aperture R_ acc = 0.1 pc. 
  Despite the strong gradient in resolution at R =0.1–1 kpc, we do not see any significant artifacts owing to dynamical splitting and merging of gas elements, except for a short initial transient that we discard in the analysis.  We have experimented extensively with increasing the mass resolution slowly as the simulation evolves and also refining directly down to the highest resolution level; while the former minimizes the perturbations introduced, both methods yield similar long term behavior.  
  We have also experimented with different radial functional forms for χ_ ref(R), but the flexibility of this particle splitting technique allows for hyper-refinement with any arbitrary geometry (see  for hyper-refinement of magnetohydrodynamic jets propagating in a cluster environment).
   We find that χ_ ref∝ R^2 provides a convenient mapping between standard resolution at 1 kpc and χ_ ref^ min such that we reach uniform resolution at ≲100 pc.
  
  
  
  Table <ref> summarizes the main properties of our fiducial simulations.  The initial conditions are taken at each of the redshift snapshots corresponding to the , , and  conditions, where we remove all black hole particles other than the one located at the center of galaxy A4.  The actual black hole mass is set to M_ BH = 10^8  such that we can compare the inflow rate under different conditions for the same M_ BH (we test different M_ BH values in additional simulations; see <ref>).  
  We use the same gravitational softenings for dark matter, ϵ_ DM = 57 pc, but use fully-adaptive softenings for gas (matching the hydrodynamic solver) with an arbitrarily small minimum allowed (but in practice minimum softenings reached are given in Table <ref>), and treat the black hole correctly as a Keplerian point-like particle (using a much smaller softening, ϵ_ BH=0.001pc, intentionally smaller than the “accretion radius”). We adopt ϵ_⋆ = R_ acc = 0.1pc for stars, motivated roughly by the softening of gas at the densities where stars form in the densest hyper-refined regions, but in Appendix <ref> we show that our results are robust relative to changes in ϵ_⋆, including adaptive and variable softenings.
  Hyper-refinement operates within R_ ref = 1 kpc in all simulations presented here and we use fiducial maximum refinement factors χ_ ref^ min = 1/1024 ( phase) and χ_ ref^ min = 1/2048 ( and  phases).  Additional simulations with lower refinement levels are employed for resolution convergence tests (<ref>).  
  We save data snapshots for every 0.01 Myr of evolution in all simulations.
  Except for the radical increase in resolution and the treatment of black hole growth, our cosmological hyper-refinement simulations implement identical physics as the original FIRE-2 simulation.
  
  
  We reach minimum timesteps ∼0.1 yr and evolve our hyper-refinement simulations for Δ t ∼ 3–20 Myr (Table <ref>).  
  The initial conditions considered have characteristic dynamical times t_ dyn≡ R/v_ c∼ 1–2  Myr→ 0.1–0.2 Myr given the circular velocity v_ c at R∼ 1  kpc→ 100 pc scales.  We are thus evolving the simulations for ∼2–20 dynamical times at the scale where hyper-refinement starts (1 kpc) and ∼20–200 dynamical times in the nuclear region (≲100 pc) within which the maximum resolution is achieved, with the inner 1–10 pc evolving for hundreds to thousands of dynamical times.
  Our simulations are thus reaching a reasonable statistical steady state in which the high resolution region matches the lower resolution conditions on larger scales in a physically meaningful way.  In addition, we find characteristic flow times t_ flow≡ M_ gas / Ṁ_ in∼ 6–22 Myr given the enclosed gas mass and inflow rate at 1 kpc (<ref>), roughly similar to our total evolution times, which emphasizes the importance of including larger scale gas flows when evolving the nuclear region over many dynamical times.
  
  
  
  
  
  
  § OVERVIEW OF SIMULATIONS
  
  
  
  
  
   §.§ From zoom-in initial conditions to hyper-refined galactic nuclei
  
  
  
  
  
  
  Figure <ref> shows the gas distribution in the central 100 pc for the three initial conditions considered, comparing the original resolution from the parent cosmological zoom-in simulation (left) and that obtained in the cosmological hyper-refinement simulations (right).  
  Each pair of simulations is shown at the same evolution time, corresponding to ∼0.6 Myr (), ∼3.8 Myr (), and ∼6.1 Myr () after reaching mass resolution m_ b∼ 30  in the hyper-refinement simulation (we reach up to one resolution level higher in the  and  conditions; Table <ref>).
  The top panels correspond to the  conditions at z=2.28, with average Σ_ gas∼ 10^10  in the inner 100 pc.  The initial conditions (top left) already capture a remarkable level of sub-structure for a massive galaxy simulated in a full cosmological context, with a clumpy, dense, circumnuclear ring and a hint of filamentary structures that appear to connect with a central gas concentration.  The hyper-refinement simulation (top right) confirms the presence of the ring-like structure in the nuclear region, which is now resolved by N_ gas∼ 2× 10^6 resolution elements 
  into a complex distribution of clumps, filaments, and spiral structures.
  
  
  The middle panels of Figure <ref> show the nuclear gas distribution for the  conditions at z=2.25, characterized by an extreme gas concentration reaching Σ_ gas > 10^11  in the central 10 pc.  The parent zoom-in simulation suggests the presence of up to four nuclear spiral arms while the ×1000 increase in mass resolution yields a variety of morphological features with more than three orders of magnitude in density contrast, including an ultra-dense eccentric disc within ∼ 10 pc that breaks up into multiple spiral arms.
  Finally, the bottom panels correspond to the  conditions at z=1.1, where the parent zoom-in simulation shows a nearly homogeneous gas distribution with average Σ_ gas≲ 10^9  in the inner 100 pc.  The cosmological hyper-refinement simulation yields a remarkably different gas distribution, with a collection of dense clumps and filaments embedded in a low density medium and the formation of a central cavity extending a few tens of pc.  
  
  
  
  
  
   §.§ Mapping gas from Mpc to pc scales
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  Figure <ref> illustrates the full dynamic range of our simulations by showing the gas distribution for the  conditions on different scales, from 1 Mpc to the central 10 pc.
  The top left panel shows the full extent of the zoom-in region, which is populated with high resolution dark matter and gas elements and embedded in a low-resolution volume of [100 ]^3 that allows us to capture tidal torques from large scale structures.  High density peaks correspond to the location of galaxies, tracing large scale filaments.  The white dashed circle in the center indicates the virial radius of the main halo, with R_ vir = 130 kpc and  M_ vir = 2.3×10^12  at z=2.28.  Zooming into the central 100 kpc, the CGM shows a complex gas distribution owing to the interplay between filamentary accretion from the IGM, cooling of hot halo gas, and galactic winds from the central galaxy and satellites <cit.>.  A merging satellite galaxy (M_⋆∼ 10^10 ) is approaching the central (M_⋆ = 5.9×10^10 ), currently on its second passage at a distance of ∼ 20 kpc and playing a minor role in driving instantaneous black hole growth (see <ref>).  
  
  
  On 10 kpc scales, the ISM shows a very irregular morphology, with gas surface densities in the range Σ_ gas∼ 10^7-10 .
  The highest densities are reached within the inner 1 kpc in a clumpy ring-like structure of radius ∼400 pc, with a secondary gas concentration forming in the inner 100 pc.  The black hole responds dynamically to the turbulent clumpy gas distribution by orbiting up to 10 pc away from the center of mass of the stellar component.  Gravitational capture of infalling clumps and filaments yields pc-scale gas disks around the black hole with changing orientation depending on the angular momentum of the infalling material.   
  
  
  Figure <ref> shows projected mass-weighted temperature distributions from 1 Mpc to 10 pc scales for the  conditions at z=2.25.  In this case, the Mpc scale view in the top left panel highlights the thin, large scale filaments of cool gas (T<10^5 K) feeding the central halo and the expanding hot gas (T>10^6 K) heated by accretion shocks and galactic winds, reaching well beyond  R_ vir≳ 130 kpc.  On 100 kpc scales, the CGM shows a prominent multi-phase structure with cold clumps and filaments embedded in a hot medium of virialized gas (T>10^6.5 K).  The projected temperature distribution in the central 10 kpc is highly irregular, with a large covering fraction of cold, star forming gas.  Dust continuum radiative transfer calculations show that these conditions produce a submm-bright phase <cit.>.
  Zooming into the inner 1 kpc, cold gas with low angular momentum falls toward the inner region while nuclear spiral arms form out of the turbulent ISM.  In the central 100 pc, the simulation captures ∼ 4 orders of magnitude variation in projected mass-weighted temperature, with prominent spiral structures of T∼10^4-5 K gas embedded in T > 10^7 K gas heated by SNe.  The bottom left panel highlights the temperature contrast in the inner 10 pc, where a cold (T∼ 1000 K), pc-scale, rotationally supported gas disk forms around the black hole.  
  See Figure <ref> for the corresponding gas mass surface density distributions in the  phase.
  
  
  
  Figure <ref> shows again projected gas surface density distributions on multiple scales, as in Figure <ref>, but in this case for the  conditions at z=1.1.  The main halo has doubled in mass since z∼ 2.3, with M_ vir = 4.4×10^12  and R_ vir∼ 250 kpc, and the central galaxy has developed a massive stellar component (M_⋆ = 2.3×10^11 ).  The CGM on 100 kpc scales shows a smoother gas distribution compared to earlier times, though some dense gas structures are still prominent.  With the gas fraction within 0.1 R_ vir falling below f_ gas < 5% and SFR∼ 10 , the galaxy forms a thin, rotationally supported disk in the inner ∼ 2 kpc.  A low-density contrast, bipolar structure appears to extend a few kpc above and below the plane of the disk, likely driven by an earlier episode of galactic winds.  Zooming into the central 1 kpc, the disk structure is dominated by star-forming regions with Σ_ gas > 10^9  embedded in a low density medium, with a compact spiral structure superimposed.
  In the central 100 pc, gas consumption by star formation and ejection by stellar feedback form a cavity that extends a few tens of pc and starves the black hole of fuel <cit.>.  Gravitational torques from the stellar component drive gas clumps and filaments from the inner boundary of the circumnuclear disk down to the central 10 pc, occasionally feeding the black hole (see <ref>). 
  
  
  
  
  
   §.§ Galaxy radial structure and kinematics
  
  
  Figure <ref> shows the stellar, gas, and star formation rate surface densities as a function of cylindrical radial distance R_ cyl for the three conditions simulated.  We define a cylindrical coordinate system (centered on the black hole) independently for each time and radial bin, where the z-axis is aligned with the total angular momentum inside of the 3D radial distance R = R_ cyl to the black hole, also performed independently for the gas and stellar components.  Surface densities Σ_⋆, Σ_ gas, and Σ_ SFR are then computed within these local cylindrical radial bins (logarithmically spaced) by integrating the mass/SFR within ±100 pc along the vertical direction relative to the local cylindrical plane.  This allows us to capture spatial and temporal variations in the structure and kinematics of the galaxy, but our results are not sensitive to the exact choice of local radial bins, their vertical extension relative to the local galaxy plane, or the use of spherical rather than cylindrical radial bins. 
  Note that the black hole and therefore the adopted reference frame can move relative to the center of mass of the galaxy  (Figure <ref>), with subsequent analysis reflecting the instantaneous conditions driving black hole growth. 
  We compute median surface density profiles considering the full time evolution of each simulation.  In addition, we extract the 10–90% range of variation including only timesteps with non-zero values in each radial bin, which allows us to study infrequent events such as the transport of gas clumps across the central cavity that forms in the  phase.
  
  
  
  
  
  The  conditions at z=2.28 (green) show a significant increase in stellar mass surface density Σ_⋆∼ 10^8 → 10^11  in the range R_ cyl = 1 → 0.1 kpc, with a mild increase to Σ_⋆∼ 10^11.7  at smaller scales.  The initial stellar distribution (dashed line) is very similar to the median Σ_⋆ found during the ∼ 3 Myr time evolution in the hyper-refinement simulation (solid line), but it does not extend to within the inner 1 pc.  
  New stars formed out of ultra-high resolution gas populate the inner ∼1 pc region, but the stellar distribution on larger scales does not change significantly over the short ∼3 Myr simulation period.
  The shaded region indicates the 10–90% range for each radial bin, reflecting the build up of stellar mass within 1 pc, reaching Σ_⋆∼ 10^12.5 , but also variability in Σ_⋆ owing to black hole dynamics relative to the stellar distribution. 
  The  phase at z=2.25 (orange), ∼40 Myr later, shows very similar Σ_⋆ at >100 pc but significantly higher values in smaller scales relative to the  conditions.  The initial  conditions had already reached Σ_⋆∼ 10^13  within ∼10 pc, increasing to even more extreme densities as new ultra-high resolution stars form in the inner ∼1 pc.
  Moving forward to z=1.1, the  conditions (blue) show a more extended stellar distribution owing to size growth and expansion of the inner region.  Within ∼100 pc, Σ_⋆ is roughly ×10 lower compared to the  conditions but the original stellar distribution retained the central core with Σ_⋆∼ 10^13 , further increasing during the ∼20 Myr hyper-refinement simulation period.
  The extreme stellar densities reached in the  and  conditions are likely unphysically large owing to the lack of AGN feedback.    
  
  
  The middle and bottom panels of Figure <ref> show the radial gas and SFR distributions. 
  The highest gas surface densities on galactic scales (>100 pc) are reached in the  conditions, with Σ_ gas≳ 10^9  in the central 1 kpc and Σ_ SFR∼ 10 →1000  increasing inward.  Within ∼1 pc around the black hole, the  phase maintains a median Σ_ gas≳ 10^10  over the ∼3 Myr simulation period, with variations of over two orders of magnitude, and Σ_ SFR∼ 10^5–10^7 . 
  The  phase shows lower Σ_ gas on galactic scales but more steady and significantly larger gas and SFR surface densities (∼0.5 dex) at <10 pc, with median Σ_ gas∼ 10^10.5  and Σ_ SFR∼ 10^6  in the central 1 pc.
  The  conditions show median Σ_ gas≳ 10^8 on galactic scales but a rapid decline at <100 pc owing to the formation of a low-density central cavity.  Nonetheless, gas transported down to <1 pc due to instabilities in the outer edge of the cavity can reach Σ_ gas > 10^9 , with star formation happening all the way to the black hole accretion radius.
  
  
  
  
  
  
  Figure <ref> investigates the kinematic and thermal properties of gas from R_ cyl = 0.1 pc–10 kpc for our  (green),  (orange), and   (blue) simulations.  The top two panels show the average (mass-weighted) temperature (T) and sound speed (c_ s) within each cylindrical radial bin, extending only ±10 pc along the vertical direction relative to the local plane.  As in Figure <ref>, solid lines show median values for the full hyper-refinement simulation time and shaded regions indicate the 10–90% range variation in each radial bin. 
  The  conditions are dominated by cool gas (T≲ 10^5 K) in the entire range R_ cyl = 0.1 pc–10 kpc owing to the large amount of dense gas and thus rapid cooling times.  With average c_ s∼ 10 , thermal pressure support does not play a relevant dynamical role.  Nonetheless, heating by SNe following intense star formation can sometimes increase the temperature to T > 10^6 K in the inner ∼ 100 pc; note that the black hole can also move from cold to hot dominated regions as it orbits within the central 10 pc.  The thermal properties in the  conditions are rather different, with average T∼ 10^6 K on galaxy scales and T≳ 10^7 K consistently in the central ∼10–100 pc, where there is a stark contrast between the temperatures reached by the cold and hot ISM phases (see Figure <ref>).  Despite this, gas within ∼1 pc of the black hole is predominantly cold.  By the time we reach the  conditions at z=1.1, the galaxy has already developed a hot halo and we find T ∼ 10^6.5 K at ∼10 kpc scales, dropping to T ∼ 10^4 K within the ∼ kpc scale gas disk.  At <100 pc, the low-density central cavity seen in Figure <ref> reaches T ∼ 10^8 K and c_ s∼ 1000 , while gas that penetrates the cavity down to the black hole accretion radius is predominantly cold. 
  
  
  
  
  
  
  The three bottom panels of Figure <ref> show the radial dependence of the vertical velocity dispersion (σ_ z), the radial velocity dispersion (σ_ r), and the azimuthal velocity (v_ϕ) in cylindrical coordinates for the gas component.   
  Velocity dispersions are defined as σ_ z^2 ≡⟨ v_ z^2 - ⟨ v_ z⟩^2⟩ and σ_ r^2 ≡⟨ v_ r^2 - ⟨ v_ r⟩^2⟩.
  The  conditions show remarkably flat velocity dispersion profiles, with σ_ z∼ 100  in the entire range R_ cyl = 0.1 pc–10 kpc and slightly higher dispersion σ_ r in the radial direction <cit.>.  The galaxy-scale azimuthal velocity reaches v_ϕ∼ 400  at ∼ 100 pc–1 kpc.  In the central 10 pc, the azimuthal velocity increases to >1000  with v_ϕ∝ R^-1/2, as expected for Keplerian rotation within the black hole radius of influence.       
  With increased stellar and gas mass surface densities, the  conditions show even higher levels of turbulence.  In this case, both σ_ z and σ_ r  increase from ∼ 200  on galaxy scales to ∼ 500  at ∼ 0.1 pc.  Nonetheless, rotational motion still dominates with v_ϕ > 600  at R_ cyl<100 pc.
  At later times, the rotational velocity in the  conditions reaches v_ϕ > 800  on kpc scales owing to the build up of a massive stellar component.  The gas distribution shows a rotationally supported disk with ∼ 50  turbulent motions.  The cool gas that penetrates the inner hot cavity (blue shaded region within ∼10 pc) follows a roughly Keplerian rotation curve with v_ϕ∝ R^-1/2.
  
  Overall, despite the high level of turbulent motions in the  and  conditions and the role of thermal pressure in the  central cavity, rotational support appears to dominate in the full range of scales analyzed, with 
  c_ s < σ_ z∼σ_ r < v_ϕ.   
  The  conditions reach extreme stellar densities (Figure <ref>), higher than the maximum stellar surface density predicted by analytic models due to the failure of stellar feedback <cit.>.
  In addition, the large v_ϕ velocities reached in the  conditions are comparable only to the most extreme compact galaxies observed <cit.>, suggesting the need for AGN feedback or some other additional mechanism to suppress nuclear star formation <cit.>.
  We investigate next if large inflow rates down to sub-pc scales can occur, a prerequisite for efficient “QSO-mode” black hole feedback.
  
  
  
  
  
  
  
  § MULTI-SCALE GAS INFLOWS
  
  
  
  
  
   §.§ Explicit accretion at <0.1 pc
  
  
  
  Figure <ref> shows Ṁ_ BH as a function of time measured at 0.1 pc in hyper-refinement simulations for the  (z = 2.28; green),  (z = 2.25; orange), and  (z = 1.10; blue) conditions, considering in all cases a black hole with initial mass M_ BH^ ini = 10^8 . 
  For context, the middle panel shows the black hole accretion history from z = 10 → 1 obtained by post-processing the original cosmological zoom-in simulation adopting the “free-fall" sub-grid accretion estimator Ṁ_ BH = α  M_ gas/t_ dyn, where t_ dyn = (R_0^3/GM_ tot)^1/2 is the dynamical time within R_0∼100 pc and M_ gas and M_ tot are the gas mass and total mass within R_0.  We define α≡ 5×10^-4 such that M_ BH/M_ bulge∼ 0.002 at z=1, which results in accretion rates roughly consistent with the gravitational torque accretion estimator used in <cit.>, but the characteristic accretion history obtained is independent of this choice.
  In the absence of AGN feedback, the sub-grid accretion model predicts short accretion episodes that can reach or even exceed the Eddington rate at early times (z > 2.3) and a transition to more steady accretion at late times, resembling the evolution in nuclear gas surface density shown in Figure <ref>. 
  
  
  
  
  
  
  The  conditions at z = 2.28 target the epoch when the bursty-to-steady transition is happening.  The top left panel of Figure <ref> (green) shows high level of variability in accretion when measured at 0.1 pc, with Ṁ_ BH = 0.001–10  and average ⟨Ṁ_ BH⟩ = 0.6  during a period of ∼3 Myr.  For M_ BH^ ini = 10^8 , this implies an average Eddington ratio ⟨λ_ Edd⟩≈ 0.3, where λ_ Edd≡Ṁ_ BH/Ṁ_ Edd and Ṁ_ Edd≈ 2.2  is the Eddington rate. 
  The  phase at z = 2.25 (orange panel; top right) coincides with the peak of black hole accretion according to the free-fall sub-grid accretion model, occurring ∼40 Myr after the  conditions.  In this case, the hyper-refinement simulation predicts ⟨Ṁ_ BH⟩≈ 6   with very little variability during a period of ∼4 Myr, corresponding to ⟨λ_ Edd⟩≈ 2.6.
  This implies that quasi-steady super-Eddington feeding rates onto the accretion disk are possible in the absence of AGN feedback.
  The  conditions at z = 1.1 (blue panel; bottom) show significantly lower accretion rate, with ⟨Ṁ_ BH⟩ = 0.01  during a period of ∼20 Myr, corresponding to ⟨λ_ Edd⟩≈ 0.005.  The formation of a low-density central cavity filled with T∼ 10^7 K gas (Figure <ref>) inhibits black hole growth during ∼75% of the time (with Ṁ_ BH = 0), while the instantaneous accretion rate ranges from Ṁ_ BH = 0.001–0.1  when cold gas clumps and filaments destabilize from the inner edge of the circumnuclear disk and fall toward the black hole (Figure <ref>).  The finite mass of gas resolution elements (m_ g∼ 15 ) implies an upper limit of roughly Ṁ_ BH < 10^-6  during the inactive periods in the  phase.
  
  
  Overall, our hyper-refinement simulations predict qualitatively distinct accretion phases with properties that are not trivial to infer even from the ∼100 pc physical conditions available in state-of-the-art cosmological zoom-in simulations.  Remarkably, we show that it is possible to feed the central 0.1 pc at a rate sufficient to power a luminous QSO (L_ bol≈ 3×10^46 erg s^-1 for Ṁ_ BH = 6  and ϵ_ rad=0.1) during the epoch of peak nuclear gas density in the central galaxy of a massive halo (M_ vir∼10^12.5 ) at z∼2.  At the same time, our simulations predict strong variability in the conditions preceding the QSO phase and low-duty cycle and lower luminosity for the more common conditions prevalent at later times.
  
  
  
  
  
  
  
  
   §.§ Star formation–inflow connection
  
  
  
  Figure <ref> investigates the connection between gas inflow rate and SFR over five orders of magnitude in spatial scales.  For each hyper-refinement simulation, we compute the net mass flow rate across spherical shells and show the time-averaged radial profile as either solid or dotted lines for net inflow and outflow rate, respectively.  
  Note that large-scale, coherent galactic outflow events driven by stellar feedback are common at higher redshift but absent in the conditions analyzed here, where inflows and outflows can coexist in the same radial domain.
  Colored dashed lines show the time-averaged enclosed SFR as a function of radial distance, computed from the instantaneous SFR of individual gas elements, and shaded regions indicate the 10–90% variation in each radial bin.
  
  In the  conditions (top), inflows dominate over outflows in the entire radial range R = 0.1 pc–10 kpc, with the peak inflow rate reaching Ṁ_ in∼ 475  at ∼1 kpc and a monotonic decrease to Ṁ_ in≲ 1  at ∼0.1 pc.  The enclosed SFR shows significant variability occurring during only ∼3 Myr, indicated by the green shaded region, with the global peak of star formation reaching ∼1000 .  As expected, the relative variation in enclosed SFR increases at lower radii, owing to the shorter dynamical times and the non-trivial motion of the black hole relative to the surrounding gas distribution.
  Despite this, the time-averaged enclosed SFR profile roughly matches the net inflow rate in the entire range 0.1 pc–10 kpc. 
  
  The  conditions (middle) display similar large scale inflow rates as the  simulation, but show net outflow in the intermediate range ∼50 pc–1 kpc (coincident with the maximum temperatures reached) and significantly higher net inflow rates in the inner region.  The enclosed SFR is again roughly similar to the net inflow rate on ≳1 kpc and ∼10 pc scales, but decreases rapidly in the inner 1 pc.   As expected, the global inflow rate and SFR are significantly lower in the  conditions (bottom), but the inflow rate profile again roughly follows the enclosed SFR distribution.  In this case, outflows dominate in the range ∼10–100 pc owing to the formation of a hot, low-density central cavity, but cold gas clumps driven inward by instabilities yield a net inflow rate of ∼0.01  at 0.1 pc while feeding star formation at a rate ∼100 times lower.  
  Overall, these results show that star formation can consume as much gas as provided by inflows at all scales, with some indication for comparatively lower star formation relative to the inflow rate in the central ∼1 pc.
  
  
  Figure <ref> shows more directly the relation between black hole accretion and radially integrated SFR on different scales, where we plot the time-averaged ratio Ṁ_ BH/SFR as a function of radial distance.  Ṁ_ BH can greatly exceed the enclosed SFR for the smallest scales probed here, reaching Ṁ_ BH∼ 100×SFR at 0.1 pc for the  and  conditions.  In all cases, the SFR within ∼10 pc already exceeds the accretion rate, and Ṁ_ BH/SFR decreases for the global (∼1–10 kpc) star formation rate down to 1/1000 in the  and  conditions and 1/50 for the  phase.  
  These accretion to global SFR ratios have important implications for black hole–galaxy scaling relations (<ref>).  
  
  
  
  
  
  
  
  
  
   §.§ Cool gas inflow and angular momentum
  
  
  Figure <ref> shows the fractional contribution of cool gas (T < 10^4.5 K) to the total inflow rate as a function of radial distance (solid lines).  In the  conditions, cool gas comprises almost the totality of inflowing gas within 10 kpc.  For the  conditions, cool gas contributes >90% of the inflow rate within 1 kpc, dropping to ∼85% at 10 kpc.  Similarly, the cool gas component represents ∼90% of the inflow rate within 1 kpc for the  conditions, while in this case hot gas clearly dominates the inflow rate at 10 kpc.   
  As expected, hot gas (T > 10^4.5 K) contributes a higher fraction of the mass outflow rate (dashed lines), which clearly dominates over the cool component at R>1 kpc for the  and  conditions and for the range 10 pc–100 pc in the  phase. 
  Thus, while the mass-weighted temperature can exceed 10^6 K within the central 100 pc (top panel of Figure <ref>), fueling of the central black hole relies primarily on the supply of cool gas in the nuclear region for all three conditions analyzed, with thermal pressure support playing a minor role.
  
  
  
  
  
  
  
  
  
  Figure <ref> shows the radial profile of the specific angular momentum of cool gas, j_ gas, normalized by the specific angular momentum assuming circular motion, j_ circ = √(G M_ enc(R) R), where M_ enc(R) is the total enclosed mass within R.  We compute j_ gas as the magnitude of the angular momentum of cool gas in each spherical radial shell divided by its cool gas mass content, averaging over time for each simulation.  In all cases, j_ gas varies by ∼4 orders of magnitude from kpc to sub-pc scales, roughly tracing j_ circ.
  The  conditions show nearly circular motion within ∼1 kpc, roughly the size of the thin gas disk shown in Figure <ref>, but also show strong variations of j_ gas at <100 pc owing to the formation of a hot central cavity where cold gas clumps and filaments penetrate with ∼25% duty cycle. 
  The  conditions show the lowest j_ gas relative to j_ circ at 10 pc–1 kpc, but all cases show substantial angular momentum support, with j_ gas∼ 0.8 j_ circ at <10 pc, indicating the need for torques to remove angular momentum and drive gas inflows down to smaller scales.
  
  
  
  
  
  
  Figure <ref> shows the angle θ between the angular momentum vector of cool gas in spherical radial shells and that of the angular momentum of cool gas in a thin spherical radial shell centered at R=1 kpc.  
  The  conditions show very good alignment across scales, with average θ∼ 0^∘ at 0.1–1 kpc corresponding to a well defined, rotationally supported disk.  On smaller scales (<100 pc), cold gas clumps penetrating the central cavity tend to maintain the larger scale angular momentum direction but can be occasionally torqued into other directions all the way to θ > 90^∘.   
  The gas rich, turbulent conditions in the  phase show a very different situation, with a persistent misalignment of θ∼ 70^∘ between the galaxy-scale disk and the nuclear gas disk at <100 pc.  The strong variability in the inner ∼10 pc manifests as rapid changes in angular momentum direction, covering the full range 0^∘ < θ < 180^∘.   
  The  conditions show less time variability but an increasing misalignment from kpc scales all the way to counter-rotation on pc-scales.
  Thus, despite the similarity between the radial profiles of j_ gas and j_ circ, the angular momentum direction can change significantly across scales.  The pc-scale gas disk ultimately responsible for feeding the black hole can form in any direction relative to the galaxy scale disk.  
  
  
  
  
  
   §.§ Gravitational torques and pressure gradients
  
  
  
  
  
  
  
  Figure <ref> illustrates the magnitude and time variability of the total torque exerted on the cool gas component (T < 10^4.5 K) for the  conditions.  
  We compute the average specific torque in each spherical radial bin as:
  
      τ_j = 1/M_ bin ∑_i m_i r_i× a_j( r_i),
   
  where M_ bin = ∑ m_i, the sum over i refers to the cool gas elements with mass m_i inside the radial bin, and a_j( r_i) is the acceleration exerted by different components on gas element i located at r_i.  We compute separately the torque due to gravitational forces from gas, stars, or dark matter, where a_j = ∇Φ_j( r_i) and Φ_j is the gravitational potential for each component, and the torque from pressure forces, where a_j = ∇ P( r_i) /  ρ( r_i) and P and ρ are the gas pressure and density, respectively[Torques from gravitational forces and pressure gradients are computed in post-processing using M. Grudić's pykdgrav and meshoid Python repositories available at <https://github.com/mikegrudic>.].  
  We then compute the total torque per unit mass in the direction of the angular momentum L in each radial bin, τ_ tot = [ τ_ grav + τ_ hydro ]_L̂, where τ_ grav = τ_ gas +  τ_ star + τ_ DM.
  
  
  The dotted lines of different colors in Figure <ref> show the radial dependence of the specific torque τ_ tot at three randomly chosen times during the  phase.  We find very strong instantaneous specific torques on the cool gas component that can exceed 1000 km^2 s^-2, with significant time variability and sign changes over the entire radial range 0.1 pc–1 kpc. The black solid line indicates the time-averaged specific torque τ_ tot, which shows that there are also persistent sign changes for the radial profile in a time-average sense.  The transition between negative to positive average specific torque at ∼ 20 pc roughly corresponds to the transition seen in Figure <ref> from inflow dominated dynamics to outflow dominated for the  conditions.
  Qualitatively similar results are found for the  and  conditions.
  
  
  
  
  
  
  
  Figure <ref> compares the contribution of gravitational forces from other gas (blue), stars (red), and dark matter (black), as well as pressure forces (green) to the total specific torque on the cool gas component as a function of radial distance.  Note that fractional contributions can exceed unity (|τ_ j / τ_ tot| > 1) when different components contribute to τ_ tot with opposite sign.
  The top panel shows that gravitational torques from the stars dominate the total torque for most radii in the  conditions.  Dark matter contributes the least to gravitational torques, followed by gas self-gravity which can become comparable to the stellar torques at galactic scales (>500 pc).  Torques from pressure forces are typically 1–2 orders of magnitude lower than stellar torques but increase in the inner 1 pc to the extent that pressure gradients can become the dominant mechanism of angular momentum transport at ∼ 0.1 pc.  However, local viscous stresses driven by MHD turbulence, not captured here, are believed to dominate angular momentum transport in the accretion disk around the central BH <cit.>.  The  conditions (middle panel) show qualitatively similar results, with stellar torques clearly dominating in the central 100 pc and, in this case, pressure gradients becoming important at larger scales.  The  conditions (bottom panel) show very large scatter in the fractional contribution of different torque components owing to the formation of a low density central cavity.  Nonetheless, stellar torques are again the dominant component. 
  
  
  
  Figure <ref> shows the fractional contribution of stars at different radii to the total stellar gravitational torque.  We compute the stellar torque on the gas located in each radial bin as in Equation <ref> but considering only the gravitational potential produced by stars located in an independent set of radial annuli.  
  The bottom panel shows a very clear connection between the radial distribution of gas and the location of the stars that dominate the gravitational torque for the  conditions.  The torque in the central 1 pc is dominated by stars within 1 pc, the torque on gas at 1 pc < R < 10 pc is dominated by stars in the same radial bin, and similarly for radial annuli on larger scales.  Stars located at R>100 pc contribute only ∼ 10% of the total torque at 10 pc and <1% at 1 pc, suggesting that the instantaneous gas inflow rate at <1 pc is decoupled from the larger scale stellar distribution.
  The middle panel shows qualitatively similar results for the  conditions, with the main difference being that the torque within 1 pc < R < 10 pc is dominated by stars at 10 pc < R < 100 pc.
  The  conditions (top panel) show enhanced coupling of scales due to strong large scale non-axisymmetries and lower stellar surface densities in the central (<10 pc) region.  While the overall trend is similar, the total torque at R<1 pc contains significant contributions from stars in the entire radial range 0.1 pc–1 kpc, suggesting that the instantaneous accretion rate in this case is sensitive to the large scale stellar distribution.  This result would be difficult to capture without our hyper-refinement technique enabling us to cover a large dynamic range in scales and should be considered in the development of improved sub-grid accretion models.
  
  
  
  
  
  Figure <ref> illustrates the non-axisymmetries in the stellar component that yield the multi-scale gravitational torques quantified in Figure <ref> (top panel) for the  conditions.  We make stellar mass surface density maps (Σ_⋆) on different scales, oriented face-on relative to the angular momentum of gas within 100 pc, and compute the corresponding axisymmetric component (Σ_⋆, sym) by azimuthally-averaging Σ_⋆ in cylindrical radial bins around the black hole.  We show the ratio Σ_⋆ / Σ_⋆, sym from 1 kpc (left) to 10 pc (right), which highlights the strong non-axisymmetries present in all scales.   Spiral arms with prominent stellar clumps are visible on sub-kpc scales along with a smooth asymmetric component that extends down to ∼100 pc scales, where a dense, eccentric stellar disk dominates the gravitational potential.  The non-axisymmetric component reaches order unity, with Σ_⋆ / Σ_⋆, sym∼ 0.5–2 in the strong non-linear regime.   
  On <10 pc scales (comparable to the orbital motion of the black hole), the stellar distribution becomes less coherent and non-axisymmetries of order Σ_⋆ / Σ_⋆, sym∼ 0.8–1.2 appear and vary on short timescales, while a pc-scale stellar disk forms around the black hole with frequent changes of orientation. 
  Remarkably, in the  conditions, the full range of stellar non-axisymmetries from ∼kpc to ∼pc scales contributes significantly to the total torque on sub-pc scale gas immediately responsible for black hole fueling.   
  
  
  
  
  
  
  
  
  
  
   §.§ Black hole mass dependence
  
  
  
  
  
  The inflow rates presented throughout the paper correspond to the fiducial simulations for a black hole with initial mass M_ BH = 10^8  accreting under different conditions.  Figure <ref> explores the accretion rate dependence on M_ BH for the  (top),  (middle), and  (bottom) conditions.  In each case, we replicate the fiducial simulation with identical resolution and initial conditions but vary the initial black hole mass in the range M_ BH = 10^6–10^9 .  
  Changing M_ BH has a minor impact in the gas and stellar distributions in the initial conditions prior to hyper-refinement. 
  All simulations are evolved for the same physical time except for M_ BH = 10^9  in the  and  conditions, which only follow ∼1 Myr of evolution owing to their significantly higher computational expense as a result of shorter dynamical times.
  
  The inflow/outflow rates at >100 pc are roughly similar in all simulations regardless of M_ BH, as expected given that the dynamical influence of the black hole is negligible on galactic scales.  Nonetheless, we see variations in inflow rate due to stochastic effects, including the movement of the black hole (which we use as reference frame) relative to the center of mass of the galaxy.  The  conditions show very similar inflow rates at <10 pc for M_ BH = 10^6-7  while black holes with M_ BH = 10^8-9  show significantly higher inflow rates.
  The highest average inflow rates at ∼1–100 pc occur for M_ BH = 10^9 , but sub-pc accretion rates are higher for M_ BH = 10^8 .
  This suggests an overall trend for higher Ṁ_ BH with increasing M_ BH but subject to significant stochastic effects.  
  The  and  phases are consistent with no dependence of gas inflow rate on black hole mass at all scales except for M_ BH = 10^6 , with order of magnitude lower Ṁ_ BH.  Given the mass resolution for the (preexisting) stellar and dark matter components (m_⋆ = 3.3 × 10^4  and m_ DM = 1.7 × 10^5 ), it is possible that we are not fully capturing the dynamics and gravitational capture of gas by the M_ BH = 10^6  black hole.
  Overall, while changing M_ BH may affect other properties such as gas density and relative velocity,
  these results are roughly consistent with the weak M_ BH dependence in the gravitational torque accretion model of <cit.>, and in stark contrast with the strong ∼M_ BH^2 dependence in Bondi accretion.
  
  
  
  
  
  
  
  
  § DISCUSSION
  
  
  We have presented a new technique to model black hole growth at sub-pc resolution in a full cosmological galaxy formation context.  We achieve an unprecedented dynamic range by augmenting the Lagrangian MFM hydrodynamics method in GIZMO <cit.> with a new dynamic, hyper-Lagrangian refinement technique that increases the resolution dramatically approaching the central black hole.  
  The success of this technique relies on key physics and algorithmic implementations in FIRE-2 <cit.> and our explicit sub-pc scale treatment of black holes:
  
  
  
  
    * Adaptive gravitational softenings:  We solve gravitational forces for identical gas mass distribution as hydrodynamic forces, decreasing the gravitational softenings dynamically as required by the increased resolution in the hyper-refinement region.
  
  
  
    * Adaptive star formation criteria:  Only gas that is locally self-gravitating and Jeans unstable is allowed to form stars, which is crucial in galactic nuclei where more standard density/molecular thresholds can overestimate the SFR owing to high mean densities even for Jeans stable/unbound gas.
  
  
  
    * Predictive star formation efficiency (SFE):  Gas satisfying the star formation criteria is not artificially forced to have a low SFE per free-fall time, which is rather an outcome of local self-regulation by stellar feedback <cit.>.  This is crucial to predict whether gas can inflow faster than it can turn into stars without imposing artificial constraints in the sub-grid star formation model. 
  
  
  
    * Resolved multi-phase interstellar medium:  We do not include a pressure floor nor impose an effective equation of state for ISM gas, which would prevent us from modeling stellar feedback and turbulence directly or predicting the relative efficiencies of local fragmentation and radial inflow.
  
  
  
    * Adaptive stellar feedback:  The adopted prescriptions are explicitly designed and tested to be consistent across mass resolution in the full range m_ b∼ 1–10^5  without re-calibration or tuning of free parameters <cit.>. This is crucial to maintain physical consistency and model stellar feedback robustly 
  as we transition from standard resolution at >1 kpc to the maximum hyper-refinement resolution level at <100 pc, where individual SNe are well resolved.
  
  
  
    * Time-resolved stellar evolution:  All stellar feedback quantities depend on the actual stellar ages rather than e.g. scaling with SFR, which is crucial in our hyper-refinement simulations with nuclear dynamical times vastly shorter than stellar evolution times.
  
  
  
    * Cooling physics for large dynamic range:  Our cooling module includes dust-gas thermal coupling, cosmic ray-gas coupling, and a treatment for optically-thick cooling which is crucial to capture the thermodynamic state of gas at the densities reached in the nuclear disk (n_ H > 10^7 cm^-3).  The same cooling physics has been used for individual molecular cloud-scale simulations <cit.> reaching 10^-5 pc resolution.
  
  
  
    * Resolved gas gravitational capture:  We avoid many of the uncertainties suffered by sub-grid accretion parameterizations by resolving the gravitational capture of gas by the central black hole at <0.1 pc.
  
  
  
    * Resolved black hole dynamics:  We follow the response of the black hole to the time varying gravitational potential without including artificial drag forces or repositioning algorithms, which is crucial to capture the complex dynamics in galactic nuclei that lead to the development of non-axisymmetries, missalignements, and a variety of other phenomena.      
  
  
  
  
  We stress that accurately following such a large dynamic range is not simply a matter of “brute force” resolution increase: many reasonable numerical implementations of the physics applicable on galaxy scales ≳ kpc and resolution ∼10^5 M_⊙ would break/diverge or give un-physical results if naively applied to simulations with four orders-of-magnitude smaller-scale resolution.
  At our highest resolution, we reach timesteps as low as months (Table <ref>), which makes this technique computationally infeasible for cosmological integrations down to z=0.  Crucially, we circumvent this limitation by exploiting existing FIRE-2 simulations to identify interesting phases for re-simulation at ultra-high resolution.
  As a first application, we have focused on gas transport and black hole growth in a massive galaxy (M_⋆∼ 6–20 × 10^10 ) before, during, and after its peak of nuclear gas density at z∼ 2.5→1.  
  
  
  
  
  
  
   §.§ Comparison to previous methodologies
  
  
  
  <cit.> carried out adaptive mesh refinement (AMR) simulations of a Milky Way-mass progenitor at z=3–6, performing one of the few studies of gas transport down to sub-pc scales in a cosmological galaxy formation context.  They included prescriptions for star formation and stellar feedback but neglected their effects when reaching the maximum resolution in the nuclear region, missing also a treatment for optically-thick cooling.
  While their extremely gas rich conditions are rather different, the diversity of transient features, high time variability, and strong misalignment across scales in our simulations are qualitatively consistent with their findings. 
  
  
  Smoothed particle hydrodynamics (SPH) simulations with static particle splitting have also been used to study the development of multi-scale non-axisymmetries in idealized systems <cit.>.
  These valuable studies are nonetheless limited by the use of a pressurized ISM via an effective equation of state, lack of low-temperature cooling, constant density threshold for star formation, low sub-grid star formation efficiency, or stellar feedback prescriptions designed for lower resolution.
  Some of our key results are in good agreement with earlier idealized models, including the dominant role of stellar torques and the weak dependence of inflow rate on M_ BH <cit.>, confirmed in more recent simulations with stellar feedback and multiphase cooling <cit.>.
  
  
  Complementary approaches are becoming possible with the implementation of super-Lagrangian refinement techniques in other cosmological simulations codes. 
  <cit.> implemented a dynamic refinement scheme in the moving-mesh code AREPO <cit.> to increase the resolution dynamically around the black hole, which has been used to improve black hole feedback implementations in cosmological simulations <cit.> and idealized simulations of galaxies <cit.> and clusters <cit.>.
  <cit.> presented a new refinement scheme to resolve the black hole radius of influence in RAMSES <cit.> and studied the growth of black hole seeds in isolated cooling halos.  
  These novel techniques could also provide sub-pc resolution in a full cosmological context, but exploiting the increased dynamic range while maintaining physical consistency requires including ISM/stellar physics such as those included here, capable of extrapolating to very different resolution and density regimes.
    
  
  
  
  
   §.§ Reproducing luminous QSO inflows
  
  
  The observed radiative output of QSOs with bolometric luminosities L_ bol∼ 10^46–10^47 erg s^-1<cit.> imply black hole accretion rates Ṁ_ BH∼ 2–20  assuming a radiative efficiency ϵ_ rad≡ L_ bol / Ṁ_ BH c^2 ≈ 0.1 <cit.>, and up to an order of magnitude higher for the most extreme QSOs known <cit.>.  Reproducing such high inflow rates from galactic scales down to the black hole accretion disk represents a significant challenge for models.  
  Cold flows, major mergers, minor mergers, and secular processes are generically identified as AGN triggering mechanisms in cosmological and galaxy scale simulations <cit.>, but predicted inflow rates have been often limited by either resolution, sub-grid ISM physics, idealized initial conditions, or uncertainties in black hole accretion parameterization.
  For the first time in a full cosmological context, we have shown that gas inflow rates of up to ∼10  can indeed occur at sub-pc scales in a host halo with M_ vir∼ 10^12.5  at z∼2, sufficient to power a luminous QSO.
  
  
  
  Our simulated  phase corresponds to very particular conditions selected at the peak of nuclear gas density.
  The four FIRE-2 massive halos simulated in <cit.> all show a similar phase of extreme nuclear densities, reaching Σ_ gas > 10^10.5  and Σ_⋆ > 10^12  on ∼10 pc scales, suggesting that all massive halos undergo at least one phase conducive to a luminous QSO <cit.>.
  Constraints on QSO lifetimes from black hole demographics <cit.>, clustering <cit.>, radiation proximity effects <cit.>, gas fueling arguments <cit.>, and others vary from ∼0.1–100 Myr <cit.>, though there is indication that AGN flicker on shorter timescales <cit.>.
  In the  phase, we measure an inflow rate of ∼6  at <0.1 pc almost in steady state during the full ∼4 Myr duration of the simulation, with no sign of decline.  The full history of the host galaxy (Figure <ref>) shows nuclear Σ_ gas > 10^10.5  for ∼50 Myr, suggesting that QSO-like inflow rates could be maintained for up to a few times longer than the ∼4 Myr period simulated here.   
  These results are consistent with inferred QSO lifetimes.  We note, however, that we have neglected the effects of AGN feedback throughout this study.
  
  
  
  
  
  
   §.§ Diversity of accretion phases on cosmological timescales
  
  
  
  Hyper-refinement simulations on different regimes (Figure <ref>) suggest that the  conditions may indeed correspond to a special period in the host galaxy's lifetime.
  We find strikingly different conditions a mere ∼40 Myr preceding the  phase.  While the gas and stellar surface density profiles at >100 pc are roughly similar, the  conditions show extreme variability in sub-pc scale inflow rate.  During the ∼3 Myr duration of the simulation, the accretion rate varies by four orders of magnitude from Ṁ_ BH = 0.001–10 , corresponding to L_ bol∼ 6×10^42-46 erg s^-1.  With a factor of ten lower average accretion rate, ⟨Ṁ_ BH⟩ = 0.6 , order of magnitude changes in  accretion can occur in only ∼0.01 Myr, in stark contrast with the steady inflow rate in the  phase.  
  This transition likely originates from the significantly deeper stellar potential in the  phase, providing stronger torques, enhanced retention of gas against turbulent motions and stellar feedback, and stabilizing the black hole location in the inner ∼10 pc, compared to the significantly more dynamic  conditions with a flat stellar surface density profile in the central ∼100 pc.
  
  
  Another important qualitative change in accretion properties occurs at late times (z ∼ 1), when a thin, rotationally supported gas disk forms, with lower surface density and star formation activity  <cit.>.  These  conditions show significantly lower average accretion rate, ⟨Ṁ_ BH⟩ = 0.01 , in this case characterized by ∼2 Myr-long active phases with Ṁ_ BH = 0.005–0.1  followed by longer periods of inactivity owing to the formation of a hot, low-density cavity in the central ∼100 pc that inhibits black hole growth.  Gas consumption by star formation and stellar feedback are likely responsible for the formation of the cavity <cit.>, and possibly an inner Lindblad-like resonance <cit.>, only penetrated by cold gas clumps and filaments falling toward the black hole when destabilized out of the inner edge of the circumnuclear disk (Figure <ref>).  These conditions are likely representative of more frequent, low-luminosity AGN in massive galaxies after their peak in star formation activity, with their central black holes accreting predominantly at low-Eddington ratios in the radiatively inefficient, kinetic regime <cit.>.
  
  
  Interestingly, the  phase occurs at the time when the host galaxy is transitioning from bursty star formation driving prominent galactic winds at early times to a more steady mode of star formation that promotes the formation of extended disks at late times.  This transition has been identified in the FIRE simulations for massive (M_⋆ > 10^10 ) galaxies at different redshifts  <cit.> and coincides with the virialization of the inner halo <cit.>.  The early bursty epoch represents yet another qualitatively distinct regime for black hole growth, where sub-grid accretion calculations predict short (∼1 Myr) accretion episodes that can reach or even exceed the Eddington rate followed by longer periods of inactivity owing to stellar feedback continuously evacuating gas from galactic nuclei (Figure <ref>, middle panel; ; ; Byrne et al. in prep.).  These early conditions should be investigated in future hyper-refinement simulations.
  
  
  
  
  
   §.§ Multi-scale torques and AGN triggering
  
  
  
  Our simulations show a diverse range of morphologies across scales: from large scale filaments to nested spirals, bars, ring-like features, irregular streams, and clumps (Figures <ref>–<ref>). 
  The cascade of non-axisymmetric gravitational instabilities seen on different scales resembles the classic “bars-within-bars” scenario <cit.>, generalized to a diverse range of morphologies as found in previous idealized models <cit.> and some cosmological simulations <cit.>.
  Non-axisymmetric structures with a variety of morphologies are indeed often observed in galactic nuclei <cit.> and connected to kinematic evidence for radial gas flows.
  
  
  Despite the diversity in physical conditions and accretion properties, a common feature identified for the , , and  phases is that gravitational torques from the stars dominate gas angular momentum transport down to our ∼0.1 pc resolution limit (Figure <ref>), in agreement with <cit.>.  
  While gas shows strong non-axisymmetries on all scales, gravitational self-torquing can become comparable to the stellar torques only on galactic (∼ kpc) scales for our most gas-rich  conditions.  Hydrodynamic torques by pressure gradients are stronger than gas self-torquing but still sub-dominant compared to stellar torques.
  Other mechanisms for angular momentum transport may be more important in very gas-rich systems with efficient local fragmentation, including scattering of dense gas clumps and turbulence driven by gravitational instability and stellar feedback <cit.>.  However, comparing our  and  conditions suggests that gravitational torques from a dominant asymmetric stellar nuclear structure is the most efficient mechanism to drive sustained large inflow rates.
  
  
  We have explicitly shown that it is generally local torques from nearby stars that are most relevant for driving gas inflows (Figure <ref>), as assumed in models of local angular momentum transport <cit.>.  Our best example is the  phase, where the total torque on gas at a given scale is clearly dominated by stars living in the same radial domain, implying that the instantaneous gas inflow rate at ∼ 0.1pc is decoupled from the larger scale stellar distribution.  However, the unique dynamic range of our simulations demonstrates that this is not always the case: the net torque on sub-pc scale gas in the  phase contains significant contributions from stars in the entire radial range 0.1 pc–1 kpc, which cannot be explicitly captured in idealized nuclear scale simulations or local transport models.
  In either case, the characteristic time variability of stellar non-axisymmetries varies widely on different scales (Figure <ref>) and the net torque on gas can change sign both radially, in a time average sense, and temporally, in a given spatial scale <cit.>, with important implications for observations connecting large scale torques and AGN activity.
  The small scale (<10 pc) motion of the black hole relative to the stellar center of mass in the  phase plays a minor role in the observed non-axisymmetries at >100 pc, but more off-centered black holes in low mass hosts and/or higher redshift <cit.> could significantly affect the resulting large scale gravitational torques.
  
  
  The non-trivial dependence of sub-pc inflow rates on multi-scale non-axisymmetries and their time variability can explain why observations provide inconclusive evidence for a connection between large scale bars <cit.>, clumpy disks <cit.>, or nuclear non-axisymmetries <cit.> and instantaneous AGN activity <cit.>.
  For instance, our simulated  phase shows time-averaged positive torque and correspondingly net outflow rate on ∼100 pc scales (Figures <ref> and <ref>) while sustaining a large inflow rate of ∼10  on ∼10 pc scales (of which  ∼6  reach <0.1 pc) owing to strong time-averaged negative torquing in the inner region.  
  
  
  Observations also appear inconclusive regarding the relative roles of galaxy mergers <cit.> and secular processes <cit.> driving black hole growth, 
  with idealized and cosmological simulations sometimes producing contrasting results <cit.>.
  A massive satellite galaxy (M_⋆∼ 10^10 ) is approaching the central galaxy in its second passage toward finale coalescence at the time of the  and  phases, located ∼20 kpc and ∼15 kpc away, respectively. 
  These accretion phases occur within only ∼40 Myr of each other and produce accretion rates that vary by four orders of magnitude while the stellar distribution on ∼kpc scales does not change appreciably, highlighting the difficulty of connecting global morphological features with instantaneous accretion rate.  
  The merging galaxy contributes less than ∼1 % of the total stellar torque on nuclear scale gas (<100 pc) during the  and  phases (Figure <ref>), suggesting a minor role feeding the central black hole, but the overall impact of the ongoing merger on longer timescales (e.g. increasing the density of the central ∼1 kpc) should be investigated in future work.  The lower nuclear activity in the  conditions is clearly unrelated to mergers and driven by secular processes.
  
  
  
  
  
  
   §.§ The SFR–AGN connection
  
  
  In addition to the self-gravity of stars and gas, our simulations show that it is crucial to include star formation and stellar evolution in gas transport calculations <cit.>.  Stellar feedback can produce galaxy-scale perturbations in the form of galactic winds and coherent re-accretion events <cit.>, 
   but also regulates star formation locally and generally impacts the dynamics and multi-phase structure of the ISM <cit.>.  This has important implications for AGN variability and duty cycle as seen for the  and  conditions (Figure <ref>), where the accretion rate becomes significantly suppressed when the black hole is embedded in the hot gas phase.
  Stellar mass return can become an important source of fueling at late times in gas poor galaxies <cit.>, 
  but star formation itself can decrease substantially the gas supply available for black hole growth: in our simulations, star formation can consume as much gas as is provided by inflows in the full range of scales (1 pc–10 kpc), with some indication for comparatively lower star formation on pc scales (Figure <ref>; see also ).  Recently formed stars (<5 Myr) can contribute ∼10–20% of the total stellar torque on gas at <10 pc for gas-rich conditions ( and  phases), further emphasizing the importance of including star formation down to the smallest scales. 
  
  
  Our results also have important implications for the connection between global SFR and instantaneous black hole accretion.  Positive average correlations between global SFR and AGN luminosity have been reported for a number of observations <cit.>, while others suggest that they hold only for bulge-dominated galaxies <cit.> or the highest luminosity systems <cit.>. 
  Other studies suggest little or no connection between average SFR and black hole accretion <cit.> or link luminous AGN activity with a suppression of star formation <cit.>, possibly only for galaxies above the star formation main sequence <cit.>.
  Part of the discrepancy may originate from biases and selection effects <cit.>, but different characteristic timescales for star formation and AGN activity are likely crucial as well <cit.>.  
  
  
  Our simulations demonstrate that the instantaneous accretion rate (measured at <0.1 pc) can indeed vary by orders of magnitude for a given host galaxy SFR <cit.>. 
  We find that the average accretion rate roughly follows a linear proportionality ⟨Ṁ_ BH⟩∝⟨ SFR⟩ for star formation within small radii, but time variability decouples AGN activity and global star formation.  
  These results are consistent with observations in the local universe indicating correlations between nuclear star formation and instantaneous black hole accretion in individual systems but no clear connection to global star formation <cit.>.
  
  
  We find time-averaged accretion to global star formation ratios Ṁ_ BH /  SFR∼ 1/1000 for the  and  conditions, in agreement with expectations from the M_ BH/M_ bulge ratios measured in the local universe <cit.> and similar to the typical Ṁ_ BH /  SFR ratios observed for AGN in star forming galaxies <cit.>.  
  In our  phase, the central black hole is overgrowing its host galaxy, with Ṁ_ BH /  SFR∼ 1/50, also in good agreement with observations of luminous QSOs and AGN-dominated sources <cit.>.
  This suggests that a common gas supply for star formation and black hole growth regulated by gravitational torques may play a key role driving the black hole–galaxy connection <cit.>.
  
  
  
  
  
   §.§ From galaxy scales to the accretion disk
  
  
  
  The inflow rate at 0.1 pc corresponds roughly to the scale below which the non self-gravitating accretion disk is expected to form according to analytic models <cit.> and reverberation mapping observations <cit.>: 10^3–10^4 R_ s≈ 0.1 pc for the Schwarzschild radius R_ s of black holes with M_ BH = 10^8-9 .  
  Our results have interesting implications for accretion disk models, providing realistic outer boundary conditions that differ significantly from the idealized, torus-like, equilibrium conditions often employed as initial conditions <cit.>.  We have shown that radial flows are present at ∼0.1 pc and the outer disk feeding rate can vary by orders of magnitude on ∼3 Myr (Figure <ref>), which would imply drastic changes in the geometric and radiative properties of the accretion disk, from super-Eddington fueling conditions <cit.> to low-Eddington, hot accretion flows <cit.>, on timescales comparable to the standard α-disk viscous time <cit.>.  
  
  
  
  
  Gas on sub-pc scales is preferentially rotationally supported but the specific angular momentum direction can change on timescales <0.1 Myr.  Sub-pc scale disks can thus be significantly misaligned with the galaxy-scale disk, in agreement with observations of megamaser disks <cit.> and jets <cit.>. 
  Nuclear scale disks are also often decoupled from the angular momentum of the host galaxy <cit.>, in agreement with our results <cit.>.
  Strikingly, we find persistent, increasingly misaligned orientation across scales in the  conditions, with near spin flip in the pc-scale disk orientation relative to kpc scales (Figure <ref>).
  The changing orientation of sub-pc scale disks can have drastic implications for the coupling efficiency of black hole feedback <cit.> and the spin evolution of massive black holes <cit.>, which itself has strong implications for black hole mergers <cit.> and the origin of radio jets <cit.>.
  
  
  
  
  
   §.§ Black hole accretion parameterizations
  
  
  Most modern cosmological simulations employ sub-grid black hole accretion prescriptions based on the spherical <cit.> model[A notable exception is the hybrid accretion prescription in the SIMBA simulation <cit.>, where cold gas accretes following the gravitational torque accretion model <cit.> and hot gas accretes following the Bondi parameterization.] (e.g. Horizon-AGN, Eagle, and IllustrisTNG).
  Despite the success of these models at reproducing global galaxy and black hole observables, previous idealized simulations have shown that Bondi accretion can fail to reproduce gas inflow rates by orders of magnitude under a variety of conditions <cit.>, limiting the predictive power of cosmological simulations.
  Our results indicate that the inflow at sub-pc scales is dominated by cool gas with significant rotational support, and it is primarily driven by stellar gravitational torques.  
  This is inconsistent with the assumption in Bondi accretion of spherically symmetric, pressure-supported gas in hydrostatic balance.
  
  If we apply Bondi accretion to the  phase (arguably the most favorable conditions for Bondi, with higher presence of hot, pressure-supported gas compared to the other phases), we obtain Ṁ_ Bondi = 4π  G^2  M_ BH^2 ρ  /  c_ s^3  > 10^7  for M_ BH = 10^8  and the gas density ρ and sound speed c_ s found during the active phases on pc scales, in striking contrast with the range Ṁ_ BH∼ 0.01–0.1  that we measure in the hyper-refinement simulation (Figure <ref>).  If we consider “turbulent Bondi-Hoyle”, replacing c_ s^3→ (c_ s^2+σ_ z^2)^3/2 where σ_ z is the vertical velocity dispersion, we obtain Ṁ_ Bondi,turb > 10^4 , still many orders of magnitude higher than the explicit inflow rate <cit.>.  If we further include the gas azimuthal velocity v_ϕ, replacing[Including the azimuthal velocity in this way is not well motivated theoretically.] c_ s^3→ (c_ s^2+σ_ z^2+v_ϕ^2)^3/2 similar to the prescription in Eagle <cit.>, the accretion rate is still overestimated by more than two orders of magnitude.  Considering the black hole–gas relative velocity has a minimal impact compared to σ_ z and v_ϕ, while adding a density-dependent “boost" factor <cit.> only increases these discrepancies.
  In addition, Bondi-based models predict a decrease in gas inflow rate with increasing radial distance by factors >10^3 in the range 0.1–100 pc, in clear contrast with the opposite radial trend shown in Figure <ref>, and do not capture the rather different accretion conditions in the , , and  phases.  
  In practice, Bondi-based models in cosmological simulations do not yield such unrealistically high accretion rates owing to self-regulation by AGN feedback and imposing the Eddington limit, but this illustrates the inapplicability of current black hole fueling models.
  
  In contrast, a simple free-fall accretion estimator evaluated at ∼100 pc (with overall normalization such that M_ BH/M_ bulge∼ 0.002 at z=1) reproduces the sub-pc inflow rates to within an order of magnitude (Figure <ref>).  While the free-fall estimator over/under-predicts the low/high accretion phases[A trend of overpredicted inflow rate during low-accretion phases () and underpredicted inflow rate during high-accretion phases ( and ) is expected: the constant normalization factor in the free-fall estimator reproduces the average behavior but fails to capture the full variability.  In the more complex gravitational torque model of <cit.>, the normalization depends on the amplitude of non-axisymmetric perturbations, which is higher in the  and  phases and lower in the  phase, in the right direction toward reproducing the sub-pc inflow rates.], it reproduces the correct relative differences between the , , and  phases, unlike Bondi accretion.
  We defer a detailed comparison to accretion parameterizations to future work, but note that a key element of sub-grid accretion prescriptions is the intrinsic dependence on black hole mass, Ṁ_ BH∝ M_ BH^p <cit.>: strong M_ BH dependence (p>1) requires strong feedback self-regulation, as is the case in Bondi accretion, while it is not a necessary condition for weak M_ BH dependence (p<1).
  The sub-pc accretion rates predicted in our simulations indicate a weak M_ BH dependence for a variety of conditions, 
  in good agreement with the idealized simulations of <cit.> and their gravitational torque accretion model.  
  In addition to being crucial for the origin of black hole–galaxy scaling relations <cit.>, the weak M_ BH dependence also has important implications for super-Eddington feeding.  
  In the absence of AGN feedback, we show that sustained sub-pc fueling rates can reach up to ×3 and ×10 higher than the Eddington limit for M_ BH = 10^8  and M_ BH = 10^7 , respectively, in the extreme conditions of the  phase.
  The Eddington limit is often imposed in cosmological simulations but can be greatly exceeded according to recent accretion disk models <cit.>, provided that sufficiently high inflow rates can be driven from larger scales as our simulations predict.  
  
  
  
  
  
  
   §.§ Caveats and future prospects
  
  
  We have neglected AGN feedback in order to study how gas is transported down to the black hole accretion disk in the first place.  Reported accretion rates should thus be regarded as upper limits for how efficiently stellar torques can drive gas inflows across scales.  A variety of idealized, high-resolution simulations have shown that AGN feedback can open a central cavity and thus regulate the growth of the central black hole <cit.>. 
  However, the coupling efficiency of AGN feedback in the crucial 0.1–10 pc regime and thus the extent to which self-regulation occurs remains uncertain under highly non-isotropic gas distributions and extreme inflow conditions such as those seen in the  phase.  Future cosmological hyper-refinement simulations will address the role of self-regulation by incorporating a novel implementation of fast, accretion-driven winds <cit.> and other black hole feedback channels.
  
   
  On the other hand, we have neglected accretion of stars by the central black hole, which could represent an additional source of growth <cit.>.  Future work should consider the incidence of tidal disruption events under the extreme nuclear stellar densities reached in the  and  condition, though these are probably unphysically large owing to the lack of AGN feedback in the original FIRE-2 simulation employed for re-simulation <cit.>.  Synthetic observations including dust radiative transfer show that the effective size and stellar mass surface density are in general agreement with observational scaling relations at z>2 but that the stars are too compact at lower redshift <cit.>, comparable only to the most extreme compact galaxies observed <cit.>.  Encouragingly, our results show that increasingly large inflow rates fuel the central black hole at the time that AGN feedback is becoming most needed.
  
  
  Despite the unprecedented dynamic range with full galaxy formation physics, our simulations are still ∼3 orders of magnitude from actually resolving accretion down to the innermost stable circular orbit.
  Nonetheless, our hyper-refinement technique opens up the possibility of implementing black hole feedback exactly at the scale where accretion disk simulations predict the mass, momentum, and energy output (and their angular dependence) from winds and radiation <cit.>, as well as studying their coupling efficiency and impact from sub-pc scales to the CGM.
  Additional physics missing in our calculations but readily available in GIZMO include magnetohydrodynamics, conduction, viscosity, and cosmic ray feedback <cit.>, whose effects have not been explored yet in a galaxy formation context at the resolution achieved here.
  
  
  Our simulations have been limited to exploring only three different galactic nuclei conditions corresponding to a massive galaxy (M_⋆∼ 6–20 × 10^10 ) before, during, and after its peak of nuclear gas density at z∼ 2.5→1.  The , , and  phases simulated predict strikingly different black hole accretion properties, which motivates expanding this work to a much broader set of host galaxy conditions.  Future cosmological hyper-refinement simulations should model black hole growth and feedback in the full range of environments, from active dwarf galaxies <cit.> and nearby Seyferts with obscuring “torus” structures <cit.> to the extreme conditions triggering the highest redshift QSOs <cit.>.  
  
  
  
  
  
  We thank Ena Choi, Romeel Davé, Robert Feldmann, John Forbes, Shy Genel, Melanie Habouzit, Yan-Fei Jiang, Yuan Li, Xiangcheng Ma, Eve Ostriker, Roxana Pop, Julissa Rojas-Sandoval, Matthew Smith, Tjitske Starkenburg, Kung-Yi Su, Paul Torrey, Rainer Weinberger, and Sarah Wellons for many insightful discussions and suggestions during the development of this work.  
  We thank the referee for a detailed and constructive review that helped improve the paper.
  We acknowledge outstanding support by the Scientific Computing Core group and the Center for Computational Astrophysics at the Flatiron Institute as part of the SMAUG project, which are supported by the Simons Foundation.
  DAA was supported in part by NSF grant AST-2009687.
  EQ was supported in part by a Simons Investigator Award from the Simons Foundation and by NSF grant AST-1715070.
  Support for PFH was provided by NSF Collaborative Research Grants 1715847 & 1911233, NSF CAREER grant 1455342, NASA grants 80NSSC18K0562, JPL 1589742. 
  CAFG was supported by NSF through grants AST-1517491, AST-1715216, and CAREER award AST-1652522; by NASA through grant 17-ATP17-0067; and by a Cottrell Scholar Award and a Scialog Award from the Research Corporation for Science Advancement.
  GLB was supported in part by the NSF throughs grants OAC-1835509 and AST-2006176, as well as support from the STScI under NASA contract NAS5-26555.
  DK was supported by NSF grant AST-1715101 and the Cottrell Scholar Award from the Research Corporation for Science Advancement.
  The simulations were run on Flatiron Institute's research computing facilities (Gordon-Simons, Popeye, and Iron compute clusters), supported by the Simons Foundation, and XSEDE allocation TG-AST160048, supported by NSF grant ACI-1053575.  Additional numerical calculations were run on the Caltech compute cluster “Wheeler,” allocations FTA-Hopkins supported by the NSF and TACC, and NASA HEC SMD-16-7592. 
  
  
  
  
  
  
  
  
  
  
  § 
  
  
  
  
   §.§ Resolution convergence
  
   
  Figure <ref> investigates the convergence properties of predicted accretion rates relative to mass resolution for the  and  conditions.  We show cumulative black hole growth as a function of time in simulations starting from identical initial conditions but implementing dynamic hyper-refinement of gas to different levels, where the color scale indicates the factor by which the mass resolution is increased relative to the initial conditions (see Figure <ref> for the radial profile of gas mass resolution in each level).  All other simulation parameters, including gravitational softenings and accretion aperture, are kept identical.
  There is no clear direct relationship between change in resolution and variation in accretion rate, but we find that lower resolution simulations tend to overpredict black hole growth by up to one order of magnitude relative to our highest resolution simulations.  Encouragingly, there is indication of convergence for the highest resolution levels, when the gas mass resolution in the nuclear region satisfies m_ g≲ 100 .
  
  
  
  
  
  
  
  It is interesting to note that the average black hole growth rate, measured as explicit gravitational capture of gas at <0.1 pc, varies by less than one order of magnitude in  simulations with mass resolution changing by ×1000.  The weak dependence of  conditions to mass resolution can be attributed to a number of factors:
  
  
  
    * The initial black hole (M_ BH = 10^8 ) is already very massive compared to individual gas particles even at  the lowest resolution (m_ g∼10^4 ), implying that the black hole radius of influence is well resolved and gas particles can indeed be gravitationally bound to the black hole in a resolved sense.
  
  
  
    * The total, time-integrated black hole growth (∼10^6  at the highest resolution) is also significantly larger than a single gas particle mass even at the lowest resolution, still allowing for a particle-by-particle representation of accretion over the episode.
  
  
  
    * The gas inflow rate is roughly constant within a radial distance R ≲ 50 pc (Figure <ref>) which is large enough to be resolved even at the lowest resolution, and the amount of gas flowing into this central region is roughly independent of resolution over the tested range.
  
  
  
  
  
  
  
  The accretion aperture R_ acc=0.1 pc is not resolved at low resolution, where individual particles represent a much larger volume, but particles can still be accreted when their center of mass happens to be within R_ acc.  The probability of finding a low resolution particle within R_ acc is low, but a single accretion event contributes the equivalent of many high resolution particles.  Since the radial inflow is roughly independent of R below a scale which is also resolved at our lowest resolution, the time-integrated black hole growth becomes roughly independent of resolution.
  These conditions can clearly break down for a number of reasons, including lower mass black holes, lower accretion rates, or strong radial dependence of inflow rates.  
  For example, the  conditions show more variation with resolution owing to the lower accretion rate and total integrated black hole growth.  In this case, even the highest resolution levels are not expected to fully converge during the Δ t ∼ 5 Myr period for which simulations with different resolution levels are available, which is not long enough to average over the stochastic fueling conditions and low duty cycle driven by the formation of a hot, low density cavity in the central ∼100 pc (see Figure <ref>).   
  
  Figure <ref> shows an additional resolution convergence test for the  phase, in this case comparing the full radial profile of gas inflow rate, outflow rate, and enclosed SFR (similar to Figure <ref> but separating the inflow/outflow components) for the two simulations with the highest resolution available (corresponding to m_ g∼ 30–60 ) during the same evolution time Δ t ∼ 1 Myr.  Note that 1 Myr corresponds to ∼10–100 dynamical times at 10–1 pc for the  phase, reaching a reasonable statistical steady state for the much larger and more steady accretion rates.  We find roughly similar time-averaged radial profiles for both resolution levels over the full dynamic range 0.1 pc–10 kpc, with a slight trend for higher rates at <10 pc in the lower resolution simulation.  While these measurements can be affected by stochasticity in black hole dynamics (with the black hole location used as center of reference), the sub-pc inflow rates at our highest resolution levels converge to within ∼50%.
  Overall, these tests show that our results are robust to changes in resolution and motivate the development of hybrid accretion models for full cosmological simulations, switching between (1) a continuous sub-grid accretion model when the inflow rate and black hole radius of influence are unresolved to (2) direct gravitational capture of individual gas particles when they are resolved.
  
   
  
  
  
   §.§ Gravitational softenings and stellar splitting
  
  
  We test the robustness of our results relative to the treatment of gravitational softenings and the mass resolution of the stellar component by comparing our fiducial model described in <ref> with simulations implementing:
  
  
    * Larger stellar gravitational softenings, which we increase from ϵ_⋆ = 0.1 pc in the fiducial model to the same stellar softening employed in the original FIRE-2 simulation, ϵ_⋆ =7 pc.
  
  
  
  
  
  
  
  
  
    * Adaptive gravitational softenings for collisionless particles, which are computed as the mean (kernel-averaged) interparticle spacing within the kernel as for the gas component <cit.>.  For star particles, ϵ_⋆∝ρ_ b^-1/3, where ρ_ b is the total baryonic (gas+stars) mass density while for the dark matter the adaptive softening is determined from the dark matter particle neighbor distribution, ϵ_ DM∝ρ_ DM^-1/3.  The minimum adaptive softenings for stars and dark matter are ϵ_⋆,  min = 0.1 pc and ϵ_ DM, min = 57 pc, respectively.
  
  
  
    * Variable gravitational softenings for star particles, with ϵ_⋆ = 7 pc × (m_⋆/ 3e4  M_⊙)^1/3 (and upper limit ϵ_⋆,  max = 7 pc) such that star particles represent the same stellar mass density regardless of the particle mass m_⋆.  This assigns the same stellar softening as the original cosmological zoom-in simulation for pre-existing stars at the initial conditions (ϵ_⋆ =7 pc) while providing higher force resolution for star particles forming out of higher resolution gas in the nuclear region (ϵ_⋆,  min∼ 0.6 pc for m_⋆ = 20 ).
  
  
  
    * Dynamic stellar splitting, where we increase the mass resolution of the stellar component by splitting star particles progressively according to their distance to the central black hole.  We define a mass resolution profile as in Figure <ref> with maximum refinement factor χ_ ref = 1/64, which yields m_⋆,  min≈ 500  in the central ∼10 pc.  The splitting procedure for star particles is identical to that of the gas component except that split stars retain the same velocity as the parent star particle.
  Note that this is a simplified treatment of stellar splitting, which should involve interpolation in the 6D phase space that describes collisionless systems.  In this case, we use fixed softenings with ϵ_⋆ = 0.1 pc as in the fiducial model.
  
  
  
  Figure <ref> shows the cumulative black hole growth predicted for the  phase in simulations testing different implementation details.  While the instantaneous accretion rate can vary by up to two orders of magnitude between different models, owing to stochastic variability in the very dynamic  conditions, the cumulative black hole growth follows a similar trend regardless of the exact values of gravitational softenings or numerical implementation details.  
  Analysis of gravitational torques (as in Figures <ref>, <ref>, and <ref>) shows no significant differences between simulations except for a decrease in the contribution of stars within 1 pc to the total gravitational torque on gas at <1 pc when using larger stellar gravitational softenings (ϵ_⋆ =7 pc), as expected.  While the integrated black hole growth is not significantly affected by the larger softenings (possibly due to the very strong large scale torques but also stochasticity in the  conditions), these results suggest that it is more important to achieve high gravitational force resolution than specifically hyper-refining the stellar component, with the overall results robust to numerical details.
  
   
  
  
  
  
  
  § 
  
  
  
   §.§ Multi-scale gas in the  phase
  
  
  
  Figure <ref> shows the multi-scale gas mass surface density distribution for the  conditions, for direct comparison to the gas density distributions in the  and  phases (Figures <ref> and <ref>) as well as the temperature distribution for the  conditions (Figure <ref>).  Despite the similar high level of turbulence, the  phase shows significantly more centrally concentrated gas density distribution on scales <1 kpc compared to the  conditions, reaching Σ_ gas > 10^11  in the central 10 pc and feeding the central black hole in steady state through an extremely dense and cold pc-scale gas disk. 
  
  
  
  
  
  
  
  
  
  
  
  
  
  aasjournal
  
  
  
  