Galaxy-induced transformation of dark matter haloes

Abstract

We use N-body/gasdynamical cosmological simulations to examine the effect of the assembly of a central galaxy on the shape and mass profile of its surrounding dark matter halo. Two series of simulations are compared; one that follows only the evolution of the dark matter component of individual haloes in the proper Λcold dark matter (ΛCDM) cosmological context, and a second series where a baryonic component is added and followed hydrodynamically. The simulations with baryons include radiative cooling but neglect the formation of stars and their feedback. The efficient, unimpeded cooling that results leads most baryons to collect at the halo centre in a centrifugally supported disc which, due to angular momentum losses, is too small and too massive when compared with typical spiral galaxies. This admittedly unrealistic model allows us, nevertheless, to gauge the maximum effect that galaxies may have in transforming their surrounding dark haloes. We find, in agreement with earlier work, that the shape of the halo becomes more axisymmetric: post galaxy assembly, haloes are transformed from triaxial into essentially oblate systems, with well-aligned isopotential contours of roughly constant flattening (〈c/a〉∼ 0.85). Haloes always contract as a result of galaxy assembly, but the effect is substantially less pronounced than predicted by the traditional ‘adiabatic-contraction’ hypothesis. The reduced contraction helps to reconcile ΛCDM haloes with constraints on the dark matter content inside the solar circle and should alleviate the longstanding difficulty of matching simultaneously the scaling properties of galaxy discs and the galaxy luminosity function. The halo contraction we report is also less pronounced than found in earlier simulations, a disagreement which suggests that halo contraction is not solely a function of the initial and final distribution of baryons. Not only how much baryonic mass has been deposited at the centre of a halo matters, but also the mode of its deposition. Although simple formulae might work in particular cases where galaxies form nearly adiabatically, in general it might prove impossible to predict the halo response to galaxy formation without a detailed understanding of a galaxy's detailed assembly history.

Galaxy: disc, Galaxy: formation, Galaxy: kinematics and dynamics, Galaxy: structure

1 INTRODUCTION

Over the past couple of decades, cosmological N-body simulations have established a number of important results regarding the structure of cold dark matter (CDM) haloes. In particular, broad consensus holds regarding the overall shape of the halo mass distribution and its radial mass profile. CDM haloes are distinctly triaxial systems supported by anisotropic velocity dispersion tensors (Frenk et al. 1988; Jing & Suto 2002; Allgood et al. 2006; Hayashi, Navarro & Springel 2007; Bett et al. 2007), a result that reflects the highly aspherical nature of gravitationally amplified density fluctuations as they become non-linear and assemble into individual haloes. Despite their seemingly chaotic assembly, the spherically averaged mass profile of CDM haloes is nearly ‘universal’, and can be approximated fairly accurately, regardless of halo mass or redshift, by scaling a simple formula (see e.g. Navarro, Frenk & White 1996b, 1997; Navarro et al. 2004, 2010).

One major limitation is that these results have been obtained from simulations that neglect the presence of the baryonic, luminous component of galaxy systems. Although this may be a suitable approximation in extremely dark matter-dominated systems, it is expected to fail in regions where baryons contribute a substantial fraction of the mass, as is the case in the luminous regions of normal galaxies like the Milky Way (MW). The response of the dark halo to the assembly of a galaxy is thus a crucial ingredient of models that attempt to interpret observational constraints in terms of the prevailing CDM paradigm.

There has been extensive work, both numerical and analytical, on the modification of the mass profile of a dark halo induced by the assembly of a central galaxy. Early results suggested that simple analytical prescriptions based on the conservation of adiabatic invariants gave an accurate description of the halo response. Following the early work by Barnes & White (1984), Blumenthal et al. (1986) devised a simple formula to link the dark mass profiles before and after the assembly of a galaxy. Given the initial, spherically symmetric enclosed mass profiles of the dark matter, Mⁱ_dm(r), and baryons, Mⁱ_b(r), one may derive the final dark mass profile, M^f_dm(r), once the final baryonic mass profile, M^f_b(r), is specified. The model assumes that dark matter particles move on circular orbits before and after the contraction, and that their initial, r_i, and final, r_f, radii are related by the condition:

1

where M_dm=M_dm(r_f) =M_dm(r_i) is the dark mass enclosed by each dark matter particle (i.e. no shell crossing).

Despite its simplicity, and the crude approximation on which it is based, early N-body work (e.g. Barnes & White 1984; Blumenthal et al. 1986; Jesseit, Naab & Burkert 2002) found reasonable agreement between equation (1) and the results of simulations, and helped to establish the ‘adiabatic-contraction’ formulation as the default procedure when considering the halo response to the formation of a central galaxy.

Semi-analytic models of galaxy formation have adopted this formulation in order to link the (observed) dynamical properties of luminous galaxies to the (predicted) properties of the dark haloes that surround them (see e.g. Mo, Mao & White 1998; Cole et al. 2000; Firmani & Avila-Reese 2000; Croton et al. 2006; Dutton et al. 2007; Somerville et al. 2008). This work has highlighted a number of potential problems when attempting to reconcile the predicted properties of galaxies in the ΛCDM cosmogony with observed scaling relations.

One particularly important challenge has been the inability of semi-analytic galaxy formation models to match simultaneously the zero-point of the Tully–Fisher (TF) relation and the galaxy luminosity function. Successful models require that the rotation speed of discs be of the order of the virial velocity of the haloes they inhabit (see e.g. Somerville & Primack 1999; Croton et al. 2006). However, models that include adiabatic contraction typically predict disc rotation speeds well in excess of the virial velocity.

As discussed most recently by, for example, Dutton et al. (2007), Gnedin et al. (2007) and Dutton, van den Bosch & Courteau (2008), potential solutions to this problem include: (i) revising the adiabatic hypothesis so as to reduce (or even reverse!) the halo contraction; (ii) adopting lighter stellar mass-to-light ratios in order to minimize the contribution of baryons and to allow for more dark mass enclosed within disc galaxies or (iii) modifying the cosmological parameters so as to reduce halo concentrations. We shall concentrate our analysis here on option (i), although we note that all three possibilities are probably equally important and should be treated on equal footing when addressing these issues in semi-analytic models of galaxy formation.

The possibility that the actual response of the halo to the assembly of the disc might differ substantially from the predictions of the adiabatic-contraction formalism has been noted before. Barnes (1987) and Sellwood (1999) have remarked that the adiabatic-contraction hypothesis might lead to substantial overestimation of the halo compression. As discussed by Sellwood & McGaugh (2005) and Choi et al. (2006), such deviations are likely to depend strongly on the orbital structure of a halo, performing best when most particles have large tangential motions, but increasingly poorly as the tangential velocity anisotropy is reduced.

It should also be pointed out that the studies mentioned above were carried out by perturbing simple models with a potential term designed to imitate the growth of an assembling galaxy. These studies, therefore, miss the hierarchical nature of the assembly of a galaxy and of its surrounding halo. In particular, the many mergers and accretion events expected in such formation scenario may proceed on relatively short time-scales, at odds with the assumptions on which the adiabatic response of the halo is based.

The study of Gnedin et al. (2004), on the other hand, uses the proper cosmological context, but focuses on simulations of galaxy clusters, rather than the galaxy-sized systems of interest for the issues discussed above. In this respect, our study is similar to that of Gustafsson, Fairbairn & Sommer-Larsen (2006), who studied four simulations of galaxy-sized haloes in the ΛCDM scenario. Their conclusions are in good agreement with the results we discuss here.

The assembly of a central galaxy also modifies the three-dimensional shape of the surrounding dark matter halo. This was already noted in the first simulations to include, in addition to a dark matter halo, a dissipative baryonic component, such as the early work of Katz & Gunn (1991) and Katz & White (1993). Dubinski (1994) studied this further by growing adiabatically a central mass concentration inside a triaxial dark matter halo and confirmed that the steepening of the potential leads to much rounder halo shapes. This result was confirmed by Kazantzidis et al. (2004), who also noticed the effect in their cluster simulations. We extend this body of work by focusing, as in Hayashi et al. (2007), on the shape of the gravitational potential rather than on the axial ratios of the inertia tensor, as well as on its dependence on radius.

The plan of this paper is as follows. Section 2 describes the numerical simulations and Section 3 presents the main results. After a brief general description of the evolution in Section 3.1, Section 3.2 describes the main properties of the simulated central galaxies; Section 3.3 discusses the modifications they induce on the shape of the halo gravitational potential; while Section 3.4 compares the mass profiles of dark haloes before and after the inclusion of a dissipative baryonic component. We apply these results to the MW in Section 3.5 and compare with earlier work in Section 3.6. We end with a brief summary of our main conclusions in Section 4.

2 NUMERICAL SIMULATIONS

2.1 Cosmological parameters

We adopt cosmological parameters consistent with the combined analysis of the 2dfGRS (Colless et al. 2001) and the first-year Wilkinson Microwave Anisotropy Probe (WMAP) data (Spergel et al. 2003): a present-day value of the Hubble constant of H₀= 70 km s⁻¹ Mpc⁻¹; a scale-free initial density perturbation spectrum with no tilt and normalized by the linear rms mass fluctuations on 8 h⁻¹ Mpc spheres, σ₈= 0.9. The matter-energy content of the Universe is expressed in units of the present-day critical density for closure, and contains a dominant cosmological constant term, Ω_Λ= 0.7, as well as contributions to the matter content, Ω_M= 0.3, from CDM, Ω_CDM= 0.259, and baryons, Ω_b= 0.041.

2.2 Code and initial conditions

We have performed a suite of 13 numerical simulations of the formation of galaxy-sized haloes in the ΛCDM cosmogony. Each simulation follows the evolution of a relatively small region of the Universe, excised from a 432³-particle simulation of a large 50 h⁻¹ Mpc periodic box (Reed et al. 2003), and resimulated at much higher mass and spatial resolution. Each ‘zoomed-in’ re-simulation follows the formation of a single galaxy-sized halo of mass ∼ 10¹² M_☉ and its immediate surroundings. The simulations include the tidal fields of the parent simulation, and follow the coupled evolution of gas and dark matter. The hydrodynamical evolution of the gaseous component is followed using the smooth particle hydrodynamics (SPH) technique. This re-simulation technique follows standard practice, as described, for example, by Power et al. (2003). All our simulations were performed using gasoline, a parallel N-body/SPH code described in detail by Wadsley, Stadel & Quinn (2004).

2.3 Halo selection and resimulation

The 13 haloes were selected from a list of all ∼ 10¹² M_☉ haloes in the parent simulation, with a mild bias to avoid objects that have undergone a major merger after z= 1 or that have, at z= 0, unusually low (λ < 0.03) spin parameter. The mass resolution of the resimulations is such that haloes are represented with ∼ 40 000 dark matter particles within their virial radius1 at z= 0. Two sets of simulations are performed for each halo, one where only a dark matter component is followed, and another where dark matter and baryons are included. Dark matter-only (DMO) simulations assume that the total matter content of the Universe is in CDM, Ω_M=Ω_CDM= 0.3.

Simulations including a baryonic component (‘DM+B’) include equal number of gas and dark matter particles, with masses modified to satisfy our choice of Ω_b= 0.041 and Ω_CDM= 0.259. The gaseous component is evolved using SPH including, besides the self-gravity of gas and dark matter, pressure gradients and shocks. Baryons are allowed to cool radiatively according to the cooling function of a gas of primordial composition down to a temperature of 10⁴ K, below which cooling is disabled. No star formation or feedback is included. As we discuss below, this choice maximizes cooling and favours the collection of most baryons at the centre of each dark halo. Although unrealistic as a model for galaxy formation, this choice has the virtue of simplicity and also allows us to examine the halo response when the effect of the baryonic component is maximal.

Pairwise gravitational interactions are softened adopting a spline softening length kept fixed in comoving coordinates, ε= 2.5 h⁻¹ kpc. We test for numerical resolution effects by simulating each halo with eight times fewer particles and softening lengths twice as large. In the interest of brevity, we do not present detailed results from this low-resolution series here, but we have used these runs to verify numerical convergence and to understand the reliability of our conclusions. In particular, we observe good convergence in the mass profiles down to the inner convergence radius, r_conv, where collisional relaxation time-scales become comparable to the age of the Universe (Power et al. 2003). Our simulations support the conclusion of Power et al. (2003) that, outside of this radius, estimates of the circular velocity converge to better than 10 per cent accuracy.

Furthermore, we have also tested for numerical convergence by increasing the number of particles in one of the simulated haloes. Halo S02h is equivalent to S02, but with ∼ 3.4 times more particles: at z= 0 this halo has ∼ 140 000 particles per component within the virial radius. A softening length of ε= 1.67 h⁻¹ kpc has been adopted in this case. We explicitly show the convergence in the mass profiles for the three resimulations of this halo (S02h, S02, S02l) in Fig. 7. Table 1 lists the main properties of the simulated haloes.

Figure 7

Enclosed mass profile of various components of halo S02h. The ‘DMO’ profile is shown by the thick solid curve, after scaling masses by (1−f_bar), so that, within r_vir, the total dark mass of the DMO and DM+B runs is comparable. The thin line shows a Navarro–Frenk–White halo fit to the DMO profile. Other colours and line types correspond to the various components of the DM+B run, as specified in the figure labels. The thick dashed blue curve shows the dark mass profile for the DM+B run. The assembly of most baryons into a central galaxy leads to a contraction of the dark mass profile. The profile predicted by the ‘adiabatic-contraction’ formula (equation 1) overpredicts the response of the halo. The modified adiabatic-contraction model of Gnedin et al. (2004), ‘Contra’, also overestimates the halo response. The thinner ‘DM+B’ curves correspond to two further resimulations of lower numerical resolution (S02 and S02l; see Table 1). The crosses indicate the convergence radius, r_conv, of each run, as defined by Power et al. (2003). Note the excellent agreement in mass profiles outside the convergence radius of each run.

Table 1

Main parameters of simulated haloes. M_vir is the total (dark matter plus gas) mass inside the virial radius, r_vir. N_DMO and N_DM are the number of dark matter particles inside the virial radius r_vir in the ‘DMO’ and ‘dark matter plus baryons’ simulations, respectively. Columns 6 and 7 list the number of gas particles, N_gas, within the virial radius and within the radius used to define the central galaxy, r_glx= 10 kpc. The specific angular momentum of the dark matter within r_vir and of the baryons within r_glx is also listed. Last two columns show the gravitational softening ε and the convergence radius r_conv.

Label	Mvir (1012M☉)	rvir (kpc)	NDMO (< rvir)	NDM (< rvir)	Ngas (< rvir)	Ngas (< rglx)	jDMO (kpc km s−1)	jDM (kpc km s−1)	jgas (kpc km s−1)	ε (kpc)	rconv (kpc)
S01	1.22	277.84	53 009	52 686	54 734	41 008	1847.7	2162.8	462.3	3.57	6.91
S02	0.94	254.18	40 847	40 687	39 775	31 097	1063.3	1088.3	332.8	3.57	6.89
S02h	0.95	255.89	139 748	139 729	140 618	105 000	1095.7	1200.3	359.1	2.38	4.08
S02l	0.93	253.63	5058	5083	4716	4020	989.9	1066.1	289.1	7.14	17.08
S03	1.51	298.24	64 656	65 376	66 416	38 266	4410.1	2883.4	438.1	3.57	6.73
S04	0.85	246.02	37 419	36 432	39 007	30 767	1443.1	1626.7	508.8	3.57	6.97
S05	0.99	258.64	42 564	42 334	45 270	30 450	2505.0	2543.2	173.3	3.57	6.96
S06	1.22	278.05	53 509	52 652	55 674	31 392	2736.0	3056.0	289.4	3.57	7.51
S07	0.86	246.95	37 350	37 305	36 524	26 034	1332.2	1385.1	310.0	3.57	6.83
S08	0.84	245.39	35 902	36 224	38 221	26 049	1221.3	1292.5	277.6	3.57	7.23
S09	0.84	245.23	36 619	36 259	37 441	29 142	850.5	929.5	531.8	3.57	6.80
S10	0.97	257.68	42 004	42 074	43 469	35 034	1976.6	2231.8	430.7	3.57	6.60
S11	0.88	249.53	40 181	38 190	39 551	23 640	1430.3	1676.8	289.6	3.57	7.13
S12	0.90	251.25	39 715	39 341	38 140	28 898	715.7	764.0	381.5	3.57	6.70
S13	1.08	266.72	46 863	46 337	50 223	37 270	1470.7	1607.5	261.3	3.57	6.99

3 RESULTS

3.1 General evolution

Because our simulations neglect the formation of stars and their feedback, the evolution of the baryonic component is characterized by the rapid cooling and collapse of baryons at the centre of the early collapsing progenitors of the final halo. As a result, the main mode of galaxy assembly are mergers: up to 70 per cent of all baryons in the central galaxy were accreted in the form of dense, cold, gaseous clumps that sink to the centre through dynamical friction.

As shown in earlier work (see e.g. Navarro & Benz 1991; Navarro & White 1994; Navarro et al. 1997), this mode of assembly leads to very efficient accretion of baryons into the central galaxy and to the transfer of much of their angular momentum to the surrounding halo. Because baryons are assumed to remain in gaseous form at all times, those that can cool are forced to settle into centrifugally supported discs at the centre of their surrounding dark haloes. Thus, we expect the central galaxies in our simulations to be discs that contain a large fraction of all baryons within the virial radius, and to have angular momenta well below the angular momentum content of the system as a whole.

Fig. 1 illustrates this for the case of our highest resolution halo, S02h (see Table 1). The four panels in this figure show the dark matter (black) and baryonic (coloured) components, zooming by consecutive factors of 3 towards the central galaxy. Cold gas (T < 10^4.5 K) is shown in green, whereas hot gas (T > 10^4.5 K) is shown in magenta. Note that most of the cold gas inhabits the centre of the main halo and its substructures, where it forms easily identifiable thin, massive discs. All cold baryons associated with the central disc are contained within a sphere of radius r_glx= 10 kpc, which we shall use hereafter to define the central galaxy in all runs.

Figure 1

Dark matter (black) and gas (coloured) particles for halo S02h. The colour of a gas particle reflects its temperature, T < 10^4.5 K (green) or T > 10^4.5 K (magenta). Each panel zooms into the centre of the system by consecutive factors of 3 in radius (see axis labels). The circle in the top-left panel shows the virial radius. The circle in the bottom-right panel shows the radius, r_glx, used to define the central galaxy. Note that the colder (green) gas component inhabits the centre of the main halo and of its substructures, where it forms thin, centrifugally supported discs. The hotter component (magenta) is distributed more or less uniformly across the main halo and makes up only 16 per cent of all the baryons within the virial radius. The system has been rotated so that the z-axis is aligned with the angular momentum of all cold gas particles within r_glx. Note that the gaseous disc is mildly warped and that, as a result, the inner disc is slightly tilted relative to the spin axis of the outer disc particles. The latter dominate the appearance of the disc in the top two panels, and lead to a slight but noticeable tilt of the central disc relative to the coordinate axes.

3.2 Mass, size and angular momentum of central galaxies

The central galaxy contains most of the baryons within the virial radius of the system (74 per cent in the case of S02h; see Table 1). This is shown in Fig. 2, where we plot, for all our simulations, the baryonic mass of the central galaxy, M_disc, versus the radius, r₉₀, that contains 90 per cent of its mass, scaled to the virial mass and radius of the system, respectively. The vertical dashed line indicates the universal baryon fraction in the simulations, f_bar=Ω_bar/Ω_M= 0.137. Typically, 70 to 80 per cent of all baryons within r_vir are found in fairly small central discs that are fully contained within a radius of the order of ∼ 3 per cent of the virial radius.

Figure 2

Radius containing 90 per cent of the baryonic mass, r₉₀, versus the total baryonic mass of the central galaxy, M_disc, in our simulated haloes (filled circles), expressed in units of the halo virial values. The dashed vertical line indicates the universal baryon fraction, f_bar=Ω_b/Ω_M. The central galaxy typically makes up about 70–80 per cent of all baryons within r_vir. The solid line shows the corresponding values for the MW, varying the virial mass of its halo. Symbols along the line indicate a few different values of the virial velocity, V_vir= 220, 150, 110 and 92 km s⁻¹. Note that the simulated central galaxies are substantially more massive and smaller than spirals like the MW. This maximizes the effect of baryons in transforming the shape and mass profile of the dark halo.

We may compare this with a typical spiral galaxy like the MW, where the mass and size of the baryonic component can be estimated accurately. Assuming that the disc is exponential with total mass 4.5 × 10¹⁰ M_☉ and radial scalelength R^MW_d= 2.5 kpc and that the mass of the bulge is M_bulge= 4.5 × 10⁹ M_☉, we estimate that 90 per cent of the MW baryons are confined within r^MW₉₀= 9.4 kpc.

In order to compare this with the simulation results shown in Fig. 2 we need to make an assumption about the virial mass of the MW halo. For a virial velocity of the order of the rotation speed at the solar circle (220 km s⁻¹), the baryonic component makes up only ∼ 1 per cent of the virial mass, or roughly 7 per cent of all baryons within r_vir. This is about an order of magnitude less than the fraction of baryons that collect in the central discs of the simulated galaxies. Only if the virial velocity is as low as 110 km s⁻¹ do the MW bulge+disc contain a fraction of available baryons comparable to that seen in the simulations (i.e. ∼ 75 per cent). However, in that case r^MW₉₀∼ 0.06 r_vir, a factor of ∼ 4 times larger than the simulated discs.

The reason why simulated gaseous discs are so much smaller than typical spirals is that baryons have lost a large fraction of their angular momentum to the surrounding halo. This is illustrated in Fig. 3. where we plot the angular momentum of the gaseous disc (in units of that of the system as a whole, which is dominated by the dark matter component), J_d≡J_disc/J_vir=M_discj_disc/M_virj_vir, versus the mass fraction, m_d≡M_disc/M_vir. These are the parameters commonly adopted in semi-analytic models of disc formation, such as those of Mo et al. (1998).

Figure 3

Angular momentum and mass of the central galaxy in our simulations, expressed in units of the virial values. The curve labelled j_disc=j_vir corresponds to central galaxies with the same specific angular momentum as the surrounding halo (a common assumption of semi-analytic models). Those on the curve labelled j_disc/j_vir=m_d/f_bar have specific angular momenta (in units of the halo's) that scale in proportion to the fraction of baryons within the virial radius that have collected in the central galaxy. Note that the simulated central galaxies have even lower angular momenta than in the latter assumption, highlighting the large angular momentum losses that accompany the assembly of the central galaxy, and explaining the small sizes of the discs shown in Fig. 2.

Since (i) baryons acquire during the expansion phase as much angular momentum as the dark matter, and (ii) it is unlikely that M_disc/M_vir will exceed the universal baryon fraction, we do not expect the specific angular momentum of the disc, j_disc, to exceed that of the system as a whole, j_vir. Therefore, central galaxies are unlikely to populate the shaded areas of Fig. 3.

Most semi-analytic work assumes, for simplicity, that j_disc=j_vir or, equivalently, that J_d=m_d, regardless of the value of m_d. On the other hand, Navarro & Steinmetz (2000) argue that it is unlikely that a galaxy where m_d << f_bar would have the same specific angular momentum as the whole system, and propose that j_disc, as a fraction of j_vir, should be comparable to the fraction of baryons that make up the central galaxy; i.e. j_disc/j_vir=M_disc/(f_barM_vir) =m_d/f_bar. This is illustrated in Fig. 3 by the lower diagonal line.

Our simulated galaxies are well below both lines. This indicates that baryons have specific angular momenta much lower than their surrounding haloes and explains their small sizes compared to typical spirals. These are therefore unrealistic models of spiral galaxy formation, but the combination of large mass and small size allows us to probe the response of the dark halo in the case where the deepening of the potential well due to the central galaxy is maximal.

3.3 Halo shape

As anticipated in Section 1, the halo responds to the presence of the central galaxy by becoming significantly more spherical. This is illustrated in Fig. 4, where we compare the shape of isodensity and isopotential contours for the ‘dark matter only’ (DMO) and ‘dark matter plus baryons’ runs of halo S02h. Panels on the left correspond to the DMO run, and those on the right to the DM+B run. Particles are coloured according to the local value of their gravitational potential or local density (binned in a logarithmic scale), computed using only the dark matter particles.

Figure 4

Dark matter particles corresponding, at z= 0, to the two resimulations of halo S02h, ‘DMO’ (left-hand panels) and ‘dark matter plus baryons’ (right-hand panels). Particles are coloured according to the value of the gravitational potential (binned in logarithmic units, top panels) or the local density (bottom panels). Note that the halo responds to the assembly of the central galaxy (shown in black in the right-hand panels) by becoming noticeably more spherical. As in Fig. 1, projections are chosen so that the angular momentum axis of the central disc coincides with the z-axis. The central disc is shown with black dots to illustrate the orientation of the disc relative to the halo shape. The disc rotation axis is reasonably well aligned with the minor axis of the halo, although a slight tilt between the two (probably associated with the weak disc warp seen in Fig. 1) is noticeable. Particles are plotted in sequence of increasing potential depth/density; as a result, those in the central regions overlap and hide outer particles that may be projected near the centre.

Note that using the gravitational potential leads to rounder isocontours (as expected due to the importance of the monopole term in the potential) as well as to much more stable estimates of the shape of the halo, since it is much less affected by the presence of substructures and other transient fluctuations in the mass distribution.

The isopotential contours are well approximated by ellipsoidal surfaces, and we use the axial ratios of such ellipsoids to measure, as a function or radius, the change in shape of halo S02h. We show the result in Fig. 5. As discussed by Hayashi et al. (2007), the radial dependence of the axial ratios may be approximated by the formula

2

Here, α parametrizes the central value of the axial ratio, (b/a)₀ or (c/a)₀, by 10^−2α; r_α indicates the characteristic radius at which the axial ratio increases significantly from its central value and γ regulates the sharpness of the transition.

Figure 5

‘Radial’ dependence of axial ratios b/a (open symbols) and c/a (filled symbols) for our high-resolution run S02h. Halo shapes are measured by fitting 3D ellipsoids to the position of dark matter particles in narrow logarithmic bins of the gravitational potential. The potential is computed using only dark matter particles, so it actually corresponds to the contribution of the dark matter component to the overall potential. Red (lower) symbols correspond to the DMO run; blue (upper) symbols to the DM+B run. Note that the central galaxy turns a rather triaxial, nearly prolate halo into an axisymmetric, nearly oblate one. Curves are fits using equation (2). Fit parameters for this halo and the mean value computed over the sample are listed in Table 2.

The presence of the central galaxy turns the halo from a triaxial, nearly prolate system into an axisymmetric, nearly oblate system where the axial ratio is nearly constant with radius. As a result, equation (2) is not adequate for the nearly featureless radial dependence of the axial ratios in the DM+B runs. We fit the latter with a simple constant value. Best-fitting parameters are listed in Table 2.

Table 2

Best-fitting parameters for isopotential contour profile shapes. Columns 1–3 (4–6) correspond to axial ratio b/a (c/a) for DMO runs (equation 2). Column 7 (8) corresponds to axial ratio b/a (c/a) for Dark Matter+Baryons (DM+B) runs. We quote values for the high-resolution run S02h shown in Fig. 5 and also the average for all 13 haloes (shown as thick lines in Fig. 6).

	b/a (DMO)			c/a (DMO)			b/a (DM+B)	c/a (DM+B)
	α	γ	rα (rvir)	α	γ	rα (rvir)
S02h	0.113	1.767	0.227	0.137	2.102	0.323	0.976	0.857
Average	0.081	1.383	0.139	0.103	1.427	0.265	0.948	0.857

The profile fits for all simulations are compiled in Fig. 6, and illustrate a few important points. As a result of the assembly of the central galaxy: (i) haloes become nearly oblate; (ii) axial ratios are roughly independent of radius and (iii) halo shapes are affected well beyond the size of the central galaxy, and nearly as far out as the virial radius.

Figure 6

‘Radial’ dependence of axial ratios b/a (upper panels) and c/a (lower panels) for all simulated haloes. The curves for the ‘DMO’ runs are computed by fitting equation (2) to the radial dependence of the axial ratios of particles along isopotential contours. The curves for the ‘dark matter + baryons’ (DM+B) runs are the measured axial ratios of particles along isopotential contours. The contours in all cases refer to the potential contributed solely by the dark matter component. Note that, as a result of the assembly of the central galaxy, (i) haloes become nearly oblate and (ii) axial ratios become approximately independent of radius. The effect on the shape of the halo extends almost to the virial radius, far beyond the actual size of the central galaxy. Thick lines correspond to mean values computed over the sample. Shaded areas indicate regions where the radius is smaller than the softening length.

The results shown here are in good agreement with those reported in earlier work, although a direct comparison is hindered by the fact that many studies used galaxy-cluster simulations to illustrate the effect, or simulations where the assembly of the central galaxy was not modelled self-consistently. Furthermore, most prior work focused on changes in the shape of the isodensity contours or the inertia tensor of a halo rather than on the shape of the gravitational potential, as we do here.

With these caveats in mind, our results are in good agreement with these earlier results. Dubinski (1994), for example, report modifications of up to Δ(b/a) ∼ 0.2 and Δ(c/a) ∼ 0.1 by growing a massive perturber at the centre of a triaxial halo. Kazantzidis et al. (2004) find that the central parts of their simulated galaxy clusters change by Δ(b/a) ∼ 0.4 and Δ(c/a) ∼ 0.2 and report that the changes are radially dependent: greatest near the centre and negligible at or around the virial radius. Debattista et al. (2008) remark that even as far away as half the virial radius (well outside the radius where most baryons reside) the axial ratios might still be modified by ∼ 0.2 or more.

Our results are even more dramatic than this. As one may appreciate directly from Fig. 6, haloes are rendered essentially oblate (〈b/a〉∼ 1, 〈c/a〉∼ 0.85) by the central galaxy in our simulations. More realistic galaxy formation simulations may lead to less extreme transformation, and it remains an interesting question whether typical galaxies truly make their haloes axisymmetric or are just able to reduce their triaxiality by some limited amount. The quantitative answer most likely depends on the mass and concentration of the baryonic component relative to its surrounding halo: low-surface-brightness galaxies may very well be surrounded by haloes which are on average more triaxial than those of their high-surface-brightness counterparts. We are currently working on simulations aimed at addressing these issues in detail.

3.4 Mass profile and contraction

Fig. 7 shows the enclosed mass profile of the DMO and DM+B runs corresponding to halo S02h. DMO dark masses have been scaled by (1−f_bar) so that the total dark mass in both the DMO and DM+B runs are the same. The DMO mass profile is shown by the thick red solid curve and, as expected, it is well approximated by the NFW (Navarro et al. 1996b, 1997) formula (thin solid line). The central disc that assembles in the DM+B run leads to a contraction of the dark halo: there is more dark mass in the DM+B run (thick dashed line) at all radii compared with the DMO run. This is a result common to all our simulations; in no case do we see the halo ‘expand’ as a result of the assembly of the central galaxy.

It is important to verify that these conclusions are not unduly affected by numerical resolution. To this aim, we compare in Fig. 7 the results of three DM+B resimulations of the same halo (S02, S02h, S02l; see Table 1). The dark and baryonic mass profiles of the three simulations are indicated with lines of the same type but varying thickness; the thickest and thinnest correspond to the highest (S02h) and lowest (S02l) numerical resolution case, respectively. The crosses along each curve mark the convergence radius, r_conv, of each simulation, as defined by Power et al. (2003). There is clearly very good convergence (typically better than 10 per cent) for all profiles outside r_conv, indicating that our results are reliable outside this radius.

The halo contraction seen in Fig. 7 is not as pronounced as what would be expected from the ‘adiabatic-contraction’ model discussed in Section 1. The adiabatic-contraction prediction (equation 1) is shown by the thick dotted line in Fig. 7. The discrepancy between model and numerical results is not small. At ∼ 5 kpc, a radius roughly twice the gravitational softening and that encloses more than 1300 dark particles, the adiabatic-contraction formula predicts ∼ 2.5 times more dark mass than found in our simulations. Even at a radius of 10 kpc, equation (1) overestimates the halo contraction by more than ∼ 50 per cent in mass. The modified contraction proposed by Gnedin et al. (2004) (see curve labelled ‘Contra’ in Fig. 7) fares better, but it still overpredicts the halo contraction relative to the results of our simulations.

A similar result to the one shown in Fig. 7 for halo S02 applies to all of our simulated haloes. We illustrate this in Fig. 8, where we plot the ratio between radii that contain a given number of dark matter particles in the DMO run (r_i) and the DM+B run (r_f) versus the ratio between M_i, the total mass within r_i in the DMO run, and M_f, the total mass within r_f in the DM+B run. With this choice, the adiabatic-contraction prediction is simply r_f/r_i=M_i/M_f, which is traced by the 1:1 line in Fig. 8.

Figure 8

Dark halo response to the assembly of a central galaxy. The ordinate shows the ratio, r_f/r_i, between the radius containing a given amount of dark mass in the DMO and DM+B runs, respectively, after scaling DMO masses by 1−f_bar. The smaller r_f/r_i the stronger the halo contraction. The x-axis shows the ratio between the total mass, M_i, contained within r_i (in the DMO run) and M_f, that enclosed within r_f (in the DM+B run). The adiabatic-contraction formula (equation 1) predicts that r_f/r_i=M_i/M_f; this is indicated by the 1:1 line in the figure. Numerical results are shown for individual haloes, from the radius that contains 1000 dark particles outwards. The adiabatic-contraction formula overestimates the halo response. Gnedin et al. (2004)'s modified adiabatic contraction (‘Contra’) does better but still overpredicts the halo response at most radii, and especially near the centre. Symbols with error bars trace the mean and rms of all haloes in our series. The thick upper curve is a fit using equation (3). Crosses indicate the convergence radius, r_conv, and open circles the gravitational softening, ε. Halo S02h (shown in Fig. 7) is highlighted with a thicker line.

The numerical simulations are shown by the upper curves in this figure. Near the centre, the baryons dominate over the DMO mass profile, and therefore M_i/M_f≪ 1. Further out, the contribution of baryons to the total enclosed mass decreases, approaching the universal baryon fraction at the virial radius, where M_i/M_f tends to unity. The innermost radius plotted for each profile indicates the location of the radius that contains 1000 dark matter particles. The gravitational softening is shown by open circles and the ‘convergence radius’, r_conv, by crosses.

The deviations from the adiabatic contraction and the ‘Contra’ predictions occur in regions that are well resolved numerically, and therefore are not the result of numerical artefact. The solid circles with error bars in Fig. 8 trace the median and rms scatter of all our simulation results. The ratio r_f/r_i tends to a constant for M_i/M_f≪ 1, suggesting that the halo response approaches saturation in regions where baryons dominate. A simple formula captures the average behaviour well (thick upper solid line),

3

with a= 0.3 and n= 2. As an application, we shall use this expression below to explore what constraints this implies for the mass and concentration of the halo of the MW. We caution, however, that equation (3) applies, strictly speaking, only to galaxies whose formation histories resemble the rather extreme assembly paths of simulated galaxies in our series.

3.5 Application to the Milky Way

The MW offers a case study for the results described above. Because the mass and radial distribution of the baryonic component, as well as the rotation speed of the local standard of rest (LSR), are relatively well known, the total dark mass contained within the solar circle is firmly constrained. Adopting the same quantities for the bulge and disc of the MW adopted in Section 3.2, we find that baryons contribute (in quadrature) ∼ 171 km s⁻¹ to the circular velocity of the LSR, which we assume to be 220 km s⁻¹. This implies that the dark mass of the MW within the solar circle (R_☉= 8 kpc) is ∼ 3.5 × 10¹⁰ M_☉. In order to allow for the possibility that some of this dark mass may be baryonic, we shall treat this as a formal upper limit in the analysis that follows.

The second constraint comes from the total baryonic mass of the MW, under the plausible assumption that the mass of the central galaxy cannot exceed the total mass of baryons within the virial radius, ≈f_barM_vir. Thus, the minimum virial mass allowed for the MW halo by this constraint is 3.64 × 10¹¹ M_☉, which corresponds to a virial velocity of ∼ 92 km s⁻¹.

Are these constraints compatible with ΛCDM haloes? The top-left panel of Fig. 9 shows the dark mass expected within 8 kpc versus halo virial velocity. Each dot in this panel corresponds to an NFW halo with concentration drawn at random from the mass–concentration relation (including scatter) derived by Neto et al. (2007) from the Millennium Simulation (Springel et al. 2005). As expected, the mass within 8 kpc increases with the virial velocity (mass) of the halo, modulated by the fairly large scatter in concentration at given halo mass.

Figure 9

Dark mass within 8 kpc versus virial velocity for haloes drawn at random from the ΛCDM mass–concentration relation (including scatter) of Neto et al. (2007). The DMO panel (top left) shows the results neglecting the effects of baryons and assuming an NFW profile for the dark haloes. The black ‘wedge’ highlights the region of parameter space compatible with observations of the MW. The panel labelled ‘AdiabCont’ shows the result of contracting each halo adopting the adiabatic-contraction formula (equation 1). Very few haloes would be consistent with the MW if this formula holds. The bottom-left panels show the result of applying our simulation results for the halo contraction, as captured by equation (3). The bottom-right panel repeats the same exercise, but after reducing all concentrations by 20 per cent in order to mimic the expected concentration change resulting from adopting the latest cosmological parameters from the 5-year analysis of WMAP satellite data (Duffy et al. 2008).

The black ‘wedge’ in this panel indicates the region allowed by the MW constraints discussed above. As discussed by Eke, Navarro & Steinmetz (2001), most ΛCDM haloes satisfy the constraints, and would be consistent with the MW if they were somehow able to avoid contraction. Note as well that the larger the MW halo virial velocity the lower the allowed concentration. For V_vir= 220 km s⁻¹ (the value required by semi-analytic models to match the TF relation and the luminosity function; see Section 1) only haloes with c_vir < 9.5 would be consistent with the MW. For comparison, the average concentration of V_vir= 220 km s⁻¹ haloes is 〈c_vir〉∼ 10.4 and its (lognormal) dispersion is σ_{log c}= 2.9, so this condition effectively excludes only 57 per cent of such haloes from the allowed pool.

The situation is rather different if haloes are adiabatically contracted (panel labelled ‘AdiabCont’ in Fig. 9). The contraction substantially increases the dark mass contained within 8 kpc, so that very few haloes satisfy the MW constraints. For example, essentially no halo with V_vir≈ 220 km s⁻¹ could host a galaxy like the MW, and the few that could would need to have c_vir < 3.2, many sigma away from the average concentration of haloes of that mass. The halo of the MW would need to be a very special halo of unusually low concentration if adiabatic contraction holds.

Adopting the halo contraction of our numerical simulations improves matters. This may be seen in the bottom-left panel of Fig. 9, where we have contracted each halo using equation (3). The range of allowed halo masses and concentrations is broader; even some haloes with V_vir∼ 220 km s⁻¹ could be consistent with the MW, provided that c_vir < 6.5. This is significantly lower than the average at that mass, but conditions are much less restrictive if the halo of the MW is less massive; for example, half of all V_vir∼ 130 km s⁻¹ haloes satisfy the MW constraints. The possibility that the virial velocity of the MW is significantly lower than 220 km s⁻¹ has indeed been advocated by a number of recent observational studies (Smith et al. 2007; Sales et al. 2007; Xue et al. 2008).

Lowering the average concentration at given mass also helps. This is shown in the bottom-right panel of Fig. 9, where we have repeated the exercise, but lowering all concentrations by 20 per cent. According to Duffy et al. (2008), this is approximately the change in average concentration when modifying the cosmological parameters from the values we adopt here (which are consistent with the first-year analysis of the WMAP satellite data) to those favoured by the latest 5-year WMAP data analysis. With this revision, 18 per cent of V_vir∼ 220 km s⁻¹ haloes and roughly half of all V_vir∼ 160 km s⁻¹ haloes would be consistent with the MW.

3.6 Comparison with earlier work on halo contraction

The failure of the adiabatic-contraction formalism to match the results of our simulations suggests that the response of the halo is more complex than what can be captured with a simple model for the contraction of spherical shells. This is confirmed by comparing our results with earlier work.

Our results disagree not only with the traditional adiabatic-contraction formula, but also with the modified formalism proposed by Gnedin et al. (2004) (shown with thin dash–dotted green lines in Fig. 8). These authors calibrated their results with numerical simulations not dissimilar to ours, except for the mass scale – they used mainly galaxy cluster haloes. Halo contraction is known to depend on the velocity distribution of the halo (see e.g. Valluri et al. 2010). However, cluster-sized haloes show the same mild radial anisotropy in their velocity profiles as galaxy-sized haloes, so it is unlikely that differences in orbital structure are the root cause of this result.

A clue to the reason for the discrepancy is provided by the cases shown with solid lines in Fig. 10. These correspond to two gasdynamical simulations of galaxy formation that differ from the ones presented in this paper mainly in their inclusion of star formation and feedback effects (Steinmetz & Navarro 2002). One of those simulations (labelled ‘S’ in Fig. 10) follows a halo with a quiet dynamical formation history and results in a galaxy with a prominent disc component (Abadi et al. 2003a,b). The other (labelled ‘E’ in Fig. 10) corresponds to halo with an active merging history that results at z= 0 in a spheroid-dominated ‘elliptical’ galaxy (see e.g. Meza et al. 2003).

Figure 10

Same as Fig. 8, but for the haloes of two galaxies simulated including star formation and feedback. One of the simulations form a galaxy with a prominent stellar disc (curve labelled ‘S’). This galaxy has a particularly quiet dynamical formation history and has been analysed in detail by Abadi et al. (2003a). Another (curve labelled ‘E’) corresponds to a galaxy that has undergone a late major merger so that most stars are in a spheroid resembling an elliptical galaxy (see e.g. Meza et al. 2003). The ‘S’ halo contraction is in reasonable agreement with the prediction of ‘Contra’ (Gnedin et al. 2004). The ‘E’ halo contraction, on the other hand, is much closer to the results presented in this paper (shown by the filled circles and thick upper line). This suggests that the merger history of a galaxy plays an important role in the resulting halo contraction. Crosses and open circles indicate r_conv and ε, as in Fig. 8.

Interestingly, the halo contraction differs significantly in each simulation. The ‘elliptical’ galaxy halo contracts much less than the ‘spiral’ galaxy halo. Barring numerical artefact, these results seem to suggest that the halo response does not depend solely on the initial and final distribution of baryons. Not only how much baryonic mass has been deposited at the centre of a halo matters, but also on the mode of its deposition. It is certainly plausible that central galaxies assembled through merging of dense subclumps may lead to different halo response than galaxies assembled through a smooth flow of baryons to the centre. Mergers of baryonic subsystems may in principle pump energy into the dark halo (through dynamical friction), altering its central structure and softening its contraction. This possibility has been argued before (see e.g. El-Zant, Shlosman & Hoffman 2001; Ma & Boylan-Kolchin 2004; Mo & Mao 2004; Jardel & Sellwood 2009; Pedrosa, Tissera & Scannapieco 2010), and seems to find favour in our simulations, where mergers between massive baryonic clumps are frequent and plentiful and halo contraction is less strong than reported in earlier work.

Recent work has also speculated that mergers may lead to halo expansion and that this would help to reconcile the properties of disc galaxies with ΛCDM haloes (Dutton et al. 2007). Our haloes always contract, and it seems safe to conclude that our results reflect the maximum effect of mergers on the halo response. We conclude therefore that mergers alone are unlikely to result in halo expansion. If such expansion is truly needed to reconcile disc properties with ΛCDM haloes, it should come as a consequence of other processes not considered here, such as feedback-driven winds that may remove a substantial fraction of baryons from the central galaxy (see e.g. Navarro, Eke & Frenk 1996a; Babul & Ferguson 1996).

4 SUMMARY

We have used a suite of cosmological N-body/gasdynamical simulations to examine the modifications to the dark halo structure that result from the assembly of a central galaxy. The formation of 13 ΛCDM haloes is simulated twice, with and without a baryonic (gaseous) component. For simplicity, the gasdynamic simulations include radiative cooling but neglect star formation and feedback. This favours the formation of massive central baryonic discs at the centre of the early collapsing progenitors of the final halo. As these systems merge, the gaseous clumps merge and re-form a disc at the centre of the remnant.

The merger process leads to the transfer of a large fraction of the angular momentum from the baryons to the halo. At z= 0, the simulated central galaxies are too massive and too small to be consistent with observed spiral galaxies. Although unrealistic as a disc galaxy formation model, these simulations allow us to probe the dark halo response in the interesting case where the deepening of the potential well resulting from the formation of the central galaxy is maximized. Our main conclusions may be summarized as follows.

• Dark haloes become significantly more axisymmetric as a result of the assembly of the central galaxy. The triaxial, nearly prolate systems that form in the absence of a baryonic component are transformed into essentially oblate systems with a roughly constant axial isopotential ratio 〈c/a〉≈ 0.85. This confirms hints of earlier work in the sense that more massive/concentrated central galaxies lead to rounder haloes. Since our simulations result in unrealistically massive galaxies, the question of how triaxial the typical halo around a normal galaxy truly is remains pending.

• Haloes always contract in response to the formation of the central galaxy. The ‘adiabatic-contraction’ formalism overestimates the halo contraction in our simulations. The discrepancy increases towards the centre, where the effect of baryons is larger. A simple empirical formula (equation 3) describes well our numerical results, and it provides a useful estimate of the importance of halo contraction in the extreme case where a concentrated galaxy forms through a sequence of major mergers.

• The halo contraction in our simulations is less pronounced than found in the earlier numerical work (see e.g. Gnedin et al. 2004), and suggests that the response of a halo does not depend solely on the final mass and radial distribution of baryons in the central galaxy, but also on the mode of their assembly.

• We apply these results to the MW, where accurate estimates of the mass of baryons and dark matter inside the solar circle exist. These allow us to probe the range of halo virial mass and concentration consistent with the constraints. Only haloes of unusually low mass and concentration would match such constraints if ‘adiabatic contraction’ holds. This restriction is less severe if one uses the halo contraction reported here (equation 3): although few ΛCDM haloes of virial velocity ∼ 220 km s⁻¹ would be eligible hosts for the MW, the situation improves for lower virial velocities. Haloes of average concentration with virial velocity as large as V_vir= 135 km s⁻¹ would be consistent with the MW constraints.

The dark halo response to the formation of a galaxy seems inextricably linked to the full coupled evolution of baryons and dark matter and may vary from system to system. Progress in our understanding of the distribution of dark matter in a baryon-dominated system will thus likely be made in steps with our understanding of the particular assembly history of each individual system.

1

We define the virial radius, r_vir, of a system as the radius of a sphere of mean density Δ_vir(z) times the critical density for closure. This definition implicitly defines the virial mass, M_vir, as that enclosed within r_vir, and the virial velocity, V_vir, as the circular velocity measured at r_vir. Quantities characterizing a system will be measured within r_vir, unless otherwise specified. The virial density contrast, Δ_vir(z), is given by Δ_vir(z) = 18π²+ 82f(z) − 39f(z)², where f(z) =[Ω₀(1 + z)³/(Ω₀(1 + z)³+ Ω_Λ))]− 1 and Ω₀=Ω_CDM+ Ω_b (Bryan & Norman 1998). Δ_vir≈ 100 at z= 0.

We thank James Wadsley, Joachim Stadel and Tom Quinn for allowing us to use their excellent gasoline code for the numerical simulations reported here. We also thank Oleg Gnedin for making publicly available and advice on his excellent ‘Contra’ code. An anonymous referee's report helped us to improve on a previous version of this paper.