Global torques and stochasticity as the drivers of massive black hole pairing in the young Universe

The forthcoming Laser Interferometer Space Antenna (LISA) will probe the population of coalescing massive black hole (MBH) binaries up to the onset of structure formation. Here we simulate the galactic-scale pairing of $\sim10^6 M_\odot$ MBHs in a typical, non-clumpy main-sequence galaxy embedded in a cosmological environment at $z = 7-6$. In order to increase our statistical sample, we adopt a strategy that allows us to follow the evolution of six secondary MBHs concomitantly. We find that the magnitude of the dynamical-friction induced torques is significantly smaller than that of the large-scale, stochastic gravitational torques arising from the perturbed and morphologically evolving galactic disc, suggesting that the standard dynamical friction treatment is inadequate for realistic galaxies at high redshift. The dynamical evolution of MBHs is very stochastic, and a variation in the initial orbital phase can lead to a drastically different time-scale for the inspiral. Most remarkably, the development of a galactic bar in the host system either significantly accelerates the inspiral by dragging a secondary MBH into the centre, or ultimately hinders the orbital decay by scattering the MBH in the galaxy outskirts. The latter occurs more rarely, suggesting that galactic bars overall promote MBH inspiral and binary coalescence. The orbital decay time can be an order of magnitude shorter than what would be predicted relying on dynamical friction alone. The stochasticity, and the important role of global torques, have crucial implications for the rates of MBH coalescences in the early Universe: both have to be accounted for when making predictions for the upcoming LISA observatory.


INTRODUCTION
The coalescence of massive black hole (MBH) binaries will constitute one of the loudest gravitational-wave (GW) events in the low-frequency band of the forthcoming spaceborne observatory Laser Interferometer Space Antenna (LISA; Amaro-Seoane et al. 2017; Barack et al. 2019). LISA will detect mergers between 10 4 -10 7 M MBHs up to and above z ∼ 20, thus probing the galaxy and MBH clustering up to the onset of structure formation.
The evolution of MBH pairs formed as a result of a E-mail: elisa.bortolas@uzh.ch galaxy merger has been first pictured in the pioneering work by Begelman et al. (1980): in simple terms, dynamical friction (DF; Chandrasekhar 1943) brings the MBHs down to the scale at which they form a bound system; shortly afterwards, the binary energy and angular momentum get drained via repeated three-body scatterings with interacting stars (e.g. Saslaw et al. 1974) or via its interaction with gas; finally, at roughly mpc scales, GWs start dominating the evolution and lead to coalescence (e.g. Thorne & Braginskii 1976).
The details of this rather simple picture have been largely investigated in the past few decades. A number of pieces of literature focused on the relatively small scale evolution (100-0.01 pc), and in particular on the effects of gas (e.g. Armitage & Natarajan 2002;Escala et al. 2005;Dotti et al. 2007;Mayer et al. 2007;Goicovic et al. 2017) and on the efficiency of three-body stellar scatterings in ensuring the GW phase is eventually reached (e.g. Milosavljević & Merritt 2003;Khan et al. 2011;Vasiliev et al. 2015;Bortolas et al. 2016;Gualandris et al. 2017;Bortolas et al. 2018a,b). The majority of these studies hint to an effective inspiral in both stellar and gaseous environments, with characteristic time-scales ranging from ∼10 Myr to a few Gyr, depending on the properties of the background environment.
The larger-scale, DF-driven inspiral phase has been long regarded as relatively easy to model: a vast number of works confirmed to a good extent the well-known Chandrasekhar predictions for the decay time-scale (e.g. Boylan-Kolchin et al. 2008), highlighting a prompt inspiral for galaxy collisions with mass ratios above 1:5 (e.g. Taffoni et al. 2003), and a less efficient pairing otherwise, as the DF-induced deceleration scales with the inverse of the perturber's mass. Moreover, a smaller intruder galaxy is likely to get severely tidally-and ram-pressure-stripped prior to reaching the main galaxy's interiors (e.g. Callegari et al. 2011;Van Wassenhove et al. 2014;).
However, nearly all studies on the topic have been carried out in very idealized frameworks: either two MBHs are placed in an equilibrium (galaxy) model (e.g. Gualandris et al. 2017;Tamburello et al. 2017), or two galaxies and their embedded MBHs are allowed to come together in isolation (e.g. Capelo et al. 2015;Bortolas et al. 2018b). This idealized approach has the advantage of requiring reasonable computational resources while allowing to study the MBH pairing down to sub-kpc and sometimes even sub-pc scale (e.g. Khan et al. 2018). On the other hand, this strategy discards the cosmological context and its implications, e.g. the possibility of distorted, non-symmetric galaxy morphologies, galaxy fly-bys, morphological transformations, and so forth. Accounting for the larger-scale environment is important when addressing the MBHs' path to coalescence in the present-day Universe, and becomes crucial in the infant Universe that LISA is going to probe.
In fact, even in the absence of a complete cosmological scenario, a series of studies on the pairing of MBHs in star-forming galaxies at z ∼ 1-2 highlighted that the formation of numerous massive star-forming clumps within the galaxy could scatter the inspiralling MBH orbit and render its evolution stochastic at the scale of the circumnuclear disc (Fiacconi et al. 2013;del Valle et al. 2015;Souza Lima et al. 2017) or even the entire galactic disc (Roškar et al. 2015;Tamburello et al. 2017). Khan et al. (2016) were the first to address the MBH pairing from cosmological distances all the way down to the hardening and coalescence phase. However, they focused on a very massive galaxy at z ≈ 3 [this would become a brightest cluster galaxy (BCG) by z = 0, see Feldmann & Mayer 2015]; this relatively rare object, which underwent quenching after a merger-driven starburst, rapidly exhausted its gas reservoir, allowing the authors to turn the calculation into a pure N -body one, thus discarding the role of hydrodynamics. They found a surprisingly short coalescence time-scale (∼10 Myr), about two orders of magnitude smaller than what was previously obtained in idealized merger simulations at z ≈ 0 (e.g. Khan et al. 2011).
More recently, attempts to account for the cosmological framework including hydrodynamics started to come forward (e.g. Tremmel et al. 2018a,b;Pfister et al. 2019). These studies hint to a more complex evolution compared to what would be predicted via a simple application of the theory of DF (e.g. Chandrasekhar 1943;Ostriker 1999;Colpi et al. 1999) and suggest that the time-scale for the pairing in these more realistic environments is much less predictable, especially for the relatively low-mass (<10 5 M ) MBHs that experience a weaker DF (e.g. Bellovary et al. 2019). However, the mass and spatial resolution limits in the aforementioned studies often require to introduce DF as an external drag force in order to properly follow the orbital decay of MBHs at sub-galactic scales.
Motivated by these considerations, here we follow the kpc-scale pairing of LISA MBHs in a typical, non-clumpy main-sequence star-forming galaxy at z = 7-6 using unprecedented mass and force resolution.

METHODS
We study the evolution of the pairing of MBHs embedded in a star-forming main-sequence galaxy at z 7. Our initial conditions are taken from the z = 7.32 snapshot of the high-resolution cosmological zoom-in simulation PonosHydro (PH; presented in Fiacconi et al. 2017) -the precursor run -to which we add the MBHs, as described in detail in Section 2.2.1. We run our simulations using the smoothed particle hydrodynamics, N -body code gasoline2 (Stadel 2001;Wadsley et al. 2004Wadsley et al. , 2017, using a setup that is virtually identical to the parameters set for the PH run in Fiacconi et al. (2017). The only different parameter between our runs and the one presented in Fiacconi et al. (2017) is the maximum allowed timestep, which in our case is ≈ 0.2 Myr, whereas it was ≈ 16 Myr in the precursor run. This choice allows us to follow the orbital evolution of the MBHs in much greater detail.
In the following sections, we describe the precursor and current runs. Unless otherwise stated, all quantities are expressed in physical units, the time coordinate has to be intended as the time since the Big Bang, and the adopted Λ-CDM cosmology is consistent with the results of the Wilkinson Microwave Anisotropy Probe 7/9 years (Ωm,0 = 0.272, ΩΛ,0 = 0.728, Ω b,0 = 0.0455, σ8 = 0.807, ns = 0.961, and H0 = 70.2 km s −1 Mpc −1 ; Komatsu et al. 2011;Hinshaw et al. 2013).

Precursor run
PH is a high-resolution, hydrodynamics, N -body version of the dark matter-only zoom-in cosmological simulation PonosV from the Ponos suite by Fiacconi et al. (2016), which was run in a box of 85.5 cMpc using the code pkdgrav3 (Potter et al. 2017).
The initial mass of gas and (high-resolution) dark matter particles in PH is 883.4 and 4397.6 M , respectively.
The physical gravitational softening of those particles is, for z ≤ 9, 47.9 pc for baryons and 81.8 pc for dark matter. 1 In the precursor run, gas can cool radiatively, via both non-equilibrium cooling of atomic primordial species and photoionization equilibrium cooling of metal species (Shen et al. 2013), assuming a redshift-dependent photoionizing UV background (Haardt & Madau 2012).
Star formation and stellar feedback are modelled following the recipes of Stinson et al. (2006). More specifically, gas particles are allowed to form stars only if the local gas density and overdensity are higher than 10 mH cm −3 and 2.64, respectively, the local gas temperature is lower than 10 4 K, and the flow is convergent and locally Jeans-unstable. If all these criteria are met, then gas particles are stochastically selected such that, on average, dM * /dt = * Mgas/tdyn, where M * and Mgas are the mass of stars and gas involved, respectively, * = 0.05 is the star formation efficiency, and t dyn is the local dynamical time. The newly formed stellar particles have a mass of 353.4 M : being much more massive than real individual stars, they represent a stellar population, which is described by a Kroupa (2001) initial mass function.
Stars with masses between 1 and 8 M release part of their mass as stellar winds, with the returned gas retaining the metallicity of those stars. Stars with masses between 8 and 40 M , depending on their mass-dependent lifetime (Raiteri et al. 1996), can explode as Type II supernovae: they inject into the surrounding gas a given amount of total, iron, and oxygen mass, dependent on the progenitor stellar mass (Woosley & Weaver 1995;Raiteri et al. 1996), and energy (10 51 erg), according to the 'blastwave model' of Stinson et al. (2006), in which cooling is temporarily disabled. Type Ia supernovae are also modelled and inject 10 51 erg and a fixed amount of mass and metals (Thielemann et al. 1986;Raiteri et al. 1996).
In addition to the standard SPH formulation, we include thermal energy and metal turbulent diffusion (with a diffusion coefficient C = 0.05; Wadsley et al. 2008;Shen et al. 2010), and a pressure floor (Roškar et al. 2015), to have a minimum number of resolution elements that must resolve the Jeans length and ensure the correct fragmentation behaviour of the gas. The gas becomes therefore non-ideal when the pressure floor is reached, since the temperature can cool below the ideal-gas equivalent of the pressure floor (see, e.g. Capelo et al. 2018).
For more details about the PH simulation and its dark matter-only counterpart (PonosV), we refer the reader to Fiacconi et al. (2017) and Fiacconi et al. (2016), respectively.

Current run
The simulation presented in this paper is relatively expensive and the human time needed for the integration can get to several months even with the best current present-day facilities. On the other hand, the result of a single run has to be taken with caution, as its statistical significance cannot be addressed.
For these reasons, here we adopt a workaround that allows us to study the pairing of multiple MBH pairs within a 1 For z > 9, these quantities need to be multiplied by 10/(1 + z). 7.26 7.08 6.91 6.75 6.59 6.45 6.31 6.18 z Figure 1. Angle θ between the gas angular momentum within the central 0.5 kpc and its value in the initial conditions as a function of time. The orientation of the gas angular momentum within this radius is used as a proxy for the galaxy orientation throughout the paper. single run: we place six secondary MBHs in the main galaxy and we study the pairing of each of those with the (unique) primary MBH. Our idea is justified by the different massscales in our system: that is, the mass of the typical particle in our system (m particle 4000 M in the high-resolution region) is much smaller than the mass we assign to the MBHs (∼10 6 M ) which in turn is much smaller than the mass of the host galaxy (Msystem ≈ 5 × 10 9 M within the central 1.25 kpc in the initial conditions). For this, unless the MBHs reach distances small enough to scatter each other, we can assume them to evolve independently. In Section 2.2.2, we describe in detail how we check for spurious scatterings.

Initialization of the MBHs
In what follows, we refer to the 'centre of the galaxy' as the centre of mass of all particles (excluding the MBHs) in the galaxy's innermost 0.5 kpc, iteratively recomputed until a tolerance of 0.1 pc is achieved; the velocity of the galaxy is defined as the mass-weighted velocity of the same particles. Analogously, the galaxy angular momentum (and galaxy plane) is obtained by computing the gas angular momentum within the central 0.5 kpc. These quantities are recomputed for each snapshot in the simulation and used as the reference frame for the analysis; in particular, we choose the z axis to be parallel to the galaxy angular momentum.
For reference, Figure 1 shows the evolution of the angle between the galaxy angular momentum and the same quantity computed at the beginning of the integration. This angle never exceeds 20 degrees, suggesting that the galaxy plane does not significantly oscillate about its initial position. The presented reference frame is also adopted for initializing the secondary MBHs. In fact, these objects respond to the bulk motion of the galaxy. We set the six secondary MBHs in two triangular configurations as shown in Figure 2: three of them (labelled as 1-3) have an initial separation from the centre of 0.75 kpc, whereas the other three (labelled as 4-6) are placed at a radius of 1.25 kpc. The MBHs initially have a tangential velocity equal to the velocity of a circular orbit given by the enclosed mass within their radius as if the galaxy were spherically symmetric: such velocity is equal to 116 and 128 km s −1 for MBHs initially at 0.75 and 1.25 kpc from the centre, respectively. The initial angular momentum direction of the MBHs coincides with the galactic one.
The primary MBH is not placed at the galaxy centre but in the densest region near it, found via the shrinking sphere mechanism, as this is a more likely region for its formation; its initial velocity is equal to the centre of mass velocity of particles in a sphere of 0.2 kpc around it. Note that the primary MBH and the centre of the galaxy are initially separated by a distance of 0.29 kpc, which gradually shrinks to zero in ∼200 Myr (the evolution of this separation in the run is shown in Figure 4).
The primary MBH has a mass M1 = 4 × 10 6 M i.e. 1/500 of the total mass in stars within the galaxy virial radius, in relatively good agreement with what suggested in Trakhtenbrot & Netzer (2010) and Lupi et al. (2019). The secondaries are assigned a mass M2 = 10 6 M ; they are inserted in the simulation without any host galaxy remnant, in the assumption that all the mass of the parent system that brought them into the main halo was stripped. Such stripping event is more likely for mergers with mass ratio 1:10 (e.g. Callegari et al. 2011;Van Wassenhove et al. 2014). However, we keep the mass ratio between the primary and the secondaries M2/M1 = 0.25, as this could mimic the presence of a stellar compact nucleus around the secondary, that sur-vived stripping; in addition, a large M2 enhances the effect of DF that is not completely resolved in our runs, as better discussed in Section 3.
The MBHs have a softening equal to 10 pc. MBH accretion and feedback are switched off in the present run, as this allows us to better control the effect of galactic structures on the orbital evolution, net of radiative effects; note that MBH feedback is suspected to have a small effect on the overall galactic structure if the MBH mass is smaller than ∼10 7 M (e.g. Anglés-Alcázar et al. 2017).

Check for MBH scatterings
The mass enclosed within a sphere of 0.1 kpc around each secondary MBH at the beginning of the integration is larger than 2M2, implying that the MBHs need to get closer than 0.1 kpc to scatter each other. Whenever this happens, we find the smallest distance attained by the two MBHs and compute the mass enclosed within a sphere of radius equal to their minimum distance (excluding the mass of the MBHs themselves). If this mass is smaller than 2M2, we take one of the two MBHs out of the simulation and we restart it before the scattering occurred.
Spurious scattering occurs only once in the run: MBHs 1 and 5 encounter each other at t ≈ 797 Myr; the former is thus removed. MBHs 3 and 4 also get very close to each other at t ≈ 795 Myr, but they do not meet the conditions for one of them to be excluded.
Note that, whenever a secondary gets close enough to the primary MBH, it enters a regime in which we are unable to resolve the subsequent, small-scale inspiral. In particular, we consider the large-scale inspiral to be completed when the periapsis reaches 0.1 kpc; we remove the MBH shortly afterwards. The orbit of the primary is not going to be strongly affected by the inspiralling secondaries, as there is a dense nucleus around it whose mass at the closest approach with each secondary is always larger than 4 × 10 6 M .
Another source of spurious scattering could be constituted by contaminating (low-resolution) dark matter particles; we checked that the contaminant fraction is always zero within the innermost 15 kpc, i.e. within a radius much larger than our region of interest.

RESULTS
The properties of the main galaxy are broadly detailed in Fiacconi et al. (2017): the galaxy harbours gas in a hot and turbulent state, and it lacks massive star forming clumps. It displays an evolving, geometrically thick stellar and gaseous disc whose typical stellar (gaseous) scale-length, scale-height, and total mass are, respectively, about 0.5 kpc, 0.2 kpc, and 5×10 8 M (0.8 kpc, 0.3 kpc, and ∼5×10 8 M ), whereas the galaxy's virial mass is ∼ 10 11 M ; these numbers are estimated at the beginning of the integration but can be referenced at least on an order-of-magnitude level for the whole integration.   7.26 7.08 6.91 6.75 6.59 6.45 6.31 6.18 z 7.26 7.08 6.91 6.75 6.59 6.45 6.31 6.18 z Figure 4. Orbital evolution of the MBHs as a function of time. The left-hand panels refer to the primary MBH (MBH 0) and to the MBHs with an initial separation from the centre of 0.75 kpc (MBHs 1-3), whereas the right-hand panels show the evolution of the MBHs initially displaced 1.25 kpc from the centre (MBHs 4-6). The top row shows the angle between the MBHs orbital plane and the galaxy plane; the central row shows the separation between each MBH and the galaxy centre; and the bottom row shows the distance of each secondary MBH from the primary. The hatched area in the bottom panels marks the region in which the orbital evolution of the secondaries cannot be resolved in our run, as their separation from the primary drops below 0.1 kpc. Figure 3 shows the position of the MBHs in the galaxy in three different snapshots of the simulation, whereas Figure 4 shows their orbital evolution as a function of time.

Black hole orbital evolution
The primary MBH initially oscillates around the centre of the galaxy, as its orbital decay is not completed yet. In fact, the nuclear overdensity in which the primary has been placed is the remnant of a lighter galaxy that merged with the main galaxy some ≈ 50 Myr before the start of the present simulation. The intruder's dense nucleus overtook the nucleus of the main galaxy in a nuclear coup fashion (Van Wassenhove et al. 2014), and it is on its way to the centre of the galaxy.
We now focus on the orbital-decay time of the secondary MBHs (last column of Table 1). Since MBH 1 is prematurely removed from the simulation, we can only conclude that its inspiral time would have been longer than ≈ 60 Myr. MBHs 2 and 3 take ≈ 50 and ≈ 110 Myr, respectively, to reach a separation smaller than 0.1 kpc from the primary. This factor >2 in the inspiral times suggests a certain degree of stochasticity in the orbital evolution, especially considering that both MBHs started at the same radius (0.75 kpc).
The effect of stochasticity is even more prominent in the orbital evolution of MBHs with an initial separation of 1.25 kpc from the centre: among those, MBHs 4 and 5 complete their inspiral, respectively, in 140 and 160 Myr (i.e. their inspiral time differs by only ∼15 per cent); on the other hand, MBH 6 never attains a separation smaller than 0.5 kpc from the centre, and its final radius (and angular momentum magnitude) is doubled if compared to the initial one; its inspiral time is likely going to be much longer than the simulated time of ≈ 210 Myr. We anticipate here that the overall orbital evolution of MBHs 4-6 (and perhaps even MBH 3) is significantly affected by the development of a strong bar whose dynamical effect can be inferred already at time ≈ 840 Myr, as discussed below; the presence of the bar becomes obvious even by eye-inspection after t = 880 Myr (see Figure 3).

Dynamical friction versus gravitational torques
In this section, we address the main mechanisms responsible for the MBHs orbital evolution. In particular, we compare the global gravitational torques 2 acting on each MBH to the DF-induced torques, whereas the effect of the bar alone is detailed in Section 3.3. The global gravitational torques acting on each MBH are obtained by computing the total gravitational acceleration atot due to all particles in the system, 3 including particles in the immediate vicinity of the MBH; this implies that atot contains the component of the acceleration associated to DF. The global torque acting on each MBH is then computed as τ tot = r × atot, where r is the MBH position vector and × is the cross product.
The DF-induced torque is obtained via the semianalytical approach detailed below.

Dynamical friction estimate
Our simulation features three different components that contribute to DF, i.e. gas, stars, and dark matter. The gasinduced deceleration can be estimated as (Ostriker 1999) where ρg is the gas density, M = |vg|/cs is the Mach number, cs is the speed of sound, vg is the local relative velocity between the MBH and the surrounding gas, and G is the gravitational constant. Note that this prescription is valid only if M > 1, whereas, for subsonic motion, we replace the natural logarithm in Equation (2) with M 3 /3. The value of the minimum and maximum impact parameters, bmin and bmax, are discussed below.
The DF-induced deceleration due to the two nondissipative components, i.e. stars and dark matter, can be written as (Chandrasekhar 1943) where the subscript x refers to quantities computed with respect to the stellar or dark matter background, ρx is the associated local density, vx is the relative velocity between the MBH and the surrounding medium x, and σx is the onedimensional velocity dispersion. The value of the maximum impact parameter, bmax, is set equal to 0.6 kpc for the stellar and gaseous component; this value corresponds to a few disc scale-heights, as the underlying assumption of quasi-homogeneous density fails at larger scales. The value of bmax is chosen to be 3 kpc for dark matter, i.e. of the order of the scale radius of a Navarro et al. (1996) profile whose fit was performed at the beginning of the integration. For each component, the minimum impact parameter is computed as where σ bg is the velocity dispersion of each background species.
In order to constrain the effects of DF, we compute the relevant quantities in Equations (2)-(4) within a sphere of radius R = 0.15 kpc around each secondary MBH. The Mach number is obtained by mass-averaging the speed of   sound of all particles in the sphere. This quantity is computed as cs = γakBT /(µmp), where T is the gas temperature, γa = 5/3 is the adiabatic index, kB is the Boltzmann constant, µ ≈ 0.6 is the mean molecular weight, and mp is the proton mass. The relative velocity (vg, vx) is obtained as the mass-weighted mean of the relative velocity between the MBH and the particles of each background component within the sphere; the associated velocity dispersion is estimated as i σ 2 i /3, where i runs over the three spatial directions.
Finally, the DF-induced torque can be obtained via for each component bg of the background (gas, stars, and dark matter), where r is again the MBH position vector. In order to assess the validity of our computation, we checked that the DF-torque given by the sum of the three background components does not change significantly (i) when varying the radius R of the sphere over which the relevant quantities in Equations (2) instead of taking the cross product. In all cases, we obtain a similar evolution for the DF-induced torque and we recover the results described in the next sections. Figure 5 shows the time evolution of the quantities used in Equations (2)-(3) to estimate DF; the bottom panels show the DF-induced torque due to each background component. Gas gives the dominant contribution for most of the time (84,55,73,79,89, and 50 per cent of the time for MBHs 1-6, respectively), followed by stars, whereas the dark matter contribution is almost negligible. The only exception is MBH 6, for which dark matter dominates DF 20 per cent of the simulation time: this is due to the fact that the MBH gets ejected in the outskirt of the system.
The significance of DF is very sensitive to the relative velocity between the MBH and the surrounding medium (v bg ). If, for simplicity, we focus on the gas component only, we can appreciate that |vg| is in the range that allows for both M > 1 and |vg| < 100 km s −1 for a large fraction of time (at least before the bar starts to have a significant dynamical impact, at t ≈ 840 Myr). The relatively weak DF in the very first Myr of the evolution is due to our set-up: the initial relative velocity between the MBH and the gas is close to zero by construction, hence the DF torque magnitude is low, as the motion is initially subsonic. Figure 6 compares the global gravitational torque experienced by each MBH to the DF torque obtained by the sum of the three contributing background species; in what follows, we will always refer to the DF-induced torque as the sum of these components. Beware that the global torque is computed from all particles, including the closest ones, thus it intrinsically includes the DF acting in the simulation. It is important to note that the DF might be underresolved in the simulation owing to its limited resolution. To verify this, we compare the DF computed either adopting the physical minimum impact parameter bmin given in Equation (4) or the maximum between the physical bmin and minimum separation that can be resolved in the simulation; the physical bmin ranges from a few to ∼20 pc for gas, and it is of the order of 1 pc for stars and dark matter; the minimum resolved separation in the simulation amounts to nearly 10 pc. 5 of the interacting particle and the one of the MBH (10 pc), as this is what happens in the integrator. We find that we under-resolve DF by ≈ 10-20 per cent only, as shown in Figure 6. In the remaining of the paper, we will always refer to DF computed adopting the physical bmin (Equation 4). Figure 6 clearly shows that the global torque experienced by each MBH is typically much larger than that from the DF; on average, the ratio between the global and DF torque amounts to 10-100 for MBHs 1-5 and to ∼300 for MBH 6. This is true in the case we consider the magnitude of the global gravitational torque (dark blue), and also when we compute the torque component along the MBH angular momentum direction (cyan); the latter happens to be very small for very short time ranges, and it is intrinsically more noisy, as it oscillates between positive and negative values. This fact suggests that the DF is not the main culprit driving the orbital evolution of the MBHs, which seems to be triggered by the inhomogeneities in the large-scale matter distribution, such as the bar and the spiral structures that develop in the galaxy (see Figure 3).

Dynamical friction and global torques: a direct comparison
Such consideration is further supported by the analysis of the local torque τ local . We define this quantity as the torque originating from a sphere around each MBH with twice the radius enclosing 2M2. We find that the magnitude of the local torque (also shown in Figure 6) is on average a factor of 10-20 smaller than the global one, meaning that the torque acting on each MBH is not dominated by the MBH DF wake. The local torque is sometimes notably larger than the DF torque; to explain this, note that τ local is estimated within a sphere of typical radius 100-150 pc; however, such radius can become ≈ 400 pc as the MBH wanders in a low-density region (note in particular MBH 6); when this happens, the sphere's radius is close to or larger than the disc scale-height, and the large density difference within the sphere may dominate the torque estimate. Also note that we are showing the modulus of the local torque and not its component along the MBH angular momentum.
A legitimate question is whether the effect of global 5 In fact, this value is computed by performing a mass-weighted average of the softening.  torques averages out over a full orbit. To investigate this aspect, we computed the torque as the time derivative of each MBH's angular momentum. 6 This latter approach would not be justified if we were to directly compare which effect is more relevant at each instant of time (as in Figure 6, as the derivative of the angular momentum implicitly assumes we are integrating the global torque over the frequency of our snapshots). However, computing the torque as the time derivative of the angular momentum is a better strategy when one wants to average the torque over a certain time-scale. Figure 7 shows the time evolution of the global gravitational torque component along the MBH angular momentum and the associated moving average (computed over a full azimuthal oscillation); the plot also shows the DFinduced torque, both instantaneous and averaged over a full azimuthal oscillation. While the global gravitational torque averages to a value close to zero prior to the bar formation, the same is not true after t ≈ 840 Myr; this is the time at which a stable barred structure starts to appear, and coincides with an enhancement in the torque experienced by each MBH. The figure suggests that the final orbital evolution of MBHs 4-6, and most probably even 3, has to be attributed to the presence of the bar, that produces a significant torque and determines the final fate of the MBHs.

Bar analysis and bar-induced torque
In order to better understand the development of the bar that dominates the dynamical evolution of the MBH (Figure 7), we characterize the stability of the disc in terms of the Toomre parameter, where Σ is the gas or stellar surface density, κ = 3(Vt/R) 2 is the epicyclic frequency defined via the azimuthal velocity Vt and the cylindric radius R (this approach is analogous to what has been done in e.g. Fiacconi et al. 2017;Ceverino et al. 2010); w is a constant equal to π for gas and to 3.36 for stars; finally, V is equal to the radial velocity dispersion σr,stars for stars, to the speed of sound cs for gas and, if the gas is turbulent, it is better approximated by V = c 2 s + σ 2 r,g where σr,g is the radial velocity dispersion of gas. We also compute the total Toomre parameter, evaluated as where Qstars (Qgas) is the stellar (non-turbulent gaseous) Toomre parameter (Romeo & Falstad 2013;Fiacconi et al. 2017). The Toomre parameter as a function of radius is shown in Figure 8 for three different snapshots. Q is computed in a disc of radius 2 kpc and half-height 0.8 kpc; we did not consider stars belonging to the small nucleus around  The Toomre parameter Q associated to gas is large, suggesting that gas is not in a gravito-turbulent state, as already broadly detailed in Fiacconi et al. (2017). The total and stellar Toomre parameters, instead, are only slightly larger than one: this hints to a nearly bar-unstable disc that can likely be triggered toward the formation of a large scale structure. In fact, a necessary, but not sufficient, condition for the formation of a bar is that Q < 2-3 (e.g. Binney & Tremaine 2008); whether the bar forms or not depends on other structural parameters of the galactic disc and halo, as well as on numerical details such as the force and time resolutions (e.g. Debattista et al. 2006;Dubinski et al. 2009).
The bar that develops in our simulation is clearly visible both in the stellar and in the gaseous projected density profiles. This suggests that the triggering for the bar could have been external, as gaseous bars are not likely to form in absence of an external perturbation (Mayer & Wadsley 2004). However, a much deeper analysis would be required to address the genesis of the bar: this aspect will be explored in a forthcoming paper.
In order to give a quantitative representation of the bar structure that is growing in the system, we perform a Fourier decomposition of the stellar surface density distribution Σ * , following in detail the process described in Zana et al. (2019). In particular, we compute the ratio between the second and the zeroth components of the Fourier development and the angular phase of the overdensity, as 7 Stars belonging to the nucleus have been selected as the stars closer than 0.3 kpc from the primary MBH whose circularity |jz/jc| < 0.5; here jz is the angular momentum in the z direction and jc the angular momentum of a circular orbit with the same radius. The angular momenta and radii here are computed in a reference frame centred on the position and velocity of the primary MBH. where the summation is carried over all the j-th stellar particles within an annulus of width 2 pc, height 1 kpc and centred at the cylindrical radius R. θj is the azimuthal angle of the j-th particle of mass mj. When the stellar disc hosts a barred overdensity over a certain radial range, the parameter A2(R) reaches a peak near the outer edge of the structure and the value of A2,max ≡ max[A2(R)] provides an estimate for the bar strength. Another outcome of the Fourier decomposition which is significant for our analysis is the parameter A2(< R), evaluated through Equation (8), but including all the particles that are enclosed within the radius R. In this study we compute the length of the bar by checking the fluctuation of the phase Φ(R). If R peak is the radius where the function A2(< R) has its local maximum, we define the bar extent, RΦ, as the radius R of the last annulus where the condition |Φ(R peak ) − Φ(R)| < arcsin(0.15) is valid, thus the phase of the stellar overdensity starts to vary significantly.
The time evolution of the parameters A2,max, and RΦ is shown in Figure 9. Given the chaotic nature of the high-redshift disc, with developing substructures and disturbances from the external environment, the assessment of the bar parameters is not straightforward. In order to avoid the noise due to the unsteady environment, we do not take into account any deviation from axisymmetry with A2,max(R) < 0.15. For this reason both the strength and length of the bar show significant ranges where it is impossible to retrieve their values unless a much deeper analysis is carried out, which is beyond the purpose of this work.
The bar seed could be already growing around t = 750 Myr but it is only after about 840 Myr that a stable structure is in place [A2,max(R) > 0.15] -this epoch is marked in Figure 7 in order to relate the raise of the torque magnitude with the presence of the bar. The bar rapidly reaches a strength A2 max 0.5 and a length of about RΦ 0.5 kpc, and maintains these values till the end of the run.

Stochasticity generators
Our analysis so far clearly suggests that the stochastic orbital evolution of the MBHs is to be attributed to the global torques arising from the whole host system, rather than DF alone. In the previous sections, we highlighted that the presence of the bar significantly contributes to the global torques. The torques that cannot be attributed to the bar itself (or the torques acting prior to the bar formation) have to be associated to some other perturbations in the potential.
We checked where the global gravitational torques mainly come from, at each snapshot of the simulation: most of them originate from the bar and spiral structures, whereas the impact of the central overdensity in which the primary MBH is embedded (which oscillates around the galactic centre for most of the evolution) seems not to have a major impact, unless the separation between the inspiralling and the primary MBH is closer than ∼0.2 kpc. We note that galactic fly-bys may also give a relevant contribution to the overall torque.
We should also mention that, to some level, DF and global torques are two sides of the same coin. Stochasticity is at least partially due to the fact that the MBHs are not evolving in a smooth background with relatively fixed properties; instead, the turbulent status of the medium implies that, at each time, the MBH finds itself in a background with a different sound speed, density, temperature; as a consequence, the significance of the DF itself is intrinsically stochastic. In short, gravitational inhomogeneities are sources of stochasticity as they induce random net torques, but the same inhomogeneities also randomise the significance of DF.
A legitimate question is to which degree the significant stochastic torques acting in the presently discussed highredshift galaxy are in place at later epochs, and especially in the present day Universe. On the one hand, it is clear that high-redshift systems are more chaotic and violent environments compared to their present-day counterparts (e.g. Table 1. Comparison between the inspiral time-scale inferred from DF and the effective inspiral time in the simulation. For each MBH (first column), the table shows the magnitude of the initial angular momentum j 0 (second column), the time T 0 to complete the first full azimuthal oscillation (third column), the average DF torque over T 0 (fourth column), the time-scale associated to the DF (fifth column), and the effective time for the inspiral (sixth column). The reported time-scales suggest that the DF is not the main responsible for the orbital decay. Mortlock et al. 2013), thus the overall significance of their stochastic torques is probably much larger than today. On the other hand, in this study we see most of the torque regulating the MBHs orbital evolution coming from the bar and associated spirals. Such structures are abundant in the local Universe (Kelvin et al. 2018 -a prominent example is our own Milky Way, Portail et al. 2017) and presumably at all times. It is thus reasonable to expect that the stochastic orbital evolution of MBHs, driven by bars and spiral structures, occurs at all cosmic epochs.

Orbital evolution and DF time-scale
The time-scale for the MBH orbital decay is the key parameter that encodes the efficiency of the MBH pairing. Table 1 compares the effective decay time of each MBH to the decay computed as j0/τDF,0, where j0 is the modulus of the MBHs initial angular momentum and τDF,0 is the magnitude of the torque associated to the DF, averaged over the first full azimuthal oscillation T0 (also in the table). The inspiral time one would infer relying on the DF only is a factor 4-15 larger than the actual inspiral time for the MBHs that manage to reach the centre of the system. 8 The final angular momentum of MBH 6 is twice its initial, and the DF it experiences by the end of the run is ∼100 times smaller than the initial one: its inspiral time is thus most likely going to exceed j0/τDF,0 by a large amount, unless global torques play in its favour in the subsequent evolution. In this respect, it is worth mentioning that a galaxy merger with mass ratio ∼1:5 is about to happen (≈ 50 Myr after the end of our integration the intruder enters the disc of the main galaxy) and the bar gets severely trimmed when this happens (see also Guedes et al. 2013), thus it may not influence the further orbital evolution of MBH 6.

Bars in the young Universe
What is clear from our analysis is that, once the bar forms, it heavily affects the orbital evolution and it is the most prominent torque source, together with the spirals. Given this, it would be important to infer the fraction of strong bars in the early Universe that LISA is going to probe. The fraction of barred galaxies amounts to > 30 per cent in the present day Universe (Consolandi 2016, and references therein). Although it is much harder to obtain similar statistics for the high-redshift Universe, numerical experiments suggest that the fraction of galaxies harbouring a bar is significant (≈ 40 per cent; Algorry et al. 2017;Rosas-Guevara et al. 2020) at z ∼ 1, thus it is reasonable to expect that such fraction would remain significant at larger redshift (z ∼ 7-6).
In general, high-redshift galaxies are expected to suffer much more frequent dynamical interactions due to external perturbations compared to low-redshift ones. Although external perturbations seem not to be critical for the formation of a stellar bar, as highlighted by recent cosmological zoomin simulations, some specific combinations of impact parameter and relative velocity between the main galaxy and an external perturber could boost the bar formation in potentials that are already bar-unstable (see Zana et al. 2018a for a recent discussion). Note that high-redshift bars could be frequently killed or trimmed by the frequent merger events and fly-by passages, therefore resulting short-lived. Nevertheless, they could re-appear multiple times if the external perturbations do not deeply modify the galactic potential, as found in Zana et al. (2018b). Therefore, concerning our case, the reason for the sudden weakening of the structure for a few Myr around t = 875 and 900 Myr (see Figure 9) could be due to some internal or external perturbations to the disc, that are strong enough to undermine the structural integrity of the bar, but not enough to ultimately alter the potential well.

Small-scale decay via stellar hardening
After the MBHs reach a pericentre of less than 100 pc, their dynamics cannot be resolved in our simulation. However, the fact that the primary MBH is surrounded by a dense nucleus dominated by stars allows us to infer an upper limit for the decay time down to the GW emission phase.
The MBHs need to get close enough to form a hard MBH binary, i.e. a binary whose semimajor axis is smaller than a h = GM2/(4σ 2 ) (Merritt 2013) which, in our case is ∼0.1 pc if σ is computed from stars in a sphere of 150 pc around the primary. Stars dominate the density profile at any time within 150 pc from the primary (the stellar mass fraction here is 60-80 per cent). Therefore, it is reasonable to consider stars as the sole contributors to the DF (the presence of a gaseous circumnuclear disc would probably boost the inspiral, Mayer et al. 2007, although see Souza Lima et al. 2017Lima et al. , 2020. In the nucleus, the stellar motion is isotropic to a good approximation, so that the average relative speed between the MBH and the medium coincides (10) here re is both the scale radius of the central stellar nucleus (whose properties are discussed below) and the separation at which the small-scale DF inspiral starts; re here is also used as the maximum impact parameter entering in the Coulomb logarithm, while the minimum one is set to 0.1 pc. The velocity dispersion around the primary MBH increases from about 50 to 150 km s −1 throughout the simulation, meaning that the DF-driven decay is of order ∼100 Myr. Afterwards, the binary separation could be reduced by three-body scatterings. The time-scale to reach the final coalescence can be estimated following Sesana & Khan (2015): where ρ infl and σ infl are the stellar density and velocity dispersion, respectively, within the binary influence radius r infl (i.e. the sphere containing twice the binary mass in stars), H ≈ 15 is a parameter obtained via scattering experiments (Quinlan 1996;Sesana et al. 2006), c is the speed of light in vacuum, and f (e) = (1 − e 2 ) −7/2 (1 + 73/24e 2 + 37/96e 4 ) is the eccentricity (e) enhancement function for GW emission (Peters 1964). The length-scale a * /gw is the typical separation at which the binary shrinking transits from being dominated by three-body scatterings to being mainly driven by GWs and practically represents the separation at which a MBH binary spends most part of the time during the stellar hardening evolution. Figure 10 shows the density profile of stars around the  (11). The first column shows the inner density slope γ of the assumed stellar density profile; the second, third, and fourth columns show the decay time-scale for MBH binaries whose eccentricity in the transition between stellar hardening and GW-decay is respectively equal to 0, 0.5, and 0.9. Details on the assumed profiles can be found in the text and in Figure 10.
primary MBH near the end of the integration (t = 920 Myr). We fit it using a Dehnen (1993) density profile, and obtain a total mass Mt = 4.3 × 10 9 M , a scale radius r0 = 0.1 kpc, and an inner density slope γ ≈ 0; the fitted profile is also shown in Figure 10. If one estimates the velocity dispersion via σ 2 infl = G(M1 + M2)/[(1 + γ)r infl ] (Alexander 2005), the decay time-scale obtained via Equation (11) would amount to the values shown in the first line of Table 2. We can note that only a large binary eccentricity 9 at the onset of the GW phase guarantees a decay time shorter than the previous evolutionary time-scales.
In fact, the limited resolution of the simulation does not allow us to properly resolve the density profile in the innermost 100 pc, which is likely steeper than what we found. For this, we extrapolated the inner density slope by using a Dehnen (1993) density profile with Mt obtained from the previous fit, but imposing an inner slope γ = 1 and 1.5, and ensuring in both cases that the radius containing 0.1Mt is the same in the fitted and extrapolated profiles. The extrapolated profiles are shown in Figure 10. The decay time-scales obtained via Equation (11) for the extrapolated profiles are also listed in Table 2. The decay time-scale is in these cases of the order of or shorter than the decay time of the previous stages, meaning that the coalescence would supposedly happen promptly once the hard binary is formed. Nonetheless, the hardening and inspiral time-scales we inferred in Table 2 are definitely longer than the ≈ 15 Myr found by Khan et al. (2016), who addressed the pairing of ∼ 10 8 M MBHs in a massive early type galaxy at z ≈ 3, which would become a BCG by z ≈ 0. The discrepancy in the time-scales has to be attributed to the large central density in their merging systems, which is nearly 10 5 M pc −3 at 10 pc (their figure spiral, even though in some cases they may ultimately hinder it; (iv) the effective decay time for 4 out of 5 inspiralling MBHs is a factor 4-15 shorter than what would be inferred relying on the DF alone; (v) the whole pairing process, including the final shrinking due to stellar interactions, occurs on time-scales not much shorter than the age of the Universe at z ≈ 6, and the final coalescence is likely going to occur at z = 5-3.5 (except for the case in which the bar ejects the inspiralling MBH).
Our simulation has the important limitation that MBH accretion and feedback were not taken into account. However, it is reasonable to expect that feedback would render DF even less efficient by evacuating the gas mass around an MBH (see, e.g. Park & Bogdanović 2017;Souza Lima et al. 2017; del Valle & Volonteri 2018 for a smaller scale analysis of this aspect). Each inspiralling MBH may also enhance its mass by accreting gas along the decay: in this respect, it is worth remembering that the MBH masses we adopted in the current run are relatively large for z 6, thus we do not expect this to represent a major issue. We also note that accretion on to MBHs may be used in future high-redshift surveys to sample offset or dual active galactic nuclei (e.g. De Rosa et al. 2019), in order to estimate how many MBHs could have been ejected by non-axisymmetric structures such as the bar described in this work.
A further simplistic assumption we made is that MBHs are initially placed on circular orbits aligned with the disc plane. Note however that the stochastic torques randomize the MBH orbits relatively quickly in the simulation, rendering them eccentric and at least mildly off-plane. 11 The possibility of an erratic orbital decay is not completely new to the study of MBH inspiral: as already mentioned, giant star-forming clumps have been proposed as further sources of stochasticity (Fiacconi et al. 2013;Roškar et al. 2015;Tamburello et al. 2017;Souza Lima et al. 2017); however, it is unclear whether the blobs found in observations of galaxies at z ∼ 1-2, that have been long interpreted as giant clumps, are instead under-resolved smaller object (e.g. Ivison et al. 2020). The existence of bars and spiral structures is obvious even in the low-redshift Universe, thus we believe the source of stochasticity in the current study to be more general, and potentially acting at all redshifts.
To conclude, our results clearly indicate that the standard (Chandrasekhar 1943 andOstriker 1999) prescriptions that model DF as a local drag force in analytical and semianalytical models can be greatly inadequate once the cosmological environment is taken into account, and especially in the young Universe that LISA is going to probe. A large effort is thus needed in this direction, in order to infer realistic time-scales for the inspiral in the presence of bars, spiral arms, and perturbations from external galaxies, as such disturbances may ultimately determine the final fate of an inspiralling MBH.