## Abstract

We describe techniques for incorporating feedback from star formation and black hole (BH) accretion into simulations of isolated and merging galaxies. At present, the details of these processes cannot be resolved in simulations on galactic scales. Our basic approach therefore involves forming coarse-grained representations of the properties of the interstellar medium (ISM) and BH accretion starting from basic physical assumptions, so that the impact of these effects can be included on resolved scales. We illustrate our method using a multiphase description of star-forming gas. Feedback from star formation pressurizes highly overdense gas, altering its effective equation of state (EOS). We show that this allows the construction of stable galaxy models with much larger gas fractions than possible in earlier numerical work. We extend the model by including a treatment of gas accretion onto central supermassive BHs in galaxies. Assuming thermal coupling of a small fraction of the bolometric luminosity of accreting BHs to the surrounding gas, we show how this feedback regulates the growth of BHs. In gas-rich mergers of galaxies, we observe a complex interplay between starbursts and central active galactic nuclei (AGN) activity when the tidal interaction triggers intense nuclear inflows of gas. Once an accreting supermassive BH has grown to a critical size, feedback terminates its further growth and expels gas from the central region in a powerful quasar-driven wind. Our simulation methodology is therefore able to address the coupled processes of gas dynamics, star formation and BH accretion during the formation of galaxies.

## Introduction

It is now recognized that galaxy collisions and mergers are relevant to a wide range of phenomena associated with both ordinary and active galaxies. The collisionless dynamics of this process is relatively well understood. The seminal work of Toomre & Toomre (1972) using restricted *N*-body methods demonstrated that gravitational tidal forces between discs interacting transiently can account for the narrow tails and bridges commonly seen in morphologically peculiar galaxies. Subsequent work using self-consistent models (Barnes 1988, 1992; Hernquist 1992, 1993c) verified these conclusions and made it possible to extend the calculations to investigate mergers and the detailed structure of the remnants left behind. These studies showed that in many respects the remnants closely resemble elliptical galaxies, as proposed by Toomre (1977).

However, it is clear that pure stellar dynamics are not adequate to explain the extreme variety of objects linked to galaxy encounters. This possibility was already anticipated by Toomre & Toomre (1972) who suggested that collisions could ‘bring *deep* into a galaxy a fairly *sudden* supply of fresh fuel’, leading to a period of violent, enhanced star formation. The significance of non-gravitational processes in galaxy interactions was confirmed in the 1980s when it was recognized that the brightest infrared sources seen with the *IRAS* satellite are invariably associated with peculiar galaxies (Sanders et al. 1988; Melnick & Mirabel 1990). A natural interpretation of these systems is that the infrared emission is powered by a starburst, triggered through mergers.

There are also numerous studies indicating that quasars, radio galaxies and active galactic nuclei (AGN) are preferentially found in tidally disturbed objects (for reviews, see e.g. Barnes & Hernquist 1992; Jogee 2004). Indeed, much circumstantial evidence suggests an evolutionary connection between merger-induced starbursts and quasar activity (Sanders et al. 1988). In addition, in recent years, strong dynamical evidence has been accumulated indicating that supermassive black holes (BHs) reside at the centre of most galaxies (Kormendy & Richstone 1995; Magorrian et al. 1998; Kormendy & Gebhardt 2001). Moreover, a remarkable connection between supermassive BHs and the properties of their host galaxies has been discovered: the masses of central supermassive BHs are tightly correlated with the stellar velocity dispersion of their host galaxy bulges (Ferrarese & Merritt 2000; Gebhardt et al. 2000; Tremaine et al. 2002), as well as with the mass of the spheroidal component (Magorrian et al. 1998; McLure & Dunlop 2002; Marconi & Hunt 2003). This strong link suggests a fundamental connection between the growth of BHs and the formation of stellar spheroids in galaxy haloes (Kauffmann & Haehnelt 2000; Di Matteo et al. 2003; Volonteri, Haardt & Madau 2003; Wyithe & Loeb 2003; Di Matteo et al. 2004; Granato et al. 2004; Haiman, Ciotti & Ostriker 2004, and references therein).

A better understanding of galaxy interactions and mergers will, therefore, require simultaneous accounting of the physics of interstellar gas, star formation, the growth of supermassive BHs in galactic nuclei, and various forms of feedback associated with both massive stars and AGN. As we discuss in detail below, it is not immediately obvious how these effects should be incorporated into simulations, primarily because we presently lack a sufficiently developed theory of star formation, but also because the current generation of computer models cannot resolve the complex structure of star-forming gas on the scales of whole galaxies. For these reasons, efforts to study non-gravitational processes in galaxy mergers have to rely on strong simplifications, which have often been very restrictive in previous work.

Negroponte & White (1983) used a sticky-particle treatment to show that galaxy mergers can concentrate gas in the inner regions of a remnant. Later, Hernquist (1989) and Barnes & Hernquist (1991, 1996) examined the fate of gas in minor and major mergers with a smoothed particle hydrodynamics (SPH) algorithm, neglecting star formation and feedback and approximating the structure of the interstellar medium (ISM) with an isothermal equation of state (EOS). The first serious attempts to model star formation in galaxy mergers were made by Hernquist & Mihos (1995) and Mihos & Hernquist (1996), who also used an SPH method with an isothermal gas and included a minimal form of kinetic feedback from massive stars (Mihos & Hernquist 1994b). Gerritsen & Icke (1997) and Springel (2000) began to generalize these calculations by implementing simple approaches to allow for departures from an isothermal EOS.

While the earlier calculations yielded many successes by providing an explicit theoretical link between galaxy mergers and starbursts, many related phenomena remain poorly understood, partly because of the approximations made in handling the consequences of feedback. In this paper, we introduce a new methodology to incorporate star formation and BH growth into galaxy-scale simulations using a subresolution approach. Starting from a ‘microscopic’ physical theory for, e.g., the ISM, we form a ‘macroscopic’ coarse-grained representation that captures the impact of star formation and BH growth on resolved scales. This strategy is not tied to a particular underlying model and can be used to graft any well-specified theoretical description onto the simulations. As a specific example, in what follows we will illustrate our technique using a subresolution multiphase model for star-forming gas developed by Springel & Hernquist (2003a, hereafter SH03). As an additional component, we will add a treatment of gas accretion onto supermassive BHs, modelled as collisionless ‘sink’ particles, and we will treat feedback processes associated with the accretion.

Note that the use of a subresolution approach to include unresolved effects in simulations is a general concept also used, for example, in *N*-body representations of collisionless fluids, where six-dimensional phase space is discretized into fluid elements that are computationally realized as particles. Our treatment is analogous to this, except that we endow the particles with internal degrees of freedom that characterize, for example, the thermodynamic state of star-forming gas. These internal degrees of freedom evolve in a manner consistent with the basic theory from which they are derived, enabling us to couple the ‘microphysics’ of the ISM to the larger scale dynamics of the galaxies.

This paper is organized as follows. We first discuss the construction of compound galaxy models used in our study in Section 2. We then summarize our microscopic descriptions of star formation and BH growth in Sections 3 and 4, respectively, and discuss aspects of our numerical treatment in Section 5. We use our models to examine the impact of various forms of feedback on isolated discs in Section 6 and major mergers in Section 7. Finally, we conclude in Section 8.

## Model Galaxies

Each galaxy in our study consists of a dark matter halo, a rotationally supported disc of gas and stars, and a central (optional) bulge. The parameters describing each component are independent, so that a wide range of morphological types can be specified. The models are constructed in a manner similar to the approach described by Hernquist (1993a) and Springel (2000), but with a number of refinements. In particular, we use halo density profiles motivated by cosmological simulations and include pressure gradients in setting the equilibrium structure of the gas. The former consideration is important for the development of tidal features in interacting systems (Dubinski, Mihos & Hernquist 1996; Mihos, Dubinski & Hernquist 1998; Dubinski, Mihos & Hernquist 1999; Springel & White 1999), while the latter improvement is significant when the gas fraction of the disc is large, or if the gas has a relatively stiff EOS.

### Density profiles

#### Haloes

In the earlier work of Hernquist (1993a), haloes were described by truncated isothermal spheres with constant density cores. However, *N*-body simulations such as those of Dubinski & Carlberg (1991) and Navarro, Frenk & White (1996, hereafter NFW) indicate that cosmological haloes have density profiles that are centrally peaked and drop at large radii more rapidly with radius than isothermal spheres.

In view of this, we model the dark matter mass distribution with a Hernquist (1990) profile, i.e.

with cumulative mass profile*M*(<

*r*) =

*M*

_{dm}

*r*

^{2}/(

*r*+

*a*)

^{2}. There are several motivations for this choice. In its inner parts, the shape of this profile agrees with the NFW fitting formula, but due to its faster decline in the outer parts, the total mass converges, allowing the construction of isolated haloes without the need for an ad hoc truncation. Furthermore, in many situations it is convenient to work with components having analytical distribution functions, as is the case for the Hernquist (1990) profile, but not for the NFW model.

To make contact with common descriptions of haloes in cosmological simulations, we associate the Hernquist profile with a corresponding NFW halo with the same dark matter mass within *r*_{200}, the radius at which the mean enclosed dark matter density is 200 times the critical density. We also require that the inner density profiles are equal (i.e. ρ_{dm}=ρ_{NFW} for *r*«*r*_{200}), which implies a relation between *a* and the scalelength *r*_{s} of the NFW profile. The latter is often given in terms of a concentration index *c*, conventionally defined as *c*=*r*_{200}/*r*_{s}, where *r*_{s} is the scalelength of the NFW halo. We then have the relation

*c*= 10, this gives

*a*≃ 1.73

*r*

_{s}. The Hernquist profile then contains 75 per cent of its total mass within

*r*

_{200}and 99 per cent within 1.1

*r*

_{200}.

In Fig. 1, we compare the density profile of an NFW halo with a Hernquist model of the same concentration. In the inner regions, the two mass distributions match closely, while at large radii the density profiles asymptote to ρ_{NFW}∝*r*^{−3} and ρ_{Hern}∝*r*^{−4}, respectively. At present, observations indicate that either form can give a good description of the density profiles of actual haloes when extended beyond the virial radius (Rines et al. 2000, 2002, 2003, 2004).

#### Discs and bulges

We model disc components of gas and stars with an exponential surface density profile of scalelength *h*; i.e.

*M*

_{d}= (

*M*

_{gas}+

*M*

_{⋆}) =

*m*

_{d}

*M*

_{tot}, where

*m*

_{d}is dimensionless and

*M*

_{tot}is the total mass of the galaxy.

The bulge is taken to be spherical, for simplicity. We also model it with a Hernquist profile:

We treat the bulge scalelength*b*as a free parameter that we parametrize in units of the disc scalelength and specify the bulge mass as a fraction

*m*

_{b}of the total mass, i.e.

*M*

_{b}=

*m*

_{b}

*M*

_{tot}.

We set the disc scalelength *h* by relating it to the angular momentum of the disc. Following Mo, Mao & White (1998) and Springel & White (1999), we first write the total halo angular momentum of the halo as

*J*

_{d}=

*m*

_{d}

*J*, corresponding to conservation of specific angular momentum of the material that forms the disc. Assuming that the disc is centrifugally supported, this then implies a one-to-one relation between λ and

*h*. In order to obtain this connection, we compute the disc angular momentum assuming strict centrifugal support of the disc and negligible disc thickness compared with its scalelength. Then we can write where the circular velocity is given by Here

*y*=

*R*/(2

*h*), and the

*I*and

_{n}*K*are Bessel functions. Note that we will later drop the thin-disc approximation for the potential of the disc in favour of an accurate representation of a thick disc. However, for the conversion of the free parameter λ into a disc scalelength

_{n}*h*, the accuracy of equation (9) is sufficient.

We specify the vertical mass distribution of the stars in the disc by giving it the profile of an isothermal sheet with radially constant scalelength *z*_{0}. The three-dimensional (3D) stellar density in the disc is hence given by

*z*

_{0}as a free parameter that effectively determines the ‘temperature’ of the disc and set the velocity distribution of the stars such that this scaleheight is self-consistently maintained in the full 3D potential of the galaxy model.

However, a similar freedom is not available for the gas, because once cooling and star formation processes are taken into account, we cannot choose the temperature freely. Rather, in a broad class of models the gas will stay close to an (effective) EOS of the form *P*=*P*(ρ), implying a tight relation of gas pressure and gas density. Examples for such situations include isothermal gas, or the subresolution multiphase model by SH03. For a given surface density, the vertical structure of the gas disc then arises as a result of self-gravity and the pressure given by this EOS, leaving no freedom for prescribing a certain vertical profile. In this situation, the vertical structure has to be computed self-consistently, as discussed next.

#### Vertical structure of the gas disc

We assume that the vertical structure of the gas disc in our axisymmetric galaxy models is governed by hydrostatic equilibrium, i.e.

where F is the total gravitational potential of all mass components. This equation can be rewritten as where γ= (*d*ln

*P*)/(

*d*ln ρ) is the local polytropic index of the EOS. Note that for a given potential F, the solution of this equation is determined by the integral constraint where S

_{gas}(

*r*) is the surface mass density we prescribed in equation (3). Assuming a given potential for the moment, we solve equations (12) and (13) by integrating the differential equation for a chosen central density value in the

*z*= 0 mid-plane. This chosen starting value is then adjusted until we recover the desired surface density for the integrated vertical structure of the gas layer. We carry out this process as a function of radius, using a fine logarithmic grid of points in the

*R*—

*z*plane to represent the axisymmetric gas density distribution.

### Evaluating the potential

The above invokes the problem of determining the potential and the resulting gas distribution self-consistently. We address this problem with an iterative method. Starting with a guess for an initial potential, we compute the implied vertical gas structure as described above, treating the potential as fixed. Then, we recompute the potential for the updated gas distribution and repeat the whole process. In doing this, the potential and the gas distribution converge quickly, within a few iterations, to a self-consistent solution.

However, in doing this, we also need a method to accurately compute the potential and the resulting force field for the radially varying density stratification of the gas. Because we want to be able to set up quite massive gas discs that are in equilibrium right from the start, we cannot utilize simple thin-disc approximations, but rather need to compute the full potential accurately. The generality of the gas distributions used here makes this a non-trivial exercise.

We use a somewhat contrived, yet extremely flexible solution for this task. We set up a discretized mass distribution that represents the desired disc components accurately and then compute the potential numerically with a hierarchical multipole expansion based on a tree code. For a given (current) mass model of the galaxy, we use for this purpose a suitably distorted grid of particles of typically 2048 × 64 × 64 particles per component (the different numbers refer to radial, azimuthal and vertical directions, respectively). While we find that this gives sufficient accuracy for our purposes, we note that the number of particles can be easily made much larger to generate a still smoother potential, if needed. Unlike in a normal *N*-body code, here we are not interested in evaluating the forces or the values of the potential for all these particles. We only use them as markers for the mass and instead compute the potential at a fixed set of spatial coordinates. The work required for this potential computation then scales only with log(*N*). The primary cost of making *N* large lies therefore in the memory consumed on the computer and not in the CPU time. The potential and force field of the spherical mass distributions used for the dark matter halo and the bulge are known analytically, so we simply add them to the result obtained with the tree for the disc components.

### Velocity structure

Once the full density distribution has been determined, we compute approximations for the velocity structure of the collisionless and gaseous components. For the dark matter and stars in the bulge, we assume that the velocity distribution function only depends on energy *E* and the *L*_{z} component of the angular momentum. Then, mixed second-order moments of the velocity distribution vanish, 〈*v*_{R}*v*_{z}〉=〈*v*_{z}*v*_{f}〉=〈*v*_{R}*v*_{f}〉= 0, as well as the first moments in radial and vertical directions, 〈*v _{R}*〉=〈

*v*

_{z}〉= 0. The velocity distribution can then be approximated as a triaxial Gaussian with axes aligned with the axisymmetric coordinate system. The non-vanishing second moments can be obtained with the Jeans equations. We have

In the azimuthal direction, there can be a mean streaming component 〈*v*_{f}〉, which is not determined by the Jeans equations. Once it is specified, the dispersion of the Gaussian velocity distribution in the azimuthal direction is given by

*v*

_{f}〉 to zero, meaning that the bulge is assumed to have no net rotation. For the dark matter halo, we assume that it has the same specific angular momentum as the disc material. For simplicity, we impart this angular momentum by making 〈

*v*

_{f}〉 a fixed fraction of the circular velocity, i.e. 〈

*v*

_{f}〉=

*f*

_{s}

*v*

_{c}(Springel & White 1999). Note that we then have

*f*

_{s}« 1, i.e. the net streaming of the resulting dark halo is quite small.

The velocity structure of the stellar disc can in principle be much more complicated, and it will in general have a distribution function that does not only depend on *E* and *L*_{z}. For simplicity, we however continue to approximate the velocity distribution with a triaxial Gaussian that is aligned with the axisymmetric coordinated vectors.

As Hernquist (1993a) points out, observationally there is good evidence that 〈*v*^{2}_{R}〉 is proportional to 〈*v*^{2}_{z}〉 for the stellar disc. We hence assume σ^{2}_{R}=〈*v*^{2}_{R}〉=*f _{R}*〈

*v*

^{2}

_{z}〉. Note that larger values of

*f*will increase the Toomre

_{R}*Q*parameter, making the disc more stable against axisymmetric perturbations. In most of our default models, we actually set

*f*= 1 which corresponds to the approximations made for the dark halo and the bulge. Note that the Toomre

_{R}*Q*also varies in response to the vertical dispersion 〈

*v*

^{2}

_{z}〉, which sensitively depends on the assumed disc thickness.

To specify the mean streaming, we employ the epicyclic approximation, which relates the radial and azimuthal dispersions in the stellar disc:

where Using equation (15) for 〈*v*

^{2}

_{f}〉, we then set the streaming velocity for the stellar disc as This completes the specification of the velocity structure of the collisionless components.

For the gas, we deal with a single valued velocity field, where only the azimuthal streaming velocity has to be specified. The latter is determined by the radial balance between gravity on one hand, and centrifugal and pressure support on the other hand. Specifically, we here have

Note that the gas temperature at each coordinate is given by the EOS based on the self-consistent mass distribution derived above.We carry out the integrations over the force field needed in the Jeans equations with the help of a fine logarithmic grid in the *R*—*z* plane. This grid is also used to tabulate the final velocity dispersion fields. The differentiations needed in equation (15), for example, are approximated by finite-differencing off this grid. Once all of the density and velocity distributions have been computed, we initialize particle coordinates and velocities by drawing randomly from the respective distributions. Values for the velocity structure at individual particle coordinates are obtained by bilinear interpolation of the *R*—*z* grid.

Kazantzidis, Magorrian & Moore (2004) recently examined the accuracy of initializing spherically symmetric dark matter haloes with a Gaussian velocity distribution with dispersion derived from the Jeans equations. They found that the innermost parts of haloes are not precisely in equilibrium when this approximation is applied. As a result, the core relaxes within the first few crossing times to a density profile, which lies slightly below the input profile in the centre, such that the central cusp becomes softer. Using a more accurate computation of the velocity distribution function that takes higher order moments into account, they also demonstrated that this behaviour can in principle be avoided.

In Fig. 2, we show an example of this effect. To illustrate this in a manner unaffected by two-body relaxation, we have set up a pure dark matter halo with four million particles using our code for constructing compound galaxies. When evolved forward in time, the density profile in the very centre fluctuates below the desired input value during the first crossing times. Afterwards, the density profile shows no further secular evolution and stays extremely stable for times of order the Hubble time. The magnitude of the initial relaxation effect appears to be consistent with what Kazantzidis et al. (2004) found. However, we argue that this perturbation is acceptably small for our purposes. This is because the effect on the central dark matter profile is really quite moderate and occurs in a radial region that is already dominated by baryonic physics when a disc component is included as well. In fact, our disc galaxy models show remarkably little secular relaxation when started from initial conditions constructed in the above fashion, even for high gas fractions. We here benefit from the included treatment of gas pressure forces and the improved potential computation compared with previous methods (Hernquist 1993a; Springel 2000).

### Galaxy model parameters

In summary, we specify a disc galaxy model with the following parameters.

- (i)
The total mass is given in terms of a ‘virial velocity’,

*v*_{200}. We set*M*_{tot}=*v*^{3}_{200}/(10*G H*_{0}). - (ii)
The total mass of the disc is given in terms of a dimensionless fraction

*m*_{d}of this mass, i.e.*M*_{disc}=*m*_{d}*M*_{tot}. Likewise, the mass of the bulge is given in terms of a dimensionless fraction*m*_{b}of the total, i.e.*M*_{bulge}=*m*_{b}*M*_{tot}. The rest of the total mass is in the dark matter halo. - (iii)
A fraction

*f*_{gas}of the disc is assumed to be in gaseous form, the rest in stellar form. The bulge is taken to be completely stellar. - (iv)
The stellar disc is assumed to have an exponential profile with radially constant vertical scaleheight

*z*_{0}. We specify the latter in units of the radial disc scalelength and typically use*z*_{0}≃ 0.1–0.2*h*. - (v)
Similarly, the bulge scalelength is specified as a fraction of the radial disc scalelength.

- (vi)
A value for the spin parameter λ of the halo is selected. Its value effectively determines the radial scalelength of the disc.

In Fig. 3, we show rotation curves for two galaxy models that we use in many of the numerical simulations in this study. The two models have a mass comparable to the Milky Way, *M*_{tot}= 0.98 × 10^{12}*h*^{−1} M_{⊙}, corresponding to a virial velocity *v*_{200}= 160 km s^{−1}. Their disc mass is specified by *m*_{d}= 0.041. The model that includes a bulge (shown on the right) has a rotation curve very similar to the one used by Hernquist (1993a). Note that the model on the left, with its slowly rising rotation curve due to the absence of a bulge (a similar model was used, for example, by Mihos & Hernquist 1994b), is more susceptible to tidal perturbations.

## Star Formation and Stellar Feedback

Current state of the art simulations lack the dynamic range needed to follow the detailed structure of star-forming gas on scales characterizing entire galaxies. For example, in simulations like those presented in the following sections of this paper, tens or hundreds of thousands of particles are used to characterize the gas in individual disc galaxies, so that each particle represents roughly ≳10^{5} M_{⊙} of interstellar material for typical disc gas fractions. However, in order to resolve bound structures in SPH, an object must contain at least a number of particles comparable to the number averaged over when forming smoothed quantities. Here, we set this parameter typically to be ≃64 neighbouring particles, so that the simulations cannot resolve bound objects with masses smaller than about ∼6 × 10^{6} M_{⊙}. Simulations with total particle numbers 100 times larger than those presented here are feasible, but even then we would not be able to resolve self-gravitating structures with masses ≲10^{5} M_{⊙}. Furthermore, such calculations would be limited to a spatial resolution of order ∼20 pc. These mass and spatial scales are barely sufficient to resolve individual giant molecular clouds and are inadequate to describe the details of the star-forming gas contained within them.

More seriously, there is at present no complete theory for star formation, so the important physics is not well understood. Even if the simulations could include the entire dynamic range needed to resolve the essentials of star formation on galactic scales, the description of this process would therefore still have to be parametrized using bulk averaged quantities. For these reasons, we are motivated to adopt an approach in which we discard detailed information about star-forming gas in favour of a procedure that captures generic consequences of star formation and associated feedback effects on resolved scales. In this manner, we do not attempt to characterize all aspects of the physics but instead seek to identify general trends.

As an illustration of our procedure, we employ the subresolution multiphase model for star-forming gas developed by SH03. In this picture, the ISM is pictured to consist of cold clouds where stars form, embedded in a hot, pressure-confining phase. The gas is assumed to be thermally unstable to the onset of a two-phase medium at densities above a threshold ρ_{th}. The fraction of mass in the two phases at higher densities is set by star formation and feedback, evaporation of the cold clouds through thermal conduction, and the growth of clouds owing to radiative cooling. The rates at which the mass in the two phases evolves are given by equations (5) and (6) in SH03. The thermal budget of the cold clouds and hot gas is determined by the same effects and evolves according to equations (8) and (9) in SH03.

Our coarse-grained representation of this model reduces largely to two principle ingredients: (i) a star formation prescription or ‘law’; and (ii) an effective EOS for the ISM. For the former, we adopt a rate motivated by observations (Kennicutt 1989, 1998):

where dρ_{*}/d

*t*is the star formation rate, β depends on the initial mass function (IMF) and is the mass fraction of stars that die immediately as supernovae, ρ

_{c}is the density of cold clouds, and

*t*

_{*}is a characteristic time-scale. For a Salpeter type IMF, we have β≈ 0.1. In our multiphase picture, ρ

_{c}≈ρ, where ρ is the total gas density (see Fig. 1 in SH03). A good match to observations is obtained if

*t*

_{*}is assumed to be proportional to the local dynamical time where

*t*

^{0}

_{*}is a constant parameter.

The effective EOS of the gas in this model is given by

where ρ_{h}and

*u*

_{h}are the density and specific thermal energy of the hot phase, with similar quantities defined for the cold clouds, while γ= 5/3 is the adiabatic index of the gas. In detail, the EOS depends on the values of the parameters in our model.

As described by SH03, the parameter values are constrained by the nature of the cooling curve for the gas and by requiring that the EOS be a continuous function of density. SH03 also showed that if star formation is rapid compared with adiabatic heating or cooling owing to the motion of the gas, then the multiphase picture leads to a cycle of self-regulated star formation where, in equilibrium, the growth of cold clouds is balanced by their evaporation from supernova feedback. Under these conditions and for a given IMF, the EOS is determined by three parameters: the constant appearing in the star formation law, *t*^{0}_{*}, a normalization of the cloud evaporation rate, *A*_{0}, and a supernova ‘temperature’, *T*_{SN}, that reflects the heating rate from a population of supernovae for the adopted IMF. The requirement that the gas be thermally unstable at densities above ρ_{th} fixes the ratio *T*_{SN}/*A*_{0}≈ 10^{5}, while in the self-regulated regime, the condition that the EOS be continuous implies that ρ_{th} depends mainly on the ratios *T*_{SN}/*A*_{0} and *T*_{SN}/*t*^{0}_{*} (see equation 23 of SH03).

In Fig. 4, we show an example of the form of our effective EOS for the parameter choices *T*_{SN}= 10^{8}, *A*_{0}= 1000 and *t*^{0}_{*}= 2.1 Gyr. For an isolated disc galaxy, these values reproduce the observed correlation between star formation and gas density (see figs 2–3 of SH03). For densities higher than ρ_{th}, the EOS shown in Fig. 4 is fitted to an accuracy of 1 per cent by (Robertson et al. 2004)

*P*

_{eff}, would increase linearly with density, as indicated by the dashed line for densities higher than ρ

_{th}. However, when we account for supernova feedback in our multiphase model, the EOS in star-forming gas is stiffer (solid curve in Fig. 4). Suppose a region at densities lower than ρ

_{th}is compressed. If the gas were isothermal, the effective pressure would increase proportional to the density. However, in our multiphase model, once ρ > ρ

_{th}, the gas is thermally unstable and stars can form. As the density rises further, the star formation rate increases even faster, as does the rate of energy injection into the gas from supernova. If this feedback energy is retained by the gas, its effective EOS will be stiffer than isothermal.

The concept of an ‘effective EOS’ thus incorporates the impact of local feedback in regulating the thermodynamic properties of the ISM. Using this concept, we can describe the dynamics of star-forming gas on galactic scales and account for the consequences of stellar feedback on large scales. For example, the dynamical stability of a pure gas disc to axisymmetric perturbations is determined by the condition (Toomre 1964)

where*c*

_{s}is the sound speed, κ is the epicyclic frequency and S is the gas surface density. This relation represents a competition between the stabilizing influences of pressure (

*c*

_{s}) and angular momentum (κ) against the destabilizing effect of gravity (S). If we consider a small region in the disc, the appropriate EOS with which to compute

*c*

_{s}is not one that captures the detailed structure of the ISM, but an effective EOS on the scale of the instability. For parameters typical of galactic discs, the fastest growing mode has a wavelength that is comparable to the size of the disc (Binney & Tremaine 1987) and so stability is determined by the response of the ISM to compressions on this scale.

In the absence of a complete theory for the ISM, the choices for the star formation law and effective EOS are not unique. For example, Barnes (2004) has recently proposed a star formation prescription that is based on the rate at which gas generates entropy in shocks. His numerical tests indicate that this formalism may yield a better match to observations of interacting galaxies than a purely density-dependent star formation law.

In our approach, we will consider the macroscopic star formation law and effective EOS as elements of our description that can in principle be varied independently to isolate the large-scale physics that regulates star formation in colliding and merging galaxies. In particular, we will also consider simulations where we ‘soften’ the EOS by linearly interpolating between that for the standard implementation of our multiphase model and an isothermal EOS. We will denote this softening by a parameter *q*_{EOS} so that *q*_{EOS}= 1 will correspond to, for example, the ‘stiff’ EOS shown by the solid curve in Fig. 4 and *q*_{EOS}= 0 will correspond to the ‘soft’, isothermal EOS, indicated by the dashed curve in Fig. 4.

While we do not expect that the effective EOS employed here is a precise description of real star-forming gas, the approach that we adopt is flexible and can be applied to any well-specified theory for the ISM. In particular, Springel & Hernquist (2003b) have shown that our coarse-graining procedure leads to a numerically converged estimate for the cosmic star formation history of the Universe (Hernquist & Springel 2003) that agrees well with low redshift observations.

The hybrid multiphase model of SH03 also includes a phenomenological representation of galactic winds. In this picture, some of the energy available from supernova feedback is tapped to drive an outflow if the star formation rate exceeds a critical threshold. For simplicity, we do not include this additional mode of feedback in the present study.

## Black Hole Accretion and agn Feedback

It is now widely believed that BH growth and associated feedback energy from this process may be important for a variety of phenomena related to the evolution of ordinary galaxies as well as unusual behaviour in quasars, starbursts and AGN. It has been suggested that the *M*_{BH}—σ relation may arise if strong outflows are produced in response to a major phase of BH accretion, which via their interaction with the surrounding gas would inhibit any further accretion and hence BH growth (Silk & Rees 1998; Fabian 1999; King 2003; Wyithe & Loeb 2003). Indeed, X-ray observations of a number of quasars (mostly broad absorption line systems) reveal significant absorption, implying large outflows with a kinetic power corresponding to a significant fraction of the AGN bolometric luminosity (Chartas, Brandt & Gallagher 2003; Crenshaw, Kraemer & George 2003; Pounds et al. 2003). In the case of radio-loud quasi-stellar objects (QSOs), there is also evidence that up to half of the total power is injected in the form of jets (Rawlings & Saunders 1991; Tavecchio et al. 2000). Inevitably, such outflows must have a strong impact on the host galaxies. Models of interacting and merging galaxies should, therefore, account for feedback from both star formation and accretion onto central, supermassive BHs. Note that BH accretion is expected to be responsible for the majority of the BH mass growth and to provide sufficient energy supply for driving associated outflows.

As with star formation, current numerical simulations cannot resolve the properties of the accretion flow around nuclear BHs on galactic scales. For example, consider a BH of mass *M*_{BH} accreting spherically from a stationary, uniform distribution of gas whose sound speed at infinity is *c*_{8}. The gravitational radius of influence of the BH is then (Bondi 1952)

*M*

_{BH}= 10

^{7}M

_{⊙}in a gas with a sound speed

*c*

_{8}= 30 km s

^{−1}, this is numerically The Schwarzschild radius of a BH of this mass is Our understanding of the nature of accretion onto supermassive BHs is sufficiently poor that it is not clear what range in spatial scales would be required to obtain an accurate description of the impact of BH growth and feedback on galactic scales. However, given the relatively poor resolution that can be achieved in simulations like those presented here, it is clear that they cannot represent the full complexities of this process in any detail. By analogy with out treatment of star formation and supernova feedback, we are hence led to adopt a coarse-graining procedure.

In what follows, we use an effective subresolution model to characterize the growth of supermassive BHs in galactic nuclei and the consequences of feedback from accretion on surrounding gas. Technically, we represent BHs by collisionless particles that can grow in mass by accreting gas from their environments. A fraction of the radiative energy released by the accreted material will be assumed to couple thermally to nearby gas, and influence its motion and thermodynamic state.

Like our procedure for coarse graining the ISM, our method is flexible and can be applied to any model for BH accretion. As a starting point, we relate the (unresolved) accretion onto the BH to the large-scale (resolved) gas distribution using a Bondi—Hoyle—Lyttleton parametrization (Hoyle & Lyttleton 1939; Bondi & Hoyle 1944; Bondi 1952). In this description, the accretion rate onto the BH is given by

where ρ and*c*

_{s}are the density and sound speed of the gas, respectively, α is a dimensionless parameter, and

*v*is the velocity of the BH relative to the gas. We will also assume that the accretion is limited to the Eddington rate where

*m*

_{p}is the proton mass, σ

_{T}is the Thomson cross-section and e

_{r}is the radiative efficiency. The latter is related to the radiated luminosity,

*L*

_{r}, and accretion rate, , by i.e. it simply gives the mass to energy conversion efficiency set by the amount of energy that can be extracted from the innermost stable orbit of an accretion disc around a BH. For the rest of this study, we adopt a fixed value of e

_{r}= 0.1, which is the mean value for radiatively efficient Shakura & Sunyaev (1973) accretion onto a Schwarzschild BH. We ignore the possibility of radiatively inefficient accretion phases. The accretion rate is then

We will assume that some fraction e_{f} of the radiated luminosity *L*_{r} can couple thermally (and isotropically) to surrounding gas in the form of feedback energy, viz.

_{f}∼ 0.05, so that ∼0.5 per cent of the accreted rest mass energy is available as feedback. This value fixes the normalization of the

*M*

_{BH}—σ relation and brings it into agreement with current observations (Di Matteo, Springel & Hernquist 2005).

## Numerical Approach

Nothing in our formulation of supernova and BH feedback is tied to a particular numerical scheme for solving the dynamical equations for the gas and collisionless matter. Our subresolution models and coarse-graining procedure could be incorporated into either particle- or grid-based methods. In our simulations, we employ an *N*-body algorithm for the collisionless material, which, in our case, includes dark matter, stars and BHs. For the hydrodynamics, we employ an SPH code, which represents fluids elements by particles (Gingold & Monaghan 1977; Lucy 1977; Monaghan 1992). In SPH, fluid properties at a given point are estimated by local kernel-averaging over neighbouring particles and smoothed versions of the equations of hydrodynamics are solved for the evolution of the fluid.

The particular code we use is an improved and updated version of gadget (Springel, Yoshida & White 2001). The code implements a TreeSPH algorithm (Hernquist & Katz 1989) where gravitational forces are computed with a hierarchical tree method (Barnes & Hut 1986) and SPH is used for the hydrodynamics. In the new version gadget2 employed here, the hydrodynamical equations are solved using a fully conservative technique (Springel & Hernquist 2002), which maintains strict entropy and energy conservation even when smoothing lengths vary adaptively (Hernquist 1993b).

Besides gravity and hydrodynamics, the code follows radiative cooling processes of an optically thin primordial plasma of helium and hydrogen in the presence of a UV background (Katz, Weinberg & Hernquist 1996). Star formation and the dynamics of the highly overdense gas are treated with the multiphase model described earlier. Independent star particles are stochastically spawned out of the gas phase (SH03), thereby avoiding artificial coupling between gaseous and collisionless matter.

BHs are represented by collisionless ‘sink’ particles, which only feel gravitational forces. For these particles, we compute estimates for the gas density in their local environments, in the same fashion as it is done for normal SPH particles. Similarly, we also determine the average gas temperature in the local SPH smoothing environment around BH particles, as well as the gas bulk velocity relative to the BHs. These quantities are then used to estimate the BH accretion rates, based on the equations specified in the previous section.

To implement the actual accretion, we follow a similar stochastic approach as it is applied for regular star formation. To this end, we compute for each gas particle *j* around a BH a probability

*t*is the time-step, ρ is the gas density estimated at the position of the BH and

*w*is the kernel weight of the gas particle relative to the BH. We then draw random numbers

_{j}*x*uniformly distributed in the interval [0, 1] and compare them with the

_{j}*p*. For

_{j}*x*<

_{j}*p*, the gas particle is absorbed by the BH, including its momentum. On average, this procedure ensures that the BH particle accretes the right amount of gas consistent with the estimated accretion rate .

_{j}However, because the accretion rate depends sensitively on *M*_{BH}, this procedure would be quite noisy under conditions of poor resolution, where accreting a single gas particle can change the BH particle mass by a substantial factor. We circumvent this problem by giving the sink particle an internal degree of freedom that describes the BH mass in a smooth fashion. The value of this variable represents the BH mass for the case of ideal resolution where the gas mass is not discretized. The real dynamical mass of the sink particle tracks this internal mass closely, albeit with stochastic fluctuations around it. For this reason, we use the internal mass for computing the accretion rates, while the growth of the hole is followed both in terms of the internal mass and the actual dynamical mass. The latter increases in discrete steps when whole gas particles are absorbed, but follows the internal variable in the mean, with fluctuations that become progressively smaller for better resolution. Using this procedure, we can reliably follow the growth of BHs in terms of their internal mass variable even in haloes with just a few hundred particles, which is particularly important for cosmological simulations, where the earliest generations of galaxies are typically not very well resolved. Note that we conserve momentum when gas particles are absorbed by BH sink particles. If the BH is moving relative to the gas and has a high accretion rate, this can effectively act like a friction force, which reduces the relative motion.

Together with the accretion, we compute the rate of feedback associated with the BH growth. This energy is added kernel-weighted to the thermal reservoir of the gas in the local environment around the BH. The shape of the kernel is of the usual SPH form and we choose its size *h* in analogy to the smoothing length determination used for gas particles, i.e. the kernel size *h* is determined based on the implicit solution of the equation (4π/3) *h*^{3}ρ=*M*_{sph} (Springel & Hernquist 2002), where ρ is the kernel estimate of the gas density at the position of the BH and *M*_{sph} is the mass of ∼ 64 SPH gas particles.

We still need to specify how this energy feedback is included in our multiphase model for the star-forming gas. In the latter, we do not directly enforce the effective EOS expected for equilibrium conditions between supernova heating and radiative cooling. Rather, we keep following the thermal energy content (or equivalently the specific entropy, which is what we actually use) of gas particles as independent thermodynamic variable at all densities. If gas is dense enough to be in the multiphase state, we employ the framework of the simplified treatment of the multiphase model described in SH03 (sections 2 and 5.2), assuming that the specific thermal energy content decays on the time-scale τ_{h} given by equation (11) of SH03 back to the specific energy expected based on the effective EOS. If the local cooling time τ_{cool} (computed assuming that all gas is in the ‘hot phase’) is shorter than this relaxation time-scale, we instead use normal radiative cooling to reduce the excess energy. In this way, we can allow for local deviations of the gas pressure from the effective multiphase model, while treating the decay back to the effective EOS consistently with the expected radiative cooling losses. This formalism also allows an inclusion of the phenomenological softening of the EOS explored in the present paper. We note that if the local cooling rate of the gas is not high enough to radiate away all of the energy injected by the BH on a short time-scale, an increase of the gas sound speed occurs, which can then throttle the accretion in the Bondi-dominated regime.

In the final stages of a galaxy merger, the cores of the galaxies coalesce to form a single dark matter halo and a single stellar system. Presumably, this also means that a central binary system of two supermassive BHs is formed, which may subsequently harden and eventually lead to physical collision and merger of the BHs themselves. However, it is a controversial question how long it would take to harden the BH binary by stellar-dynamical (Makino & Funato 2004) or hydrodynamical processes (Escala et al. 2004) until finally gravitational wave emission becomes important, leading to quick orbital decay. If BH merger processes were inefficient, it should be a common situation that a third supermassive BH is brought in by the next galaxy merger. This BH could then interact with the binary and lead to sling-shot ejection of one of the holes, with the remaining pair being ejected in the opposite direction. As this may seriously impair the growth of BHs to the large masses that are observed, the existence of supermassive BHs may be taken as circumstantial evidence that BH binaries probably do merge. We therefore assume that two BH particles merge instantly if they come within the smoothing lengths of each other and their relative velocities are at the same time smaller than the local sound speed. Note that as with the accretion itself, we lack the dynamic range in simulations of whole galaxies to study the hardening process of BH binaries directly.

## Isolated disc Galaxies

We have run a series of simulations of isolated galaxies using the methods described in Section 2 to initialize them. In the following, we examine their stability using models with and without bulges, having the rotation curves shown in Fig. 3. The structure of these models was motivated by properties of galaxies like the Milky Way and to enable us to compare our results with earlier work, such as Hernquist & Mihos (1995) and Mihos & Hernquist (1996). Here, we focus on the effect of varying the stiffness of the EOS and the gas fraction in the disc. A significant advantage of our formalism is that we are now able to construct stable equilibrium models with a larger gas fraction than possible in previous studies.

### Stability of models without BHs

In order to facilitate comparison between different cases, in this section we fix the parameters associated with our multiphase model and vary only the disc gas fraction and the stiffness of the effective EOS. For definiteness, we take *t*^{0}_{*}= 8.4 Gyr, *A*_{0}= 4000 and *T*_{SN}= 4 × 10^{8}. This choice for *t*^{0}_{*} is a factor of 4 larger than that employed by SH03 to match the Kennicutt Law. However, we find that for our isolated galaxies, *t*^{0}_{*}= 2.1 Gyr yields global star formation rates ≳4 M_{⊙} yr^{−1} for gas fractions of 10 per cent and structural properties similar to the Milky Way. In our Galaxy, the inferred average star formation rate is only ∼1 M_{⊙} yr^{−1}. Increasing *t*^{0}_{*} by a factor of 4 gives therefore better agreement with the long gas consumption time-scale inferred for the Galaxy. This also simplifies a comparison with Hernquist & Mihos (1995) and Mihos & Hernquist (1996), who chose parameters in their star formation prescription to give global star formation rates ∼1 M_{⊙} yr^{−1} for galaxies like the Milky Way with 10 per cent of their disc mass in gas.

Once *t*^{0}_{*} is fixed, we select *A*_{0} and *T*_{SN} as above so that the critical density for the gas to be thermally unstable, ρ_{th}, implies a critical projected gas surface density for star formation to occur that is similar to observations (Kennicutt 1989, 1998). Because ρ_{th} in our formalism depends mainly on the ratios *T*_{SN}/*A*_{0} and *T*_{SN}/*t*^{0}_{*}, a factor of 4 lengthening of *t*^{0}_{*} relative to SH03 requires a factor of 4 increase in *T*_{SN} and *A*_{0} for a critical gas surface density ∼ 10 M_{⊙} pc^{−2} (see figs 2–3 of SH03). As we indicated earlier, the effective EOS depends primarily on ratios between these three parameters. Scaling *t*^{0}_{*}, *T*_{SN} and *A*_{0} by the same factor implies that the EOS for our isolated discs is the same as the ‘stiff’ case in Fig. 4, only the gas consumption time-scale is larger. In the stability analysis shown below, we hold these parameters fixed, so that the threshold density for star formation is the same, but we artificially vary the stiffness of the EOS by linearly interpolating between the stiff and isothermal limits in Fig. 4. EOS softenings of *q*_{EOS}= 1 and *q*_{EOS}= 0 refer to the stiff and isothermal cases, respectively. We also vary the fraction *f*_{gas} of the disc mass in gas, but keep the total disc mass constant.

As a typical example, Fig. 5 shows the time evolution of the projected gas density in an isolated galaxy model with a bulge. The gas fraction is *f*_{gas}= 0.1 and the EOS softening parameter is *q*_{EOS}= 0.5. Initially, the particles randomly sample the exponential surface density distribution, so there are no coherent structures present in the disc. However, a steady spiral pattern develops within a fraction of a rotation period owing to the action of swing amplification (Toomre 1981). This behaviour is similar to that reported in previous numerical studies of disc structure. For the parameter choices in this example, the disc is stable for many rotation periods and does not evolve strongly. The gas is converted into stars at a relatively low rate according to our star formation prescription but, otherwise, the disc, bulge and halo do not evolve.

However, strong, unstable evolution is possible if the gas fraction is large and the EOS is too soft. In Fig. 6, we compare the distribution of gas in model galaxies without bulges after approximately one rotation period of evolution. The gas fraction is increased from *f*_{gas}= 0.1 to *f*_{gas}= 0.99 along rows (left to right), while the EOS is softened from *q*_{EOS}= 1 (top) to *q*_{EOS}= 0.05 (bottom) along columns. Some interesting features and trends are apparent from Fig. 6.

First, for a relatively stiff EOS, *q*_{EOS}≳ 0.3, increasing the gas fraction actually suppresses the amplitude of spiral structure in the disc. This is because in our full multiphase model, the effective temperature of the gas is *T*_{eff}≳ 10^{5} K over most of the disc for *q*= 1 and so the gas is dynamically hotter than the old stars. Increasing the gas fraction in this case makes the disc more stable and less susceptible to amplifying non-axisymmetric patterns. Thus, in the top row of Fig. 6, the spiral structure that is apparent for *f*_{gas}= 0.1 (top left) is mostly washed out for *f*_{gas}= 0.99 (top right). This tendency does not occur if the EOS is soft because then the effective temperature is such that the gas is dynamically colder than the old stars. For *q*_{EOS}= 0.125 (bottom two rows), this results in unstable behaviour if the gas fraction is large.

Secondly, the set of models shown in Fig. 6 clearly delineates stable cases from unstable ones. For our full, stiff EOS (*q*_{EOS}= 1) even pure gas discs are stable (top row). This is also true for intermediate cases with *q*_{EOS}= 0.5. However, if the effective temperature of the gas is lower than that of the stars, very gas-rich discs are unstable and fragment within a rotation period. Thus, for example, with a nearly isothermal EOS (*q*_{EOS}= 0.05, bottom row), discs with a gas fraction *f*_{gas}≳ 0.3 are already unstable (bottom row).

The previous studies by Hernquist & Mihos (1995) and Mihos & Hernquist (1996) treated the ISM as an isothermal gas and included only weak feedback in the form of a small input of random kinetic energy. These authors were unable to construct stable models with large gas fractions owing to disc fragmentation like that shown in Fig. 6. Supernova feedback in our multiphase model pressurizes the gas making the effective EOS stiffer, stabilizing the disc.

A similar separation between stable and unstable behaviour is seen in our galaxy models that include bulges. We illustrate this in another manner in Figs 7 and 8 where we show the evolution of the global star formation rates in the discs for models without and with bulges, respectively.

The behaviour in Figs 7 and 8 clearly separates stable and unstable models. In unstable discs, fragmentation leads to a runaway growth in density perturbations, yielding unsteady evolution in the global star formation rate. Thus, all our galaxies are stable for stiff EOSs with *q*_{EOS}= 0.5. Models with a bulge are stable for all gas fractions with *q*_{EOS}= 0.25, however those without a bulge (Fig. 7) are unstable with this EOS for *f*_{gas}≳ 0.6. Softer EOSs yield unstable behaviour for galaxies without bulges for *f*_{gas}≳ 0.4. For models with bulges, instability sets in at *f*_{gas} > 0.8 for *q*_{EOS}= 0.125 and *f*_{gas}≳ 0.6 for *q*_{EOS}= 0.05. The scalings with *q*_{EOS} and *f*_{gas} are broadly consistent with the Toomre criterion.

### Models with BHs

We have also run models similar to those described in Section 6.1, but where we added a sink particle to represent a BH and allowed it to accrete gas from the galaxy, increasing its mass and adding a source of feedback energy in addition to supernovae. The examples we discuss in this section are not intended to provide a physical picture for the growth of supermassive BHs in isolated galaxies because we have made no attempt to tie the origin of the BHs to galaxy formation. Nevertheless, we can use these simulations to identify some generic features of our description of the interaction between BHs and the ISM that are relevant for dynamically evolving situations, like galaxy mergers.

The coupling between BH accretion and surrounding gas has several implications. First, heating from AGN feedback energy can drive outflows from galaxies if the accretion rate is sufficiently high. This is shown in Fig. 9, where the effect of BH feedback is illustrated for an isolated galaxy with a bulge. For the low accretion rates relevant to the times indicated, the deposition of thermal energy into the gas in the vicinity of the BH can drive a weak wind perpendicular to the plane of the disc, as indicated by the velocity vectors in this figure. These outflows are reminiscent of those in our supernova-driven winds (figs 5–8 in SH03). More powerful outflows can be produced during mergers when gas is driven into the centre of the remnant, fuelling strong nuclear accretion.

As implied in Section 4, our description of BH accretion yields several distinct phases of evolution, depending on the properties of the gas near the BH. For accretion rates lower than Eddington, the BH mass will evolve according to

where*M*

_{0}is the initial BH mass and This result is valid only when the properties of the gas near the BH and the relative velocity do not evolve significantly with time. When the accretion rate is higher than Eddington, the BH grows exponentially with time, where the Salpeter time, depends only on physical constants and on the radiative efficiency. For a radiative efficiency of 10 per cent,

*t*

_{S}≈ 4.5 × 10

^{7}yr.

In the absence of feedback effects, a low-mass BH would grow first according to equation (36). Once the BH becomes sufficiently massive, given the properties of the surrounding gas, it then enters a period of exponential growth, described by equation (38). Note that the first, Bondi-limited growth can be very slow if the initial BH seed mass is small, or if a small value for the coefficient α is selected. We typically start our calculations with BH seed masses of order 10^{5} M_{⊙}, placed with zero velocity at the centre of a galaxy model, and set α large enough to allow a hole in a gas-rich environment to reach the Eddington regime in about ∼ 0.5 Gyr.

We illustrate this behaviour in Fig. 10, which shows the growth of BHs in isolated, gas-rich galaxies, with and without the impact of accretion feedback on the gas. If feedback is neglected, the BH grows in a manner consistent with the expectations noted above, first according to equation (36) and later exponentially. When the effects of feedback are included, however, the growth is self-regulated and leads to a saturated final state, which depends on the initial gas content of the galaxy. The self-regulation occurs because, as the accretion rate increases, the amount of energy available to affect nearby gas also grows. If this material is then heated and dispersed, the supply of gas to grow the BH will decline as accretion is shut off.

## Major Mergers

The approach we have adopted for handling various forms of feedback makes it possible to explore physical effects that were not accessible previously. For example, Hernquist & Mihos (1995) and Mihos & Hernquist (1996) were able to model galaxies with only small gas fractions because their treatment of feedback was insufficient to stabilize highly gas-rich discs. With our method, we can examine mergers between galaxies with a very large gas content, as might be relevant to objects at high redshifts. Also, by including processes related to BH growth and accretion, we can study the interplay between stellar and AGN feedback that is likely to be important for understanding the nature of ultraluminous infrared galaxies (ULIRGs).

In what follows, we illustrate some of the consequences of our models for stellar and AGN feedback by looking at major mergers of disc galaxies. A representative case is given in Fig. 11, which shows the distribution of gas during the merger of two equal-mass disc galaxies from a prograde, parabolic orbit. As the discs pass by one another for the first time (*t*= 0.2–0.3), strong tidal forces yield extended tails and bridges. Gas is shocked at the interface between the two galaxies and the tidal response drives gas into the central regions of the two galaxies. When the galaxies finally coalesce (*t*= 1.2–1.5), much of the remaining gas is either converted into stars during an intense burst, or shock heated to temperatures characteristic of the virial temperature of the remnant.

### Mergers without black holes

#### Comparison with previous work

One of our goals is to see to what extent the strength of stellar feedback influences starbursts during mergers. As a starting point, we compare our results to those obtained in earlier studies, in particular, to the work of Hernquist & Mihos (1995) and Mihos & Hernquist (1996). These authors employed a highly simplified treatment of star formation by modelling the ISM as an isothermal gas with a temperature *T*_{ISM}= 10^{4} K, and including a weak form of supernova feedback by injecting small amounts of kinetic energy into the gas. The amplitude of this feedback was adjusted so that isolated discs forming stars quiescently had gas profiles similar to those observed.

We have modified our simulation code to mimic the star formation and feedback algorithm of Mihos & Hernquist (1994b) in an attempt to reproduce their results. In Fig. 12, we show the evolution of the star formation rate in prograde mergers of equal-mass disc galaxies with and without bulges. We have selected the parameters and galaxy models to match the results of Mihos & Hernquist (1996) as closely as possible. The results in Fig. 12 can be compared directly to, for example, fig. 5(a) in Mihos & Hernquist (1996).

The morphologies of the merging galaxies roughly follow the time sequence in Fig. 11, with only a modest dependence on their structural properties. When the galaxies first pass by one another, relatively strong starbursts are triggered in each disc in the case when bulges are not present. Galaxies with bulges experience much weaker starbursts during first passage. Mihos & Hernquist (1996) showed that this difference arises from the relative ease with which the discs can amplify *m*= 2 disturbances in response to tidal forcing. During final coalescence, a much stronger starburst is excited in the models with bulges because much of the gas in the bulgeless galaxies has already been consumed by this time. These trends agree well with those identified by Mihos & Hernquist (1996).

The amplitudes of the starbursts shown in Fig. 12 are comparable to, but somewhat stronger than, those shown in Fig. 5(a) of Mihos & Hernquist (1996) and considerably stronger than those found by Barnes (2004, see, for example, his fig. 5). By varying the parameters of our simulations, we have found that the evolution of the star formation rate in mergers with only relatively weak feedback is sensitive to numerical resolution. If the star formation rate is tied to gas density, the amplitudes of merger-induced starbursts depend on the compressibility of the gas, which is influenced by both the stiffness of the EOS, as well as the dynamic range in resolution of the numerical algorithm. We have verified that differences in resolution are mainly responsible for the residual discrepancies between the star formation rates shown in Fig. 12, and those reported by Mihos & Hernquist (1996) and Barnes (2004).

#### Multiphase treatment of ISM

As a next step, we have performed a grid of merger simulations using the multiphase treatment of supernova feedback described in Section 3. We have varied the structural properties of the galaxies, the orbit of the mergers and the parameters in our multiphase model. Below, we limit the discussion mainly to the consequences of variations in the fraction of gas in each galaxy and the strength of feedback to highlight the new types of outcomes that are possible.

In Fig. 13, we show the evolution of the star formation rate from one of our simulations, a prograde, major merger of two disc galaxies without bulges. In this case, each disc initially consisted of *f*_{gas}= 99.9 per cent gas and only 0.1 per cent of old stars. The EOS was relatively stiff with a softening parameter of *q*_{EOS}= 0.5. (That is, the EOS was linearly interpolated to lie midway between the curves in Fig. 4 for ρ > ρ_{th}.) Even with this extreme gas content, the model discs are stable in isolation, owing to the pressurization of the gas from supernova feedback. The parameters of this simulation are deliberately chosen to be extreme to simplify the discussion, but qualitatively similar results are obtained in more general situations.

In the example shown in Fig. 13, the star formation rate peaks at roughly 500 M_{⊙} yr^{−1}, considerably higher than in earlier simulations of equal-mass mergers of discs, owing to the large gas content of each galaxy. Star formation rates at this level are similar to those inferred for systems at high redshift such as Lyman-break galaxies and SCUBA sources, suggesting that some of these objects may result from mergers of gas-rich discs.

Fig. 13 also demonstrates that while the general conclusions of Mihos & Hernquist (1996) are obtained under broader circumstances, some of their results clearly depend on the strength of feedback from star formation. In models with weak supernova feedback, as in Fig. 12, galaxies without bulges experience stronger starbursts at first passage than during final coalescence. The evolution shown in Fig. 13 is rather different, with the strongest starburst accompanying the final merger. While the merging galaxies in this example still develop a bar-like structure in response to tidal forcing, the gas is sufficiently hot dynamically that it is not strongly compressed during this phase of evolution, limiting the strength of the initial starburst. The more intense starburst during the final merger results from strong shock compression of the gas and gravitational torquing when the galaxies coalesce, in a manner reminiscent of the mergers with bulges examined by Barnes & Hernquist (1991, 1996) and Mihos & Hernquist (1996).

An interesting aspect of the simulation used in Fig. 13 concerns the structure of the merger remnant. A considerable amount of gas is left behind following the merger. Much of it settles into an extended disc, surrounding a relatively compact, rapidly rotating remnant of newly formed stars. Overall, the remnant is morphologically closer to a spiral galaxy with a bulge than to an elliptical. Further implications and a detailed analysis of this simulation are presented in Springel & Hernquist (2005).

### Mergers with black holes

An even greater diversity of outcomes results during a merger if the impact of central, supermassive BHs is included in addition to feedback from star formation. As we discuss in Di Matteo et al. (2005) and Springel, Di Matteo & Hernquist (2005), mergers of galaxies that involve star formation, supernova feedback, BH growth and accretion, and AGN feedback yield a generic time history in which the tidal interaction drives gas into the centre of the remnant, fuelling a starburst and triggering the rapid growth of the central BH(s). Feedback from accretion onto the BH(s) influences the thermodynamic properties of the surrounding gas, eventually expelling it from the inner regions of the remnant, terminating the episode of rapid star formation and self-limiting the growth of the BH(s). The starburst and AGN activity are slightly offset in time owing to the detailed form of the response of the gas to the feedback. Thus, while we expect starbursts and AGN activity to be correlated in this picture, the remnant will be primarily seen in different stages as the merger progresses and will evolve from one type of object into the other, depending on the observed wavelength regime.

In Fig. 14, we show three snapshots from a simulation that includes our model for the growth of BHs. In this example, the galaxies included bulges and the discs were 10 per cent gas. Following the first passage (shown on the left), but before the galaxies coalesce, the discs are distorted by their mutual tidal interaction, but only a relatively weak starburst is triggered, because the bulges stabilize the discs against a strong *m*= 2 response. The BHs do not accrete significantly during this phase of evolution and so neither galaxy would be observable as an AGN. When the galaxies merge but before the gas has been consumed (middle), tidal forces trigger a nuclear starburst and fuel rapid growth of the BHs. The peak in the starburst occurs near the time indicated in Fig. 14, but the peak in the BH growth is offset, owing to the delayed action of AGN feedback on the gas. In principle, the system should now be both a starburst and an AGN, but it is likely that the quasar activity would be obscured by surrounding gas and dust, at least during the beginning of the burst. Only during the final stages of evolution could the remnant be visible as an AGN, when outflows remove the dense layers of gas around it.

As the remnant settles into a relaxed state, star formation and BH accretion are rapidly quenched by the expulsion of gas from the centre, as a consequence of AGN feedback. The remnant would no longer be observable as either a starburst or an AGN and will age rapidly, quickly resembling an evolved red galaxy (Springel et al. 2005). At this point, the BH growth ceases and the BH mass attained is consistent with those expected from the observed *M*_{BH}—σ relation (Di Matteo et al. 2005).

The interplay between BH growth and the physics of the ISM in simulations like that shown in Fig. 14 has significant implications for the observed properties of these systems. In Fig. 15, we show the evolution of the intrinsic star formation rates in simulations of galaxy mergers with and without bulges and including or neglecting central BHs. This figure shows that the presence of the BHs and the feedback energy derived from their growth can dramatically alter the nature of the starburst resulting from a merger. The peak amplitude of the starburst in the merger that included bulges is lowered by about a factor of 5, owing to AGN feedback. We note, however, that the observed emission from such an object would include both the radiation from the starburst as well as the energy output from BH accretion. If most of the AGN emission is reprocessed into optical or infrared radiation by surrounding gas and dust before it escapes from the inner regions of the remnant, it could be incorrectly attributed to the starburst, rather than being associated with a buried quasar. A more complete understanding of the implications of these results will require a detailed modelling of the spectral energy distribution from the combination of the starburst and the obscured AGN.

We have varied the orbital geometry and structural properties of the merging galaxies to see if the results described above are generic. While the detailed outcome does depend on, for example, the orbit, many of the overall properties are relatively unaffected. An example is shown in Fig. 16 where the evolution of the total BH mass is given for six different orbital configurations. The top three panels show prograde—prograde, prograde—retrograde and retrograde—retrograde mergers (left to right), while the bottom three panels show cases where the disc orientations were chosen randomly. While the detailed response varies from model to model, the final total BH mass is insensitive to the orbital parameters. In all cases, the final BH mass is ∼2–3 × 10^{7} M_{⊙}. This value is consistent with the observed correlation between BH and spheroid mass for the remnants produced in these mergers, as shown in more detail in Di Matteo et al. (2005). The simulation results are in accord with the view that BH growth is regulated by AGN feedback and the dynamical response of gas to this supply of energy.

We note that in our model the AGN feedback energy is deposited isotropically in the gas around the BH. While this appears to be a natural approximation given that we cannot resolve the detailed accretion physics close to the hole, observations show that BH outflows are typically bipolar (although not necessarily strongly collimated) and often result from powerful jets. It is possible that our isotropic coupling enhances the effects of the feedback, but this effect should be compensated in part by our comparatively conservative assumption for e_{f}. We also caution that the physical nature of the coupling of AGN feedback to the surrounding gas is not fully understood yet. In particular, there may be important feedback mechanisms other than thermal heating. For example, radiation pressure may play a very important role as well (Murray, Quataert & Thompson 2005).

The presence of central BHs can also have a significant impact on the stellar density structure of merger remnants. Mihos & Hernquist (1994a) found that their remnants typically had centrally peaked luminosity profiles from the compact stellar cores left behind by starbursts. In our new simulations, the remnants have smoother luminosity profiles into their centres, in much better agreement with observations of elliptical galaxies (Springel, Di Matteo & Hernquist, in preparation). This outcome is driven by a combination of the consumption of dense gas by the BHs, the dynamical heating of the gas as the BHs spiral together and merge, and the dispersal of gas by AGN feedback. Our results thus indicate that BH growth may be an integral part of spheroid formation.

## Conclusions

In this paper, we have introduced a new methodology for simultaneously modelling star formation and BH accretion in hydrodynamical simulations of isolated and merging galaxies. Our approach uses the concept of forming ‘macroscopic’, coarse-grained representations of subresolution physics, allowing us to study the effects of these processes on larger, resolved scales. To illustrate the principal effects of our methods and of the physics we describe, we have focused on simulations of isolated and merging galaxies. Throughout, we have concentrated on the methodological aspects of our work, rather than on an in-depth analysis of the physical results of our simulations. We provide the latter in related papers (e.g. Di Matteo et al. 2005; Springel et al. 2005).

For our simulations, we have constructed compound galaxy models that consist of a dark matter halo, a rotationally supported disc of gas and stars, and a central bulge, with independent parameters describing each of the structural components. Our approach improves upon previous methods by accounting for gas pressure forces, as well as featuring a flexible numerical approach for computing the gravitational potential for non-trivial, radially varying mass distributions in the disc. As a result, we are able to construct very stable disc galaxy models even for high gas fractions.

In simulations of isolated galaxies, we have emphasized the importance of the treatment of feedback processes in the dense star-forming ISM for the stability of disc galaxies. Simulations that use a simple isothermal EOS, or a weak form of kinetic feedback, have gas discs that are very susceptible to axisymmetric perturbations. As a result, galaxies with a large gas surface mass density cannot be evolved in a stable fashion for a long time in such models. Our multiphase model for the ISM encapsulates the effects of local supernova feedback in the form of an EOS, which is comparatively stiff. When this description for the ISM is chosen, the dense gas is stabilized by supernovae pressurizing the gas, so that much larger gas fractions than in previous studies become possible.

We used simulations of isolated galaxy models to examine the properties of the various phases of accretion possible in our model for the growth of BHs. In particular, we showed how the BH growth of a small seed BH accelerates in the Bondi regime and eventually reaches Eddington-limited, exponential growth. However, the growth process can be slowed down and regulated by the feedback energy associated with the accretion. We assumed a simple thermal coupling of a small fraction of the bolometric luminosity of the hole to the surrounding gas; this can increase the local sound speed of the gas and thereby reduce the Bondi accretion rate, or in extreme cases of high accretion, generate a powerful, pressure-driven quasar wind.

Using a number of simulations of major galaxy mergers, both with and without BHs, we have explored a few basic consequences of our methods. When we model the star-forming gas with an isothermal EOS and weak feedback, our results are in good agreement with previous work by Mihos & Hernquist (1996). However, our multiphase treatment of the ISM allows us to investigate galaxy mergers with much higher gas fractions than possible before. Here, new types of outcomes are possible for major mergers. For example, very high star formation rates of several hundred solar masses per year can be reached in gas-rich mergers. In merger simulations with BHs, we have shown that the presence of accreting BHs can dramatically alter the dynamics of the merger. The feedback energy associated with the growth of the hole owing to the tidally triggered inflow of gas has a direct effect on the strength of the nuclear starburst. During the final galaxy coalescence, the BH can expel the gas from the centre in a powerful outflow, quenching the starburst on a short time-scale. The resulting elliptical remnant is then gas poor and reddens quickly.

Already, these exploratory results show clearly that adding an accreting supermassive BH in the centre can have a decisive influence on galaxies. In fact, we think this indicates that the build-up of supermassive BHs and massive galaxies is an intertwined process, which requires a self-consistent treatment to be meaningfully addressed. Our new numerical methodology represents the first attempt to do carry out this demanding task in cosmological hydrodynamical simulations, opening up the possibility to jointly study the formation of galaxies and the build-up of supermassive BHs hosted by them. Our methods thus provide a powerful tool to gain new insights into the process of hierarchical galaxy formation.

### Acknowledgments

This work was supported in part by NSF grants ACI 96-19019, AST 00-71019, AST 02-06299 and AST 03-07690, and NASA ATP grants NAG5-12140, NAG5-13292 and NAG5-13381. The simulations were performed at the Centre for Parallel Astrophysical Computing at the Harvard-Smithsonian Centre for Astrophysics.