Eddington limit, for 100 per cent of the age of the universe to reach their observed masses at z > 7. But observations at all lower redshifts, and theoretical estimates of the effect of SNe and BH feedback and BH dynamics all argue for much lower duty cycles (see e.g. Johnson & Bromm 2007;Whalen et al. 2008;Alvarez, Wise & Abel 2009;Milosavljević, Couch & Bromm 2009;Habouzit, Volonteri & Dubois 2017. An obvious possible solution is to form more massive seeds: It has been proposed that primordial gas at high-z could experience inefficient cooling and fragmentation, producing massive Population III stars (Bromm & Larson 2004) which could collapse to BH seeds as large as ∼ 100 M (e.g. Madau & Rees 2001;Li et al. 2007;Volonteri 2012;Hirano et al. 2014) or even hypermassive quasi-stars which could leave seeds as large as ∼ 10 4 -10 5 M (e.g. Bromm & Loeb 2003;Hosokawa, Omukai & Yorke 2012;Hosokawa et al. 2013;Hirano et al. 2017;Inayoshi, Li & Haiman 2018), or directly collapsing to BHs as massive as 10 5 M (Lodato & Natarajan 2006, 2007). Yet several authors argue that this requires vanishingly improbable conditions (see e.g. Corbett Moran, Grudić & Hopkins 2018 and discussions in sections 5.2 and 5.3 from Inayoshi et al. 2020). But even these most-optimistic models only reduce the time-scales by a logarithmic factor [as time-scales scale as log (M SMBH /M seed )]: even in these models, a phase of highly super-Eddington accretion -either resulting from runaway gas capture in high-gas-density regions (e.g. Alexander & Natarajan 2014;Madau, Haardt & Dotti 2014;Lupi et al. 2016;Pezzulli, Valiante & Schneider 2016;Regan et al. 2019;Natarajan 2021), or runaway mergers of massive stars (e.g. Portegies Zwart et al. 2004;Devecchi & Volonteri 2009;Katz, Sijacki & Haehnelt 2015;Reinoso et al. 2018) or of other seed BHs (e.g. Davies, Miller & Bellovary 2011;Lupi et al. 2014) at the centre of a common potential minimum undergoing dynamical relaxation -is likely needed to explain SMBHs at z > 7 (Haiman 2004;Kawashima et al. 2012;Pacucci, Volonteri & Ferrara 2015;Inayoshi, Haiman & Ostriker 2016;Ryu et al. 2016;Takeo et al. 2019).

However, in the past two decades, many independent studies (e.g. focused on galaxy mergers (Governato, Colpi & Maraschi 1994;Volonteri & Perna 2005;Callegari et al. 2009;Bellovary et al. 2010;Tremmel et al. 2018b), dwarf galaxy evolution (Tamfal et al. 2018;Bellovary et al. 2019;Boldrini, Mohayaee & Silk 2020) and/or BH growth/dynamics (Callegari et al. 2011;Volonteri et al. 2016Volonteri et al. , 2020;;Biernacki, Teyssier & Bleuler 2017;Anglés-Alcázar et al. 2017b;Pfister et al. 2019; Barausse et al. 2020;C ¸atmabacak et al. 2020)) have pointed out that all of these models face a different and potentially even more severe challenge: what we refer to as the 'sinking problem.' In brief: essentially all of the rapid/efficient accretion models require that BHs sink 'efficiently' and remain tightly bound to the galaxy centre or potential minimum, where densities are on average highest. This usually requires a well-defined and stable dense central region in a relatively massive galaxy at lower redshift (z ࣠ 4) (Tremmel et al. 2017(Tremmel et al. , 2019;;Ricarte et al. 2019), but it may not be possible dynamically for even 'high' mass seeds in realistic turbulent, clumpy, high-redshift (z 7) galaxies which undergo frequent dynamical perturbations (from e.g. mergers and 'bursty' star formation and stellar feedback) and lack such central regions, especially in the short time-scales available. Observationally, SMBHs are seen in the Galactic Centre for most massive quasistellar objects (QSOs) (including those at high-z where imaging is possible, e.g. Venemans et al. 2017;Bañados et al. 2019;Decarli et al. 2019;Novak et al. 2019;Wang et al. 2019, and almost all massive galaxies comparable to QSO hosts at low redshifts, see e.g. Ferrarese & Merritt 2000;Gebhardt et al. 2000;Tremaine et al. 2002;Graham et al. 2011;Beifiori et al. 2012). But in spatially resolvable low-z dwarf galaxies where star formation is known to be 'bursty' (Weisz et al. 2014;Sparre et al. 2017;Faucher-Giguère 2018;Velázquez et al. 2021) and there is no well-defined dynamical centre (see e.g. Kallivayalil et al. 2013), AGNs are extremely rare and those identified are randomly scattered in position around the galaxy (Mezcua & Domínguez Sánchez 2020;Reines et al. 2020). As numerical simulations of high-z galaxies have improved in both numerical resolution and incorporating the physics of star formation and stellar feedback in a turbulent, multiphase ISM (Anglés-Alcázar et al. 2014, 2017b;Kretschmer & Teyssier 2020), most models have converged toward the prediction that high-z galaxies are clumpy, bursty, chaotic, and dynamically unrelaxed systems (even more so than most local dwarfs, e.g. Muratov et al. 2015;Oklopčić et al. 2017;Ma et al. 2018a;Tamfal et al. 2018;Kim et al. 2019;Meng & Gnedin 2020), in agreement with deep observations with the Hubble Space Telescope (HST) (Elmegreen et al. 2007;Overzier et al. 2010;Swinbank et al. 2010). Although there is some evidence for rotation in some hosts as noted by, e.g. Decarli et al. (2018), Venemans et al. (2019), they usually exhibit very large dispersion with σ ∼ v, consistent with the simulations analyzed in Ma et al. (2017), which does not challenge the conclusion. But in almost all models for rapid BH growth at near-Eddington or super-Eddington rates at z 7, the most optimistic assumption possible is usually made: namely that the BH remains 'anchored' to the local potential minimum at the centre of some well-ordered galaxy (e.g. Li et al. 2007). To accrete gas, the BH must first capture it from the surroundings, and dimensional estimates for the 'capture rate' drop highly superlinearly and extremely rapidly if the BH or background medium are moving relative to one another and/or if the BH lies outside of the galactic density maximum (Hoyle & Lyttleton 1939). Models like runaway stellar mergers or BH-BH seed mergers for rapid growth fundamentally depend on the idea that both the 'main seed' BH and all other stars/seeds are anchored to and sinking rapidly towards a common dynamical centre (Portegies Zwart & McMillan 2002;Gürkan, Freitag & Rasio 2004;González et al. 2021;Shi, Grudić & Hopkins 2021) Historically, the 'sinking' of BH seeds in high-z galaxies has largely been studied by assuming (1) seeds form at the centres of their protogalaxies (rather than where stars form or at local density maxima), (2) galaxies are smooth objects with well-defined dynamical centres and centrally peaked density profiles (i.e. bulge + disc or isothermal sphere models, rather than messy, non-relaxed systems), and (3) that BH and merging galaxy orbits decay according to dynamical friction (DF), which is a statistical accumulative effect caused by successive two-body gravity encounters, effectively acting like a 'drag force' proportional to the BH/merging galaxy mass, in which the traditional Chandrasekhar (1943, hereafterC43) DF formula (assuming a homogeneous, infinite, idealized background medium) is applied. In this paper, we therefore revisit the 'sinking' and 'retention' problems for seed BHs in early galaxies. We use high-resolution cosmological simulations which include the crucial physics described above, combined with both direct ('live') N-body integration of seed BH trajectories and semi-analytic orbit integration in post-processing, to follow a wide range of possible BH seed populations with different formation properties and locations. In post-processing, we apply a modified DF estimator developed in a companion paper (Ma et al., in preparation), which is more flexible, accurate, and computationally efficient. In Section 2, we describe our numerical simulations and the semi-analytic post-processing method.

The plan of this paper is as follows: In Section 3, we present the results from simulations and semi-analytical integration of sample orbits, and show that seed BHs are generally not able to sink efficiently or be retained even at high seed masses. In Section 4, we discuss possible solutions to this problem, but also use our simulations to highlight how these solutions encounter still other problems. We summarize in Section 5.

Throughout, we assume a standard flat Lambda cold dark matter cosmology with m = 0.31, = 1m , b = 0.046, and H 0 = 68 km s -foot_0 Mpc -1 (e.g. Planck Collaboration VI 2020).

The simulations we study are re-simulations of the high-redshift (z > 5) galaxies presented in Ma et al. (2018aMa et al. ( , 2018bMa et al. ( , 2019) ) based on the Feedback In Realistic Environments (FIRE; Hopkins et al. 2014Hopkins et al. , 2018) ) project. 1 Specifically, we re-simulate the cosmological zoomin simulations centred around the galaxies 'z9m12a' and 'z5m12b'. Each of these represents a galaxy which has reached a halo mass 10 12 M , a stellar mass > 10 10 M , and a star formation rate 150 M yr -1 by redshifts z ∼ 9 and 5, respectively. As discussed in Ma et al. (2019), these are chosen to be plausible analogues to the observed hosts of the highest redshift, brightest QSOs. We note that while there are many other well-resolved galaxies in each cosmological zoom-in volume, we follow the most massive galaxy as it is the best candidate for a QSO host (but our conclusions about failure of BHs to 'sink' are even stronger in lower mass galaxies).

The simulations are run with an identical version of the GIZMOfoot_2 code (Hopkins 2015) to their original versions in Ma et al. (2018b). We use the mesh-less finite-mass (MFM) mode for solving hydrodynamic equations, with the identical FIRE-2 implementation of star formation and stellar feedback. The detailed baryonic physics included are all described extensively in Hopkins et al. (2018), but briefly summarized here. Gas cooling includes a variety of processes (molecular, atomic, fine structure, recombination, dust, free-free, Compton, etc.) accounting for 11 separately tracked species (H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe), following the meta-galactic ultraviolet background from Faucher-Giguère et al. ( 2009) with selfshielding. Stars are formed on the free-fall time from gas which is locally self-gravitating, molecular/self-shielded, denser than n > 1000 cm -3 , and Jeans-unstable following Hopkins, Narayanan & Murray (2013). Each star particle, once-formed, represents an IMFsampled population of known mass, age and metallicity, and we explicitly account for stellar mass-loss (from OB and AGB outflows), core-collapse and Ia supernovae, and radiative feedback (in the forms of photoionization and photoelectric heating, and single and multiple-scattering radiation pressure), with rates tabulated from standard stellar evolution models (Leitherer et al. 1999).

The only difference between our simulations and those in Ma et al. (2018b) is that we re-run them including a 'live' model for the formation of a broad spectrum of BH seeds, which are allowed to follow the full N-body dynamics. We emphasize that we do not artificially 'force' BHs to follow the potential minimum or decay their orbits via any prescriptions of sub-grid DF, as in some cosmological simulations (e.g. Hopkins et al. 2005Hopkins et al. , 2006Hopkins et al. , 2008a;;Springel, Di Matteo & Hernquist 2005;Sijacki et al. 2015;Anglés-Alcázar et al. 2017a).

We form BH seeds as follows: Whenever gas meets all the star formation criteria above and is about to be transformed into a star particle, it is assigned a probability of instead becoming a BH seed. Instead of setting the probability as an adjustable constant as in, e.g. Bellovary et al. (2011), it is weighted so that BH seeds form preferentially at the lowest metallicities (Tremmel et al. 2017) and highest surface densities/gravitational acceleration scales: specifically, we adopt p ∝ exp (-Z/0.01 Z ) [1exp (-/ 0 )], where ∼ M/R 2 is integrated to infinity with the Sobolev estimator from Hopkins et al. (2018) and 0 = 1 g cm -2 , with 0.01 Z = 1.4 × 10 -4 . The metallicity weighting is motivated to be consistent with our current understanding of seed BH formation models, all requiring low metallicities. For instance, Pop III stars and direct collapse models require low-metallicity primordial gas, while models of runaway mergers in star clusters strongly favour low-metallicity due to the lower mass-loss of massive stars in such environments (González et al. 2021). The value of 0 is specifically chosen because it is the density where analytic models (Fall, Krumholz & Matzner 2010) and numerical simulations (Geen, Soler & Hennebelle 2017;Grudić et al. 2018;Kim et al. 2019) of individual star formation and BH growth have shown robustly that stellar feedback fails to 'blow out' gas from the region efficiently, leading to runaway collapse/accretion. Exceeding this limit is required in many (but not all) models for massive BH seeds, either to prevent extended accretion disks from being destroyed by radiation from the accreting proto-quasi-star in direct collapse models, or as a necessary requirement to form superdense star clusters, which are the essential prerequisite for star cluster-based IMBH formation models (e.g. runaway merging) to initiate rapid growth (see Grudić et al. 2018;e.g. Shi et al. 2021). The normalization of p is chosen to form the maximum number of seeds before they begin to represent an appreciable fraction of the total galaxy mass and therefore perturb the dynamics. If the particle is selected to become a BH seed, then we draw a BH seed mass uniformly in log M from M = 10 3 -10 7 M .

Because we wish to only study the dynamics of BH seeds, we ignore BH accretion or feedback. These will be studied in future work.

Our 'default' simulations have an approximately constant baryonic mass resolution of m i ∼ 7000 M and a five times higher DM resolution. This is sufficient to explicitly resolve N-body DF and other effects on the more massive seeds ( 10 5 M ) we simulate: Depending on the details of the gravity scheme, one generally achieves this for seed masses M (10-100) m i . 3 To assess the effects of resolution on the dynamics of lower mass BH seeds, we briefly re-simulate one of our galaxies after applying a super-Lagrangian (AMR-like) refinement step (e.g. Anglés-Alcázar et al. 2021), to run with 800-M baryonic resolution 4 , and measure whether there is any significant difference in the 'sinking rate' of seeds at any BH mass after 100 Myr. We find no measurable difference. There is a simple reason why the detailed numerical accuracy of the DF forces on such low-mass seeds has little effect: The actual DF time for low-mass seeds (with e.g. M 10 5 M ) is far longer than the Hubble time at these (high) redshifts, so DF plays an essentially negligible role in their dynamics on a galactic scale.

Several authors who have implemented DF as a sub-grid routine (e.g. Pfister et al. 2019) pointed out that sub-grid corrections could make a difference in the seed BH orbits. This may also be an issue for the accuracy of direct simulations, especially for low-mass seed dynamics. It is therefore useful to check the validity of our 3 We enable the additional improvements to the gravitational time-step criteria, tidal force treatment, tree-opening, and integration accuracy detailed in Guszejnov et al. (2020), Grudić & Hopkins (2020) where they were developed for simulations of star formation which require accurate evolution of stellar binaries and multiples, and set the force softening of the BH seeds to a very small value (10 -3 pc) to represent real sink particles while using adaptive force softening for all other types to represent a smooth background. Detailed studies have shown that using adaptive softening as we do to ensure a smooth background force and with the more strict time-step and integration accuracy criteria used here, DF-like forces can be accurately captured for BHs with masses 10 times the background particle mass, while with less accurate integration often used in cosmological simulations which do not intend to resolve few-body effects, the pre-factor is more like ∼100 (van den Bosch et al. 1999;Colpi et al. 2007;Boylan-Kolchin, Ma & Quataert 2008;Hopkins et al. 2018;Pfister et al. 2019;Barausse et al. 2020;Boldrini et al. 2020). 4 Since the gravitational acceleration for BHs we study is strongly dominated by baryonic masses near the Galactic Centre (we confirm the N-body forces from dark matter are sub-dominant by order-of-magnitude or more), we did not refine the dark matter mass/force resolution in these re-simulations, as it is largely irrelevant to our conclusions. The image shows the clumpy structure of high redshift galaxies. Right-hand panel: The BH particles in this simulation at this particular snapshot, ranging from 10 3 to 10 7 M , covering a wide range of possible masses from different seed BH formation scenarios. BHs appear mostly randomly distributed in the galaxy, but with enhanced clustering near the Galactic Centre. However, we do not see significant seed-BH mass dependence, and the apparent Galactic Centre clustering simply reflects the overall concentration of mass (the galaxy effective radius here is ∼ kpc).

simulations with some alternative approach. Hence, we implement a semi-analytic analysis for the dynamics of BH seeds in postprocessing, both as a check of our direct numerical simulations, and a way to gain analytic insight and explore even larger parameter spaces prohibited by the resolution and computational expense of our simulations. In post-processing, we can create an arbitrary sample of BH seeds at any desired time, and evolve them in time-independent potentials taken directly from the numerical simulations, allowing us to map the dynamics in detail.

To do so, we re-calculate the trajectories of 100 BH 'test particles,' taking background potentials from the simulations and adding an analytic DF force explicitly in post-processing, during which we apply a newly developed DF estimator that is discussed in a companion paper (Ma et al. in preparation). We approximate the N-body dynamics of a seed of mass M with an acceleration a M = a ext + a df , where a ext is the 'normal' external gravitational acceleration on a test particle (computed identically to how the forces are computed in-code, for the adaptively force-softened potential from all N-body particles in the simulation). Then a df is the 'DF force' -the next-order (non-linear) term which represents the drag force arising from deflection of bodies by M. Specifically, we adopt the following expressions which can be directly computed from the simulation data (either on the fly or in post-processing): Here, a ext and a df are defined as a sum over all N-body particles i, with N-body masses m i , relative position

is the impact parameter. S i (r i ) is the usual dimensionless forcesoftening kernel to prevent numerical divergences, defined as

We refer interested readers in our expression for a df to the companion paper (Ma et al. in preparation). But briefly, our expression reproduces exactly the classical (C43)

V in cases consistent with the assumptions of C43, i.e. when the background distribution function is spatially homogeneous (constant density and velocity), timeinvariant, Maxwellian, and single-component. But it allows more naturally for cases which violate these conditions. Our expression also removes the ambiguity of the C43 expression in estimating a number of ill-defined continuum quantities, when applied to discrete simulation N-body data (e.g. how and on what scales to evaluate ρ, σ , andV; what value of to use). Usually, α i 1 such that a df ∝ α -1 i ∝ M, which means as expected that the DF acceleration is the largest for the most massive BHs, and potentially negligible for small BHs.

Here, we present the results from direct simulations, focusing on the clustering behaviour of BH particles. In Fig. 1, we show a projected image of the galaxy 'z9m12a' at redshift z = 10.4, as a typical high-redshift snapshot in our simulations. The left-hand panel shows the total non-BH mass (i.e. including dark matter, gas, and stars) density distribution, with the Galactic Centre located at the origin. The image shows the extremely clumpy appearance of typical highz galaxies, with multiple local density maxima near the Galactic Centre, consistent with both other simulations and observations. In the right-hand panel, we overplot the positions of BH particles near the Galactic Centre. The colour labels their masses, ranging from 10 3 to 10 7 M , which cover a wide range of seed BH masses from different formation scenarios. There is no significant position dependence upon mass for BH particles in the galaxy, with some mild clustering near the Galactic Centre. No significant mass dependence is observed.

To analyse the sinking problem of seed BHs, we show the magnitudes of galactocentric distance r and velocity v of BH We define a BH particle being trapped and efficiently sinking if it is located within <0.5 kpc from the Galactic Centre with a speed less than 10 km s -1 (shaded area). The colours label the mass of each BH particle. From our simulations, there are no BH particles trapped in this manner, nor any significant dependence on their masses of their positions and velocities.

particles selected from nine different snapshots in Fig. 2. Specifically, the BH particles are selected from snapshots in 'z5m12b' at z = 9.0, 7.7, 7.0, 5.9, and 5.0, and snapshots from 'z9m12a' at z = 10.9, 10.4, 9.9, and 9.5. Although snapshots at later redshifts contain BH particles that are already present at earlier redshifts in the same galaxy, the different snapshots are well separated in time such that the positions and velocities of these BH particles can be considered to be statistically independent. If a BH particle is located within 0.5 kpc from the Galactic Centre with a (relative) velocity less than 10 km s -1 (Fig. 2 shaded area), we consider it to have 'efficiently' undergone sinking and trapping in the Galactic Centre. Fig. 2 suggests that none of our BH particles in the mass range of 10 3 -10 7 M has achieved this at the redshift they are observed. There is also no clear dependence of BH positions and velocities on their masses, indicating their dynamics is basically independent of their masses if BH masses are below 10 7 M , i.e. the dynamics is dominated by the mass-independent external gravity, while the mass-dependent DF plays a negligible role.

Here, we present the results from semi-analytic post-processing, with our new DF estimator, to cover a wider range of BH masses. Specifically, we select snapshots from 'z5m12b' at z = 5.0, 7.0, 9.0, and 'z9m12a' at z = 9.5. In each snapshot, we place 100 test particles to integrate their dynamics, whose initial parameters are generated in the following way: The masses are randomly selected from 100 to 10 10 M (uniformly sampling log of mass), while the initial positions and velocities are chosen randomly from star particles in the corresponding snapshot, which is not only a convenient sampling method, but physically motivated since we would expect seed BHs are mostly born in similar locations to star clusters. With such sampling, are also able to study the effects of initial galactocentric distances/velocities on sinking (so we stress that our conclusions are completely independent of how we perform this sampling). In post-processing, we ignore the dynamics of background particles, i.e. we apply a time-independent gravity potential, as we would expect the static background to represent random sample of typical chaotic high-z galaxies, not an accurate reflection of some certain galaxy. The assumption of a static but realistically clumpy mass distribution allows us to gain insight into the effects of spatial inhomogeneities in the gravitational potential expected in typical, chaotic high-z galaxies. However, the orbits that we calculate in this way are not necessarily fully realistic since they neglect the time dependence of the potential. We note, though, that time dependence of the potential seems unlikely to accelerate sinking relative to a static-potential calculation; if anything time-dependence of the potential could further contribute to keeping seeds away from the Galactic Centre. The external gravity and DF are calculated by equation (1). Essentially, the difference between our 'live' dynamics simulations and these postprocessing calculations allows us to see how the time-dependence of the potential alters (in aggregate) the dynamics of sinking BH seeds.

To further see how the 'clumpiness' of the potential alters the BH dynamics, we re-run our semi-analytic orbit integration in a 'spherically smoothed' version of the potential. In these calculations, we take the exact same spherically averaged mass profile from the full simulation snapshot studied above, ρ(r) ≡ dM enc (< r)/4π r 2 dr in narrow radial annuli dr, and then use this as the background potential for our orbit integration. So, by definition, this has the same spherically averaged M enc (< r) and circular velocity profile, but no substructure.

In Fig. 3, we show several sample orbits for test particles of different masses in the z = 7.0 snapshot of 'z5m12b' overlaid on its mass density distribution. The orbits in the original snapshots are shown in the upper panel, while in the lower panel we show the trajectories integrated from the spherically smoothed version of this snapshot, with the same test-particle initial conditions. The thin lines show the trajectories and the black cross shows the final positions of test particles. The test particles follow chaotic orbits in the clumpy snapshot with no significant dynamical centre (as we would expect for a high-z galaxy). It appears that for the most massive test particles M 10 8 M , their velocities significantly decrease within a Hubble time at z = 7 (∼ 1 Gyr), and their final positions lie within the very central region of the galaxy. But there is no significant sinking for low-mass test particles. In the smooth galaxy, the particles behave similarly, yet it takes a shorter interaction time for the most massive test particles to sink. The velocity evolution of one particular test particle of 8.7 × 10 7 M is shown in Fig. 4, and it is shown that the velocity decay time-scale is about one order of magnitude shorter in the smooth galaxy compared to the clumpy galaxy. This suggests that the clumpy nature of early galaxies may increase the sinking time of seed BHs by an order of magnitude, by introducing chaotic dynamics to their orbits. In Fig. 5, we show the initial and final positions of all test particles we integrate in this particular snapshot, and its spherically smoothed version. We also show their initial and final velocity magnitudes as a function of mass in the lower panel. In the clumpy galaxy, while the test particles are randomly distributed in the galaxy initially, those with M 10 8 M show clustering behaviour near the centre after the integration, and their speeds decay to less than a few kilometers per second, indicating that they sink to the Galactic Centre after the integration. The remaining low-mass particles remain scattered around, with no significant decay of their speeds. The smooth potential reduces the minimum sinking Each panel is 8 kpc across in spatial scale. We find that in the high-z galaxy, the most massive test particles do sink to the Galactic Centre within a Hubble time at z = 7 (∼ 1 Gyr), while the low-mass seeds are simply experiencing chaotic orbits. In the smooth galaxy, the sinking behaviour is not very different for these five samples, yet for the massive seeds which are able to sink, their sinking time reduces drastically. This suggests that clumpy galactic backgrounds generally inhibit the sinking of massive seeds.

The evolution of the magnitude of the BH velocity as a function of interaction time for our integration of an 8.7 × 10 7 M test particle. We see that in both the clumpy and spherical smoothed galaxy, the velocity decays within 1 Gyr. But in the smooth galaxy, the decay time is lower by about one order of magnitude than the clumpy case, suggesting again the clumpy and chaotic nature of early galaxies may drastically increase the sinking time of seed BHs. mass to ∼ 10 7 M , when test particles are integrated over an order of the Hubble time at z = 7. 5 5 There is a trend of increasing final speed with test particle mass in Fig. 5 for the smooth galaxy. This turns out to be a reflection of the different integration time of these particles: we apply a time-step control proportional to |v/a| to avoid numerical errors. The massive particles, with larger DF (larger a), hence have smaller time-steps and shorter integration time compared to the less massive ones (see also the 'interaction time' label in Fig. 3), experiencing less deceleration in the integration. This effect does not appear in the clumpy galaxy, since the lack of dynamical centres of these galaxies makes the particle It appears that the clumpy nature of early galaxies may increase the 'minimum sinking mass' by one order of magnitude. It is worth noting that the sinking massive particles in the clumpy galaxy also do not sink exactly to the same place near the centre (as they do in the smooth galaxy). This implies that a clear definition of Galactic Centre with resolution of a few hundred pc is still ambiguous for these galaxies, and has potentially major implications for the demographics of BH-BH mergers at high redshift.

In Fig. 6, we show the initial and final magnitudes of galactocentric distance r and velocity v of all our test particles across different snapshots. The coloured points show the final velocities and distances of test particles while the thin grey lines connect their final values with initial values. We define the 'sinking' region in phase space, as in Section 3.1. Since we are covering a larger mass range of test particles than what we did in direct simulations for BH particles, some of the most massive particles do efficiently sink to the 'trapped region' this time. Specifically, particles with M 10 8 M sink to the centre region of the galaxy after the integration, regardless of their initial positions and velocities. For low-mass (M ࣠ 10 8 M ) particles, their final position and velocity distributions appear to be statistically similar to their initial configurations. This confirms the robustness of our results from direct simulations, in which all BH particles are less than 10 8 M and are therefore not experiencing significant sinking.

It is also worth noting that the sinking criterion almost depends entirely on the particle mass, not on initial velocities/distances to Galactic Centre. This is in contrast to what one would naively infer from the simplest DF-time calculations which assume a smooth potential with a constant circular velocity and BHs on slowly decaying nearly circular orbits, in which case the sinking time depends explicitly on the initial distances t sink ∝ r 2 (Binney & The magnitude of initial velocities and final velocities as a function of the BH mass. We see that for the clumpy galaxy, the high mass (M 10 8 M ) test particles sink to the Galactic Centre after the integration, while the low mass particles remain randomly distributed. For the smooth galaxy, the minimum mass for sinking reduces to M 10 7 M , about one order of magnitude lower. DF and sinking are negligible for the lower-mass seeds in both cases 5 . Figure 6. The initial and final magnitudes of velocities and galactocentric distances of all our test particles across different snapshots. The coloured points show the final velocities and distances (with any final velocities less than 10 -3 km s -1 interpreted as 10 -3 km s -1 for clarity). We define a BH particle as 'trapped' as in Fig. 2. The thin grey line connects the final properties with initial properties of each particle. The colours label the mass of each particle. We can see that after our integration nearly all particles with masses 10 8 M sink to the Galactic Centre (with a significant decline of velocity and distance), yet lower mass particles are still randomly distributed. Tremaine 1987). Physically, this can be explained by three factors:

(1) for highly eccentric or radial orbits, the dependence on initial radius is much weaker, independent of the assumed density profile or details of the DF scaling (Hopkins et al. 2008b); (2) the chaotic dynamics of seed BHs in clumpy (i.e. non-smooth) galaxies effective erase the memories of their previous orbits, which makes the initial positions less important to their orbital decay; and (3) the traditional r 2 dependence of t sink depends explicitly on the implicitly assumed isothermal mass density profile of the galaxy -but more generally the DF acceleration scales as a DF ∝ ρ(r)/v 2 c . In a clumpy high-z galaxy, however, the density ρ is not necessarily falling rapidly as in an isothermal sphere (and is not a trivial smooth monotonic function of galactocentric radius), again wiping out the naively predicted rdependence of t sink .

From both direct simulations and semi-analytic post-processing calculations, we have found that seed BHs less massive than 10 8 M generally cannot sink to Galactic Centres via DF in high-z galaxies.

To have at least one seed BH positioned in the Galactic Centre, so that it could accrete to ∼ 10 9 M and provide a plausible origin for luminous high-redshift quasars, we discuss two categories of possible solutions.

The first option is to use numbers as a trade off for efficiency: Although one low-mass seed BH is not likely to sink and accrete, a large number (which we estimate quantitatively below) of low-mass seeds could possibly give an opportunity for a 'lucky one' to sink and grow. Since the dynamics of BH particles and star particles are identically solved in our simulations (both as collisionless dynamics with external gravity), and the masses of star particles are around 10 3 M , below the low-mass end where DF drag is significant, we can use the star particles in our simulation as an ensemble of test particles to estimate the fraction of stars and therefore relics (ignoring processes like kicks), which can be trapped in local clustering structures ('clumps'). We apply such analysis to two particular snapshots, namely, 'z5m12b' at z = 7.0 and 'z9m12a' at z = 10.4.

We are only interested in clumps broadly near the Galactic Centre, hence, we identify the four densest clumps within 1.6 kpc near the Galactic Centre for each snapshot respectively, as shown in the upper panels of Fig. 7. The centre of the clumps are identified as the local density maxima, and their geometrical shapes are treated as spherically symmetric with radius 100 pc enclosing almost all of the clump mass, a fair approximation as shown in Fig. 7.

The lower left-hand panel of Fig. 7 shows the enclosed stellar mass and trapped stellar mass as a function of radius around each clump. If a star particle at radius r has a maximum possible apocentric radius r max from the clump centre (using the energy and angular momentum of each to evaluate its orbit, assuming the clump is static over its orbital time-scale), we then say it is instantaneously 7. We plot the cumulative distribution of distances between the centre-of-mass of the main clump progenitor and the newly formed star particle, at the time each star particle formed. We see that at least > 80-90 per cent of star particles in these clumps form 'in situ,' at distances 1 kpc from the clump centre. Only a small fraction are formed outside the clump and later captured. Of those, almost all form in the same galaxy at distances < 5 kpc (as opposed to in satellites or different progenitor galaxies).

enclosed within r and 'trapped' within r max . The gravity potential is calculated assuming a static potential around each clump with spherical symmetry (the clumps themselves, by definition, do not have substantial substructure). We see that the stellar masses in each clump (M enclosed (|r| < 100 pc)) range from 10 7 to 10 8 M . The mass fractions of trapped stars differ for different clumps and around 30-50 per cent of stellar mass could be trapped in an ∼0.1 kpc radius of the clumps, yet this value decreases as we go deeper into the clump centre, and the clumps could eventually trap only a few per cent of enclosed star particles within ∼50 pc. For all clumps, 90 per cent of their mass is in stars (as opposed to gas or dark matter).

Some low-mass objects are trapped in the dense clumps that represent the protobulge of these galaxies. But do they actually 'sink' or get trapped dynamically, or did they simply form in situ? To track the formation history of these star particles, we show their distances to their centre-of-mass at the particular redshift when most of them are just formedfoot_5 in Fig. 8. It turns out for almost all clumps, > 80-90 per cent of the star particles which we defined as 'trapped' in these clumps are formed within 1 kpc from the clump-progenitor centre-of-mass, which means most trapped star particles are formed in situ. The only seemingly exception is clump 6, where at first glance it appears that only about ∼ 70 per cent of the trapped star particles are in situ particles, but a detailed analysis shows that the remaining particles are actually formed in another clump which merges with clump 6, which does not challenge the conclusion (though it does relate to the hypothesis discussed in Section 4.1.2). Taken together, this means that while it is possible in principle for 'lucky' low-mass objects to be 'trapped,' it is quite rare: Comparing the total stellar mass of the galaxy to the mass of stars, which form ex situ and are trapped near clump centres yields a probability of about ∼10 -5 -10 -3 (depending on how generously we define 'trapped') for a low-mass seed formed randomly in the galaxy to migrate to being 'trapped' in the central < 100 pc of a clump by z ∼ 7.

Even if this occurs, the metallicity of the star particles which undergo this processes may create new problems for seed models. While the first Pop III stars or 'direct collapse primordial clouds,' which are candidates for forming massive seed BHs, could form very early at metallicities Z 10 -5 Z , the metallicity of star particles enclosed/trapped in clumps (even restricting to the 'ex situ' stars) is generally much higher, and turns out to be the highest for the most massive clump, as shown in the lower right-hand panel of Fig. 7. This indicates that the trapped star particles in these clumps may not represent a fair sample of the ex situ seed BH particles, which are formed before the clumps themselves are formed. The earliest forming stars are actually the least likely to be trapped in such clumps: they tend to form in mini-haloes at much earlier times and therefore across many different progenitors and thus have to migrate in from the furthest distances, while the 'ex situ but trapped' stars primarily still form in situ (in the same galaxy) just at distances of ∼ 1 kpc from the clump.

For all seed BHs, either in situ or ex situ, a related problem is related to the tension between the required clump masses and their ages. In many SMBH formation mechanisms, seed BHs have a higher probability both to be initially trapped and to subsequently accrete gas rapidly in the most dense/massive clumps, but these clumps are preferentially formed later, hence providing less time for BHs to migrate and to accrete. The average age of star particles inside clumps, as shown in the lower right-hand panel of Fig. 7, is far less than the Hubble time at the redshift we examined, providing a strict constraint on duty cycle if seed BHs are indeed hyper-Eddington accreting to become SMBHs in these clumps. Nevertheless, it is worth noting that SMBH seeding prescriptions are still highly uncertain, and other mechanisms may be able to circumvent these constraints.

From the semi-analytic calculations in Section 3.2, we have found that only seed particles as massive as 10 8 M can efficiently or reliably sink to Galactic Centres in a Hubble time. Such a large mass, however, is already an SMBH. On the other hand, our analysis in the previous section has shown that dense young star clusters as massive as 10 7 -10 8 M are present near the Galactic Centre. In the previous section we also show that most trapped star particles within those clumps are already formed in situ. This suggests another possibility: While randomly formed seed BHs are generally not massive enough to decelerate individually via DF, their preferential formation in tightly bound structures with large 'effective mass' is more realistic, as clusters could scatter with other components in the galaxy and sink effectively to the Galactic Centre. Indeed, in Ma et al. (in preparation), we show that the most-massive clumps do merge efficiently as these simulations are run to lower redshift and form the 'protobulge' of the galaxy.

There have been numerous papers arguing that runaway mergers in dense globular (star) clusters are a potential way to produce intermediate mass BHs (IMBHs, with typical masses 10 2 -10 5 M , see e.g. Portegies Zwart & McMillan 2002;Gürkan et al. 2004;González et al. 2021;Shi et al. 2021), which naturally becomes a preferential way to embed massive BH seeds in dense clusters as described above. Such channels, however, suffer from other problems like large gravitational recoils that can remove the formed IMBHs from the cluster (e.g. Holley-Bockelmann et al. 2008). There are also works arguing that gas accretion in nuclear star clusters (NSCs) and starburst clusters can also build up the mass of IMBHs rapidly (Kroupa et al. 2020;Natarajan 2021), which could be another way to apply this solution here. Yet observations have put upper limits on IMBHs masses (e.g. Lützgendorf et al. 2013Lützgendorf et al. , 2015;;Kamann et al. 2016;Zocchi, Gieles & Hénault-Brunet 2017), which introduce additional constraints on these channels. It should also be noted that, while globular clusters are usually assumed to be mainly pristine clusters that formed at very high redshift in mini-haloes, hence define an 'old' population for astrophysicists in the local universe, they are not so much older than the stars at z 7. In fact, the overwhelming majority of the clusters form in situ in the galaxy as it evolves from in situ gas, not from mini-haloes merging in. This means that the metallicity and timing problems discussed in Section 4.1.1 apply to this scenario, as well.

Our conclusions are consistent with other recent works focusing on slightly different aspects of this problem. For instance, Roškar et al. (2015) and Tamburello et al. (2017) study the co-evolution of SMBH pairs, finding that galactic clumps (originated either from high-z star-forming regions or a clumpy interstellar medium created by galaxy mergers) significantly perturb their orbital evolution, which potentially delay the decay process. Tamburello et al. (2017) and Tamfal et al. (2018) also point out that SMBH/IMBH pairs are still separated by 0.1-2 kpc after ∼ 1 Gyr in their simulations, which is consistent with our findings that no well-defined Galactic Centres can be identified on sub-kpc scales under these conditions. Bortolas et al. ( 2020) simulate a 10 6 M BH in a non-clumpy galaxy embedded in a cosmological environment at z = 6-7 and they show that DF torques are usually unimportant compared to the large-scale stochastic gravitational torques in determining the BH decay, even if no clumpy structures are considered. These works support to our conclusion that the chaotic structures of high-z galaxies could drastically change the sinking time-scale (hence the minimum sinking mass), if only DF is considered. Pfister et al. (2019) presented a complementary study to ours, focusing on more idealized simulations analogous to lower redshift systems, and a smaller number of test cases, but considering in more detail many of the numerical details of 'live' sub-grid BH DF treatments (e.g. explicitly adding an analytic DF force term in lowresolution simulations). They concluded that even in idealized galaxies designed by construction with a well-defined dynamical centre and a single, massive, centrally peaked bulge (e.g. an exponentialdisc and an Hernquist bulge), lower level clumpiness in the gas (e.g. GMCs with typical masses ∼ 10 5 -10 6 M ) would drive wandering or ejection of BHs with seeds less massive than ∼ 10 5 M . They hence concluded that 10 5 M is the minimum required mass for a BH to be well stabilized in the centre of its host. Since observed star-forming clumps or complexes are much more massive at high redshifts (e.g. Tacconi et al. 2010;Förster Schreiber et al. 2011;Swinbank et al. 2011), this criterion should only move to higher masses at high-z, consistent with our findings. Further, from postprocessing cosmological simulations of massive galaxies with welldefined dynamical centres merging at z < 6, Pfister et al. (2019) also concluded that it was crucial that BHs are already well-anchored to the galaxy centres before and throughout mergers, and that the centres are well-defined and dense enough to avoid tidal disruption, in order for BHs to 'sink.' They specifically concluded that it was crucial that BHs be embedded either in a dense satellite nucleus or a massive nuclear star cluster. This is essentially identical to our 'solution 2" above. Pfister et al. (2019) also noted that in the cosmologically simulated galaxies at earlier times, when the universe is < 1 Gyr old, even with their most massive (∼ 10 5 M ) seeds, the model for DF does not help in keeping BHs in the centre, as the galaxy is so chaotic that BHs wander no matter the implementation of DF. This is again in good agreement with our conclusion.

A recent study by Trebitsch, Volonteri & Dubois (2020) provides another excellent illustration of our key conclusions, in a single case-study of a galaxy simulation with 'live' AGN accretion and feedback. While the authors found that they could produce rapid BH growth by z ∼ 6, they (1) had to impose a sub-grid DF model with an artificial superlinear density dependence (∝ ρ 3 at high densities) designed to 'anchor' BHs into high-density regions (essentially our solution 2, again); (2) still found almost no BH growth until z ࣠ 8, after the galaxy reaches M * 10 9 M and forms a dense, strongly peaked and well-defined central 'protobulge' structure, very much like the late-time-forming structures we argue are necessary for BH capture and retention; and (3) still only reach peak luminosities 10 43 erg s -1 in X-rays, about a factor of ∼10 3 -10 4 less-luminous than the most luminous QSOs observed at these redshifts (Shen et al. 2019), which makes them still challenging to form.

There are some recent studies, which might appear to be in contrast to our results at first glance. For instance, Tremmel et al. (2018a) have shown that host galaxies could aid SMBHs to shorter sinking time-scales, and the ROMULUS simulations (Tremmel et al. 2017(Tremmel et al. , 2019;;Ricarte et al. 2019) argue that it is possible to grow massive BHs by intermediate redshifts. But a closer comparison shows these simulations are consistent with all of our key conclusions. In these studies, the BHs are, as the authors note (Tremmel et al. 2018a), embedded in nuclear regions of the host galaxy, which are dense enough to avoid tidal disruption and much more massive than the BHs. The nuclear regions, with high 'effective mass', hence sink as a whole -again following our 'solution 2' above. This is effective because these studies focus on cases where the galaxies are already massive, with unambiguous massive central peaks in their density profiles at relatively low redshift (with z ࣠ 2-4, cf. fig. 4 in Tremmel et al. 2018a). Moreover, in e.g. ROMULUS, the simulations have an effective seed mass ∼ 10 6 -10foot_6 M , 7 close to our sinking mass threshold in a smooth galaxy. These demonstrate that, given enough time and a preexisting massive density peak to 'anchor' an SMBH, BHs can indeed grow following e.g. our solution 2 as speculated above. Our focus here is essentially on how the 'initial conditions' of these simulations (at earlier times and smaller mass and spatial scales) could arise. We focus on galaxies at much higher redshifts, where those dense central regions either do not exist, or have formed relatively recently (e.g. z < 9) and one wishes to form an extremely massive SMBH by z > 7, significantly shortening the available time for BH growth, especially from extremely low-mass seeds.

In this study, we explore high-resolution cosmological galaxy formation simulations to understand the dynamics of BH seeds at high-z and their implications for SMBH formation and growth. Our simulations and semi-analytic DF calculations show that BH seeds cannot efficiently 'sink' to galaxy centres and/or be retained at high redshifts unless they are extremely massive already, M > 10 8 M , i.e. already SMBHs. We show that this threshold is at least an order-of-magnitude higher than what one would expect in a spherically symmetric smooth galaxy potential, as commonly adopted in analytic or older simulation calculations which could not resolve the complex, clumpy, timedependent sub-structure of these galaxies. For smoother galaxies, this mass threshold reduces to 10 7 M , which does not change the key conclusion.

We therefore join the growing number of recent studies by different groups which have reached similar conclusions (see Biernacki et al. 2017;e.g. Anglés-Alcázar et al. 2017b;Tremmel et al. 2018b;Bellovary et al. 2019;Pfister et al. 2019;Barausse et al. 2020;Boldrini et al. 2020). All of these studies, like ours, have concluded that this 'sinking problem' for BH seeds may, in fact, be even more challenging than even other well-known challenges for explaining the formation and growth of the first SMBHs with masses 10 9 M in galaxy centres at redshifts z > 7. Our contributions to extending this previous work include: (a) studying fully cosmological simulations with higher resolution, a broader range of redshifts, a much broader spectrum of BH seed masses, and different (sometimes more detailed) explicit models for stellar feedback; (b) comparing direct cosmological simulations which only resolved N-body dynamics to semi-analytic post-processing models for DF, to verify that these conclusions are robust; and (c) extending our comparisons to the 'test particle limit' by treating all stars as possible BH seeds.

Like these other studies, we qualitatively conclude that the chaotic, rapidly time-evolving, clumpy, bursty/dynamical nature of highredshift galaxies, coupled to the very short Hubble times ( 1 Gyr) make it nearly impossible for any lower mass seeds to efficiently 'migrate' from 1 kpc scales to galaxy centres, and is far more likely to eject seeds than to retain them. Like these authors concluded, the clumpy, bursty nature of the ISM is crucial for these conclusions: so this can only been seen in simulations which resolve the cold phases of the ISM and explicitly model stellar feedback. It is also worth noting that for low-mass galaxies (the progenitors where, in most models, seeds are supposed to have formed), even at z ∼ 0, clumpiness and burstiness are ubiquitous, and it is not simply a question of dynamical perturbations but even more basically of the fact that dwarf and high-redshift galaxies do not have well-defined dynamical centres to which anything could 'sink.' This is true even for well-evolved galaxies such as the LMC today.

In fact, we show that even the extremely massive BHs ( 10 8 M ) which do 'sink' actually do not sink to the same location at sub-kpc scales, where their migration stalls. This has potentially profound implications for LISA detections of SMBH-SMBH mergers in highredshift galaxies. Essentially, the 'last parsec problem' so wellstudied in the extremely dense, smooth, well-defined bulges of z = 0 galaxies (where the Hubble time is long) becomes a 'last kiloparsec problem' in these galaxies.

Solutions to the 'sinking problem' for SMBH growth/formation generically fall into one of two categories which we discuss in detail.

(1) Either seeds form 'in situ' when the massive bulge finally forms and creates a deep central potential, or a large number of seeds form so that even the infinitesimally small fraction which have just the right orbital parameters to be 'captured' by this bulge can exist. In either case, the problem is that we show this deep central potential well does not form until quite 'late,' at redshift z ࣠ 9, from gas and stars which are already highly metal-enriched (metallicities 0.1 Z ). This would mean popular speculative BH seed formation channels like Pop III relics or 'direct collapse' from hypermassive quasi-stars could not provide the origin of the SMBHs. Moreover, the combination of the fact that this occurs late, and that the stellar IMF is 'normal' at these metallicities, means that the 'time-scale' problem is much more serious: stellar-relic BHs, if primarily growing by accretion in such massive bulges, must grow from ∼ 10 M to 10 9 M in 200 Myr -requiring sustained highly super-Eddington accretion. Alternatively (2) 'seed' BHs must have enormous 'effective' masses to form early and remain 'trapped' and/or sink efficiently to the growing galaxy centre. Of course, BHs 'born' with M BH 10 7 M would solve this, but only by bypassing any stage that could be called a 'seed' (moreover, no serious models involving standardmodel physics can produce seeds of such large mass). However, models where seeds preferentially form tightly bound in dense star cluster centres owing to physics not modelled here (e.g. runaway stellar mergers in the centre of dense, high-z massive star clusters; see Shi et al. 2021) could (if the cluster is sufficiently dense) have an 'effective' dynamical mass for our purposes of roughly the cluster itself, which could reach such large values. This suggests these regions may be promising sites for SMBH seed formation.

In future work, we will explore the role of BH accretion and feedback, and more explicitly consider models where BH seeds form in resolved star clusters, as well as a wider range of galaxy simulations. It is likely that all of the scenarios above require a sustained period of super-Eddington accretion, so we will also explore whether this requires seed BHs residing (or avoiding) certain regions within high-z galaxies. We have also neglected models where non-standard model physics (e.g. dissipative dark matter, primordial BHs) allows for new formation channels and test-body dynamics. We will also explore new applications of our numerical DF approximator, in a variety of other interesting contexts (e.g. pairing of SMBHs in massive galaxy mergers at low redshifts).