NIHAO X: Reconciling the local galaxy velocity function with Cold Dark Matter via mock Hi observations
We used 87 high resolution hydrodynamical cosmological simulations from the NIHAO suite to investigate the relation between the maximum circular velocity () of a dark matter halo in a collisionless simulation and the velocity width of the Hi gas in the same halo in the hydrodynamical simulation. These two quantities are normally used to compare theoretical and observational velocity functions and have led to a possible discrepancy between observations and predictions based on the Cold Dark Matter (CDM) model. We show that below 100 , there is clear bias between Hi based velocities and , that leads to an underestimation of the actual circular velocity of the halo. When this bias is taken into account the CDM model has no trouble in reproducing the observed velocity function and no lack of low velocity galaxies is actually present. Our simulations also reproduce the linewidth - stellar mass (Tully-Fisher) relation and Hi sizes, indicating that the Hi gas in our simulations is as extended as observed. The physical reason for the lower than expected linewidths is that, in contrast to high mass galaxies, low mass galaxies no longer have extended thin Hi rotating disks, as is commonly assumed.
keywords:cosmology: theory – dark matter – galaxies: formation – galaxies: kinematics and dynamics – galaxies: structure – methods: numerical
The Cold Dark Matter (CDM) scenario is very successful in reproducing the large scale structure of the Universe (Springel et al., 2005) and the anisotropies in the Cosmic Microwave Background (Planck Collaboration et al., 2014). In a CDM universe structure formation proceeds in a bottom-up fashion, with small structures collapsing earlier and then merging to form larger and larger haloes. This scenario predicts a larger number of low mass structures compared to more massive ones, or in other words a steeply declining halo mass function. Such a prediction seems to be at odds with observational data both around galaxies, e.g. the satellite abundance (Klypin et al., 1999; Moore et al., 1999) and in the local volume (Cole & Kaiser, 1989; Zavala et al., 2009; Trujillo-Gomez et al., 2011).
In the last years several solutions have been proposed to reconcile the observed number of satellites around the Milky Way and Andromeda galaxies with predictions from CDM (Bullock et al., 2000; Benson et al., 2002; Macciò et al., 2010; Font et al., 2011, and references therein). All these solutions are based (to a different extent) on the fact that galaxy formation becomes quite inefficient in low mass dark matter haloes, due to ionizing background, Super Novae (SN) explosions and gas removal due to ram pressure. As a result, by taking into account baryonic processes it is possible to bring into agreement the abundance of satellites in the Local Group and CDM predictions (e.g., Sawala et al., 2016).
A more persistent challenge is provided by field galaxies in the local volume and their velocity function. The galaxy velocity function (GVF) is defined as the abundance of galaxies with a given circular velocity and is conceptually similar to the galaxy luminosity function. The advantage is that, for a given cosmological model, circular velocities can nominally be predicted much more accurately than luminosities, since they are less dependent on the still poorly understood physics of galaxy formation. This makes velocity functions a very useful tool to test theoretical models (Zavala et al., 2009; Papastergis et al., 2011).
Recently Klypin et al. (2015, hereafter K15) have presented new measurements of the GVF in the Local Volume ( Mpc) and they have shown that while CDM provides very good estimates of the number of galaxies with circular velocities around and above , it fails quite dramatically at lower circular velocities, overestimating by a factor up to five the number of dwarf galaxies in the velocity range . As pointed out by K15 and other authors (e.g., Zavala et al., 2009) galaxies in this velocity range are practically insensitive to the ionization background and, by not being satellites, they are not affected by gas depletion via ram pressure or by stellar stripping. This makes the mismatch between the observed GVF and the CDM predictions a quite serious problem for the current cosmological model. Moreover, simple modifications to the CDM paradigm, for example by introducing a warm dark matter component (e.g., Schneider et al., 2014), seem also to be not able to fix this issue, as shown by K15.
Brook & Di Cintio (2015) showed that a population of galaxies with mass profiles modified by baryonic feedback (Di Cintio et al., 2014) is able to explain the GVF significantly better than a model in which a universal cuspy density profile is assumed. In this Letter we revise the effect of galaxy formation on the GVF using the hydrodynamical simulations from the NIHAO project (Wang et al., 2015). In contrast to many previous theoretical studies, which use proxies for Hi linewidths, we directly measure them. The combination of high spatial resolution and large sample size make NIHAO ideally suited to study the relation between the (rotation) velocity inferred by the kinematics of the Hi gas component and the maximum circular velocity of the halo, and thus to shed light on the reasons behind the large discrepancy between the CDM-based predicted and the observed GVF.
The NIHAO project is a suite of fully cosmological hydrodynamical “zoom-in” simulations that covers a large range of galaxy masses from dwarf galaxies to massive spirals (halo virial velocity from to ). One characteristic of NIHAO galaxies is that they have roughly the same relative resolution and hence the same number of dark matter particles within the virial radius (See Fig. 2 of Wang et al., 2015). The zoomed initial conditions were created using a modified version of grafic2 (Bertschinger, 2001; Penzo et al., 2014). Each halo is initially simulated at high resolution with DM-only, using pkdgrav (Stadel, 2001) then re-simulated with baryons. We refer to the DM-only simulations as N-body and the simulations with baryons as hydro or NIHAO.
The hydrodynamical simulations are evolved using an improved version of the SPH code gasoline (Wadsley et al., 2004; Keller et al., 2014). The code includes a subgrid model for turbulent mixing of metals and energy , heating and cooling include photo-electric heating of dust grains, ultraviolet (UV) heating and ionization and cooling due to hydrogen, helium and metals (Shen et al., 2010). The star formation and feedback modeling follows what was used in the MaGICC simulations (Stinson et al., 2013), adopting a threshold for star formation of cm. Stars can feed energy back into the ISM via blast-wave supernova (SN) feedback (Stinson et al., 2006) and via ionizing radiation from massive stars (early stellar feedback) before they turn into SN (Stinson et al., 2013). Metals are produced by type II and type Ia SN. These, along with stellar winds from asymptotic giant branch stars also return mass to the ISM. The fraction of stellar mass that results in SN and winds is determined using the Chabrier (2003) stellar initial mass function. We refer the reader to Wang et al. (2015) for a more detailed description of the code and the simulations.
2.1 NIHAO galaxies
Since our aim is to study the impact of baryons on the galaxy velocity function it is very important to use realistic simulated galaxies. In this respect, galaxies from the NIHAO project perform extremely well being consistent with a wide range of galaxy properties both in the local and distant universe: evolution of stellar to halo masses and star formation rates (Wang et al., 2015); cold and hot gas content (Stinson et al., 2015; Wang et al., 2016; Gutcke et al., 2016); stellar disk kinematics (Obreja et al., 2016); shallow inner dark matter density slopes of dwarf galaxies (Tollet et al., 2016); and resolve the too-big-to-fail problem for field galaxies (Dutton et al., 2016). They thus provide plausible, even though not unique, templates for the relation linking the dark matter circular velocity and the velocity from Hi dynamics.
2.2 Hi calculation
In calculating the neutral hydrogen Hi fraction we use the self shielding approximation described in Rahmati et al. (2013), based on full radiative transfer simulations presented in Pawlik & Schaye (2011). The self shielding approximation defines a density threshold above which the photoionization of Hi is suppressed. We calculate this threshold by interpolating between the redshift values of our fixed UV background based on Haardt & Madau (2001) and then using Table 2 in Rahmati et al. (2013). The overall effect of this self-shielding approach is to increase the amount of Hi (relative to the fiducial calculation in gasoline) bringing the simulations in better agreement with observations (Rahmati et al., 2013; Gutcke et al., 2016).
3 The observational data
In this Letter we will use the velocity function as computed by K15 using the Updated Nearby Galaxy Catalogue from Karachentsev et al. (2013) which contains redshift independent distances Mpc for 869 galaxies. The velocities have been obtained from Hi velocity widths for the large majority of the galaxies () while for galaxies without Hi line-width measurements they have assigned a line of sight velocity measurement using the average luminosity velocity relation in the K band for galaxies with measured stellar velocity dispersions. The resulting GVF is in reasonably good agreement with previous pure Hi results from HIPASS (Zwaan et al., 2010) and ALFALFA (Papastergis et al., 2011) in the velocity range 25-150 where they overlap. In order to better compare with simulations K15 used (called in K15) as observational velocity, which is defined as the width of the velocity distribution at half height () divided by two. We will use the original datapoints from K15 as the observational velocity function in the rest of this Letter.
Observed rotation velocities are usually based on wide area, 21 cm surveys. Thanks to their spectroscopic nature, Hi surveys automatically provide the spectral Hi line profile of every detected source. The velocity width () of each galaxy can be easily extracted and it is then possible to construct a velocity function (e.g., Papastergis et al., 2011).
As a first step we have replicated this observational process in our galaxies. Fig. 1 shows the velocity distribution of the Hi gas in three different galaxies with increasing mass from left to right. Since the Hi velocity depends on the galaxy inclination, we show for each galaxy three different orientations that yield the maximum (upper row), median (middle row) or minimum (lower row) velocity width. The two horizontal lines indicate two commonly used velocity definitions: the width at 50 per cent () and 20 per cent () of the maximum height. We calculate the maximum height as the average between the maximum flux on the positive and negative velocity sides. As expected, low mass galaxies tend to have a single peaked distribution, which is roughly Gaussian, while more massive ones show clear signs of rotational support and have double peaked Hi profiles.
For each galaxy the dark matter maximum circular velocity (defined as the maximum of ) of the corresponding N-body run) is indicated in the top left corner of every panel (black), together with the value of the rotational velocity of the galaxy (red), defined as half of the width at half high (). One can immediately notice that while for massive galaxies the two velocities ( and ) are quite similar, strong differences are present at lower masses, where the velocity from Hi kinematics can be even less than 20 per cent of circular velocity of the DM. This trend with mass is clearly visible in Fig. 2, where we plot the two velocities (, ) one against the other111The actual in the hydro run might be different from the one from the Nbody simulation, mainly due to (baryonic) mass loss in low mass haloes (e.g., Sawala et al., 2016). On the other hand in this Letter we are interested in comparing previous results from Nbody simulations with actual results from hydro runs and hence we are more interested a comparison of with .. Each galaxy is seen from one hundred different lines of sight, the (red) squares represent the median value among the hundred projections and the error bar is obtained by connecting the maximum and minimum values. Our results are very consistent with recent ones obtained by Vandenbroucke et al. (2016, see their figure 5) despite their use of a different code and a different feedback implementation.
Large N-body simulations such as the Millennium (Springel et al., 2005) and the Bolshoi (Klypin et al., 2011) runs can be used to compute the dark matter halo velocity distribution function, which provides the number of haloes in a given range of . Such a velocity function can be well approximated by a single power law, and Klypin et al. (2011) using the Bolshoi simulation suggested the following parameterization for a Planck cosmology:
The NIHAO suite is based on high-resolution simulations of single galaxies, for this reason we need to assign an appropriate weight for each halo in order to compute a velocity function. We proceeded as follows: we first computed the number of haloes in bins of () using their N-body counterpart (for a total of 87 objects). We then associated to each bin a weight, calculated as the ratio between our results and the global ones from K15 (eq. 1). Finally, we apply these weights to the whole sample of from the hydro simulations (meaning to all 100 different projections for each galaxy, for a total of 8700 values for ) resulting in a velocity function, which can be directly compared with the observations. We assign an error to each point based on the Poisson noise of that velocity bin.
Results are presented in Fig. 3, which is the key plot of this Letter: in black we show the observed velocity function (points) and its analytic fit via Eq. 11 in K15, in blue the dark matter halo velocity function from the Bolshoi simulations (see K15), and finally in red the velocity function as predicted by the NIHAO simulations accounting for the effects of galaxy formation. This figure clearly shows that, once the mass dependent offset between the Hi velocity and the maximum velocity of the halo is taken into account, hydrodynamical simulations based on the CDM model have no problems in reproducing the observed velocity function.
The NIHAO results seem to fall a bit short at the low velocity end when compared with the observations corrected for incompleteness (black filled circles). On the other hand they agree well with the raw data (open circles). This might indicate that the Hi gas in the NIHAO simulations is too turbulent, or is not extended enough (but see next section) or that the incompleteness correction applied by K15 was a bit too generous. Deeper future surveys such as with the Square Kilometer Array (Dewdney et al., 2009) will help in clarifying this issue.
4.1 Hi distribution in NIHAO
To validate our GVF results we need to show that NIHAO galaxies have realistic Hi velocity distributions and radial extents. To address this point we show several scaling relations. The upper left panel of Fig. 4 shows the relation between stellar mass and Hi velocity in NIHAO (red squares) and in observations as measured by Bradford et al. (2016). The velocities from Bradford et al. (2016) are corrected for inclination. To compute an equivalent quantity in the NIHAO galaxies for each simulation we find the median linewidth from the 100 projections and then divide it by , the median inclination correction for a random distribution of galaxies in the sky. We have dubbed this average velocity . The Tully-Fisher relations from the NIHAO simulations are in very good agreement with observation. The lower left panel of Fig. 4 shows the relation between the NIHAO linewidths measured at 50 per cent and 20 per cent of the peak flux (grey points for all 100 projections per galaxy and red squares for the median velocities). At velocities () the two linewidths are similar, reflecting the double peaked nature and steep wings of the line profiles. However, at lower velocities decreases following and at reaches a value consistent with that of a Gaussian . This change in line profile shape contributes to deviating from at low velocties (Brook & Shankar, 2016). The good agreement between NIHAO simulations and observations from Bradford et al. (2016) (solid black curve with 1 scatter) shows that NIHAO simulations have realistic Hi line profiles.
The NIHAO galaxies are also able to reproduce the relation between the total mass in Hi and its radial extent, as shown in the right panels of Fig. 4. In the upper one we compare the Hi mass with the scale radius of its radial distribution. This radius has been obtained by fitting an exponential profile to the face-on Hi surface density between the radii containing 10 per cent and 90 per cent of the Hi mass. In the bottom right panel we use the radius where the Hi surface density drops below 1 /pc (see Swaters, 1999). The NIHAO galaxies (in red) are able to reproduce the observed Hi distributions (black points) for both definitions of the Hi radius.
5 Discussion and Conclusions
The abundance of galaxies with a given velocity linewidth () provides a very powerful tool to test against observations different models for Dark Matter, from Cold (CDM) to Warm (WDM). The slope of this velocity function at low velocities () departs significantly from the expectations of the standard CDM model, leading to a difference in abundance of about a factor at (Papastergis et al., 2011, ALFALFA survey). Such a discrepancy is still in place when deeper observations of the local ( Mpc) Universe are used (K15, based on the survey from Karachentsev et al., 2013); moreover a simple cut in the matter power spectrum, as in WDM models, does not provide a viable solution (K15). These kinds of analyses suffer from the fact that Hi discs in low mass galaxies might not be extended enough to probe the full amplitude of the rotation curve (Brook & Shankar, 2016).
We reinvestigated this problem using the NIHAO simulation suite to study in detail the relation between the maximum circular velocity of Dark Matter haloes, , and the Hi velocity. We found that at high masses the Hi rotation curves are extended enough to fully capture the dynamical velocity of the dark matter, and hence the linewidth based velocities, , can be directly compared with velocities coming from dissipationless (N-body) simulations. At lower masses, namely , there is a clear bias between the circular velocity of the halo and the width of the Hi velocity distribution. Such a bias is due to both the Hi distribution not being extended enough to reach the peak of the circular velocity and to the departure of the Hi dynamics from a purely rotating disc, due to substantial dispersion in the vertical direction. When such a bias is taken into account, the simulated Hi velocity function is able to match the observed ones well, showing that the CDM model is not intrinsically flawed, but that attempts to realistically “observe” fully hydrodynamical simulations might alleviate apparent tensions between observations and theory predictions. In conclusions any comparisons between pure N-body simulations and observations must be taken with a grain of salt.
We thank Anatoly Klypin and Jeremy Bradford for sending us their results in electronic form. This research was carried out on the High Performance Computing resources at New York University Abu Dhabi; on the theo cluster of the Max-Planck-Institut für Astronomie and on the hydra clusters at the Rechenzentrum in Garching. XK is supported by the NSFC (No.11333008) and the strategic priority research program of CAS (No. XDB09000000).
- Benson et al. (2002) Benson A. J., Lacey C. G., Baugh C. M., Cole S., Frenk C. S., 2002, MNRAS, 333, 156
- Bertschinger (2001) Bertschinger E., 2001, ApJS, 137, 1
- Bradford et al. (2016) Bradford J. D., Geha M. C., van den Bosch F. C., 2016, ArXiv e-prints, 1602.02757,
- Brook & Di Cintio (2015) Brook C. B., Di Cintio A., 2015, MNRAS, 453, 2133
- Brook & Shankar (2016) Brook C. B., Shankar F., 2016, MNRAS, 455, 3841
- Bullock et al. (2000) Bullock J. S., Kravtsov A. V., Weinberg D. H., 2000, ApJ, 539, 517
- Chabrier (2003) Chabrier G., 2003, PASP, 115, 763
- Cole & Kaiser (1989) Cole S., Kaiser N., 1989, MNRAS, 237, 1127
- Dewdney et al. (2009) Dewdney P. E., Hall P. J., Schilizzi R. T., Lazio T. J. L. W., 2009, IEEE Proceedings, 97, 1482
- Di Cintio et al. (2014) Di Cintio A., Brook C. B., Dutton A. A., Macciò A. V., Stinson G. S., Knebe A., 2014, MNRAS, 441, 2986
- Dutton et al. (2016) Dutton A. A., Macciò A. V., Frings J., Wang L., Stinson G. S., Penzo C., Kang X., 2016, MNRAS, 457, L74
- Font et al. (2011) Font A. S., et al., 2011, MNRAS, 417, 1260
- Gutcke et al. (2016) Gutcke T. A., Stinson G. S., Macciò A. V., Wang L., Dutton A. A., 2016, ArXiv e-prints, 1602.06956,
- Haardt & Madau (2001) Haardt F., Madau P., 2001, in Neumann D. M., Tran J. T. V., eds, Clusters of Galaxies and the High Redshift Universe Observed in X-rays. (arXiv:astro-ph/0106018)
- Karachentsev et al. (2013) Karachentsev I. D., Makarov D. I., Kaisina E. I., 2013, AJ, 145, 101
- Keller et al. (2014) Keller B. W., Wadsley J., Benincasa S. M., Couchman H. M. P., 2014, MNRAS, 442, 3013
- Klypin et al. (1999) Klypin A., Kravtsov A. V., Valenzuela O., Prada F., 1999, ApJ, 522, 82
- Klypin et al. (2011) Klypin A. A., Trujillo-Gomez S., Primack J., 2011, ApJ, 740, 102
- Klypin et al. (2015) Klypin A., Karachentsev I., Makarov D., Nasonova O., 2015, MNRAS, 454, 1798
- Macciò et al. (2010) Macciò A. V., Kang X., Fontanot F., Somerville R. S., Koposov S., Monaco P., 2010, MNRAS, 402, 1995
- Moore et al. (1999) Moore B., Ghigna S., Governato F., Lake G., Quinn T., Stadel J., Tozzi P., 1999, ApJ, 524, L19
- Obreja et al. (2016) Obreja A., Stinson G. S., Dutton A. A., Macciò A. V., Wang L., Kang X., 2016, MNRAS, 459, 467
- Papastergis et al. (2011) Papastergis E., Martin A. M., Giovanelli R., Haynes M. P., 2011, ApJ, 739, 38
- Pawlik & Schaye (2011) Pawlik A. H., Schaye J., 2011, MNRAS, 412, 1943
- Penzo et al. (2014) Penzo C., Macciò A. V., Casarini L., Stinson G. S., Wadsley J., 2014, MNRAS, 442, 176
- Planck Collaboration et al. (2014) Planck Collaboration et al., 2014, A&A, 571, A16
- Rahmati et al. (2013) Rahmati A., Pawlik A. H., Raicevic M., Schaye J., 2013, MNRAS, 430, 2427
- Sawala et al. (2016) Sawala T., et al., 2016, MNRAS, 457, 1931
- Schneider et al. (2014) Schneider A., Anderhalden D., Macciò A. V., Diemand J., 2014, MNRAS, 441, L6
- Shen et al. (2010) Shen S., Wadsley J., Stinson G., 2010, MNRAS, 407, 1581
- Springel et al. (2005) Springel V., et al., 2005, Nature, 435, 629
- Stadel (2001) Stadel J. G., 2001, PhD thesis, UNIVERSITY OF WASHINGTON
- Stinson et al. (2006) Stinson G., Seth A., Katz N., Wadsley J., Governato F., Quinn T., 2006, MNRAS, 373, 1074
- Stinson et al. (2013) Stinson G. S., Brook C., Macciò A. V., Wadsley J., Quinn T. R., Couchman H. M. P., 2013, MNRAS, 428, 129
- Stinson et al. (2015) Stinson G. S., et al., 2015, MNRAS, 454, 1105
- Swaters (1999) Swaters R. A., 1999, PhD thesis, , Rijksuniversiteit Groningen, (1999)
- Tollet et al. (2016) Tollet E., et al., 2016, MNRAS, 456, 3542
- Trujillo-Gomez et al. (2011) Trujillo-Gomez S., Klypin A., Primack J., Romanowsky A. J., 2011, ApJ, 742, 16
- Vandenbroucke et al. (2016) Vandenbroucke B., Verbeke R., De Rijcke S., 2016, MNRAS, 458, 912
- Wadsley et al. (2004) Wadsley J. W., Stadel J., Quinn T., 2004, New Astron., 9, 137
- Wang et al. (2015) Wang L., Dutton A. A., Stinson G. S., Macciò A. V., Penzo C., Kang X., Keller B. W., Wadsley J., 2015, MNRAS, 454, 83
- Wang et al. (2016) Wang L., Dutton A. A., Stinson G. S., Macciò A. V., Gutcke T., Kang X., 2016, ArXiv e-prints, 1601.00967,
- Zavala et al. (2009) Zavala J., Jing Y. P., Faltenbacher A., Yepes G., Hoffman Y., Gottlöber S., Catinella B., 2009, ApJ, 700, 1779
- Zwaan et al. (2010) Zwaan M. A., Meyer M. J., Staveley-Smith L., 2010, MNRAS, 403, 1969