Accretion Disc Particle Accretion in Major Merger Simulations
A recent approach to simulating localized feedback from active galactic nuclei (AGN) by Power et al. (2011) uses an accretion disc particle to represent both the black hole and its accretion disc. We have extrapolated and adapted this approach to simulations of Milky Way-sized galaxy mergers containing black holes and explored the impact of the various parameters in this model as well as its resolution dependence. The two key parameters in the model are an effective accretion radius, which determines the radius within which gas particles are added to the accretion disc, and a viscous time-scale which determines how long it takes for material in the accretion disc to accrete on to the black hole itself. We find that there is a limited range of permitted accretion radii and viscous time-scales, with unphysical results produced outside this range. For permitted model parameters, the nuclear regions of simulations with the same resolution follow similar evolutionary paths, producing final black hole masses that are consistent within a factor of two. When comparing the resolution dependence of the model, there is a trend towards higher resolution producing slightly lower mass black holes, but values for the two resolutions studied again agree within a factor of two. We also compare these results to two other AGN feedback algorithms found in the literature. While the evolution of the systems vary, most notably the intermediate total black hole mass, the final black hole masses differ by less than a factor of five amongst all of our models, and the remnants exhibit similar structural parameters. The implication of this accretion model is that, unlike most accretion algorithms, a decoupling of the accretion rate on to the black hole and the local gas properties is permitted and obtained; this allows for black hole growth even after feedback has prevented additional accretion events on to the disc.
keywords:black hole physics – galaxies: interactions – galaxies: active – methods: numerical
Observational evidence suggests that supermassive black holes exist at the centre of all galaxies with stellar spheroids (e.g. Kormendy & Richstone, 1995; Ferrarese & Merritt, 2000). Further evidence, such as the – relationship (e.g. Silk & Rees, 1998; Ferrarese & Merritt, 2000; Gebhardt et al., 2000; Tremaine et al., 2002; King, 2003; Gültekin et al., 2009) and the – relationship (Magorrian et al., 1998; McLure & Dunlop, 2002; Marconi & Hunt, 2003), suggests that the black hole and spheroid evolution are related. Moreover, galaxies are relatively smaller today than might be naively predicted from the hierarchical model of galaxy formation, indicating that some mechanism has limited the growth of these galaxies. One favoured explanation is that during mergers gas from the merger fuels both star formation and AGN activity (e.g. Sanders et al., 1988; Scannapieco et al., 2005). The latter is likely a self-regulated process: outflows from the black hole following a strong accretion event interact with the surrounding gas, inhibiting further accretion events, and hence limiting black hole growth (e.g. Silk & Rees, 1998; Fabian, 1999; Scannapieco & Oh, 2004; King, 2005; King, 2010). Ultimately, the feedback from the increased AGN activity blows away some or possibly all the gas, truncating the star formation and leading to an elliptical galaxy (e.g. Springel et al., 2005a).
Using various models, AGN feedback has been implemented in many numerical simulations (e.g. Springel et al. 2005b; Thacker et al. 2006; Okamoto et al. 2008; Booth & Schaye 2009; Kurosawa et al. 2009; Debuhr et al. 2011). The goals of this research have been varied, from reproducing the observed relationships between the black hole and spheroid, to other factors impacted by AGN activity, such as galaxy cluster properties (e.g. Puchwein et al., 2008) or even the impact on CMB foregrounds (e.g. Scannapieco et al., 2008). In all these simulations, the accretion rate of gas on to the black hole is dependent on extrapolating the macroscopic gas properties to the microscopic scale around the black hole. This is true for viscous accretion (e.g. Debuhr et al., 2011; Hopkins & Quataert, 2011), drag accretion (Okamoto et al., 2008), and the numerical form of Bondi accretion rate (Bondi, 1952),
where is a numerical parameter used to account for the limited dynamic range in the simulations, and are the gas density and sound speed around the black hole, is the relative velocity between the gas and the black hole, and is the mass of the black hole. This accretion rate has been used in merger (e.g. Springel et al. 2005b; Di Matteo et al. 2005) and in cosmological (e.g. Sijacki et al., 2007; Booth & Schaye, 2009) simulations.
Most algorithms also limit the mass accretion to be no greater than the Eddington accretion rate,
where is the proton mass, is the Thomson cross section, and is the radiative efficiency (i.e. the mass-to-energy conversion efficiency); we set (Shakura & Sunyaev, 1973). This rate differs from the previously considered accretion rates in that it depends only on the black hole mass, and not the local gas properties.
One of the assumptions frequently used in these models is that any accreted gas is immediately transferred to the black hole. This is clearly an unrealistic assumption as material must first shed its angular momentum. Indeed some models argue that understanding this process is one of the key factors in modelling AGN feedback; for example, the viscous accretion rate in Debuhr et al. (2011) accounts for angular momentum of the gas. However, the gas in this model is nonetheless instantly accreted on to the black hole. In reality the gas would be expected to settle on to a circular orbit of radius , which is set by the angular momentum of the gas and whatever processes cause it to shed angular momentum while collapsing (cf. Hobbs et al., 2011), resulting in an accretion disc around the black hole. The gas with the lowest angular momentum would then travel through the disc and eventually be accreted on to the black hole (e.g. King, 2010).
Motivated by these ideas, Power et al. (2011) (hereafter PNK11) have implemented a two-stage accretion algorithm using an accretion disc particle (ADP), which models both the black hole and accretion disc. In the first stage, nearby gas, the bulk of which will have shed large amounts of angular momentum to reach this radius, is accreted on to the accretion disc; this accretion rate is dependent only on the relative positions of the black hole and the gas particle. In the second stage, the gas is accreted from the accretion disc on to the black hole. This method incorporates the delay between the time gas is accreted on to the disc and when it is finally accreted on to the black hole, and at the same time decouples the black hole’s accretion rate from the instantaneous gas properties around the black hole.
While the ADP model implemented in Power et al. (2011) relies upon resolving scales far smaller than those that can be resolved in merger simulations, the two stage approach, and perhaps more specifically the incorporation of a delay period before accreting on to the black hole, is an issue worthy of investigation in merger simulations which can have time-steps in the few thousand year range. But this consideration highlights the fact that resolution dependence must be considered. For example, at very low mass resolution and the associated low time resolution, concerns about accretion delays are likely less significant as the ratio of time-step and delay time get closer to unity. Therefore investigating the precise resolution dependence is important. We also emphasize that the ADP model needs additional features to be implemented in a merger simulation, but we have drawn on prior work and use methods that have been well studied in the literature.
The layout of this paper is as follows: In section 2 we discuss salient details of the simulation including the PNK11 model and how we have augmented it to allow for black hole tracking. In section 3 we compare the behaviour of the merger simulations and the sensitivity of final state diagnostics, specifically the black hole mass, as a function of the model parameters. We end with a brief review.
2 Numerical Simulations
We test the ADP method in a simulation of a major merger of two Milky Way-sized galaxies. Each galaxy contains a dark matter halo, hot gas halo, stellar disc, gas disc, a stellar bulge, and a black hole. The galaxies are initially separated by 70 kpc, and placed on a parabolic trajectory around one another similar to the merger considered in Springel et al. (2005b). The construction of the galaxies and definitions of the structural parameters are given in Appendix A, while the parameters of each component are presented in Table 4. The masses and particle numbers used are summarized in Table 1.
The simulations were run using the parallel version of Hydra (Couchman et al., 1995; Thacker & Couchman, 2006), which uses an Adaptive Particle-Particle, Particle-Mesh algorithm (Couchman, 1991) to calculate gravitational forces and the Smooth Particle Hydrodynamics method (SPH; Gingold & Monaghan, 1977; Lucy, 1977) to calculate gas forces. The star formation algorithm is described in Thacker & Couchman (2000).
To implement an AGN feedback algorithm motivated by the PNK11 model, we have had to make some modifications and additions. We have modified how feedback is returned to the local gas, and added both a black hole advection algorithm and a black hole merger algorithm. None of the issues related to black hole advection or mergers are addressed within the PNK11 paper since they consider accretion of gas on to a black hole that is embedded within an initially spherical gas cloud.
2.1 Implementing the PNK11 model: Approach to accretion
The black hole accretion algorithm is the same as in PNK11. In this method the accretion disc particle (ADP) is a collisionless sink particle containing both the black hole and its tightly bound accretion disc. The mass of the ADP is given by
where is the mass of the black hole and is the mass of the accretion disc. As in PNK11, we initialise , while we choose the initial black hole mass in each galaxy to be 10 M. The ADP has an associated ‘smoothing length’, , which is calculated at every iteration by , where a sphere with radius around the ADP includes 60 gas particles, and is the smallest resolved smoothing length in the SPH solver. The smoothing length is used to calculate the gas properties at the location of the ADP, which is used for analysis and returning feedback energy. The ADP is also given an accretion radius, , which is a free parameter of the model that determines at what radii particles are considered to have accreted on to the accretion disc. Nominally, and indeed in PNK11, this radius is expected to be on the order of a pc. While merger simulations can reach these small radii by considering arbitrarily small values of , it is worth emphasizing that the properties of the gas on these scales is not resolved. Indeed this is a generic problem with feedback models in general, namely that they use input parameters at the edge of model resolution. However, since this is an unavoidable problem, we continue on with this issue noted and examine the impact of varying the value of by up to a factor of 10.
Whenever a gas particle comes within of the ADP, it is instantly accreted on to the accretion disc, thus , where is the mass of the gas particle. The accreted particle’s mass and momentum are added to the ADP, and then the accreted particle is removed from the simulation. This capture rate is not limited in any way, and solely depends on the particles’ relative locations. Once accreted on to the accretion disc, it takes time for the gas to be transported through the disc so that it can finally be accreted on to the black hole. This time delay is on order of the disc viscous time, , which must be larger than the dynamical time at the accretion radius; PNK11 set as a free parameter, but argue that it should be 10–100 Myr in galaxy formation simulations.
The mass accretion rate on to the black hole from the accretion disc is given by,
To preserve mass continuity, the mass added to the black hole is removed from the accretion disc although the overall ADP mass remains the same. As in most numerical implementations of AGN feedback, the accretion rate on to the black hole is Eddington-limited; this moderates the growth of the black hole based upon a physical limit rather than just from the numerical accretion on to the disc.
2.2 Implementing the PNK11 model: Changes and additions for merger simulations
PNK11 use the same feedback scheme as described in Nayakshin & Power (2010). In this scheme, feedback energy is returned by adding wind particles which are radially directed away from the black hole; the wind particles have momentum , where is the velocity dispersion of the host galaxy. In our simulation, it is impractical to continually add wind particles since the particle load in the solver will climb rapidly, so we instead follow the wind prescription used in Debuhr et al. (2011). Here, the feedback rate is given by
where is the assumed infrared optical depth, and is the luminosity, where is the radiative efficiency. The momentum is returned radially, such that every gas particle within the ADP’s radius of influence, , receives an equal acceleration. We reiterate that has no explicit dependence on and is being recalculated at every iteration.
2.2.2 Black hole advection and mergers
The black hole advection algorithm is a modified version of that presented in Okamoto et al. (2008) and used in Model WT in Wurster & Thacker (2013). Here, the ADP is displaced towards the centre of mass of the sphere with radius which is centred on the ADP. The distance it is displaced is
where is the velocity of the black hole, is the time-step, and is the distance from the black hole to the centre of mass. The coefficients are based upon those in Okamoto et al. (2008), but modified for our higher resolution. Following the approach of Okamoto et al. (2008), we have selected the parameters to suppress the oscillations of the black hole particles such that the maximum amplitude of the black hole oscillation between the start of the simulation and first periapsis is less than 0.25 per cent of the core radius. This method continually displaces the black hole towards the centre of mass to counter any two-body interactions that may try to displace the black hole from the bottom of the potential well; limiting the distance preserves the possibility of the black hole remaining in a gas void. One option we considered but did not implement was to couple the ADP to the gas particle that has the lowest potential energy, is within of the ADP and satisfies , where is the relative velocity between the ADP and the gas particle and is the local sound speed; generally, this can lead to large artificial displacements and prohibits the black hole from remaining in a gas void, but with this specific accretion algorithm, it would artificially increase the accretion rate by arbitrarily moving the ADP within of a gas particle. A second option we considered was to use a tracer particle whose mass is 100–1000 times greater than any other particle; although this produces a smooth black hole trajectory, the additional mass affects the evolution time when compared to non-tracer particle models.
Finally, when the ADP’s pass within each others smoothing lengths and have a relative velocity less than the local sound speed, they are assumed to instantly merge. The merged black hole has the combined mass of the progenitor black holes, and the merged accretion disc has the combined mass of the progenitor discs, thus all masses are conserved. This is similar to the numerical merger procedure used in Springel et al. (2005b).
2.3 Parameter space and resolution dependence
To understand the impact of the two free parameters, and , we ran a suite of simulations within this parameter space and considered two separate resolutions for a total of 17 simulations. We plot where the simulations lie in the parameter space in Fig. 1, where each point corresponds to a model. Our parameter space is
In Table 2, we convert the accretion radius to physical units for both resolutions. Our model naming convention is understood as follows: ‘PNKrt’. If we do not include a resolution when referring to a model, then we are referring to both resolutions. In terms of wall-clock time, despite having a modest particle content, the fiducial resolution models still take over a month to run due to the large number of time steps required, the lack of multiple time-steps in the Hydra solver and also an slowdown that occurs as large numbers of SPH particles reach the minimum smoothing length of the solver. We thus have a limited number of fiducial resolution simulations.
|Fiducial resolution||Low resolution|
3.1 General evolution
Each of our models was evolved for 1.5 Gyr through a merger event, similar to that of Springel et al. (2005b). We found that each model followed a similar qualitative history; in Fig. 2, we plot the evolution of the gas column density of Model PNKr05t05.
From this evolution, we see four significant epochs: first periapsis, apoapsis, second periapsis and core merger, occurring at 170, 480, 880 and 990 Myr, respectively; the low resolution models reach second periapsis and core merger 25 and 6 Myr, respectively, earlier than their fiducial resolution counterparts. In all our fiducial resolution models, the black holes merge at 1.01 Gyr, and their low resolution counterparts merge 10–60 Myr later. The later merger in the low resolution models is a result of the binary black holes oscillating about one another with a greater amplitude than in the fiducial resolution models; this is from both the low resolution black holes being more massive by this epoch (see section 3.6), and there being fewer particles to induce drag on the black hole system.
The final remnant is a reformed gas disc and flattened stellar ellipsoid surrounded by a hot gas halo; the remnant is further discussed in section 3.7. The final gas discs in the fiducial resolution models have scale lengths of approximately 0.5 kpc, compared to the initial scale length of 2.46 kpc and are essentially an order of magnitude less massive than the discs in the initial conditions. Of the gas that was initially in a disc, 89 per cent of it is either accreted on to the black hole or converted into stars. With a final star formation rate of less than 0.5 M/yr, the final configuration is a red and dead elliptical, as expected (Springel et al., 2005a).
To quantify this evolution, we have plotted the total black hole mass, the accretion rates on to the black hole, and the gas density and gas temperature within of the ADP for Model PNKr05r05 in Fig. 3; the accretion rates and gas properties are geometrically averaged over both black holes and plotted in bins of 10 Myr.
For the first 300 Myr, , thus there are only 60 gas particles within this radius; after 300 Myr, , and there are several thousand particles within this radius. There is an increasing accretion rate after first periapsis, which peaks near apoapsis; there is a second peak at core merger. The gas density in the core also increases from first periapsis to apoapsis, and again at core merger. The core temperature stays relatively moderated after first periapsis. The initial decline in gas density and temperature is a result of the feedback energy carving out a void around the black hole in the initial density field which is quite symmetric. The SPH particles that contribute to the density and temperature at the location of the black hole are at , thus provide only a minimal contribution. Thus, these low densities and temperatures are essentially a result of the smoothing kernel breaking down due to a poor distribution of the particles within 2. At later times the number of particles increases and symmetry is lost which prevents any spurious temperature and density values from being calculated.
3.2 Comparison to other models in the literature
We briefly compare Model PNKr05t05 to two other models found in the literature: Model SDH, which is designed to reproduce the algorithm used in Springel et al. (2005b) and Model DQM, which is designed to reproduce the algorithm used in Debuhr et al. (2011). Model SDH uses the Bondi accretion rate given in (1), where , to calculate the rate of gas accretion on to the black hole; half a per cent of the accreted mass is converted in to feedback energy which is returned thermally by adding energy to the nearby gas particles, following the smoothing kernel. Model DQM uses a viscous accretion rate,
where is the dimensionless viscosity, is the mean gas surface density, and is the rotational angular velocity of the gas, which is based upon the multi-scale SPH simulations by Hopkins & Quataert (2010); the feedback algorithm is same as in our PNK models. We plot the total black hole mass, the accretion rates on to the black hole, and the gas density and gas temperature within of the black hole for Models SDH and DQM in Fig. 3.
We briefly compare Model PNKr05t05 to two other models found in the literature: Model SDH, which is designed to reproduce the algorithm used in Springel et al. (2005b) and Model DQM, which is designed to reproduce the algorithm used in Debuhr et al. (2011). Model SDH uses the Bondi accretion rate given in (1) to calculate the rate of gas accretion on to the black hole; the feedback energy is returned thermally by adding energy to the nearby gas particles. Model DQM uses a viscous accretion rate, which is based upon the multi-scale SPH simulations by Hopkins & Quataert (2010); the feedback algorithm is same as in our PNK models. We plot the total black hole mass, the accretion rates on to the black hole, and the gas density and gas temperature within of the black hole for Models SDH and DQM in Fig. 3.
Although all three accretion and feedback algorithms make different physical assumptions, the final total black hole masses in the fiducial (low) resolutions are equal within a factor of 1.52 (4.62). The accretion histories, however, vary considerably. As expected, Models SDH and DQM have continual accretion which starts immediately, but the accretion rate in PNKr05t05 is punctuated by periods of nearly zero accretion, especially at early times. The latter is a result of feedback from the first few accretion events creating a small and transient void around the ADP, as can be seen in the core density profile. After first periapsis, tidal torques are large enough and feedback energy is low enough to allow additional accretion events.
The core density of Model DQM is the lowest of these three models; in Model DQM, kinetic feedback builds up in the gas, and after 60 Myr, a void of is formed and maintained for the duration of the simulation. Although Model PNKr05t05 uses the same feedback algorithm, the initially lower accretion rate prevents this build up of kinetic feedback, thus this void never forms, hence the higher core density. The core temperature is highest for Model SDH, where thermal feedback energy directly increases the temperature of the gas. The core heating of Models PNK and DQM is from shock heating and star formation; the higher accretion and star formation (not shown) rates in PNKr05t05 after first periapsis results in more heating than in Model DQM.
3.3 Acceptable parameter ranges
We ran 12 low resolution simulations to probe our entire – parameter space. At our fiducial resolution we only tested five models due the time required to complete individual simulations. After analysing the results, each simulation was classified as either physical or unphysical, as is indicated by the dot type in Fig. 1. The basis for our definition of physical or unphysical relies upon a combination of structural evolution and how close the final remnant lies to the – relationship. We note that relying upon the – relationship is perhaps not a strong constraint because the masses of the seed black holes could be changed, however it remains a commonly used approach in these simulations and we thus proceed cautiously with this issue noted.
For values that are very small the accretion rate on to the black hole will be limited. The expectation in this case is that the resulting remnant will fall below the – relation. Our numerical experiments show this to be true for all the models with , with the black hole mass in PNKr02t05 being 4.5 times lower than expected, assuming the velocity dispersion remains constant. These models also suffer from sensitivity to operation ordering in the solver. At a given resolution, we would expect the absorption of particles to occur at the same time provided that the operations in the solver are executed in the same order, and hence produce the same round-off. However, when running in parallel with dynamic load-balancing this is no longer the case and round-off in accumulations can produce subtle differences around or just above machine precision. Thus, with too small of an accretion radius, pseudo-random issues due to machine precision dominate the accretion, and the results of a given simulation are not easily reproducible. On the basis of lying beneath the – relation and the exhibited numerical sensitivity, we classify all of our models as ‘unphysical’.
At the other extreme of our parameter space, there can be unphysical results if the accretion rate on to the disc is too high. If is too large, then many particles can pass within this radius and be accreted on to the disc, regardless of the amount of feedback. Likewise, a large means that the feedback energy is being returned at a lower rate, thus particles have less of an obstacle to overcome to pass within when influenced by outside forces (i.e. tidal torques from the interacting galaxies). Both result in a large accretion rate on to the disc, where the accretion disc mass can surpass the mass of the black hole (see bottom right panel of Fig. 4; the remaining three panels show the profiles for physical models, where the peak accretion disc mass is M).
This unreasonably large accretion rate on to the disc yields a large and continual accretion rate from the disc on to the black hole, which results in a vast amount of feedback energy. Once a critical amount of kinetic feedback is injected into the gas, namely that sufficient to overcome the nuclear gravitational binding energy, a large, unphysical void is formed in the galaxy. This in turn suppresses the accretion rate on to the ADP, but not from the disc on to the black hole. The accretion disc is ultimately depleted, the feedback ends and the void recollapses to begin another cycle. In situations where a large radius and large viscous time-scale are included it is possible to remove the gas from the system.
These unphysical voids often foreshadow the total disruption of the remnant. Starting at 1.04 Gyr in Model PNKr05t10, there is a catastrophic explosion, and the system is totally disrupted within 60 Myr. In other models, specifically our low resolution models, the core merger causes an increase in the accretion rate, leading to a final feedback episode that disrupts the system, although often less catastrophically than Model PNKr05t10. While the system is highly disrupted, there are often no noticeable voids, and the process is ‘gentle’ enough that the gas begins to recondense into a disc and cloud. We cautiously and liberally classify these results as physical.
In Fig. 5, we have plotted the actual and the Eddington accretion rates for six models; we also give the length of time each system is undergoing Eddington-limited accretion in the top left corner of each panel.
The unphysical model PNKr05t10 is undergoing Eddington-limited accretion for 11.7 per cent of the time, thus a (relatively) high percentage is a possible indicator of an unphysical model. However, the unphysical Model PNKr02t05 undergoes Eddington-limited accretion 6.3 per cent of the time, which is higher than the per cent of time the physical Model DQM is undergoing Eddington-limited accretion. While there is a correlation between the amount of time spent in Eddington-limited accretion and the final black hole mass for the PNK models, this is not true when comparing across different models. Further, these numbers indicate that relying upon the – relationship appears to place the amount of time in Eddington-limited accretion in a comparatively narrow band. It is thus difficult for us to draw detailed conclusions about using the amount of time spent in Eddington-limited accretion as a strong discriminant between the physical or unphysical nature of models.
Ultimately, the amount of time accreting at the Eddington limit is very closely related to the duty cycle of the quasar. In the context of models such as Small & Blandford (1992), where the population of SMBH is predicted as a function of luminosity (or mass) and time, the duty cycle directly impacts the global population because the luminosity function is a product of the number density, mass and duty cycle. Models of this type have been extensively developed by Shankar (e.g. Shankar et al., 2013; Shankar et al., 2010) and show that while Eddington accretion is likely to be more common at lower masses ( M), once the SMBH reaches 10 M, matching the luminosity function requires that growth at the Eddington rate must be considerably less frequent or radiative efficiencies must be low. While obviously not directly equivalent, matching the – relationship appears to place a very similar constraint on the allowed duty cycle.
Lastly, the gas discs in the simulation are quite stable and in the absence of any additional supply of infalling gas actually decrease their mass during the simulation by 63 per cent prior to the start of the main merger. The stability of the discs is a product both of our choosing a comparatively stable initial configuration and also that the SN feedback routine keeps the gas comparatively well supported against cold collapse. Thus the discs do not fuel the black hole through large scale fragmentation due to them becoming unstable.
3.4 Parameter sensitivity:
We have four models (two models each of fiducial and low resolution) with Myr that produce physically plausible results. In Fig. 6, we have plotted total black hole mass, the accretion rates on to the black hole, and the gas density and gas temperature within of the ADP.
As expected, the first accretion event occurs in Model PNKr10t05, and Model PNKr10t05 has a slightly higher accretion rate for most of the simulation. By increasing the radius of accretion by a factor of two (hence the volume in which a particle can be accreted, , by a factor of eight), there is only a factor of 1.35 and 1.70 increase in final total black hole mass for fiducial and low resolutions, respectively. The gas density and temperature behaviour within are very similar at both resolutions. The value of has no explicit dependence on and all models produce the same qualitative behaviour for , which settles to between first periapsis and apoapsis. We find that has negligible influence on and a minimal impact on the final black hole masses. Thus, we conclude that, as long as the parameters are in the physically acceptable range, results are not overly sensitive to the exact value of .
3.5 Parameter sensitivity:
We have five models with (two at the fiducial resolution and three at low resolution). In Fig. 7, we have again plotted total black hole mass, the accretion rates on to the black hole, and the gas density and gas temperature within of the ADP.
The total black hole masses in the fiducial resolution simulations are similar throughout the simulations, and their final masses differ by 13 per cent. At low resolution, the total black hole masses for models PNKr05t05 and PNKr05t01 match within 18 per cent at the end of the simulation; if we include PNKr05t10 (a model that produces unphysical behaviour at the fiducial resolution), then the low resolution black hole masses differ at most by a factor of 1.68.
In Fig. 8, we plot two segments of the accretion rate for one black hole in our three low resolution models with .
As expected, the first accretion event happens simultaneously in each model (top panel). This is followed by a period of Eddington accretion, which stops first for Myr and last for Myr depending upon when the rate falls below the Eddington rate. Once the Eddington accretion has stopped, the rate of feedback (where ) differs for each model, but at least at this early time the accretion rate on to the black hole can still be calculated: an exponentially decaying rate with the time dependence given by . Although PNKr05t01 accretes at the Eddington limit very slightly longer than the other two models, the faster decrease in the accretion rate means that feedback effectively stops sooner, thus giving particles more time to slow down and fall within .
At late times, each accretion event does not add enough mass to the disc to permit Eddington accretion (bottom panel); however, there is a local accretion peak followed by exponential decay. As previously discussed, the faster decay of Myr allows for rapid accretion events on to the disc; this creates large variances in the accretion rates, spanning a few orders of magnitude. The slower decay of Myr yields a moderated accretion rate, thus there is only a factor of a few between the local peak accretion rate and the minimum rate prior to the next accretion event. The time-average accretion rate yields a higher accretion rate for the lower viscous time-scales, resulting in more feedback being returned to the gas which can, at least temporarily, expel considerable amounts of gas from the system. Summarily, we find a general trend to more massive black holes with decreasing .
The torques from the interacting galaxies modify this argument by introducing disc instabilities Although the feedback rate of the Myr models can hinder gas from accreting on to the disc during a quiescent phase, there is not enough cumulative feedback to prevent a large gas flow on to the disc produced by the disc instabilities. This results in a large accretion epoch, followed by powerful outbursts of feedback energy. Thus, the tidal torques at first periapsis and core merger are strong enough to over come the feedback and cause a catastrophic (or near catastrophic) outburst, as in Model PNKr05t10.
Similar to our discussion in section 3.4, we find quantitatively similar gas behaviour within , which is a result of the ’s being similar across all models. As with , we conclude that, as long as lies with the physically acceptable range, then the simulations are comparatively insensitive to its exact value (at least to within 1/2 an order of magnitude at the fiducial resolution).
3.6 Resolution sensitivity
We have three models that are deemed to have physically plausible results at both resolutions. In Fig. 9, we have plotted total black hole mass, the accretion rates on to the black hole, and the gas density and gas temperature within of the ADP.
When we directly compare two models at different resolutions, there are significantly more differences in accretion and evolution histories than we do for when the model parameters are changed (at least for the acceptable range). As would be expected, in all cases the fiducial resolution model accretes the first particle, while the lower resolution simulations produce a larger jump in the accretion rate at early times due to the particles in those simulations being more massive. By the time the galaxies reach first periapsis at 166 Myr, the low resolution gas has a greater kinetic feedback obstacle to overcome before it can pass within of the ADP. Thus, the first low resolution accretion event immediately begins to moderate the accretion on to the disc, and this moderation continues to persist after first periapsis.
Because it has a smaller kinetic feedback obstacle to overcome, the fiducial resolution gas can more easily fall within of the ADP after first periapsis, leading to a large and essentially unmoderated accretion on to the disc. As the gas accretes from the disc on to the black hole, it modifies the environment to suppress further accretion events on to the disc; however, the disc remains massive from the previous accretion episode. Thus, there is a decrease in gas density in the core, but the accretion on to the black hole continues as the gas in the disc is continually transferred to the black hole. This major epoch of accretion on to the disc after first periapsis results in a steeper black hole growth between first periapsis and apoapsis for the fiducial resolution models.
For all three models, the fiducial resolution models have the expected higher core gas densities. Since the feedback is being returned kinetically, shock heating and star formation are the primary heating mechanisms, and both scale with resolution in turn producing similar core temperatures.
Thus, based upon the two resolutions we test, we can conclude that resolution has a greater impact on the results than the values of the free parameters, and . However, these differences are predictable since there are more epochs of discrete behaviour in the fiducial resolution models than in the low resolution models, allowing for a more continuous modelling of the evolution. Also, each fiducial–low resolution pair is more similar to one another than to Models SDH or DQM of the same resolution, indicating that the model can be distinguished from other models even at low resolution.
3.7 Final states
3.7.1 Stellar remnant
A common test of numerical accretion and feedback algorithms is the – relation. For all of our physical models, we have calculated the stellar velocity dispersion around every black hole using the method described in Debuhr et al. (2011). These velocity dispersions are averaged over 1000 random lines of sight, and are plotted on the – relation in Fig. 10.
For the six models in Fig. 9, we also plot the full range of velocity dispersions, as well as those taken preferentially along the -, - and -lines of sight. The large range of is a result of the highly triaxial stellar remnant; see Fig. 11, where we have plotted a fiducial and low resolution stellar remnant. The fiducial resolution models have average values very near the observed – relation, and a range that nearly stays within the one-sigma scatter. The low resolution stellar remnants are more elliptical than their fiducial resolution counterparts, thus they have a greater range of ’s. While some of the average ’s fall outside of the one-sigma range, they all obtain ’s near the observed – relation if the line of sight is near the plane of the disc. To verify these results, we have recalculated the velocity dispersions of the fiducial resolution models based upon the gravitational potential of the black hole, and found these values to be consistent within five per cent of the average values reported above.
3.7.2 Gas properties
In Figs. 12 and 13, we plot the gas column density of ten remnants at 1.5 Gyr. These show, respectively, the inner 100 and 20 kpc of six PNK models (three models as each resolution) and Models DQM and SDH.
As discussed in sections 3.4 and 3.5, there is much qualitative similarity amongst the three fiducial remnants of the PNK models, as well as Models DQM and SDH. In all cases, there is a condensing gas cloud with a reformed disc.
The radius, surface density profiles and total gas mass of the reformed discs in the PNK models are only slightly dependent on the free parameters; the surface density profile is plotted in the top panel of Fig. 14 and the gas disc masses are given in Table 3.
|Model||Disc mass ( M)|
After core merger, PNKr10t05 has a slightly larger accretion rate than PNKr05t05 due to its larger accretion radius, this leads to more feedback energy and a slightly extended disc with lower surface density. Model PNKr05t01 has a higher accretion rate at core merger than PNKr05t05, which results in a small outburst at 1 Gyr, expelling some gas from the system. By 1.5 Gyr, the remaining gas is less bound leading to the slightly larger radius and lower surface density.
The fiducial resolution PNK models have different surface density profiles than Models DQM and SDH; see the bottom panel of Fig. 14 and note the difference range on the vertical axis. The disc in Model DQM has a dense core and a moderately dense torus (resulting in a high-mass disc), whereas the profile for Model SDH is a torus due to low level feedback activity carving out a region of the disc.
The three fiducial resolution PNK remnants are qualitatively more similar to one another than to their low resolution counterparts. Both PNKr05t01 and PNKr05t05 underwent a major outburst starting at 1.05 Gyr, blowing away much of the gas and preventing the reformation of the gas disc. In the subsequent few 100 Myrs, the gas begins to recondense into the cloud presented here. In PNKr10t05, this major outburst never occurs, thus allowing the disc to reform.
We have implemented the accretion disc particle (ADP) method of Power et al. (2011) into a major merger simulation of two Milky Way-sized galaxies. We ran five fiducial resolution simulations and twelve low resolution simulations varying the free parameters and . Our primary conclusions are as follows:
For accretion radii that are too small (i.e. ), the final black hole mass is far smaller than predicted from the – relationship. Thus all of our models are classified as unphysical.
For accretion radii that are too large (i.e. ), the accretion rate on to the accretion disc (hence on to the black hole) is unreasonably large. The resulting feedback is enough to catastrophically disrupt the system. We thus classify all of our models as unphysical.
For large viscous time-scales (i.e. Myr), feedback from an accretion event persists long enough to hinder secular accretion. Tidal forces from the interacting galaxies can overcome the low amount of feedback and funnel considerable amounts of gas on to the accretion disc. Depending on resolution and accretion radius, this short accretion epoch can be large enough such that its resulting feedback energy can catastrophically disrupt the system. Thus a few of our Myr models are classified as unphysical.
The values of and , assuming they were in the allowed parameter space, had minimal affect on the gas properties within of the black hole. This result was expected since had no explicit dependence on or .
The exact value of has only a minimal affect on the resulting system, assuming that it is in the allowed parameter space. By doubling the accretion radius from to , the final black hole mass increases by only a factor of 1.35 (1.70) for our fiducial (low) resolution simulations.
As we decrease the value of the final black hole mass increases; the final range of black hole masses spans a factor of 1.13 (1.68) for our fiducial (low) resolution models. Thus, the exact value of has only a minimal affect on the resulting (physical) system.
Decreasing the resolution increases the final black hole mass; for any given fiducial–low resolution pair, the final black hole mass differs by at most a factor of 1.90. The fiducial resolution models experience steeper black hole growth between first periapsis and apoapsis, and the lower resolution models are more prone to major outburst events shortly after core merger. Given the parameters tested here, resolution has the largest impact on the outcome of the model. We understand that these differences are unavoidable and are a result of a single accretion event in the low resolution models returning 7.6 times more energy to the gas than the fiducial resolution models.
In AGN feedback models where the accretion rate is dependent on the gas properties near the black hole, the accretion and energy feedback rates are never zero, and span only a few orders of magnitude; in our simulations, the accretion and feedback begins immediately. In Model DQM, an initially high accretion rate leads to an immediate and dramatic modification of the black hole’s environment. In Model SDH, the accretion rate does not become significant until apoapsis, at which point the black hole undergoes a major accretion epoch, with its mass increasing by a factor of 20.5 in 470 Myr.
With the ADP model, the initial accretion rate is zero or very small. The minimal modification of the ADP’s environment by first periapsis allows a large accretion event on to the disc, which results in a large accretion rate on to the black hole, leading to a very rapid black hole growth; in the case of Model PNKr05t05, the black hole mass increases by a factor of 19.6 in 200 Myr. Although these three models, Models SDH, DQM and PNK, have very different algorithms and different evolution histories, there are similar remnants in all three cases.
The ADP algorithm decouples the accretion rate on to the black hole from the gas properties around it and allows a high accretion rate even after feedback has prevented additional accretion events on to the disc. Yet to have confidence in this conceptual model a precise observational strategy for determining parameters is necessary (essentially looking at the correlation, or lack of, between activity and the nuclear gas environment). Recent observational work (Wild et al., 2010) has unveiled how star formation activity and black hole growth appear tied together in spheroids. These results undoubtedly give useful hints on how black hole growth proceeds outside of rapid accretion phases, but fueling of black hole via cold gas accretion remains the biggest uncertainty in our models currently. To date, studies of nuclear CO morphology such as the NUGA project (e.g. García-Burillo et al. 2005; van der Laan et al. 2011) have found no obvious morphological links between local AGN activity and mid-scale CO morphology. But such studies are obviously limited by resolution concerns, as well as probing an entirely different AGN luminosity regime. However, ALMA in its full configuration with the largest baselines will provide resolutions that are sufficient to study the nuclear regions the nearest active quasar systems. The prospect of examining morphology dependence in these systems is immensely exciting and should provide great insight into their evolution and how they can best be modelled.
We thank the anonymous referee for helpful comments that improved the content and clarity of the manuscript. We thank Larry Widrow for providing his NFW halo generator. JW is supported by NSERC and Saint Mary’s University. RJT is supported by a Discovery Grant from NSERC, the Canada Foundation for Innovation, the Nova Scotia Research and Innovation Trust and the Canada Research Chairs Program. Simulations were run on the CFI-NSRIT funded St Mary’s Computational Astrophysics Laboratory.
Appendix A Galaxy models
We construct our model galaxy to have similar mass components to the galaxies presented in Springel et al. (2005b). We begin by using the GalactICs package (Kuijken & Dubinski, 1995; Widrow & Dubinski, 2005; Widrow et al., 2008) to create a Milky Way-sized galaxy that consists of a stellar bulge, stellar disc, and a dark matter halo; this is done through an iterative process to produce a self consistent system. The free parameters, and the values we chose, are listed in Table 4.
|Bulge||292 km s|
|119 km s|
|Dark Matter Halo||13.6 kpc|
|330 km s|
|mean molecular weight||0.6|
The dark matter halo profile follows
where is the halo scale length, is the cutoff radius, is the central cusp strength and is a (line of sight) velocity scale that sets the mass of the halo. The truncation function, , smoothly goes from one to zero at over width .
The stellar disc has a truncated density profile that falls off approximately exponentially in and follows sech in ; the disc has radial and vertical scale heights and and truncation distances and . The radial velocity dispersion is given by , where is the central velocity dispersion and for simplicity.
The stellar bulge density profile is given by
which yields a Sérsic law for the projected density profile if , where is a free parameter. The constant is defined using , where is the depth of the gravitational potential associated with the bulge and is the radial scale parameter. The variable is adjusted such that encloses half the total projected light or mass.
The galaxy is then modified in three ways. First, ten per cent of the total stellar mass is converted into gas to create a gas disc. The gas disc is the same as the stellar disc except that it has been reflected in the plane to avoid coincidence with the star particles. The gas scale height is initially larger than physically motivated, however, cooling allows the gas to collapse into a thin disc within a few 10 Myr. This vertical collapse produces a short transient evolution of the gas accompanied by a brief increase in the star formation rate. At resolutions higher than presented here, this produces a strong ring-shaped shock which propagates outwards; one possible solution is to reduce the scale height of the gas disc, but studies in Williamson & Thacker (2012) show that this is not necessary for the resolutions presented here.
Second, a hot gas halo is added following the observationally motivated -profile (e.g. Cavaliere & Fusco-Femiano, 1976):
where is the central density, is the core radius, and is the outer slope parameter; we choose kpc and as done in Moster et al. (2011). We set by choosing the mass of the hot gas halo within 40 kpc to be equal to two per cent of the total disc mass (Rasmussen et al., 2009). To conserve total halo mass, we reduce the mass of the dark matter halo by the total mass of the hot gas halo. By assuming isotropy and hydrostatic equilibrium, the temperature profile of the hot gas halo is given by (Kaufmann et al., 2007)
where is the mean molecular weight, is the proton mass, is the Boltzmann constants, and is the total mass interior to . The halo gas is given an initial angular momentum which scales with the circular velocity, , where is the distance from the spin axis of the galaxy (Moster et al., 2011).
Finally, a black hole (sink) particle is placed at the centre of the galaxy. The galaxy has a total mass of M, 1 287 743 (168 351) particles for the fiducial (low) resolution simulations, and a Plummer softening length of pc ( pc).
- Bondi (1952) Bondi H., 1952, MNRAS, 112, 195
- Booth & Schaye (2009) Booth C. M., Schaye J., 2009, MNRAS, 398, 53
- Cavaliere & Fusco-Femiano (1976) Cavaliere A., Fusco-Femiano R., 1976, A&A, 49, 137
- Couchman (1991) Couchman H. M. P., 1991, ApJ, 368, L23
- Couchman et al. (1995) Couchman H. M. P., Thomas P. A., Pearce F. R., 1995, ApJ, 452, 797
- Debuhr et al. (2011) Debuhr J., Quataert E., Ma C.-P., 2011, MNRAS, 412, 1341
- Di Matteo et al. (2005) Di Matteo T., Springel V., Hernquist L., 2005, Nature, 433, 604
- Fabian (1999) Fabian A. C., 1999, MNRAS, 308, L39
- Ferrarese & Merritt (2000) Ferrarese L., Merritt D., 2000, ApJ, 539, L9
- García-Burillo et al. (2005) García-Burillo S., Combes F., Schinnerer E., Boone F., Hunt L. K., 2005, A&A, 441, 1011
- Gebhardt et al. (2000) Gebhardt K. et al., 2000, ApJ, 539, L13
- Gingold & Monaghan (1977) Gingold R. A., Monaghan J. J., 1977, MNRAS, 181, 375
- Gültekin et al. (2009) Gültekin K. et al., 2009, ApJ, 698, 198
- Hobbs et al. (2011) Hobbs A., Nayakshin S., Power C., King A., 2011, MNRAS, 413, 2633
- Hopkins & Quataert (2010) Hopkins P. F., Quataert E., 2010, MNRAS, 407, 1529
- Hopkins & Quataert (2011) Hopkins P. F., Quataert E., 2011, MNRAS, 415, 1027
- Kaufmann et al. (2007) Kaufmann T., Mayer L., Wadsley J., Stadel J., Moore B., 2007, MNRAS, 375, 53
- King (2003) King A., 2003, ApJ, 596, L27
- King (2005) King A., 2005, ApJ, 635, L121
- King (2010) King A. R., 2010, MNRAS, 402, 1516
- Kormendy & Richstone (1995) Kormendy J., Richstone D., 1995, ARA&A, 33, 581
- Kuijken & Dubinski (1995) Kuijken K., Dubinski J., 1995, MNRAS, 277, 1341
- Kurosawa et al. (2009) Kurosawa R., Proga D., Nagamine K., 2009, ApJ, 707, 823
- Lucy (1977) Lucy L. B., 1977, AJ, 82, 1013
- Magorrian et al. (1998) Magorrian J. et al., 1998, AJ, 115, 2285
- Marconi & Hunt (2003) Marconi A., Hunt L. K., 2003, ApJ, 589, L21
- McLure & Dunlop (2002) McLure R. J., Dunlop J. S., 2002, MNRAS, 331, 795
- Moster et al. (2011) Moster B. P., Macciò A. V., Somerville R. S., Naab T., Cox T. J., 2011, MNRAS, 415, 3750
- Nayakshin & Power (2010) Nayakshin S., Power C., 2010, MNRAS, 402, 789
- Okamoto et al. (2008) Okamoto T., Nemmen R. S., Bower R. G., 2008, MNRAS, 385, 161
- Power et al. (2011) Power C., Nayakshin S., King A., 2011, MNRAS, 412, 269, (PNK11)
- Puchwein et al. (2008) Puchwein E., Sijacki D., Springel V., 2008, ApJ, 687, L53
- Rasmussen et al. (2009) Rasmussen J., Sommer-Larsen J., Pedersen K., Toft S., Benson A., Bower R. G., Grove L. F., 2009, ApJ, 697, 79
- Sanders et al. (1988) Sanders D. B., Soifer B. T., Elias J. H., Madore B. F., Matthews K., Neugebauer G., Scoville N. Z., 1988, ApJ, 325, 74
- Scannapieco & Oh (2004) Scannapieco E., Oh S. P., 2004, ApJ, 608, 62
- Scannapieco et al. (2005) Scannapieco E., Silk J., Bouwens R., 2005, ApJ, 635, L13
- Scannapieco et al. (2008) Scannapieco E., Thacker R. J., Couchman H. M. P., 2008, ApJ, 678, 674
- Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
- Shankar et al. (2010) Shankar F., Crocce M., Miralda-Escudé J., Fosalba P., Weinberg D. H., 2010, ApJ, 718, 231
- Shankar et al. (2013) Shankar F., Weinberg D. H., Miralda-Escudé J., 2013, MNRAS, 428, 421
- Sijacki et al. (2007) Sijacki D., Springel V., Di Matteo T., Hernquist L., 2007, MNRAS, 380, 877
- Silk & Rees (1998) Silk J., Rees M. J., 1998, A&A, 331, L1
- Small & Blandford (1992) Small T. A., Blandford R. D., 1992, MNRAS, 259, 725
- Springel et al. (2005a) Springel V., Di Matteo T., Hernquist L., 2005a, ApJ, 620, L79
- Springel et al. (2005b) Springel V., Di Matteo T., Hernquist L., 2005b, MNRAS, 361, 776
- Thacker & Couchman (2000) Thacker R. J., Couchman H. M. P., 2000, ApJ, 545, 728
- Thacker & Couchman (2006) Thacker R. J., Couchman H. M. P., 2006, Computer Physics Communications, 174, 540
- Thacker et al. (2006) Thacker R. J., Scannapieco E., Couchman H. M. P., 2006, ApJ, 653, 86
- Tremaine et al. (2002) Tremaine S. et al., 2002, ApJ, 574, 740
- van der Laan et al. (2011) van der Laan T. P. R. et al., 2011, A&A, 529, A45
- Widrow & Dubinski (2005) Widrow L. M., Dubinski J., 2005, ApJ, 631, 838
- Widrow et al. (2008) Widrow L. M., Pym B., Dubinski J., 2008, ApJ, 679, 1239
- Wild et al. (2010) Wild V., Heckman T., Charlot S., 2010, MNRAS, 405, 933
- Williamson & Thacker (2012) Williamson D. J., Thacker R. J., 2012, MNRAS, 421, 2170
- Wurster & Thacker (2013) Wurster J., Thacker R. J., 2013, In Preparation