Feeding SMBHs through supersonic turbulence and ballistic accretion

Feeding SMBHs through supersonic turbulence and ballistic accretion

Alexander Hobbs, Sergei Nayakshin, Chris Power, Andrew King
Dept. of Physics & Astronomy, University of Leicester, Leicester, LE1 7RH, UK

It has long been recognised that the main obstacle to accretion of gas onto supermassive black holes (SMBHs) is large specific angular momentum. It is feared that the gas settles in a large scale disc, and that accretion would then proceed too inefficiently to explain the masses of the observed SMBHs.

Here we point out that, while the mean angular momentum in the bulge is very likely to be large, the deviations from the mean can also be significant. Indeed, cosmological simulations show that velocity and angular momentum fields of gas flows onto galaxies are very complex. Furthermore, inside bulges the gas velocity distribution can be further randomised by the velocity kicks due to feedback from star formation. We perform hydrodynamical simulations of gaseous rotating shells infalling onto an SMBH, attempting to quantify the importance of velocity dispersion in the gas at relatively large distances from the black hole. We implement this dispersion by means of a supersonic turbulent velocity spectrum. We find that, while in the purely rotating case the circularisation process leads to efficient mixing of gas with different angular momentum, resulting in a low accretion rate, the inclusion of turbulence increases this accretion rate by up to several orders of magnitude. We show that this can be understood based on the notion of “ballistic” accretion, whereby dense filaments, created by convergent turbulent flows, travel through the ambient gas largely unaffected by hydrodynamical drag. This prevents the efficient gas mixing that was found in the simulations without turbulence, and allows a fraction of gas to impact the innermost boundary of the simulations directly. Using the ballistic approximation, we derive a simple analytical formula that captures the numerical results to within a factor of a few. Rescaling our results to astrophysical bulges, we argue that this “ballistic” mode of accretion could provide the SMBHs with a sufficient supply of fuel without the need to channel the gas via large-scale discs or bars. We therefore argue that star formation in bulges can be a strong catalyst for SMBH accretion.

11footnotetext: E-mail: alexander.hobbs@astro.le.ac.uk

1 Introduction

The growth of supermassive black holes (SMBHs) at the centre of galaxy bulges is known to be correlated with observable properties of the host spheroid. The best established correlations are between the SMBH mass, , and the velocity dispersion, , i.e., the relation (Gebhardt et al., 2000; Ferrarese & Merritt, 2000), and the relation (Magorrian et al., 1998; Häring & Rix, 2004), where is the mass of the bulge. Understanding the formation and growth mechanisms of SMBHs is therefore believed to be important in determining the evolution of the larger host systems.

High luminosities of high redshift quasars at suggest SMBH masses of up to (Kurk et al., 2007), implying accretion onto the black holes close to or even above the Eddington limit (King & Pringle, 2007) for up to 1 Gyr. Such a strong inflow is most likely the result of major mergers between galaxies, which have been shown via hydrodynamical simulations to drive gas inward through gravitational tidal fields (Barnes & Hernquist, 1991). Simulations have also shown that galaxy mergers within large dark matter halos at may be able to drive gas in to the centre of the galaxy at a sufficient rate to grow a seed BH of to an SMBH of in the required time i.e., by (Li et al., 2007). Recent cosmological simulations of dark matter halos that include gas physics and black hole feedback processes are able to build up the required mass by via extended Eddington-limited accretion (e.g., Sijacki et al., 2009).

However, an ab-initio treatment of all the relevant physical processes (gravity, star formation, magneto-hydrodynamics, radiative feedback etc.) over the full range of scales required is still beyond current numerical capabilities. Therefore, little is known about the mechanism(s) that connect the large-scale flow ( kpc) to the small-scale accretion flow in the immediate vicinity of the SMBH ( pc). Current cosmological simulations require sub-resolution prescriptions to encapsulate the accretion process (e.g., Springel et al., 2005; Booth & Schaye, 2009; Sijacki et al., 2009). These models commonly use Bondi-Hoyle accretion (Hoyle & Lyttleton, 1939; Bondi & Hoyle, 1944; Bondi, 1952), corrected upwards by orders of magnitude when the resolution limitations cause an under-prediction of the desired SMBH growth rate. Recently, DeBuhr et al. (2009) have performed hydrodynamical simulations of major galaxy mergers employing a sub-grid model for viscous angular momentum transport to dictate the accretion rate, the results of which suggest that the Bondi-Hoyle prescription actually over-predicts accretion onto the SMBH by up to two orders of magnitude. The applicability of the Bondi-Hoyle solution to the SMBH accretion rate is therefore not clear.

Looking at the problem from the other end, where dynamical and viscous timescales are short compared with cosmological timescales, it is expected that geometrically thin and optically thick accretion discs (Shakura & Sunyaev, 1973) form due to a non-negligible amount of angular momentum in the gas. While it is possible to study such discs at a far higher resolution and with more physics, the approach is obviously limited by the fact that the disc is commonly assumed to just exist a priori. How these discs operate on pc scales is also not clear as the discs are susceptible to self-gravity and hence may collapse into stars (Paczyński, 1978; Kolykhalov & Sunyaev, 1980; Collin & Zahn, 1999; Goodman, 2003; Nayakshin & Cuadra, 2005).

There is thus a genuine gap between these two scales and approaches. The purpose of this paper is to try to bridge this gap to some degree with the help of numerical hydrodynamical simulations of gas accreting onto a SMBH immersed in a static, spherical bulge potential. In particular, we study the ‘intermediate-scale’ flow, from the inner 100 pc of a galactic bulge (just below the resolution limit for some of the better-resolved cosmological simulations) to the inner parsec, where the gas would be expected to settle into a rotationally-supported feature. Instead of assuming the presence of the disc, we allow it to form self-consistently from the larger scale flow. For simplicity of analysis of the results, we start with a spherically-symmetric shell of gas of constant density, where the gas is isothermal throughout the simulation. The first non-trivial element in our study is the rotation of the shell, which should force the gas to settle into a rotating disc.

While this turns out to be an interesting system in itself, forming a narrow ring rather than a disc, we feel motivated to consider a more complex initial velocity field in the shell. First of all, cosmological simulations show the presence of cold streams (e.g., Dekel et al., 2009) immersed into a hotter, lower density medium. These complex gaseous structures must also propagate to smaller scales within the bulge. Secondly, the abundance of gas required to fuel the SMBH should also fuel star formation in the host. Massive stars deposit large amounts of energy and momentum in the surrounding gas (Leitherer et al., 1992), which is probably one of the dominant sources of turbulence in the interstellar medium (e.g., McKee & Ostriker, 1977; Mac Low & Klessen, 2004). In general, where we can observe it, it is evident that turbulence is ubiquitous in astrophysical flows over a wide range of scales and systems (for a review see e.g., Elmegreen & Scalo, 2004). Strong supersonic turbulence is a key ingredient in all modern simulations of star formation (e.g., Krumholz & Bonnell, 2007; Bate, 2009).

Therefore, we must expect that, far from the SMBH, while the specific angular momentum of gas, , is too large to be captured by the SMBH within the central parsec (or less), the dispersion in is also large. It is not unreasonable to expect that there will be some gas counter-rotating with respect to the mean rotation of the flow. To test the importance of these ideas for accretion of gas onto SMBHs, we draw on the machinery developed for turbulent flows in the star formation field by adding a turbulent velocity spectrum to the gas in the initial shell. By varying the normalisation of the turbulent velocities we essentially control the dispersion in the initial gas velocities.

We find that the role of the turbulence and in general any disordered supersonic velocity field is two-fold. First of all, it broadens the distribution of specific angular momentum, setting some gas on lower angular momentum orbits. Secondly, as is well known (e.g., McKee & Ostriker, 2007), turbulence creates convergent flows that compresses gas to high densities. We find that the dynamics of such regions can be described reasonably well by the ballistic approximation in which hydrodynamical drag and shocks are of minor importance for the high density “bullets” moving through the lower density background fluid.

These two effects increase the accretion rate onto the SMBH by up to several orders of magnitude, largely alleviating the angular momentum barrier problem. Note that several authors have already suggested that star formation feedback and supersonic turbulence inside accretion discs (Wang et al., 2009) can promote SMBH accretion by amplifying angular momentum transfer (Collin & Zahn, 2008; Kawakatu & Wada, 2008; Chen et al., 2009). An astrophysical conclusion from our simulations is that star formation in the host can actually promote AGN accretion.

The paper is arranged as follows. In Section 2 we outline the computational method used, and in Section 3 we discuss our initial conditions. Section 4 presents an overview of the gas dynamics. Sections 5 & 6 present the results for the no turbulence and turbulent cases, respectively, and Section 7 details the fate of the gas that makes it inside the accretion radius. Sections 8 & 9 comprise our discussion and conclusion respectively.

2 Computational method

We use the three-dimensional smoothed particle hydrodynamic (SPH)/N-body code GADGET-3, an updated version of the code presented in Springel (2005). Smoothing lengths in the gas are fully adaptive down to a minimum smoothing length of pc, which is much smaller than the scales that we resolve in the simulation. The code employs the conservative formulation of SPH as outlined in Springel & Hernquist (2002), with the smoothing lengths defined to ensure a fixed mass within the smoothing kernel, for . Each simulation below starts with particles, with an individual SPH particle possessing a fixed mass of . The Monaghan-Balsara form of artificial viscosity is employed (Monaghan & Gingold, 1983; Balsara, 1995) with and .

The calculations are performed in a static isothermal potential with a central constant density core to avoid divergence in the gravitational force at small radii. We also include the black hole as a static Keplerian potential. The mass enclosed within radius in this model is:


where , , kpc, and pc and . The one dimensional velocity dispersion of this potential is km s.

As we concentrate on the hydrodynamics of gas accreting on the SMBH, we make two simplifying assumptions which avoid further non-trivial physics (to be explored in our future work). First, the position of the SMBH is held fixed during the simulations here. This is a reasonable approximation as all of our initial conditions are either exactly spherically or azimuthally symmetric, or have these symmetries when averaged over the entire simulation volume (e.g., turbulence is assumed to be isotropic in this case). Allowing the black hole to move self-consistently would require relaxing the static potential assumption.

Secondly, gravitational forces between the particles are switched off to avoid complications that might arise from gas self-gravity. While this could be viewed as a shortcoming of our work, we believe that inclusion of self-gravity, ensuing star formation and stellar feedback would only strengthen our results and conclusions. The gravitational contraction of gas clouds would create even higher density contrasts and thus make ballistic trajectories even more likely. Star formation feedback would drive its own turbulence and hence amplify the effects of the turbulence that we seed.

One further simplification is that the gas is isothermal throughout the entirety of the simulations. This is a fair assumption as cooling times are short for the high densities we are considering. Accretion of gas onto the SMBH is modelled with the “accretion radius” approach - particles that come within an accretion radius of pc are removed from the simulation and we track the total mass and each component of the net angular momentum vector of the accreted particles. This contrasts with the Bondi-Hoyle accretion formulations that are frequently used in cosmological simulations (e.g., Sijacki et al., 2009; Booth & Schaye, 2009). We believe the accretion radius approach, frequently used in the star formation field (e.g., Bate et al., 1995), is essential for the problem at hand as it prevents unphysical accretion of SPH particles with too large an angular momentum, whereas the Bondi-Hoyle formulation does not.

The units of length, mass, time and velocity used in the simulations are, respectively, cm (1 kpc), g (), Myr, and km s.

3 Initial conditions

The starting condition for all our simulations is that of a uniform density, spherically symmetric gaseous shell centered on the black hole. The inner and outer radii of the shell are and kpc, respectively, for most of the simulations (see Table 1 in the Appendix). In principle one can expect the outer radius of the shell to be much larger in a realistic bulge with effective radius of a few kpc, but we are forced to limit the dynamic range of the simulations due to computational resources. Furthermore, we believe we understand how our results scale with the outer radius of the shell (cf. Section 6.2.4), and hence the dynamic range limitation does not actually influence our conclusions.

The total mass of the shell is . To minimise initial inhomogeneities we cut the shell from a relaxed, glass-like configuration. The initial velocity field is composed of two parts: net rotation and a seeded turbulent spectrum. The net rotation is described by the azimuthal velocity component, , where is a parameter of the simulation (see Table 1 in the Appendix). In all our runs, is below the circular velocity at the shell radii, meaning that our initial conditions are not in equilibrium. We stress that this is deliberate, as our investigation here is the formation of a disc from infalling gas, rather than from already rotationally-supported gas.

Our approach to setting up the turbulent velocity field follows that of Dubinski, Narayan & Phillips (1995). We assume a Kolmogorov power spectrum,


where is the wavenumber. The key assumption here is that the velocity field is homogeneous and incompressible, and so we can define in terms of a vector potential such that . This is useful because the components of can then be described by a Gaussian random field with an associated power spectrum


This is a steep power-law, implying that the variance in will diverge sharply as decreases, and so we follow Dubinski, Narayan & Phillips (1995) and introduce a small scale cut-off . Equation 4 can then be written as


where is a constant that sets the normalisation of the velocities. For our purposes we set it equal to unity and normalise the velocity field once the statistical realisation has been generated. Physically the small scale cut-off can be interpreted as the scale as the largest scale on which the turbulence is likely to be driven.

Our approach to generating the statistical realisation of the velocity field is straightforward. First we sample the vector potential in Fourier space, drawing the amplitudes of the components of at each point from a Rayleigh distribution with a variance given by and assigning phase angles that are uniformly distributed between and . Second we take the curl of ,


to obtain the Fourier components of the velocity field. Finally we take the Fourier transform of this to obtain the velocity field in real space. We use a periodic cubic grid of dimension when generating the statistical realisation of the velocity field and we use tricubic interpolation to estimate the components of the velocity field at the position of each SPH particle.

The initial parameters for each run are given in Table 1 in the Appendix.

4 Overview of gas dynamics and main results

Before we embark on a quantitative study of the simulation results, we present several snapshots from the simulations that illustrate graphically the nature of the gas flow. In particular, here we discuss the overall gas dynamics for two simulations that typify the extremes of behavior that we find - S30 and S35 (cf. Table 1). Both simulations have an initial rotation velocity which results in a mean circularisation radius of (see §5.1). Simulation S30 has no imposed turbulence. Simulation S35 has turbulence characterised by , implying that gas turbulent motions are of the order of the velocity dispersion in the bulge.

4.1 A shell with no initial turbulence

Our no turbulence initial condition acts as a base comparison for the simulations with seeded turbulence. Figure 1 shows the gas column density and the velocity field in the angle-slice projection for the simulation S30 at . For projection along the -axis, the gas column density is calculated by


where the limits of the integration are given by , and . The angle is chosen to be for Figures 1 and 3. This projection method is convenient as it permits an unobscured view into the central regions where the black hole resides.

Figure 1 shows the gas flow at an early time, both an edge-on (the left hand side panel of the figure) as well as a top-down projection of the shell (the right panel). Due to a non-zero angular momentum, gas in the shell is quickly pushed aside from the axis of symmetry, opening a cylindrical cavity. Initially gas falls closer to the centre of the potential than its circularisation radius, a consequence of highly eccentric orbits. As gas moves inside , centrifugal force exceeds gravitational force. Interaction with neighboring gas streams results in “radial” shocks in the -plane.

A further set of shocks occurs due to gas initially moving supersonically in the vertical direction. As the gas from above the plane collides with the gas falling in from below the plane, the particles are shocked and accumulate at due to symmetry, forming a disc. These two sets of shocks mix gas with different angular momentum. We shall discuss this interesting effect in greater detail below.

The overall dynamics of the simulation are thus relatively simple: angular momentum conservation and symmetry dictate the formation of a geometrically thin disc in the plane of symmetry of the shell. Accretion of gas onto the SMBH would be expected to proceed in an accretion disc mode, if at all – realistic models show that star formation time scales are much shorter than viscous times, depleting the disc of gas before it can feed the SMBH (e.g., Goodman, 2003; Nayakshin & Cuadra, 2005; Nayakshin et al., 2007).

We note that our particular choice of velocity field, namely a constant azimuthal velocity, when combined with a spherical shell, could be viewed as a somewhat peculiar initial condition. As we have mentioned, the constant condition means that the polar regions of the shell are quickly evacuated, the gas here spiralling outward in cylindrical radius to encounter gas falling in with different angular momenta, and mixing with the latter. In fact what we are modelling here is the simplest case of a flow that undergoes an angular momentum re-distribution shock. We shall go into more detail on the consequences of this in Section 5, but for now we make the point that such a flow is likely to occur when gas is accreting in a stochastic fashion from large scales, perhaps as the result of a merger. The strong mixing of the gas with different angular momenta is exactly the situation that we wish to explore here, and so we have implemented what is essentially a symmetric case of this mixing process. We acknowledge that the spherical setup is idealised, but it is also the best starting point from which to embark on a laboratory of tests where the turbulence and rotation is varied.

Figure 1: Edge-on angle-slice projection (see §4.1) of the gas flow at time in the simulation S30 (left) and top-down projection (right). The gas falls in on eccentric orbits, giving rise to a radial shock that propagates outwards as the disc forms.

Figure 2: Projected gas column density in simulation S35 at time , before gas has had a chance to accrete on the SMBH. Note the formation of multiple thin and dense filaments due to convergent turbulent velocity flows.

4.2 A shell with initial turbulence

The flow of gas in this case is far more complicated than in the case just considered. We first show in Figure 2 the full slice projection of the simulation volume, e.g., we use here (c.f. equation 6). The most outstanding feature of the figure are the long dense filaments that form due to convergent turbulent velocity flows. The density contrast between the filaments and surrounding gas is over two orders of magnitude by this time.

Figure 3 shows both the y-axis and z-axis projections of this simulation at the same time as 1 did for the no initial turbulence case, . Clearly, some of the filaments seen at the earlier time in Figure 2 survive and actually penetrate into the innermost region There are now density contrasts of as much as three orders of magnitude in regions that were completely uniform in simulation S30. We shall see that this distinction is a crucial one for the dynamics of the gas and SMBH feeding.

Figure 3: Edge-on angle-slice projection of the gas flow at time in the simulation S35 (left) and top-down projection (right). The velocity field appears far more isotropic that in Figure 1, although an imprint of the imposed net rotation can still be seen.

4.3 Turbulence and accretion

Figure 4 shows the mass accreted by the black hole versus time in the simulations S30–S37, e.g., the same rotation velocity (and thus angular momentum) but different levels of initial turbulence. This demonstrates the main result of our study. The accretion rate onto the SMBH strongly correlates with the strength of the imposed turbulence. The accretion rate increases rapidly with increasing while , but then saturates at an approximately constant level for . The main qualitative explanation of the simulations is that turbulence decreases the degree to which gas with different angular momentum mixes, and creates gas streams with small angular momentum. In particular,

  • At low , an increase in the turbulent velocity leads to greater variations in the density fluctuations that are created by the turbulent velocity flows before gas circularises (cf. Figure 2). Greater density contrasts decrease the amount of angular momentum mixing, resulting in a disc rather than a narrow ring. The inner edge of the disc lies closer to the accretion radius of the simulations and hence feeds the SMBH more efficiently.

  • At , random initial velocity fields set some gas on orbits with a vanishingly small angular momentum compared with the mean in the shell. The turbulent “kick velocity” in this case almost cancels the mean rotation for these regions. Since these regions move against the mean flow, they are also those that get strongly compressed. Reaching high densities, they continue to move on nearly ballistic trajectories, impacting the innermost region on randomly oriented orbits. Accretion in this regime is “chaotic” (King & Pringle, 2007) rather than large-scale disc-dominated.

We shall spend the rest of the paper investigating these results in more depth, suggesting and testing analytical explanations for the observed trends.

Figure 4: Accretion rate versus time for simulations S30–S37. This plot exemplifies the main result of the paper: the accretion rate on the black hole strongly increases with increasing levels of turbulence when rotation is present (see §4.3).

5 Without turbulence: why “laminar” accretion is so low

5.1 Analytical estimates based on circularisation of gas

We shall now argue that the accretion rate of the simulation S30, the one with rotation velocity and no initial turbulence, is surprisingly low compared with a straightforward and seemingly natural theoretical estimate.

Let us start by estimating the fraction of gas that should be accreted by the SMBH in our simulations. Specific angular momentum of gas determines how close to the SMBH it circularises. The angular momentum of a circular orbit is given by . Using equation 1, one can obtain a general solution for the circularisation radius, , for a given value of specific angular momentum . We give it in two extremes. If , the point-mass Keplerian value applies:


In the opposite case, , we have


which simplifies to for , when the second and the last terms on the right hand side of the equation are small.

We now make the simplest possible assumption here by suggesting that gas settles into a disc and that the distribution of gas in the disc follows the distribution of gas over the circularisation radius initially, at time . Essentially, we assume that the and the components of the initial angular momentum cancel out due to symmetry whereas the component is conserved without any exchange with neighboring cells.

To estimate the fraction of gas that will end up inside the accretion radius , we first note it is equation 7 that should be used for the circularisation radius of gas with a given -projection of specific angular momentum, . The requirement singles out a cylinder with cross sectional radius of . The intersection of this cylinder with the shell has volume . For the shell the fraction of the volume that can be accreted is then given by the ratio of this volume to the total volume of the shell, :


For example, for , and this gives , and at we have (note that in code units).

The latter analytical estimate yields accreted mass , whereas Figure 4 shows that the actually measured value to late times is . The analytical estimate thus significantly over predicts the amount of accretion.

5.2 Shock mixing of gas: ring formation and end to accretion

Figure 5 explains why our simple analytical theory did not work. Here we plot the distribution of gas in simulation S30 over the -component of the angular momentum vector of particles, , at three different times. The initial distribution is spread over a broad range of values, with a small but non neglibible fraction of gas having , e.g. the angular momentum of the circular orbit at . In our analytical estimate we assumed that this gas accretes onto the black hole. However, this is not what happens. The distribution of angular momentum at the later time, , shows a strong radial mixing of gas with different angular momentum. The angular momentum distribution narrows due to shocks and eventually becomes a highly peaked Gaussian-like ring.

Note that the deficit of gas in the innermost region in the second curve in Figure 5 is caused not by the SMBH accretion but by the shocks described above. Low angular momentum gas shocks and mixes with high angular momentum material before it has a chance to travel into the SMBH capture region, . The mixing continues to late times and the peak in the angular momentum distribution actually moves outwards. This is to be expected as the gas that fell in earlier has a smaller angular momentum and is in the path of eccentric orbits of gas falling in from greater distances.

We believe this strong mixing of different angular momentum orbits is quite a general result of initially “laminar” flows. Such flows thus initially form rings rather than discs. While our simulations deliberately omit gas self-gravity and hence star formation (§2), previous theoretical work shows that the fate of the material is then decided by whether the viscous time of the ring is shorter than the star formation consumption time scale. Most authors agree that large-scale gas discs form stars more readily than they accrete (Goodman, 2003; Nayakshin et al., 2007). This would suggests that “laminar” shells with a finite angular momentum would form stars more readily than feed the SMBH.

Figure 5: The distribution of gas in the simulation S30 over the -component of the angular momentum at three different times: , , for solid, long dashed and dashed curves respectively. Note how the angular momentum distribution becomes narrower with time (through the action of shocks).

6 Accretion with seeded turbulence: why is it efficient?

6.1 Analytical expectations

6.1.1 Weak turbulence

First we consider the case of the turbulent velocity fields with . Figures 6 and 7 show the distribution of gas over the angular momentum for three different times in the simulations S31 and S32, respectively. As in Figure 5 for the no initial turbulence run S30, the initial angular momentum distribution is broader than those at the later times. The small angular momentum tail of the initial angular momentum distribution is more pronounced for S31 and especially S32 compared with S30, which helps to explain the higher accretion rates measured in these simulations.

The main effect, however, is the reduction in the shock mixing of gas. Consider a “small” angular momentum part of the distribution, to be definitive, at . This part of the distribution completely disappears at later times in simulation S30, being completely assimilated into the peak region. In contrast, in simulations S31 and S32 this part of the distribution is only reduced by a factor of 2-3 compared with the initial curves.

Apparently, the initial turbulent velocity flows create density contrasts that then propagate to smaller scales. The dynamics of a dense region are different from that of a low density region, leading to a larger spread in the angular momentum distribution at later times. This increases the amount of gas captured by the inner boundary condition.

Figure 6: The same as Figure 5 but for simulation S31, with an initial turbulent velocity . There is a greater fraction of small angular material in this simulation initially, and more of it gets retained in the “tail” to small , e.g., the inner disc, at late times.
Figure 7: The same as Figure 6 but for simulation S32. The features noted in Figure 6 are even more pronounced here.

6.1.2 Strong turbulence: ballistic accretion

Now we consider the case when . As we have seen, the turbulent accretion rate is far larger in our simulations than that for the “initially laminar” runs. We have also seen that differential and chaotic velocity flows form strong density enhancements consistent with the well known results from star formation studies (e.g., McKee & Ostriker, 2007). High density regions could in principle propagate through the mean density gas without much hydrodynamical drag, as long as their column densities are much higher than that of the surrounding gas. In this case we can approximate gas motion as ballistic and use the usual loss cone argument (e.g., page 406 in Shapiro & Teukolsky, 1983) for black hole accretion.

For a direct comparison to the simulation results, we consider a thin shell of gas that has angular momentum , and calculate the fraction of gaseous mass that has angular momentum small enough to be captured within . In doing so we assume that the energy and angular momentum of gas are both conserved as the orbits are approximately ballistic. The initial specific gas energy, , is small compared with the gravitational binding energy at , therefore we can set . An orbit that just reaches our inner boundary condition at has radial velocity . This yields the maximum specific angular momentum of gas that could still be accreted:


The case of a shell that rotates slowly compared with the local circular speed, i.e., , is most interesting, since the shell could not rotate faster than the circular speed or else the centrifugal forces exceed gravity; if it rotates at the circular speed then it is rotationally supported and would most likely form a disc rather than a spherical shell.

Assuming a monochromatic distribution

As a starting point for this estimate, we approximate the distribution of gas velocities in the shell with the turbulent velocity spectrum by a monochromatic distribution randomly and isotropically distributed in directions in the frame rotating with rotation velocity . We emphasise here that this is a simple assumption, but a sensible one, as on average one would expect the turbulent velocity field to have an isotropic character. We relate to the mean velocity amplitude in the initial turbulent velocity distribution.

In the simplest case , the fraction of solid angle that yields angular momentum smaller than is approximately . Therefore, the fraction of mass that could end up inside the accretion radius is


where is the mass of the shell.

Figure 8: Schematic geometry of an isotropic (in the rest frame of the rotating shell) wind, an analogy for the effect of turbulence on the loss-cone. The shell rotates with velocity and the central axis of the loss-cone is defined in terms of the angle .

However, this approach neglects the rotation of the shell. For a non-zero rotation velocity the loss-cone approach is valid only for , and equation 11 must be modified to take account of the orbital motion. Our derivation here is initially similar to Loeb (2004), but taken to second order when applying the small-angle approximation.

We start by considering the axis of zero angular momentum material for an isotropic, monochromatic distribution at a point in orbit of the central black hole. This is defined by


where is complementary to the angle between this axis and , as shown in Figure 8. We consider a small perturbation in angle about this axis, given by , which defines the opening angle of the loss-cone:


We expand this to second-order in the limit that , giving us


where . We therefore find for the loss-cone opening angle the expression:


For an isotropic distribution in the rest frame of the orbiting body, the fraction of the solid angle that will be accreted can be calculated via


In this case, at a given , our second-order approximation for the loss-cone solid angle is:


where . We can express this in the limit of two extremes, depending on the relative strength of the turbulence with respect to the rotation:

i) the case when :


ii) when :


Assuming a Maxwellian distribution

The treatment we have presented so far is not the whole picture, as we have assumed a monochromatic distribution for the turbulence. In reality this is not the case. The turbulent velocity distribution, over all three normally distributed components, has a Maxwellian profile, both in the initial condition and remaining approximately so in the high density regions that form. We find that the distribution for the ballistic approximation can therefore be fitted (to within error) by the expression:


To obtain the total mass within the loss-cone we therefore need to integrate this expression over the relevant solid angle:


were we should use equation 17 for the loss-cone solid angle, but replacing with . Unfortunately, performing the integration in this case is non-trivial and must be done numerically. Luckily, we can obtain a good analytical approximation by noting that the modified equation 17 is a maximum when , and falls off sharply with increasing . The majority of the accretion will therefore occur in the region of the distribution where and so we can use equation 19 in the integral instead, giving us:


where our limits are a perturbation either side of , the size of which we choose to be ***it should be noted that this choice is somewhat arbitrary but our derivation here is of course approximate. Performing the integral yields



For a thick shell we should generalise equation 23 by using the density weighted average of over the shell, . For an initially constant density profile used in our simulations and , we have , and hence


As before, we can analyse this equation in various limits:

i) when :


ii) when :


iii) when :


where we have dropped numerical factors of approximately unity.

6.2 Detailed analysis of results

Our ballistic approximation to the accretion rate yielded a result that depends on several parameters of the simulations, e.g., the “accretion” radius , the turbulent velocity , the rotation velocity and the dimensions of the shell. We shall now look at simulations covering a range in the space of these parameters to discuss specific trends and determine whether the simulations support our simple analytical theory.

6.2.1 No net angular momentum case

Figure 9 shows the accreted mass (left panel) and the accretion rate (right panel) as functions of time for simulations S00 – S07 (see Table 1). The rotation velocity for these runs is set to zero, and the turbulent velocity parameter is varied from 0 (solid curve) to 2 (blue short dashed).

The situation is rather simple for these runs. For turbulent velocity parameters much smaller than unity, the flow hardly deviates from that of a freely falling shell with a negligible pressure. Due to the constant density profile in our initial shell, this naturally leads to the accretion rate increasing as (see the dotted power law in the right panel of Figure 9) starting from the free fall time of the inner shell to that of the outer shell. As turbulence increases to , however, two effects are noteworthy. Firstly, the accretion episode starts earlier and also finishes later. This is naturally due to the spread in the initial radial velocities of the turbulent gas, with some regions starting to fall in with negative velocities, and hence arriving at the SMBH earlier, and others doing the opposite. There is also a significant reduction in the accretion rate at the highest values of turbulence. We believe this is due to the same effects already discussed in §4.3. Convergent turbulent flows create high density regions with a large range of angular momentum. These regions do not mix as readily as mean density regions and thus retain their angular momentum for longer. This reduces the amount of gas accreted. However we note that in this rather unrealistic set up (zero net angular momentum) gas accretion is still very efficient, with % of the shell being accreted by time .

Figure 9: Accreted mass vs. time (left) and accretion rate vs. time (right) for a velocity field that contains no net rotation (). The different lines represent different strengths of the mean turbulent velocity, : 0 (black solid), 0.1 (black dotted), 0.2 (black dashed), 0.3 (black dot-dashed), 0.5 (brown dot-dot-dash), 1.0 (red dashed), 1.5 (pink dotted), and 2.0 (blue dashed). Analytical line (thick red dashed) displays a slope.
Figure 10: Accreted mass vs. time (left) and accretion rate vs. time (right) for a velocity field that contains high net rotation (). Linestyles are as per Figure 9.

6.2.2 The highest rotation velocity case

In the other extreme, Figure 10 shows the accreted mass (left panel) and the accretion rate (right panel) for the simulations S50 – S57, where the rotation velocity is set to the maximum explored in the paper, e.g., . The range of the turbulent velocity parameter and the respective curve coding are the same as in Figure 9.

Comparing Figures 10 and 9, it is obvious that rotation significantly decreases the amount of gas accreted, and the accretion rate. The reduction is most severe in the case of , by about 5 orders of magnitude. The role of turbulence in simulations S50 – S57, as suggested in §4.3, is to increase the accretion rate. This occurs by enhancement of the angular momentum distribution into the low angular momentum end, and also by weakening the efficiency of shock mixing of the gas.

6.2.3 Rotation and turbulence parameter space

We now combine the results of all our numerical experiments from S00 to S57 in Table 1 by considering a single characteristic – the accreted gas mass at time in code units. All of these runs have the same initial geometrical arrangement of gas and a fixed accretion radius but differ in the strengths of the initial turbulent and rotation velocities.

The results are displayed in two ways. Figure 11 shows the accreted mass versus rotation velocity, . The simulations with the same level of mean turbulent velocity, , are connected by the different lines. In particular, the solid black curve shows the results for , and the dashed blue curve shows the results for .

It is clear that increased rotation velocity, and thus net angular momentum, always reduces the accreted mass, for any value of the turbulence parameter. This is in a general agreement with analytical expectations explored in §6.1.2.

Together with the simulation results, we plot the predicted maximum accreted mass (equation 26), which shows excellent agreement with the maximum delineated by the levels of turbulence at saturation. The analytical fit here is described specifically by


since the maximum fraction that can be accreted is unity and for we would expect the entire shell to have accreted by late times.

The increased level of turbulence acts to decrease the accretion rate in the case of very small rotation velocity but increase the accretion rate at higher rotation velocities (albeit up to saturation). This is not contradictory at all, however - the effect of turbulence in all cases is to spread the initial angular momentum distribution and prevent it from rapidly mixing into a single peak. In the case of no initial rotation this reduces the accretion rate by moving gas from zero to finite angular momentum orbits, whereas in the case of “large” initial net angular momentum the accretion rate is increased as some gas is moved to the low angular momentum orbits.

Figure 11: The accreted mass trend with rotation velocity, plotted for different levels of turbulence at (left) and the accreted mass trend with mean turbulent velocity, plotted for different levels of rotation at (right).

6.2.4 Accretion and the size of the shell

Equation 24 predicts that the fraction of accreted mass should be inversely proportional to the shell’s size, which we characterise by . To test this prediction we repeated the simulation S35 for a scaled-up version, , and a scaled-down, , shell (see runs S60 and S61 in Table 1).

To be definitive, we compare the accreted gas mass at the dynamical time at the outer radius of the shell, which corresponds to just after the initial peak in the accretion rate. This is approximately the time at which we would expect the ballistic mode to come to an end and a disc-dominated accretion mode to begin. The results are plotted in Figure 12 (left panel) as a function of .

Figure 12: Scaling of accreted mass with (left) and (right) for simulations with and . Symbols correspond to the dynamical time at the outer radius of the shell i.e., , , (left) and (right). Error bars are somewhat arbitrary, as it is not clear when the main ballistic mode of accretion ends - here they denote . Normalised analytical fits are equation 24 (red), the actual fit that would be expected for the simulation parameters, and equation 26 (green), the maximum accretion expected at saturation.

6.2.5 Trend with accretion radius

The accretion radius is a “nuisance” parameter of our study in the sense that it is introduced only to allow the simulations to run in a reasonable time. For physical reasons it would be desirable to make as small as possible. Therefore it is important to test whether the suggested analytical scaling of the results indeed holds. With this in mind, we repeated same simulation S35 for 5 different values of , spanning a range in from 0.001 to 0.005. The accreted mass from these simulations is shown in Figure 12 (right panel). The red dotted power-law gives the analytical prediction. Although similar, it is somewhat less steep than the simulation results.

6.2.6 Visualisation of the circumnuclear disc

Figure 14 compares the surface density projected along the -axis for simulations S30 (left) and S35 (right) in the inner region. Figure 14 shows the same for these two snapshots but projected along the -axis.

As was already clear from Figure 5, the gas forms a narrow ring in the simulation S30. This occurs due to efficient mixing of material with different angular momentum. This mixing is both radial (the ring is narrower) and vertical (the ring is vertically less extended).

In contrast, in the simulation S35, turbulence creates a broad distribution of angular momentum which does not mix so well. Mixing in the vertical direction is actually reasonably efficient, as the disc is thin and lies close to the plane (see the right panel in Figure 14). However the radial mixing is not so strong and the disc radial structure is very extended in S35 compared with S30.

Figure 13: Face-on projection of the disc that forms by with no turbulence (left) and mean (right). The initial rotation velocity for both was .
Figure 14: Slices of the CND (using ) formed with no turbulence and with mean , both for at t=0.418

6.2.7 Radial structure of the disc

Our analytical theory also makes a prediction for the surface density of the disc in the innermost part. Indeed, the mass captured by the inner boundary of radius scales as . This mass would likely sit in a rotationally supported disc of roughly size. The surface density in the disc should then behave as .

This analytical prediction (red power-law) is compared to the simulation results in Figure 15 which shows the surface density defined on radial shells for all the different levels of turbulence for . While the region beyound shows separate rings and highly variable surface density profiles for different simulations, the disc inside the region has a similar radial surface density profile for all simulations with . The agreement between simulations and theory is quite reasonable in that region except for runs with a small turbulent velocity . For these simulations the ballistic approximation is not appropriate and the gas does not penetrate into the region. At radii larger than 0.01 the surface density distribution is dominated by the interactions with the bulk of the initial shell as it continues to fall in. The ballistic approximation is of a limited use here because gas settling in the disc to late times (i) went through many collisions/interactions and hence the orbits are not ballistic; and (ii) has a relatively large angular momentum, so the small angle approximation we made in §5.1 is violated.

The observed dependence in the inner disc is another piece of simple but convincing evidence that the ballistic accretion approximation that we made in §6.1.2 is a reasonable one for the innermost region of the flow.

Figure 15: Azimuthally-averaged surface density versus radius for runs S50-S57 at time . Coloured curves correspond to simulations with different levels of mean turbulent velocity, as in Figure 9. The red line above the curves shows the theoretically expected (un-normalised) power law .

6.2.8 Vertical structure

Figure 14 showed edge-on projections of the disc for the simulation S30 (no turbulence) and (S35) initial turbulent velocity . It is clear that the vertical stucture of the disc is actually rather complex, especially beyond region. This is intuitively to be expected as the dynamical time is longest at larger radii, and therefore circularisation and decay of vertical velocity dispersion in the disc should take longer here. We can further quantify this by assuming that the vertical velocity dispersion is pumped at the rate at which the matter is brought in. Generation of momentum due to turbulence will therefore be described by . These motions should decay on the local dynamical time, . Assuming a steady-state allows us to estimate the vertical velocity dispersion as , where is the circular velocity at . Therefore, the expected disc aspect ratio in the outer parts (where the potential is isothermal) is .

We can calculate a dynamical, azimuthally-averaged disc scaleheight for the circumnuclear disc that forms by using the projection of the velocity dispersion along the ‘vertical’ direction for each annulus i.e., the mean angular momentum vector in each case. We therefore set,


where is the enclosed mass at , is the ‘vertical’ velocity dispersion at . Figure 16 plots ratios as functions of radius for simulations S50-S57. The solid line is for simulation S50, for which no initial turbulence is present. It is notable that the disc formed by the flow in this case is geometrically much thinner than in all the cases S51-S57. Among the turbulent cases, the general trend in is exactly as predicted above. The large scatter is explained by a strong time dependence: the disc is being stirred by strong waves that would eventually circularise the disc.

In the region inside , however, the observed dependence is flat while the analytical expectation is . We note that the vertical disc scale height is much larger than the SPH smoothing length, implying that the result inside is not due to under resolving of that region. We believe that the larger than expected vertical scale height is due to dynamical heating of the inner disc by the waves already noted at larger radii. Indeed, the inner disc is far less massive than the outer one, and any waves present in the outer region have to either dissipate or get reflected in this region.

Figure 16: Dependence of vertical disc scaleheight, , versus radius, for all levels of turbulence, with (S50-S57), at time t=0.418. Line styles of the curves are the same as in Figure 9. The red line above the curves shows the theoretically expected result, for the (un-normalised) slope .

7 Accretion inside the capture radius

While we cannot resolve the hydrodynamics of gas within the accretion radius, (1 pc for most of the simulations presented here), we track both the mass and the angular momentum of the particles that are captured. We can therefore calculate the mean angular momentum of the gas that settled inside . If the disc viscous time is comparable or longer than the duration of our simulations, which is likely (e.g., Nayakshin & Cuadra, 2005), then the direction of the mean angular momentum of the gas captured within region is related to that of the sub-parsec disc. Note that in principle the disc can be warped and then its structure is more complex (Lodato et al 2009).

Figure 17 shows the evolution of the direction of the mean specific angular momentum for the gas particles accreted inside in the simulation S30 (left panel, no initial turbulence) and S37 (right panel, ). The latter simulation is initialised with the highest levels of turbulence we considered in this paper, and therefore the effects of disc plane rotation are the strongest. The “north pole” in the diagram corresponds to a specific angular momentum vector oriented along the -axis. Clearly, there is negligible deviation from that direction for the “laminar” simulation S30, as expected. For the simulation S37, however, there is a complicated time evolution of the innermost disc direction. At one point the disc appears to be completely counter-aligned to the mean angular momentum of the shell.

Figure 17: Plots of the accreted angular momentum direction in the Aitoff projection. Both plots are with , with no turbulence (left) and (right). Crosses mark every 0.001 in code time units. Note that the title of “SMBH spin direction” refers to the spin orientation that the SMBH would have if it had accreted all of the gas within the capture radius; in reality this is of course not the case.

The accretion disc “plane wandering” at should also be seen in the larger scale flow of gas. Figure 18 shows the innermost region for the simulation S37 at time . We can clearly see here that the disc in the immediate vicinity of is also strongly displaced in terms of direction from the nominally expected one. One should also note the flows of gas along directions roughly orthogonal to the disc plane. These flows must have very different orientations of angular momentum, and it is these that cause the plane of the “sub-grid disc” to execute the non-trivial path seen in Figure 17.

Figure 18: Projected gas density and velocity field for simulation S37, the innermost region, for . Note how strongly the disc is misaligned with the -plane. At later times the misalignment decreases as more material accumulates into the disc.

8 Discussion

8.1 Summary of the results

We have studied the hydrodynamics of an initially uniform rotating gaseous shell, seeded with an initial turbulent velocity spectrum as it falls into the inner part of a spherically symmetric static potential. Particular attention was paid to the amount of gas that made it inside the “accretion radius” .

In the control runs without initial turbulence, gas settles into a nearly circular rotationally supported ring (see Figures 5, 14, 14). This was somewhat of a surprise. Intuitively we expected a wide () disc, as the initial distribution of angular momentum is broad (Figure 5). Furthermore, only a few hundred solar masses of gas were accreted, compared with the expected on the basis of a simple circularisation radius argument (§5.1). In hindsight, both of these effects are caused by the effective mixing of gas by the shocks. Low angular momentum gas gets mixed with high angular momentum gas before the former reaches the capture radius . The inner edge of the ring is a significant distance from the accretion radius. In the context of AGN accretion, such a ring would probably be strongly gravitationally unstable and form stars rather than feed the SMBH (e.g., Goodman, 2003; Nayakshin et al., 2007).

In contrast, runs in which the initial shell was seeded with turbulence show qualitatively and quantitatively different results. The most significant of these is that the accretion rates are higher by up to 3-4 orders of magnitude compared with the simulations that had no initial turbulence imposed. This is particularly striking since the net angular momentum of these simulations is the same.

It is obvious that turbulence broadens the angular momentum distribution, setting some gas on low angular momentum orbits. A less trivial conclusion from our numerical results is that turbulence also prevents the efficient gas mixing pointed out above. Physically, high density regions are created by supersonic turbulent convergent flows. These regions can then travel nearly ballistically through the average density gas. Some of this gas is able to reach the accretion radius, increasing the SMBH accretion rate compared with runs without turbulence. This parameter space trend lasts until the initial mean turbulent velocity becomes greater than the net rotation velocity, at which point the accretion trend saturates and begins to slowly decrease as turbulence is increased further.

The development of high density regions in supersonic turbulence is by now a classical result from a number of papers (for reviews see e.g., Elmegreen & Scalo, 2004; McKee & Ostriker, 2007). There is also a body of work that highlights the complicated dynamics of shocks in the multi-phase interstellar medium. For example, Bonnell et al. (2006); Dobbs (2008) suggest that shocks in a multi-phase medium generate (or rather preserve) differential motions, allowing the high density regions to slip through the shocks. These effects are similar in nature to the ballistic behavior of high density clumps found in our simulations.

Schartmann et al. (2009) simulated young nuclear star clusters injecting the energy into the surrounding gas via supernovae and winds. As in our study it was found that high density streams form. Those with small angular momentum feed the small scale “torus”. These results, and earlier results by Cuadra et al. (2006) for stellar winds in the central parsec of the Milky Way also support the viability of the “ballistic” AGN feeding mode.

We have attempted to build a simple analytical theory that would predict the accretion rate for our capture inner boundary condition. We first calculated the accretion rate in the ballistic approximation assuming that gas in the shell receives monochromatic velocity kicks in random directions in the frame rotating with velocity . This “monochromatic” theory predicts that gas capture rate (fraction) is the highest at . We then took into account in an approximate way the corrections associated with the fact that the real turbulent velocity distribution function is not monochromatic. The fraction of gas accreted within calculated in this way turned out to describe our numerical results well. Furthermore, the theory predicted that the inner disc surface density should have a power-law scaling: , which is indeed born out by the simulations.

The approximate agreement of theory and simulations in this paper implies that we can use our results for systems of astrophysical interest by an appropriate rescaling of parameters, which we attempt to do in the next section.

8.2 Can quasars and AGN be fed by ballistic accretion?

The most serious challenge to theories of quasar and AGN feeding is gravitational instability of these discs and the ensuing fragmentation of discs into stars (e.g., Paczyński, 1978; Kolykhalov & Sunyaev, 1980; Collin & Zahn, 1999; Goodman, 2003; Nayakshin et al., 2007). Observations (e.g., Paumard et al., 2006; Bartko et al., 2009) and numerical simulations (Bonnell & Rice, 2008; Hobbs & Nayakshin, 2009) suggests that this process operated in the central 0.5 pc of the Milky Way, creating young stars and probably leaving the resident SMBH – Sgr A* – without gaseous fuel to shine now (e.g., Levin & Beloborodov, 2003; Nayakshin & Cuadra, 2005).

One obvious way in which the self-gravity “catastrophe” of AGN discs can be avoided is to deposit the gas at small enough radii (e.g., King & Pringle, 2007). Namely, one can show that gaseous discs are hot enough to avoid gravitational fragmentation within the self-gravity radius, , which is of the order of pc depending on the SMBH mass and some of the disc parameters (see Goodman, 2003). Inclusion of feedback from star formation allows the disc to survive and yet maintain an interesting accretion rate on the SMBH out to a few times larger radii (Pedes and Nayakshin, in preparation).

We shall now try to estimate the plausibility of such a small radius gas deposition in the context of the “ballistic” gas deposition model we have developed. To this end we shall use the approximate scalings for the fraction of gas accreted from a rotating turbulent shell derived in this paper in §6.1.2.

The size of the shell, , enters this estimate directly. Existence of the relation implies (e.g., Häring & Rix, 2004) that the SMBH must be fed by gas located roughly the size of the bulge away from the SMBH (King, 2003, 2005), or else it is difficult to see why the quoted feedback link must exist (Nayakshin & Power, 2009). We estimate that the effective radius of the bulge, , is about kpc for velocity dispersion km/sec (see Nayakshin et al., 2009). If the rotation velocity at is a small fraction of , say km s, then the fraction of the gas that can be accreted from a turbulent shell is


where . Here we have used equation 26, assuming that the mean turbulent velocity is of the order of the rotation velocity, which is our maximum in the accretion trend (although it should be noted that going to stronger turbulence does not decrease the accreted mass by much - cf. §6.2). The relation here gives the stellar mass in the bulge as , and if we assume that a half of the bulge gas was used up in star formation (McGaugh et al., 2009), we find that the action of supernova feedback in the context of our model yields an accreted mass of . This accretion would take place over roughly the dynamical time at the outer edge of the bulge, giving us an accretion rate a few yr. In this simple estimate, then, the SMBH is indeed able to grow as massive as the observations require, assuming that our ballistic accretion mode can be sustained for the required time. Future work should test these ideas more directly, with “turbulence” driven directly by feedback from star formation in the bulge, where the bulge gas distribution has been derived from initial conditions in the gas on larger scales, e.g., the virial radius of the host dark matter halo.

8.3 A positive feedback link of star formation to AGN activity?

Observations show that starburst and AGN activity are correlated in a number of ways (e.g., Cid Fernandes et al., 2001; González Delgado et al., 2001; Farrah et al., 2003). One might think that this is entirely natural, as both phenomena require a source of gas to be present. However, star formation on kpc-scales does not have to be casually connected with SMBH activity. As an example of this, consider a case when the angular momentum of the gas is large, so that a kpc-scale disc is formed. If star formation in this disc simply consumes the gas then there is nothing left to feed the SMBH. This is what may have happened on smaller scales in the central parsec of our Galaxy million years ago (Nayakshin & Cuadra, 2005; Nayakshin et al., 2007). Naturally, if such a process can rob the SMBH of fuel at parsec scales then the situation becomes even worse at kiloparsec scales.

The opposite situation is also possible because the SMBH mass is a small fraction of the bulge mass. It is not implausible to have a sub-pc scale, non self-gravitating accretion disc that could feed a moderately bright AGN for a long ( years, say) time. The amount of gas involved in this could be , i.e., miniscule by galaxy standards. Why such activity should always be connected with a powerful starburst is not clear.

However, the inclusion of star formation feedback could alter this picture. For example, within a kpc-scale disc, feedback from stellar winds and supernovae gives velocity kicks to the surrounding gas. For a standard IMF, the total momentum in supernova shells and winds from massive stars, integrated over their life time, is about , where , and is the total mass of the stellar population (Leitherer et al., 1992; Thompson et al., 2005). If this momentum is absorbed by gas with mass then the average kick velocity for the gas is , or km s for . There is thus certainly enough momentum input to make the gas highly turbulent. The cancellation of oppositely directed momentum results in a reduction in this estimate. However, we have discussed here only the momentum input; including supernova kinetic energy could act to increase the above value (e.g., Dekel & Silk, 1986).

As we have shown in the paper, the effect of such feedback is likely to result in some small angular momentum clouds or streams that may feed the SMBH. Our model therefore predicts a positive correlation of starburst activity with SMBH activity, with the latter being somewhat offset in time more often than not (it is also possible to activate an AGN earlier if there is some low angular momentum gas in the initial shell).

Further detailed physics-based modeling of AGN feeding, feedback, and star formation, as well as comparisons to observations should shed light on the mode of quasar feeding, and also on the mode of bulge formation.

8.4 Cosmological cold stream SMBH feeding?

We have employed turbulent initial conditions as a simple and mathematically convenient way to introduce strong disordered differential motions in the gas for controlled numerical experiments in AGN feeding. We expect realistic galaxies to be always in the “turbulent” regime in the above sense because the gas contracting onto a dark matter halo is unlikely to be of uniform density or possess a uniform velocity field. In fact, large-scale numerical simulations of galaxy formation show cold gas streams that penetrate the shock-heated media of massive dark matter haloes (Dekel et al., 2009). This is physically similar to the ballistic motion of gas we observed in our simulations. While it is very unlikely that these very large-scale streams (hundreds of kpc) would be sufficiently well aimed to actually feed the SMBH, collisions of several of such streams in the centre of a galaxy would pump very strong differential motions and lead to some gas being set on small angular momentum orbits. We therefore speculate that cold large-scale gas streams may promote AGN feeding by providing energy input to turbulence/differential motions on galactic scales.

9 Conclusions

The classical challenge to AGN feeding is the large angular momentum of gas expected to result in circularised discs of kpc scales (e.g., Shlosman, Frank & Begelman 1989). Transporting gas from these scales is extremely difficult due to long viscous times and the rapid consumption of gas in star formation (e.g., Goodman, 2003; Nayakshin et al., 2007). Here we pointed out that turbulence in the bulge, e.g., driven by supernova explosions, may be a way to overcome the AGN feeding debacle by setting some gas on ballistic orbits.

To test these ideas, we have performed simulations of hydrodynamical gas flow with angular momentum, initially distributed in a spherical shell, in the static spherical potential of a bulge with a central SMBH. The duration of our simulations is a few dynamical times of the shell, allowing us to study the formation of an accretion disc and the capture of gas within a small inner boundary region. For all values of the rotation velocity, the shell was not in equilibrium with the background potential. This was done deliberately, as the goal of this project was to study the formation of a disc from an infalling gas distribution from first principles.

We found that angular momentum mixing in the runs without initial turbulence is very effictive, resulting in an initial ring rather than an extended disc. Further evolution of the ring depends on whether it can transfer angular momentum quickly enough and avoid completely collapsing into stars. If star formation consumes the ring quickly then such rings are a dead end as far as SMBH feeding is concerned.

Turbulent motions in the shell overcome the angular momentum mixing problem by creating high density regions that can travel nearly ballistically, retaining their initial angular momentum. Provided there is enough gas with small angular momentum the SMBH can be fed much more efficiently than in the “laminar” regime. The main conclusion here is that the star formation in the bulge may not simply deplete the gas (hence depriving the SMBH of its fuel) but may actually promote SMBH accretion by creating turbulent motions.

We note that turbulence in this particular project is just a mathematically convenient formalism to introduce strong disordered differential motions in the gas. We expect that realistic galaxies are always in the turbulent regime because the state of the gas on large scales is unlikely to be that of uniform density and angular momentum.

Further work with self-consistently driven turbulence from star formation and realistic initial (large-scale) boundary conditions is needed to establish the relevance of the “ballistic” accretion to feeding the SMBH, although our research suggests that it is a promising model.

10 Acknowledgments

Helpful discussions with Rashid Sunyaev, Eugene Churazov, Simon White, Walter Dehnen, Justin Read and Peter Cossins are gratefully acknowledged. Theoretical Astrophysics group at the University of Leicester is funded by a research grant from STFC. AH acknowledges an STFC studentship.


  • Balsara (1995) Balsara D. S., 1995, Journal of Computational Physics, 121, 357
  • Barnes & Hernquist (1991) Barnes J. E., Hernquist L. E., 1991, ApJL, 370, L65
  • Bartko et al. (2009) Bartko H., Martins F., Fritz T. K., et al., 2009, ApJ, 697, 1741
  • Bate (2009) Bate M. R., 2009, MNRAS, 397, 232
  • Bate et al. (1995) Bate M. R., Bonnell I. A., Price N. M., 1995, MNRAS, 277, 362
  • Bondi (1952) Bondi H., 1952, MNRAS, 112, 195
  • Bondi & Hoyle (1944) Bondi H., Hoyle F., 1944, MNRAS, 104, 273
  • Bonnell et al. (2006) Bonnell I. A., Dobbs C. L., Robitaille T. P., Pringle J. E., 2006, MNRAS, 365, 37
  • Bonnell & Rice (2008) Bonnell I. A., Rice W. K. M., 2008, Science, 321, 1060
  • Booth & Schaye (2009) Booth C. M., Schaye J., 2009, MNRAS, 398, 53
  • Chen et al. (2009) Chen Y., Wang J., Yan C., Hu C., Zhang S., 2009, ApJL, 695, L130
  • Cid Fernandes et al. (2001) Cid Fernandes R., Heckman T., Schmitt H., González Delgado R. M., Storchi-Bergmann T., 2001, ApJ, 558, 81
  • Collin & Zahn (1999) Collin S., Zahn J., 1999, A&A, 344, 433
  • Collin & Zahn (2008) Collin S., Zahn J., 2008, A&A, 477, 419
  • Cuadra et al. (2006) Cuadra J., Nayakshin S., Springel V., Di Matteo T., 2006, MNRAS, 366, 358
  • DeBuhr et al. (2009) DeBuhr J., Quataert E., Ma C., Hopkins P., 2009, ArXiv e-prints
  • Dekel et al. (2009) Dekel A., Birnboim Y., Engel G., et al., 2009, Nature, 457, 451
  • Dekel & Silk (1986) Dekel A., Silk J., 1986, ApJ, 303, 39
  • Dobbs (2008) Dobbs C. L., 2008, MNRAS, 391, 844
  • Elmegreen & Scalo (2004) Elmegreen B. G., Scalo J., 2004, ARA&A, 42, 211
  • Farrah et al. (2003) Farrah D., Afonso J., Efstathiou A., Rowan-Robinson M., Fox M., Clements D., 2003, MNRAS, 343, 585
  • Ferrarese & Merritt (2000) Ferrarese L., Merritt D., 2000, ApJL, 539, L9
  • Gebhardt et al. (2000) Gebhardt K., Bender R., Bower G., et al., 2000, ApJL, 539, L13
  • González Delgado et al. (2001) González Delgado R. M., Heckman T., Leitherer C., 2001, ApJ, 546, 845
  • Goodman (2003) Goodman J., 2003, MNRAS, 339, 937
  • Häring & Rix (2004) Häring N., Rix H.-W., 2004, ApJL, 604, L89
  • Hobbs & Nayakshin (2009) Hobbs A., Nayakshin S., 2009, MNRAS, 394, 191
  • Hoyle & Lyttleton (1939) Hoyle F., Lyttleton R. A., 1939, in Proceedings of the Cambridge Philosophical Society, vol. 34 of Proceedings of the Cambridge Philosophical Society, 405–+
  • Kawakatu & Wada (2008) Kawakatu N., Wada K., 2008, ApJ, 681, 73
  • King (2003) King A., 2003, ApJL, 596, L27
  • King (2005) King A., 2005, ApJL, 635, L121
  • King & Pringle (2007) King A. R., Pringle J. E., 2007, MNRAS, 377, L25
  • Kolykhalov & Sunyaev (1980) Kolykhalov P. I., Sunyaev R. A., 1980, Soviet Astron. Lett., 6, 357
  • Krumholz & Bonnell (2007) Krumholz M. R., Bonnell I. A., 2007, ArXiv e-prints, 712
  • Kurk et al. (2007) Kurk J. D., Walter F., Fan X., et al., 2007, ApJ, 669, 32
  • Leitherer et al. (1992) Leitherer C., Robert C., Drissen L., 1992, ApJ, 401, 596
  • Levin & Beloborodov (2003) Levin Y., Beloborodov A. M., 2003, ApJ, 590, L33
  • Li et al. (2007) Li Y., Hernquist L., Robertson B., et al., 2007, ApJ, 665, 187
  • Loeb (2004) Loeb A., 2004, MNRAS, 350, 725
  • Mac Low & Klessen (2004) Mac Low M., Klessen R. S., 2004, Reviews of Modern Physics, 76, 125
  • Magorrian et al. (1998) Magorrian J., Tremaine S., Richstone D., et al., 1998, AJ, 115, 2285
  • McGaugh et al. (2009) McGaugh S. S., Schombert J. M., de Blok W. J. G., Zagursky M. J., 2009, ArXiv e-prints
  • McKee & Ostriker (2007) McKee C. F., Ostriker E. C., 2007, ARA&A, 45, 565
  • McKee & Ostriker (1977) McKee C. F., Ostriker J. P., 1977, ApJ, 218, 148
  • Monaghan & Gingold (1983) Monaghan J. J., Gingold R. A., 1983, Journal of Computational Physics, 52, 374
  • Nayakshin & Cuadra (2005) Nayakshin S., Cuadra J., 2005, A&A, 437, 437
  • Nayakshin et al. (2007) Nayakshin S., Cuadra J., Springel V., 2007, MNRAS, 379, 21
  • Nayakshin & Power (2009) Nayakshin S., Power C., 2009, ArXiv e-prints
  • Nayakshin et al. (2009) Nayakshin S., Wilkinson M. I., King A., 2009, MNRAS, 398, L54
  • Paczyński (1978) Paczyński B., 1978, Acta Astron., 28, 91
  • Paumard et al. (2006) Paumard T., Genzel R., Martins F., et al., 2006, ApJ, 643, 1011
  • Read et al. (2010) Read J. I., Hayfield T., Agertz O., 2010, MNRAS, 405, 1513
  • Schartmann et al. (2009) Schartmann M., Meisenheimer K., Klahr H., Camenzind M., Wolf S., Henning T., 2009, MNRAS, 393, 759
  • Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  • Shapiro & Teukolsky (1983) Shapiro S. L., Teukolsky S. A., 1983, Black holes, white dwarfs, and neutron stars: The physics of compact objects
  • Sijacki et al. (2009) Sijacki D., Springel V., Haehnelt M. G., 2009, MNRAS, 400, 100
  • Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  • Springel et al. (2005) Springel V., Di Matteo T., Hernquist L., 2005, MNRAS, 361, 776
  • Springel & Hernquist (2002) Springel V., Hernquist L., 2002, MNRAS, 333, 649
  • Thompson et al. (2005) Thompson T. A., Quataert E., Murray N., 2005, ApJ, 630, 167
  • Wang et al. (2009) Wang J., Yan C., Li Y., et al., 2009, ApJL, 701, L7
S00 0.03 0.1 0.0 0.0 0.0 0.001
S01 0.03 0.1 0.0 0.1 0.001 0.0 0.001
S02 0.03 0.1 0.0 0.2 0.004 0.0 0.001
S03 0.03 0.1 0.0 0.3 0.009 0.0 0.001
S04 0.03 0.1 0.0 0.5 0.026 0.0 0.001
S05 0.03 0.1 0.0 1.0 0.103 0.0 0.001
S06 0.03 0.1 0.0 1.5 0.231 0.0 0.001
S07 0.03 0.1 0.0 2.0 0.410 0.0 0.001
S10 0.03 0.1 0.1 0.0 0.004 0.001
S11 0.03 0.1 0.1 0.1 0.001 0.004 0.001
S12 0.03 0.1 0.1 0.2 0.004 0.004 0.001
S13 0.03 0.1 0.1 0.3 0.009 0.004 0.001
S14 0.03 0.1 0.1 0.5 0.026 0.004 0.001
S15 0.03 0.1 0.1 1.0 0.103 0.004 0.001
S16 0.03 0.1 0.1 1.5 0.231 0.004 0.001
S17 0.03 0.1 0.1 2.0 0.410 0.004 0.001
S20 0.03 0.1 0.2 0.0 0.011 0.001
S21 0.03 0.1 0.2 0.1 0.001 0.011 0.001
S22 0.03 0.1 0.2 0.2 0.004 0.011 0.001
S23 0.03 0.1 0.2 0.3 0.009 0.011 0.001
S24 0.03 0.1 0.2 0.5 0.026 0.011 0.001
S25 0.03 0.1 0.2 1.0 0.103 0.011 0.001
S26 0.03 0.1 0.2 1.5 0.231 0.011 0.001
S27 0.03 0.1 0.2 2.0 0.410 0.011 0.001
S30 0.03 0.1 0.3 0.0 0.019 0.001
S31 0.03 0.1 0.3 0.1 0.001 0.019 0.001
S32 0.03 0.1 0.3 0.2 0.004 0.019 0.001
S33 0.03 0.1 0.3 0.3 0.009 0.019 0.001
S34 0.03 0.1 0.3 0.5 0.026 0.019 0.001
S35 0.03 0.1 0.3 1.0 0.103 0.019 0.001
S36 0.03 0.1 0.3 1.5 0.231 0.019 0.001
S37 0.03 0.1 0.3 2.0 0.410 0.019 0.001
S40 0.03 0.1 0.4 0.0 0.026 0.001
S41 0.03 0.1 0.4 0.1 0.001 0.026 0.001
S42 0.03 0.1 0.4 0.2 0.004 0.026 0.001
S43 0.03 0.1 0.4 0.3 0.009 0.026 0.001
S44 0.03 0.1 0.4 0.5 0.026 0.026 0.001
S45 0.03 0.1 0.4 1.0 0.103 0.026 0.001
S46 0.03 0.1 0.4 1.5 0.231 0.026 0.001
S47 0.03 0.1 0.4 2.0 0.410 0.026 0.001
S50 0.03 0.1 0.5 0.0 0.034 0.001
S51 0.03 0.1 0.5 0.1 0.001 0.034 0.001
S52 0.03 0.1 0.5 0.2 0.004 0.034 0.001
S53 0.03 0.1 0.5 0.3 0.009 0.034 0.001
S54 0.03 0.1 0.5 0.5 0.026 0.034 0.001
S55 0.03 0.1 0.5 1.0 0.103 0.034 0.001
S56 0.03 0.1 0.5 1.5 0.231 0.034 0.001
S57 0.03 0.1 0.5 2.0 0.410 0.034 0.001
S60 0.015 0.05 0.3 1.0 0.091 0.008 0.001
S61 0.06 0.2 0.3 1.0 0.119 0.041 0.001
S70 0.03 0.1 0.3 1.0 0.103 0.019 0.002
S71 0.03 0.1 0.3 1.0 0.103 0.019 0.003
S72 0.03 0.1 0.3 1.0 0.103 0.019 0.004
S73 0.03 0.1 0.3 1.0 0.103 0.019 0.005

here refers to the initial value of circularisation radius for each of the particles, calculated in the region of the potential that they start in, based solely on the rotation velocity.

value of mass accreted taken at in code units.

Table 1: Initial conditions for each simulation. The ratio with respect to the gravitational potential energy gives an indication of how well supported the shell is against the external potential through virial motions: for a virialised shell.

Appendix A An alternative velocity distribution for the shell

On the referee’s request, we have re-simulated a fiducial parameter set with an alternative velocity field. In this case we set , the specific angular momentum, to be constant, rather than the rotation velocity. The magnitude of the velocity at a given radius is then


where is a constant. This is a more conservative initial condition in the sense that this distribution would be expected to form an infinitely thin ring, since the distribution is a delta function in angular momentum space. The broadening of this distribution due to turbulence would therefore be more difficult, and any enhancement to the accretion rate that this broadening provides is therefore a significant result. We note that the other “natural” choice of velocity field for a spherical shell, namely solid body rotation, is only reasonable for a system in hydrostatic equilibrium, which our shell is not.

Figure 19: Accretion rate versus time for a set of runs with a constant angular momentum velocity field, with in code units. The other parameters for these runs were held at the fiducial values of S30, and color/linestyles is as standard for the other plots in the paper. The main result of the paper still holds: the accretion rate on the black hole strongly increases with increasing levels of turbulence when rotation is present (see §4.3).

As can be seen from Figure 19, our main conclusion is unchanged for this alternative velocity distribution. The accretion rate onto the inner parsec is increased by orders of magnitude from zero/low turbulence to high turbulence. We note the the constant angular momentum within a spherical shell is, too, a strongly idealised setup, with similar behaviour at the polar regions to that of the constant condition. However, we are again attempting to model a simple case of gas mixing due to different angular momentum at similar radii, and unfortunately there is no “realistic” single-parameter velocity field that can describe a spherical distribution with these properties.

Appendix B Convergence tests on numerical viscosity and resolution

A great deal of the accretion that we find in our simulations is the result of turbulent motions creating dense, shocked gas. Due to the artificial viscosity present in the SPH code it is possible that some of this accretion is numerical rather than physical. We have therefore performed convergence tests both on the value of the artificial viscosity parameter, (cf. Section 2) and the numerical resolution of our fiducial simulations S30 (no turbulence, ) and S35 (, ), comparing the accretion rate between the two. The results are shown in Figures 20 and 21 for the artificial viscosity and particle resolution respectively.

Figure 20: Accretion rate versus time for convergence tests on the value of the artificial viscosity parameter, . Solid linestyles refer to simulations without turbulence while dashed linestyles refer to simulations with . The rotation velocity for all runs here was . Colors are: (green), (blue), (red), (magenta) and (black; our fiducial value for the simulations in the paper). The simulations with no seeded turbulence show good convergence with a trend towards an increased accretion rate with an increasing , while the simulations with show good convergence but with no clear trend.

For Figure 20 there is clear convergence both in the runs with no turbulence and those with . This suggests that the role of the artificial viscosity is negligible in determining the accretion rate through the SMBH accretion radius, at least until the main ballistic mode of feeding ends (it is likely that once the accretion becomes dominated by an accretion disc mode, i.e., to late times, the numerical viscosity will play a greater role).

Figure 21: Accretion rate versus time for convergence tests on the numerical resolution. Once again solid linestyles are for runs with while dashed refer to runs with , and all simulations shown used . Colors refer to (green), (blue), (magenta), (red), (brown) and (black; our fiducial resolution). In the runs with no seeded turbulence there is a strong trend of decreasing accretion with increasing resolution, which shows signs of starting to converge by the time the highest resolution is reached. In the runs with there is convergence to within a factor of a few but without a clear trend in the accretion rate.

Figure 21 shows the convergence results on the number of SPH particles. In the runs with no turbulence it is not clear whether the accretion rate has converged with increasing resolution, as we were unable to run simulations with a yet higher number of particles due to computational limitations. However, the plot would seem to suggest that by the highest resolution there is some convergence to late times. It may of course simply be the case that for this particular value of all of the accretion is numerical for the runs with zero turbulence, and the accretion rate will therefore vanish in the limit of increasing resolution. We note that this leaves our conclusions unaffected, and indeed strengthens the case for turbulence enhancing the accretion rate onto the SMBH. In the runs with there is some reasonable convergence, albeit about a relatively large spread with no clear trend. A possible reason for this is the fact that the SPH method tends to suffer from an inability to converge with increasing resolution in regions where particles are distributed irregularly on the kernel scale (such as in strong shocks), as discussed by Read et al. (2010). However, this spread is not large enough to cast any doubt on our main conclusion, and furthermore is only a factor of a few for the 3 higher resolution runs (red, brown, black lines).

We can therefore be confident that the effect seen in our simulations of an enhanced accretion rate with increasing turbulence (up to saturation) through a ballistic feeding mode is physical rather than numerical.

Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description