The mass discrepancy acceleration relation in a ΛCDM context

Abstract

The mass discrepancy acceleration relation (MDAR) describes the coupling between baryons and dark matter (DM) in galaxies: the ratio of total-to-baryonic mass at a given radius anticorrelates with the acceleration due to baryons. The MDAR has been seen as a challenge to the Λ cold dark matter (ΛCDM) galaxy formation model, while it can be explained by Modified Newtonian Dynamics. In this Letter, we show that the MDAR arises in a ΛCDM cosmology once observed galaxy scaling relations are taken into account. We build semi-empirical models based on ΛCDM haloes, with and without the inclusion of baryonic effects, coupled to empirically motivated structural relations. Our models can reproduce the MDAR: specifically, a mass-dependent density profile for DM haloes can fully account for the observed MDAR shape, while a universal profile shows a discrepancy with the MDAR of dwarf galaxies with M^⋆ < 10^9.5 M_⊙, a further indication suggesting the existence of DM cores. Additionally, we reproduce slope and normalization of the baryonic Tully–Fisher relation (BTFR) with 0.17 dex scatter. These results imply that in ΛCDM (i) the MDAR is driven by structural scaling relations of galaxies and DM density profile shapes, and (ii) the baryonic fractions determined by the BTFR are consistent with those inferred from abundance-matching studies.

galaxies: evolution, cosmology: theory, dark matter

1 INTRODUCTION

Disc galaxies show a tight coupling between their observed dynamics and the distribution of baryons (e.g. Sancisi 2004; Swaters et al. 2012; Lelli, Fraternali & Verheijen 2013). This is quantified by the mass discrepancy acceleration relation (MDAR). Following McGaugh (2004, 2014) the mass discrepancy is defined as |$V_{\rm obs}^2(r)/V_{\rm bar}^2(r)$|⁠, where |$V_{\rm obs}^2$| is the total observed rotational velocity and |$V_{\rm bar}^2$| is the rotation velocity inferred from the distribution of baryons (gas and stars). For a spherical mass distribution, |$V_{\rm obs}^2(r)/V_{\rm bar}^2(r)=M_{\rm tot}(r)/M_{\rm bar}(r)$|⁠.

McGaugh (2004, 2014) finds that a very tight correlation exists between the mass discrepancy and the gravitational acceleration due to baryons |$g_{\rm bar}=V_{\rm bar}^2(r)/r$|⁠. The existence of the MDAR has been difficult to reconcile with a standard ΛCDM picture, hinting perhaps towards a fundamental law of nature (McGaugh 2014) and finding a natural explanation in Modified Newtonian Dynamics (MOND) theory (Milgrom 1983), which posits a characteristic acceleration scale. Some authors have previously pointed out that a correlation between mass discrepancy and baryonic acceleration can be expected in ΛCDM once dissipative processes of galaxy formation are properly accounted for. van den Bosch & Dalcanton (2000) used a semi-analytic model for the formation of galaxies, tuned to fit the Tully–Fisher relation, to show that a characteristic acceleration can be reproduced in ΛCDM models. Yet a more challenging exercise is to verify whether simulations and models including baryon physics can reproduce the small scatter about the observed MDAR (Walker & Loeb 2014). Recently, Santos-Santos et al. (2016) have shown that a suite of simulated galaxies from the MaGICC project (Brook et al. 2012; Stinson et al. 2013) can successfully reproduce both the MDAR and its scatter. Those simulations also reproduce observed galaxy scaling relations and give rise to the expected baryonic Tully–Fisher relation (BTFR), i.e. the relation between baryonic mass and galaxy velocity (McGaugh 2012).

In this work, we explore the issue using our current knowledge of observations and theory: we derive the MDAR and the BTFR using a semi-empirical-ΛCDM-based model and we compare them against a data set of 146 galaxies (Lelli et al., in preparation). We build a population of model galaxies whose distribution in stellar mass and other structural properties follows closely the observational sample, allowing us to explore the MDAR and BTFR. To this aim, we use empirical and theoretical scaling relations superimposed on dark matter (DM) haloes whose density profiles are predicted in ΛCDM with and without baryonic effects.

For the DM distribution we employ the universal NFW (Navarro, Frenk & White 1996) profile and the mass-dependent density profile DC14 (Di Cintio et al. 2014a,b). The latter takes into account the core formation process resulting from impulsive outflows of gas in star-forming region, which has been observed in hydrodynamical simulations at z=0 and across cosmic time (Governato et al. 2010; Teyssier et al. 2013; Di Cintio et al. 2014a; Chan et al. 2015; Oñorbe et al. 2015; Tollet et al. 2015; Trujillo-Gomez et al. 2015), and confirmed by analytic models (Pontzen & Governato 2012).

The DC14 profile has a peak in core formation efficiency for galaxies of M^⋆ ∼ 10^8.5 M_⊙, a result largely confirmed by the increased sample of galaxies in Tollet et al. (2015) as well as by Chan et al. (2015). Such profile has been applied to a variety of studies, having proven its validity in the context of velocity function of galaxies, Local Group predictions and rotation curves shapes (Brook 2015; Brook & Di Cintio 2015a,b). Once a DM halo is set, stellar and gaseous components are added using observed galaxy scaling relations. Theoretical rotation curves and mass profiles are then derived, allowing us to evaluate both the MDAR and the BTFR and to compare them with observations.

2 METHODS

2.1 Observational sample

The observed MDAR is derived from the SPARC (Spitzer Photometry & Accurate Rotation Curve) data set (Lelli et al., in preparation). For this study, we consider 146 galaxies with high-quality H i/H α rotation curves and inclination angles above 30°, in order to exclude face-on galaxies that require large inclination corrections on V_obs. The mass models will be presented in detail in Lelli et al. (in preparation). In short, the baryonic contribution to the rotation curve is calculated from the observed surface density profiles of stars and gas by solving the Poisson equation in cylindrical coordinates (Casertano 1983). The amplitude of V_star is scaled to the total stellar mass assuming a constant mass-to-light ratio of 0.5 at 3.6 μm for all galaxies (McGaugh & Schombert 2014). Similarly, the amplitude of V_gas is scaled to the total gas mass as |$\sqrt{1.33} M_{\rm H\,\small {i}}$|⁠, where |$M_{\rm H\,\small {i}}$| is the observed H i mass and the factor 1.33 takes the contribution of Helium into account. The galaxies cover a broad range in baryonic masses, from gas-dominated dwarfs with M^⋆ ∼ 10⁷ M_⊙ to star-dominated spirals with M^⋆ ∼ 10¹¹ M_⊙. It also spans a broad range in surface brightnesses, including many low-surface-brightness galaxies that are often misrepresented in optical surveys (e.g. McGaugh, Bothun & Schombert 1995). The statistical properties of the observational sample are shown in Fig. 1 as red histograms.

Figure 1.

Properties of the observational (red) and theoretical (black) sample of galaxies. From left to right we show the distribution in halo, stellar and gas mass, halo concentration, baryonic mass, gas fraction, bulge-to-total ratio, central surface brightness and stellar disc scalelength.

2.2 Theoretical model

We build ∼1500 model galaxies such that their distribution in stellar mass follows the observed one, as shown in the top-middle panel of Fig. 1. Each stellar mass bin is populated with 10× the observed number of galaxies, in order to have better statistic. The observational histograms (red) are rescaled to be compared with our theoretical sample. Each galaxy is then assigned to a DM halo using the abundance matching relation of Kravtsov, Vikhlinin & Meshscheryakov (2014) with a scatter of 0.10 dex. Such abundance matching is formally identical to the Moster, Naab & White (2013) one at the low mass end, but different for M_halo >10^11.5 M_⊙, such that a galaxy of M^⋆ ∼ 2 ×10¹¹ M_⊙ is assigned to a 10¹³ M_⊙ halo rather than to a 10¹⁴ M_⊙ one. The steeper slope of this relation at high masses, based on stellar mass function calibrations from improved photometry (Bernardi et al. 2013), shows an excellent agreement with recent hydrodynamical simulations (Tollet et al. 2015). We find that the less massive haloes implied by the Kravtsov et al. (2014) relation better reproduce the rotation curves of galaxies with M^⋆ ≳ 5× 10¹⁰ M_⊙ than the Moster et al. (2013) relation does. The resulting distribution in halo masses is shown in the top-left panel of Fig. 1 and covers the range 9.60 < log₁₀ (M_halo/ M_⊙) < 13.35.

Our model galaxies live in DM haloes that follow the concentration–mass relation from N-body simulations, assuming Planck cosmology (Dutton & Macciò 2014):  

\begin{equation} {\rm log}_{10}(c_{200})=0.905-0.101{\rm log}_{10}(M_{200}/[10^{12}h^{-1}{\rm M}_{{\odot }}]) \end{equation}

(1)

where c₂₀₀ = R₂₀₀/r_s, with a scale radius r_s, and the mass is defined as M₂₀₀ = M_halo = R₂₀₀ΔH²/2G, with Δ ∼ 200. The c − M_halo relation has an intrinsic scatter of 0.11 dex.

DM haloes are described by two different density profiles: the well-known universal, centrally-cuspy NFW one (Navarro et al. 1996), as found in collisionless simulations, and the DC14 model (Di Cintio et al. 2014a,b), derived from hydrosimulations of galaxy formation that take into account a realistic treatment of baryonic physics (Brook et al. 2012; Stinson et al. 2013). The DC14 profile is a double-power-law one, whose main parameters (namely inner, γ, and outer, β, slope and sharpness of the transition, α) depend on the integrated stellar-to-halo mass ratio of each galaxy, according to equation (3) of Di Cintio et al. (2014a): it is able to describe at the same time cuspy (γ ≥ 1) as well as cored (γ < 1) profiles, providing a useful parametrization for DM haloes whose inner structure have ben modified by baryonic processes. Once (α, β, γ) are specified via the quantity M^⋆/M_halo, the DC14 profile has the same number of free parameters as the NFW model, i.e. the total halo mass M_halo and the scale radius r_s.

Dwarf galaxies with M^⋆/M_halo < 10⁻⁴ do not have sufficient energy from their stellar component to generate powerful gas outflows that cause the central DM to expand, hence they retain the initial NFW profile. More massive galaxies in the range M^⋆/M_halo ∼ 3–5 × 10⁻³ will instead have a central distribution of DM shallower than NFW: the core formation mechanism is most efficient at such mass range. Galaxies with M^⋆/M_halo ≥ 0.05, like the Milky Way (MW), will return to have a cuspy NFW profile due to the inefficient effect of SN feedback in contrasting the deep gravitational potential. At the MW mass, a contraction is applied as found in cosmological simulations (Di Cintio et al. 2014a; Tollet et al. 2015). Since the DC14 profile has only been tested up to M_halo= 10¹² M_⊙, we assume that more massive haloes retain an NFW profile with usual concentration.

For each model galaxy, the gas mass is assigned using the empirical |$M_{\rm H\,\small {i}}-M^{\star }$| scaling relation obtained by Papastergis et al. (2012), which made use of the SDSS (Abazajian et al. 2009) and the ALFALFA (Haynes et al. 2011) survey to determine the atomic hydrogen versus stellar mass of galaxies with a scatter of ∼0.20 dex.:  

\begin{equation} {\rm log}_{10}(M_{\rm H\,\small {i}}/M^{\star })=-0.43{\rm log}_{10}(M^{\star })+3.75 \end{equation}

(2)

The total gas mass is given by M_gas = 1.33|$M_{\rm H\,\small {i}}$|⁠, neglecting the minor contribution of molecular gas as in the SPARC mass models.

Stars and gas are distributed in a disc with an exponential profile. The stellar disc scalelength, |$R_{\rm dis{\rm c}}$|⁠, is assigned using the mass–size relation obtained by the GAMA survey in the K-band for late-type galaxies at z < 0.1 (Lange et al. 2015):  

\begin{equation} R_{\rm eff} = 0.13(M^{\star })^{0.14}\left(1.0+\frac{M^{\star }}{14.03\times 10^{10}\,\mathrm{M}_{{\odot }}}\right)^{0.77} \end{equation}

(3)

where R_eff = 1.678|$R_{\rm dis{\rm c}}$| for en exponential disc. This equation is derived considering galaxies with M^⋆ >10⁹ M_⊙ and its extrapolation to lower masses does not correctly reproduce the properties of our sample. Hence, for M^⋆ <10⁹ M_⊙, we use the best-fitting relation to the dwarf galaxies in our sample: log₁₀(R_eff) = (log₁₀(M^⋆) − 7.79)/2.47. The scatter in both relations is ∼0.20 dex. The gas disc scalelength is R_gas = 2R_disc.

Galaxies with M^⋆ > 10^9.5 M_⊙ have a central bulge which follows a spherically symmetric Hernquist (1990) profile:  

\begin{equation} \rho (r) = \frac{M_{\rm bulge}R_{1/4}}{2\pi r(r+R_{1/4})^3} \end{equation}

(4)

where R_1/4 is the radius at which the bulge mass is equal to M_bulge/4. The half-mass radius, R_1/2 = |$(1+\sqrt{2})R_{1/4}$|⁠, is gauged from the R_1/2 − M_bulge scaling relation as in Gadotti (2009):  

\begin{equation} \log _{10}(R_{1/2}) = 0.30\log _{10}(M_{\rm bulge})-3.124 \end{equation}

(5)

We impose a linear bulge growth with stellar mass, M_bulge/M^⋆ = (log₁₀(M^⋆) − 9.5)/4.2, such that the bulge-to-total ratio (B/T) is ∼30 per cent at the MW mass, and further the relations between B/T versus M^⋆, R_1/2 and Hubble Type match this and other observational samples (Gadotti 2009; Moran et al. 2012).

The previous relations use a Chabrier (2003) initial mass function and units of M_⊙ and kpc. Further, we assume that their scatters are uncorrelated. In Fig. 1 the black histograms show that such semi-empirical model can reproduce very well the statistical properties of the SPARC observational sample of galaxies.

The total circular velocity of a halo following either an NFW or a DC14 density profile, embedding a galaxy with an exponential disc of gas and stars and a central stellar bulge, reads:  

\begin{equation} {V}_{\rm tot}^2(r) = { {V}_{\rm dis{\rm c},g}^2(r) + {V}_{\rm dis{\rm c},\star }^2(r) + {V}_{\rm bulge}^2(r) + {V}_{{\rm DM}}^2(r)} \end{equation}

(6)

where V_DM has been multiplied by a factor |$\sqrt{1-f_b} = \sqrt{0.842}$| to take into account the cosmic baryon fraction. The potential of a thin disc of mass |$M_{\rm dis{\rm c}}$| and scalelength |$R_{\rm dis{\rm c}}$| is (Freeman 1970):  

\begin{equation} {V}_{\rm disc}^2(r) = \frac{GM_{\rm disc}}{R_{{\rm disc}}}2y^2[I_0(y)K_0(y) - I_1(y)K_1(y)] \end{equation}

(7)

where I_n and K_n are modified Bessel functions of the first and second kind, respectively, and y = r/(2R_disc). The potential of the bulge and DM halo are computed assuming spherical symmetry. The velocity and the mass profile M(r) of each galaxy component can then be evaluated at every radius. We remark here that our model generates realistic rotation curves of galaxies at different masses, a result that will be presented in a companion paper.

3 RESULTS

We start by showing that our model is able to reproduce the BTFR, which is the relation between total baryonic mass versus galaxy rotation velocity along the flat part of the rotation curve (V_flat) (McGaugh 2005, 2012). We measure the model rotation curves at 8R_disc, such that we are encompassing the radius within which most of the baryonic mass is found. We note that the radius at which model rotation curves are measured is critical in the definition of the BTFR, resulting in different scatter and slopes at the low-mass end, issue that will be explored in a forthcoming paper.

The model galaxies are shown in Fig. 2 as black squares for a DC14 profile, with the best linear fit indicated as dotted pink line. The resulting slope, intercept and scatter about the average relation, for equally weighted points, are 3.50 ± 0.10, 2.73 ± 0.19 and 0.17 dex. Measuring at 8R_disc produces only a small difference in the BTFR between a DC14 and NFW model (thus not shown here), since this radius is large enough to have reached the point at which the two profiles are similar. The slope, intercept and scatter in the NFW case would read 3.74 ± 0.10, 2.18 ± 0.21 and 0.17 dex.

Figure 2.

BTFR derived by measuring the rotation of model galaxies at |$8R_{\rm dis{\rm c}} = 4R_{\rm gas}$| (black squares). The best fit is indicated as dotted pink line for the DC14 profile. BTFRs from the literature are shown as solid, dashed and dot–dashed line for observations, ΛCDM-based models and hydro simulations, respectively.

Overplotted in Fig. 2 are the best-fitting BTFRs from the most up-to-date observations (McGaugh 2012; Lelli, McGaugh & Schombert 2015) ΛCDM-based theoretical models (Trujillo-Gomez et al. 2011; Dutton 2012) and recent hydrodynamical simulations (Santos-Santos et al. 2016). The slope, intercept, and scatter of the BTFR represent one of the major predictions of MOND (Milgrom 1983), whereas it has been argued that they are difficult to reproduce within a ΛCDM context without fine-tuning (McGaugh 2005, 2012). We show here that a semi-empirical model, coupling observed scaling relations with ΛCDM haloes, can generate a BTFR in good agreement with the observed one. The scatter of 0.17 dex is somewhat high compared to the observational intrinsic scatter (∼0.11 dex, Lelli et al. 2015). Whether this is a challenge for ΛCDM, or is it rather pointing out the sensitivity of the BTFR to the measurements radius, is an open issue which deserves further investigation.

The mass discrepancy, M_tot(r)/M_bar(r), and the corresponding baryonic acceleration, g_bar(r), are computed for each model galaxy within the range r < 4R_gas = 8R_disc, which is a conservative proxy for the outermost radius measured in H i observations. In Fig. 3 we plot such mass discrepancy for the full sample of ∼1500 model galaxies against g_bar, using uniformly spaced bins of 100 pc (different binning criteria produce less, or more, dense plots, without changing the general pattern shown in Fig. 3). In the top two panels the MDAR is shown for a DC14 profile, colour coded according to the radius at which it is measured (left-hand panel), where blue indicates the innermost radii and red the outermost ones, and according to the stellar mass (right-hand panel), from dwarf galaxies in red to massive spirals in blue. The same colour coding scheme is used for an NFW model (bottom panels). Each galaxy is represented by a series of points. Low-mass galaxies overlaps with high mass ones to produce the MDAR relation.

Figure 3.

The MDAR of model galaxies in the mass range 6.8 ≤ log₁₀(M^⋆/ M_⊙) ≤ 11.3 for DC14 (top) and NFW (bottom) profiles, measured out to a maximum radius of r = 8R_disc. The MDAR of each galaxy is colour coded by radius (left) and stellar mass (right). The dashed black line shows the maximum baryonic contribution.

Fig. 3 illustrates that there is a clear difference between DC14 and NFW models. Low mass galaxies with 7 ≲ log₁₀(M^⋆/ M_⊙)≲9 have a shallow, core-like inner profile (γ ≲ 0.5) in the DC14 model, resulting in a lower mass discrepancy in their inner region than the corresponding cuspy NFW case. As a consequence the shape of the MDAR in the DC14 model between −10 ≲ log₁₀(g_bar/cm s⁻²) ≲ −8.5 is different than the NFW one. On the contrary, NFW low-mass galaxies have a central cusp and high values of M_tot(r)/M_bar(r) due to the increased contribution of DM to their centre.

A direct comparison with observations is offered in Fig. 4, whose three panels show the density map of the MDAR generated by our model galaxies using a DC14 profile (top and bottom panel) and a NFW one (central panel). Once again we used uniformly spaced bins of 100 pc, although different bin sizes do not alter our results. Overplotted as filled circles are the observational data referring to the 146 galaxies from the SPARC data set, considering only points with a relative error on the measured rotation curve errV/V≤0.5. Big circles indicate the most accurate data, with errV/V ≤ 0.1, while small circles represent data whose relative error is 0.1 ≤ errV/V ≤ 0.5. The dashed line is the observational fit from McGaugh (2014).

Figure 4.

Density map of the MDAR generated by model galaxies that follow the same stellar mass distribution of the observational sample. Results for the DC14 profile (top and bottom panels) and the NFW one (central panel) are shown, out to a radius r = 8R_disc and using bins of 100 pc. Filled circles represent the observed MDAR derived from the SPARC sample of 146 galaxies. In the bottom panel the DC14 heat plot is shown in black and white, with observational points colour-coded by their measured radius.

Our semi-empirical-ΛCDM-based model is able to fully account for both the shape and the observed scatter in the relation. Comparing the first two panels, it is also clear that the DC14 profile encompasses the full range of the observed MDAR, while the NFW one fails to account for several points in the low acceleration-low mass discrepancy space. This is the region where dwarf galaxies with inner cores are found in the DC14 model (see Fig. 3). To further investigate this issue we show, in the bottom panel of Fig. 4, the same DC14 heat-map with observational points colour-coded according to the radius at which they have been measured. The observational points falling in the bottom-left area of the MDAR plot come from dwarf galaxies with M^⋆ < 10^9.5 M_⊙ and mostly from their inner region, in agreement with DC14 expectations. We predict that dwarf galaxies falling in this area of the MDAR should have a central core-like distribution, which a universal NFW halo model cannot reproduce.

4 CONCLUSIONS

The MDAR, i.e. the relation between total-to-baryonic mass ratio as a function of baryonic acceleration at different radii of a galaxy, has been historically seen as a challenge to the ΛCDM paradigm (McGaugh 2014; Wu & Kroupa 2015). Such tight relation and its small scatter can be instead easily reproduced within a MOND context, in which a characteristic acceleration scale exists (Milgrom 1983).

In this Letter, we have shown that the MDAR can be reproduced using a semi-empirical model in a ΛCDM scenario, using theoretically motivated DM haloes populated with realistic disc galaxies that follow observed structural relations. We created model galaxies whose distribution in stellar, gas, baryonic mass, disc scalelength and central surface brightness closely match the observational SPARC data set (Lelli et al., in preparation), enabling accurate comparisons between theory and observation.

The MDAR generated by galaxies embedded in DM haloes that follow the mass-dependent DC14 profile (Di Cintio et al. 2014a,b) is able to account for the overall shape and scatter of the observed MDAR, from dwarfs to massive spirals. Instead, employing the universal NFW (Navarro et al. 1996) profile fails to reproduce the MDAR generated by dwarf galaxies with M^⋆ ≲ 10^9.5 M_⊙, a further evidence pointing towards the existence of cores in such objects. More accurate data of the inner kpc of dwarf galaxies could better discriminate the two models. Regardless, both models indicate that the MDAR can be reproduced in a ΛCDM context once the observed structural relations of galaxies are taken into account.

Our model can also reproduce the slope and normalization of the BTFR within the observational errors (McGaugh 2012), in line with some ΛCDM semi-analytic models (Dutton 2012; Trujillo-Gomez et al. 2011) and state-of-the-art hydrodynamical simulations (Santos-Santos et al. 2016). This implies that the global baryonic fractions from abundance-matching techniques (e.g. Kravtsov et al. 2014) are roughly consistent with those implied by the BTFR. As in previous models, however, the scatter in the BTFR is slightly higher than the observed intrinsic scatter, an issue that deserves further investigations.

This work extends the results of previous MDAR studies based on the analysis of hydrodynamical simulations (Santos-Santos et al. 2016) and semi-analytic models (van den Bosch & Dalcanton 2000), largely confirming that once the general patterns described by galactic structural relations are properly taken into account, capturing the physics of galaxy formation, the MDAR and the BTFR arise in a ΛCDM context.

ADC is supported by the DARK-Carlsberg independent fellowship programme. ADC thanks N. Amorisco and C. Brook for fruitful discussions. The work of FL is made possible through the support of a grant from the John Templeton Foundation.