© ESO, 2017

1. Introduction

For several decades the dark matter problem has been one of the main topics in astronomical research. The idea of missing or invisible mass was already proposed in the early 1930s by Jan Oort (Oort 1932) and Fritz Zwicky (Zwicky 1933) based on their observations of the motions of stars in the Milky Way disk and galaxies in the Coma cluster. Despite this early notion, the first sound evidence of the presence of dark matter only came in the 1970s, from the analysis of galaxy rotation curves by Freeman (1970), Roberts (1975), and Rubin et al. (1978). These authors found that the rotation curves of massive galaxies remain flat even at large galactocentric distances and well beyond the stellar disks. This could not be explained by the Newtonian gravity of the visible matter alone, but instead implied an additional extended halo of invisible matter. Furthermore it was found that the rotation curves of low mass and low surface brightness (LSB) galaxies show a slow, almost linear rise in the centre after subtraction of the baryonic contributions. To match this observed behaviour, empirical models of the dark matter distribution in galaxies typically have a central constant-density core. These models, such as the pseudo-isothermal sphere, can explain a wide variety of observed rotation curves (e.g. van Albada et al. 1985; Broeils 1992; de Blok et al. 2001), although they have no physical basis.

On the other hand, dark matter only simulations of the structure formation in the Universe consistently find dark matter haloes with a central cusp (e.g. Navarro et al. 1996b; Moore et al. 1999; Klypin et al. 2001; Diemand et al. 2005; Stadel et al. 2009). The dark matter density profiles derived from these simulations, however, give poor fits to the rotation curves of dwarf galaxies (e.g. de Blok et al. 2001; de Blok & Bosma 2002; Weldrake et al. 2003; Gentile et al. 2005, 2007). This cusp-core controversy has been one of the major problems of Λ cold dark matter (ΛCDM) for the past two decades. In the early years rotation curves were often derived in an overly simplistic way from poorly sampled velocity fields or even one-dimensional long-slit observations. Because both physical effects, such as non-circular motions and pressure support, and observational biases, such as beam smearing, often affect the observed kinematics in the central parts of galaxies, there has long been a discussion regarding whether the observed rotation curves are actually reliable and truly trace the gravitational potential of a galaxy (e.g. Swaters et al. 2003; Rhee et al. 2004; Spekkens et al. 2005; Valenzuela et al. 2007). Since the 1990s, however, both the quality of the observations and the analysis techniques have vastly improved and although the discussion still persists today (Pineda et al. 2017), there is growing consensus that modern rotation curves accurately trace the total gravitational potential in a galaxy, at least for properly selected systems. Despite these improvements and the substantially higher resolution of present-day dark matter only simulations, the discrepancy still persists.

An alternative solution to the cusp-core problem is that baryonic processes associated with galaxy formation and evolution also affect the dark matter halo. Various processes have been proposed in this context with different results. On the one hand, condensation of cooling gas towards the centre of a galaxy causes a further contraction of the dark matter halo and a stronger cusp (e.g. Jesseit et al. 2002; Gnedin et al. 2004). On the other hand, infalling gas clumps can transfer angular momentum to the dark matter via dynamical friction, ultimately resulting in a shallower central profile (El-Zant et al. 2001; Tonini et al. 2006; Romano-Díaz et al. 2008), but the efficiency of this mechanism is still under debate (e.g. de Blok 2010). Finally, feedback from supernovae (and AGN activity in high mass galaxies) can induce massive gas outflows that also cause the dark matter halo to expand. Navarro et al. (1996a) already investigated this scenario using a highly simplified outflow model and concluded that supernova feedback could have flattened the dark matter cusps in dwarf irregular galaxies, but is unlikely to be effective in more massive systems. In a more detailed study Read & Gilmore (2005) found that repeated bursts of star formation alternated by epochs of gas (re-)accretion can indeed gradually transform dark matter cusps into cores in simulated dwarf galaxies. This result was later confirmed and extended to somewhat higher mass galaxies by numerous studies (e.g. Governato et al. 2010; Macciò et al. 2012; Teyssier et al. 2013; Di Cintio et al. 2014b; Oñorbe et al. 2015), although the details of the conclusions sometimes differ. For example, while Governato et al. (2012) found that supernova feedback can only expand dark matter haloes in galaxies with M_∗ ≳ 10⁷M_⊙, Read et al. (2016a) concluded that cores also form in lower mass systems if star formation proceeds for long enough.

Hydrodynamical simulations therefore seem to suggest that stellar feedback effects can solve the cusp-core controversy. To really confirm this claim, however, the simulated haloes must be compared to real observed rotation curves. For this purpose, two new analytic dark matter density profiles were recently proposed. The DC14 profile was derived by Di Cintio et al. (2014a) from their simulated galaxies in the mass range 9.94 × 10⁹<M_halo/M_⊙< 7.8 × 10¹¹ and was recently tested on samples of observed rotation curves by Katz et al. (2017) and Pace (2016). Both studies concluded that the DC14 profile can indeed reproduce the observed rotation curves. Katz et al. (2017) also found that the derived dark matter halo parameters are in excellent agreement with the cosmological stellar mass-halo mass and halo mass-concentration relations. Pace (2016), on the other hand, concluded the opposite. He found halo masses significantly below the cosmological prediction for galaxies with stellar masses M_∗ ≲ 10⁹M_⊙, and a huge scatter of almost two orders of magnitude in the derived halo concentrations. The coreNFW model was derived by Read et al. (2016a) from simulations of tiny dwarf galaxies in the mass range M_halo ~ 10⁷−10⁹M_⊙. Read et al. (2017) also recently tested this profile on a set of dwarf galaxies extracted from the Little THINGS sample (Hunter et al. 2012; Iorio et al. 2017) with halo masses ranging from a few 10⁸M_⊙ to about 2 × 10¹⁰M_⊙. These authors also found good fits to the observed rotation curves and good agreement with the stellar mass-halo mass relation.

Although the physical mechanism that drives core formation is essentially the same for both halo models, they follow a somewhat different approach. In the DC14 model the stellar mass is used as a measure for the amount of supernova feedback energy that has become available and the shape of the dark matter halo is fully determined by M_∗/M_halo. In the coreNFW model, on the other hand, the core strength is regulated by the total time that the galaxy has been forming stars, while the radial extent of the core is linked to the radial distribution of the stars.

The coreNFW and DC14 models also probe different mass ranges and are in this sense complementary. While it is probably not meaningful to apply the DC14 model to the rotation curves of tiny dwarf galaxies or to extrapolate coreNFW to Milky Way-size systems, the two models should both be appropriate for halo masses of the order 5 × 10⁹M_⊙ ≲ M_halo ≲ 5 × 10¹⁰M_⊙. Therefore, if DC14 and coreNFW both correctly describe dark matter core formation, their predictions should agree in this overlapping mass range.

In this work we apply the coreNFW and DC14 models to a selection of 13 galaxy rotation curves with halo masses of 4 × 10⁹M_⊙ to 7 × 10¹⁰M_⊙, and compare their predictions. The DC14 halo is further applied to an additional 7 rotation curves with halo masses up to 9 × 10¹¹M_⊙ to investigate the agreement with the cosmological scaling relations. This paper is organized as follows: in Sect. 2 we describe the selection of our sample of rotation curves. The principle of mass modelling and the details of the two dark matter halo models are explained in Sect. 3. Our modelling strategy is described in Sect. 4 and the results are presented in Sect. 5. Finally we list our main conclusions in Sect. 6.

2. Sample selection

Our sample of rotation curves was compiled mainly from the Little THINGS (Hunter et al. 2012; Iorio et al. 2017), THINGS (Walter et al. 2008; de Blok et al. 2008), and SPARC (Lelli et al. 2016) datasets. For Little THINGS we used the publicly available rotation curve data from Iorio et al. (2017) and took the surface density and surface brightness profiles of the atomic gas and the stars from Oh et al. (2015). The latter were kindly provided to us by S.H. Oh and D. Hunter. The THINGS data, both the rotation curves and the baryonic profiles, were kindly made available by E. de Blok. Finally the SPARC data are publicly available and can be downloaded from the SPARC website¹. For all three datasets the stellar surface brightness profiles are based on observations at 3.6 μm.

For Little THINGS we selected only the galaxies that are marked as “clean dIrrs” by Read et al. (2017) and further eliminated NGC 6822 and DDO 210. For the former no rotation curve is presented by Iorio et al. (2017), while the rotation curve of the latter is highly uncertain and completely dominated by the asymmetric drift correction. This leaves 9 galaxies from the Little THINGS sample. de Blok et al. (2008) present rotation curves of 19 THINGS galaxies, from which we eliminated 10 because of poor sampling of the rising part of the rotation curve or strong non-circular motions. The gap in mass between the Little THINGS and THINGS galaxies is bridged with a set of low mass systems from the SPARC dataset. These are selected according to the following criteria: a total 3.6 μm luminosity L_3.6 ≲ 10⁹L_⊙, a rotation curve with quality label 1, a reliable distance estimate, an inclination between 40° and 80°, and little beam smearing. These criteria lead to an additional 5 galaxies. Finally we also included the rotation curve of M33 (taken from Corbelli et al. 2014), which is a galaxy previously claimed to have a strongly cusped dark matter halo (Corbelli et al. 2014; Hague & Wilkinson 2015). For this galaxy Corbelli et al. (2014) have not reported the stellar surface brightness profile, but have immediately derived the surface density profile of the stars from a pixel by pixel population synthesis analysis.

Our total sample thus comprises 24 rotation curves. For three of the THINGS galaxies in this sample, however, the rotation curves might actually be unreliable owing to a substantial bar (NGC 925), poorly constrained distance (NGC 3521), or uncertain inclination in the outer half (NGC 7793). Since Hague & Wilkinson (2014) have included these rotation curves in their analysis, we also kept them in our sample, but marked them as problematic and only showed their fits without including them in the further analysis. UGC 8490, from the SPARC dataset, was similarly marked problematic because McQuinn et al. (2015) have discovered a strong increase in its star formation rate over the past 100 Myr. This indicates that UGC 8490 might be experiencing a starburst, which may severely bias its kinematics.

Finally, Iorio et al. (2017) have remarked that their rotation curves of UGC 8508 and DDO 126 are unreliable up to a radius of 0.5 and 1.43 kpc, respectively, while Gentile et al. (2007) noted that elliptical streaming motions could be affecting the innermost data points of their rotation curve of NGC 3741, up to a radius of 1.2 kpc. Blais-Ouellette et al. (2001) additionally found that the inner seven data points of the rotation curve of NGC 3109 could be slightly underestimated because of weak beam smearing. These same data points also have suspiciously small error bars. Since only the inner parts of the rotation curves are affected, we still included these galaxies in our “good” sample, but excluded the affected data points from the fits. The rotation curves are still sampled well enough by the remaining points.

An overview of our complete sample is given in Table 1. For the galaxies from the THINGS and SPARC datasets we used distances from the Cosmicflows-2 catalogue (Tully et al. 2013). For M33 we took over the distance from Corbelli et al. (2014) and for the Little THINGS galaxies we used the distances from Iorio et al. (2017).

3. Rotation curve decomposition and halo models

In a disk galaxy the inward gravitational force that pulls a particle towards the centre is balanced by the outward centripetal acceleration from its rotation. The total gravitational potential acting on this particle is the sum of the potentials from the individual components: gas, stars, and dark matter. We can therefore write (1)Since the centripetal acceleration is proportional to the square of the circular velocity, this can be re-written as (2)where v_gas, v_∗, and v_dm are the circular velocities needed to balance the gravitational force exerted by the gas, stars, and dark matter, respectively. These are of course related to the mass distributions of the individual components. For the stars, however, the conversion between the observed luminosity and the mass is uncertain. The mass-to-light ratio Υ is therefore isolated as an unknown parameter in Eq. (2). In addition, the distribution of the gas often shows a hole in the centre. Inside this hole the gravitational pull from the gas is directed outward, giving a negative contribution to the total circular velocity. A better formulation of Eq. (2) is therefore (3)or (4)if the stellar distribution is decomposed in bulge and disk components. The value v_∗ is now the circular velocity from the stars for a mass-to-light ratio of 1.

Equations (3) and (4) form the basis for the mass modelling performed in this work. The total circular velocity v_c is measured by the rotation curve (although, see Sect. 3.3), while the gas and stellar circular velocities v_gas and v_∗ are derived from their observed surface brightness profiles. For this purpose a thin disk geometry is generally assumed for the gas. The stellar distribution is usually modelled as a thick disk with an exponential or sech² profile in the vertical direction. This leaves the stellar mass-to-light ratio(s) and dark matter contribution as the only unknowns. For the latter we use two different parameterizations: the coreNFW halo and DC14 halo. We express both these parameterizations in terms of the virial radius and virial mass. The former is defined as the radius inside which the average density of the dark matter halo is equal to Δ times the critical density of the Universe ρ_crit, where Δ and ρ_crit depend on the assumed cosmology. The virial mass is simply the enclosed mass at the virial radius, (5)For consistency with Di Cintio et al. (2014a) we use a WMAP3 cosmology (Spergel et al. 2007) with Δ = 93.6, H₀ = 73.0 km s^-1 Mpc^-1 and ρ_crit = 147.896 M_⊙ kpc^-3.

3.1. DC14

The DC14 profile is formulated by Di Cintio et al. (2014a) as a special case of the general and very flexible (α, β, γ) profile (Jaffe 1983; Hernquist 1990; Zhao 1996) (6)At small and large radii this profile follows a power law with slopes γ and β, respectively, and the sharpness of the transition between these two regimes is governed by α. This profile reduces to a simple Navarro-Frenk-White (NFW) profile for (α, β, γ) = (1, 3, 1) and the frequently used pseudo-isothermal halo is recovered when (α, β, γ) = (2, 2, 0). Hague & Wilkinson (2014, 2015) have recently used this profile in its most general form to model the dark matter haloes of M33 and a sample of THINGS galaxies.

Starting from an NFW profile, but in the general formulation of Eq. (6), Di Cintio et al. (2014a) allow the modification of the inner slope by stellar feedback by expressing the shape parameters α, β, and γ as a function of the integrated star formation efficiency M_∗/M_halo as follows: (7)with X = log₁₀ (M_∗/M_halo). These expressions are only valid for −4.1 < X < −1.3, which is the range probed by the simulations of Di Cintio et al. (2014a). At lower values of X, too few stars form to modify the dark matter halo. On the other hand, at X > −1.3 (corresponding to halo masses ≳10¹²M_⊙) processes not included in the simulations, such as AGN feedback, can start to play a role as well.

The variation of α, β, and γ as a function of X is shown in Fig. 1. The inner log slope γ first decreases with increasing X, since a higher stellar-to-halo mass ratio implies more energy input from supernova feedback. However, it reaches a minimum at X ~ −2.6 and goes back up at higher values of X. The reason for this turnover is the increasing gravitational potential of the stars, which at a certain point starts to dominate the feedback and pulls the dark matter back towards the centre. In the DC14 formalism more star formation therefore does not monotonically result in ever stronger cores.

For an NFW halo the concentration is defined as c = r_vir/r_s, where the scale radius r_s is equal to r_-2, the radius at which the slope of the density profile becomes –2. For the (α, β, γ) profile, the meaning of r_s depends on the values of α, β, and γ, with (8)Di Cintio et al. (2014a) therefore define the concentration of their dark matter haloes as (9)The concentration of the original, unmodified NFW halo can be recovered from this as (10)with again X = log₁₀ (M_∗/M_halo). It is this concentration that should be used to compare the DC14 halo from a fit to a rotation curve to, for example the mass-concentration relation.

3.2. coreNFW

A coreNFW (Read et al. 2016a) halo is essentially a NFW halo with the inner part modified by a spherically symmetric function fⁿ that models the effects of supernova feedback. Practically this modification is expressed at the level of the enclosed mass. For an ordinary NFW halo profile (Navarro et al. 1996b) (11)with concentration (12)the enclosed mass at a radius r is given by (13)The coreNFW profile is then defined as (14)with (15)The radial extent of the core is determined by the core radius r_c, which Read et al. (2016a) relate to the stellar half-mass radius as r_c = ηr_1/2, with an optimal value of 1.75 for the fitting parameter η .

The strength of the core is governed by the parameter n, which ranges between 0 < n≤ 1 and is defined as (16)Here κ is again a fitting parameter and the star formation time t_SF is the total time that the galaxy has been forming stars. The dynamical time t_dyn is the duration of 1 circular orbit at the scale radius in the unmodified NFW halo (17)The longer stars have been forming, the larger n and stronger the core. On the other hand, the bigger the original dark matter halo, the smaller n and more difficult it is to form a core. Following Read et al. (2016a), we set κ = 0.04 and choose t_SF = 14 Gyr.

3.3. Asymmetric drift correction

The gravitational attraction from the gas, stars, and dark matter is in fact not balanced solely by circular motion, but also for a small part by the internal pressure of the gas. The observed rotation velocity v_rot is therefore not exactly equal to the circular velocity v_c from Eq. (3). Instead it is given by (18)(Eqs. (4)–(227) of Binney & Tremaine 2008), where ρ and σ are the density and velocity dispersion of the gas. The asymmetric drift correction (term inside the square brackets) is usually simplified under the assumptions that the velocity dispersion is isotropic (σ_R = σ_φ), the velocity ellipsoid is aligned with the cylindrical coordinate system (), and the vertical scale height does not change much with radius. This leads to (19)where Σ is the surface density of the gas. The observationally derived radial Σσ² profile is typically rather rugged, leading to sometimes strong and unphysical fluctuations in its derivative. To avoid this, a smooth function is fitted to the profile and the derivative is determined analytically. The simplifications involved in deriving Eq. (19) limit its accuracy. As a consequence, Eq. (19) only provides an order of magnitude estimate of the correction.

For the rotation curves taken from the Little THINGS and SPARC datasets the asymmetric drift correction is already taken into account by the authors. On the other hand, de Blok et al. (2008) and Corbelli et al. (2014) did not consider asymmetric drift for their THINGS and M33 rotation curves. We therefore evaluated this correction based on Eq. (19). Depending on the shape of the Σσ² profile, we used one of the following analytic functions: (20)for a profile with a central core (Oh et al. 2011), and (21)for a profile showing a hole in the centre (Read et al. 2016b). In the inner halves of the rotation curves the derived corrections are consistently much smaller than the error bars and generally only of the order of 1 km s^-1 or less. For four galaxies the corrections become more substantial (of order 5−10 km s^-1) near the outer edge of the rotation curve. However, because these larger corrections occur far from the centre, we found that they have only little effect on our fits and do not change any of our conclusions. In addition, the agreement between the Σσ² profile and the analytic function is often not very good in these regions, making the corrections uncertain. This is illustrated for NGC 3621 in Fig. 2. Because of this uncertainty and the very limited effect the corrections have on our results, we decided to use the original, uncorrected rotation curves in our analysis.

4. Markov Chain Monte Carlo fitting

4.1. emcee

The dynamical models were fitted to the rotation curves with emcee² (Foreman-Mackey et al. 2013), which is an open-source python implementation of the affine invariant MCMC ensemble sampler from Goodman & Weare (2010). Markov Chain Monte Carlo or MCMC (e.g. Metropolis et al. 1953; Hastings 1970; Press et al. 2007) is a sampling technique that has been applied to the decomposition of rotation curves for several years (e.g. Puglielli et al. 2010; Hague & Wilkinson 2013). It is more efficient in sampling the parameter space than the fitting techniques used in earlier works and has the big advantage that it returns the full multidimensional probability distribution of all the parameters instead of only the best-fit model. In addition MCMC allows us to include physical knowledge about the parameters in the fits via so-called priors that are combined with the likelihood function.

The emcee algorithm explores the N-dimensional parameter space with different, randomly initialized walkers that each make their own MCMC chain. An initial burn-in phase is used to allow the walkers to move to the relevant high-likelihood areas of the parameter space. After this the walkers are reinitialized at their current positions and the actual MCMC chains are made. As a last step the chains of all the walkers are combined to form the final MCMC chain. For every fit we used 100 walkers, each taking 2000 steps of which the first 1000 were used as burn-in. These numbers are in line with the emcee recommendations (Foreman-Mackey et al. 2013) and the values typically used in other works (e.g. Kirichenko et al. 2015; Katz et al. 2017; Read et al. 2017), and ensured good convergence of our fits (see below). Our likelihood function is (22)For a multi-modal posterior distribution, part of the walkers can get stuck in isolated low probability modes if they are initialized randomly over the full range of the parameter space (within the imposed boundaries). This generates numerous irrelevant peaks in the retrieved posterior distribution. We therefore performed each fit in two iterations. First the walkers were initialized randomly over the full relevant range of parameter space. The different peaks in the posterior distribution were then investigated to find the mode with the highest likelihood. Next, as a second iteration, we redid the fit with the walkers now initialized in a small Gaussian ball centred on this mode and with σ equal to 1 percent of the allowed range for each parameter. The parameter values that are used in the figures below and reported in Table 2 correspond to the maximum likelihood model for each fit.

For good performance, an MCMC sampler should be run for at least a few (about 10) autocorrelation times and should have an acceptance fraction between 0.2 and 0.5 (Foreman-Mackey et al. 2013). With 1000 steps taken by each walker, the first condition was well met for all the fits. Appropriately setting the emcee proposal scale parameter to a value of 2 or 3 ensured that the second condition was also met. Finally we checked the convergence of the MCMC chains by performing each fit three times and evaluating the Gelman-Rubin eigenvalues with the GetDist³ python package. These values were well below 1 for all the fits, indicating good convergence.

4.2. Priors and parameter ranges

The coreNFW halo fits were performed with log₁₀ M_vir, c, and Υ (or Υ_d and Υ_b) as free parameters. We use the log of M_vir instead of M_vir itself as a parameter in the fits because of the large dynamical range involved. Flat priors were assumed for all free parameters. Log₁₀ M_vir and c were loosely constrained inside 8 < log ₁₀ (M_vir/M_⊙) < 14 and 1 < c < 100. The 3.6 μm mass-to-light ratio was confined to the range 0.3 < Υ_3.6 < 0.8, as motivated by the constraints from Meidt et al. (2014) and McGaugh & Schombert (2014). For M33 we allow the initial stellar mass to vary by a factor 0.758 < Υ < 1.319 based on the uncertainty that is mentioned in Sect. 6 of Corbelli et al. (2014). Following Read et al. (2017), η, κ, and t_SF were kept fixed at 1.75, 0.04, and 14 Gyr, respectively.

For the fits with the DC14 halo, we let V_vir, c_vir, and Υ (or Υ_d and Υ_b) free and again used a flat prior for each of these. Following Katz et al. (2017) we used wide ranges of 10 <V_vir/ (km s^-1) < 500 and 1 <c_vir< 100 for the first two parameters and the same range as before for the mass-to-light ratio: 0.3 < Υ_3.6< 0.8 (and 0.758 < Υ < 1.319 for M33).

Since our goal is to find models that fit the rotation curves well and yield physically acceptable dark matter haloes at the same time, we further imposed the cosmological halo mass-concentration and stellar mass-halo mass relations as log-normal priors in the fits. For the DC14 halo the shape parameters α, β, and γ are expressed as a function of log₁₀ (M_∗/M_halo), where M_halo = M_vir. Hence, α, β, and γ depend on the definition of the virial mass, which depends on the assumed cosmology. Since Di Cintio et al. (2014a) have assumed a WMAP3 cosmology, we did the same in our fits and we used the M_halo−c relation from Macciò et al. (2008) that was derived under this cosmology. The M_∗−M_halo relation from Moster et al. (2010) has also used the WMAP3 values. However, this relation was derived from abundance matching using the SDSS DR3 stellar mass function for halo masses down to ~3 × 10¹⁰M_⊙ and is an extrapolation at lower masses. As pointed out by Read et al. (2017), this extrapolation is not consistent with the newer and deeper SDSS data, which means that the Moster et al. (2010) relation is actually not reliable for halo masses below ~3 × 10¹⁰M_⊙. Indeed the stellar mass-halo mass relation from Behroozi et al. (2013), which is based on the newer SDSS data and the cosmological parameters used in the Bolshoi simulations (compatible with WMAP5 and WMAP7; Klypin et al. 2011), is much shallower and diverges significantly from the Moster et al. (2010) relation at low halo masses. Since a number of the galaxies in our sample fall in this low mass regime we opted to use the relation from Behroozi et al. (2013) rather than that from Moster et al. (2010) as a prior in our fits. For coreNFW the parameterization is independent of the assumed cosmology and depends only on the global original NFW profile (i.e. before alteration by stellar feedback) and on the stellar half-mass radius and total star formation time.

4.3. Uncertainties

Using the GetDist package, the uncertainties for the different parameters were determined from the multidimensional 68% confidence region of the full posterior distribution, as the extremal values of the projection of that region onto each parameter axis. As such the error bars give a good indication of how tight the constraints are for a certain parameter, but they should not be over-interpreted as the absolute range of good models. Indeed, if the fit quality of the best-fit model is very high, many models outside of the N-dimensional 68% confidence region often still provide an acceptable fit to the data. On the other hand, if we plot, for example, the halo mass versus its concentration, the area suggested by the two (orthogonal) error bars is often larger than the actual area to which the models from the MCMC chain are confined (i.e. the projection of the multidimensional confidence region onto the M_halo−c plane). This is illustrated in Fig. 3 for the best-fit DC14 model of UGC 7603.

4.4. χ²_red and fit quality

In the discussion of our results we express the quality of the fits to the rotation curves in terms of the reduced chi-squared statistic (). In the ideal case where the uncertainties on all the rotation curves are Gaussian and derived in a uniform way, and where all the points of a rotation curve have equal importance, this would be a good measure to compare the fit qualities for all the galaxies in our sample. In reality, however, our rotation curves are compiled from the literature with differing data quality and techniques used to estimate the error bars. In addition some rotation curves keep rising up to the last point, whereas others, for the more massive galaxies, include a large flat part. The latter is generally easier to reproduce and can have a large impact on the value of a fit, but is at the same time much less important in the analysis of core formation. For these reasons the values of our fits are only meaningful to compare the quality of different fits for the same galaxy and not to compare fits for different galaxies.

5. Results

Following the procedure outlined in Sect. 4 we have fitted DC14 halo models to each of the rotation curves in our sample and coreNFW models to the rotation curves of the Little THINGS and SPARC galaxies. The best-fit parameters and values of these fits are listed in Table 2. The results of the fits are discussed in the sections below.

5.1. DC14

Figure 4 shows the individual DC14 models for the galaxies in our “good” sample. As can be seen, the DC14 halo generally provides excellent fits to the rotation curves, confirming the recent results from both Katz et al. (2017) and Pace (2016). The only clear exception to this is the rotation curve of DDO 168, where the model overestimates the data in the inner part. However, the inner three points of the rotation curve are in fact already well accounted for by the gravitational potential of the gas alone, so any model with a non-zero contribution of the dark matter at these radii will overestimate the data. Similar arguments also hold for the very inner regions of NGC 2366 and NGC 3198.

In Fig. 5 we compare the best-fit parameter values from our models with the cosmological halo mass-concentration and stellar mass-halo mass relations. These relations were derived from dark matter-only simulations (in combination with abundance matching). To account for this in the comparison, we scale our inferred halo masses as M_vir/(1−f_b), where f_b is the Universal baryon fraction (0.176 according to WMAP3; McCarthy et al. 2007). The fits show excellent agreement with both scaling relations, although, somewhat surprisingly, our models seem to favour the M_halo−c relation from Dutton & Macciò (2014) that is based on the Planck cosmology over the Macciò et al. (2008) relation that was used as prior in the fits. The only galaxy that falls significantly outside the 2σ scatter of the Dutton & Macciò (2014) relation is NGC 3198, but models with the concentration forced inside this scatter actually still provide a good fit to the data.

Although the sample size is limited, our analysis thus seems to confirm the recent conclusion by Katz et al. (2017) that the DC14 halo is in good agreement with ΛCDM and can recover the predicted mass-concentration and stellar mass-halo mass relations in a sample of observed rotation curves. The latter might seem obvious since we use priors that “push” our fits towards these relations. The key point here, however, is that DC14 can produce dark matter haloes that simultaneously provide good fits to the rotation curves and agree with the scaling relations; this is something that, for example the NFW halo, cannot do. The priors do not “make” physical solutions according to the scaling relations, but merely act as a filter to retain only the most physical solutions if they exist. The fact that our models prefer the Dutton & Macciò (2014) relation over the Macciò et al. (2008) relation that was used as prior is a good illustration of this.

Our analysis contradicts the results from Pace (2016). We did not recover the huge range of halo concentrations that he found and we found no evidence for his claim that galaxies with M_∗ ≲ 10⁹M_⊙ often reside in less massive haloes than predicted. It should be noted here that the modelling strategy from Pace (2016) is somewhat different to our approach and that used by Katz et al. (2017). Pace (2016) used multi-nested sampling (e.g. Feroz & Hobson 2008) to fit his models and did not assume any physical priors between the parameters. The lack of priors results in posterior distributions that often contain multiple modes (or peaks). The mode with the lowest halo mass was selected as the final mode and galaxies for which the modes were too wide or not well separated were discarded. To investigate the effects of these differences, we performed a second set of fits without physical priors and used the criteria from Pace (2016) to select the final models. The results are shown in Fig. 6. CVnIdwA, UGC 8508, NGC 3741, NGC 3109, IC 2571, and NGC 2976 were discarded because their posterior distributions showed multiple blended modes; this was not the case in the original fits because the priors suppressed the additional modes. The best-fit model for WLM has shifted significantly. It is still consistent with the mass-concentration relation but now lies considerably below the stellar mass-halo mass relation. For the remaining galaxies the new models are essentially similar to the previous models or consistent within the uncertainties. The agreement with the scaling relations is still remarkably good, and although the number of galaxies is limited we find no evidence for the strong deviations that were reported by Pace (2016).

5.2. CoreNFW

For the coreNFW halo we limited the sample to the lower mass galaxies with M_halo ≲ 7 × 10¹⁰M_⊙. The individual coreNFW fits for these galaxies are shown in Fig. 7. The models again provide a decent description of the data with no clearly bad fits except for DDO 168. The fit results are compared to the cosmological scaling relations in Fig. 8. The agreement is very good and our models again seem to prefer the mass-concentration relation from Dutton & Macciò (2014) over that from Macciò et al. (2008).

The coreNFW halo was also fitted to the same rotation curves by Read et al. (2017). Their results for the individual galaxies sometimes differ significantly from what is found in this work: Read et al. (2017) have generally found somewhat lower halo and stellar masses, higher concentrations, and a better fit quality (lower ). These differences are, however, not unexpected. Indeed, while we derived the stellar and gas contributions from the measured surface density profiles from Oh et al. (2015), Read et al. (2017) used smooth exponential profiles based on Zhang et al. (2012) and Oh et al. (2015). In addition, Read et al. (2017) used the M₂₀₀ formalism while we used the virial mass M_vir, so the halo parameters given in their Table 2 should not be compared directly to our values in Table 2. For an identical dark matter halo, our virial mass and concentration should be somewhat higher than the M₂₀₀ and c₂₀₀ values from Read et al. (2017). The fact that we generally find lower concentrations probably comes from the fact that Read et al. (2017) did not use a mass-concentration prior, but instead set the boundaries for the concentration range based on the M_vir−c relation from Macciò et al. (2007) and the extremities of their M₂₀₀ range. The rotation curve of NGC 2366 that was reported by Iorio et al. (2017) also seems somewhat different from the curve that was used by (or at least shown in Fig. A2 of) Read et al. (2017).

Despite these individual differences, the main conclusions remain the same. Both works generally find acceptable fits to the data and a good agreement with the stellar mass-halo mass relation.

5.3. coreNFW versus DC14

The coreNFW halo model determines the strength of the dark matter core from the total time that the galaxy has been forming stars, while the radial extent of the core is related to the spatial distribution of the stars via the stellar half-mass radius. The stellar mass of the galaxy is not used. This is an important difference with the DC14 profile. On the one hand, DC14 uses M_∗ instead of t_SF as a measure of the amount of supernova feedback energy that has become available to form a core. On the other hand, this model also uses the additional gravitational potential due to M_∗ as a mechanism to counteract core formation. In addition, the coreNFW profile is essentially a pure NFW profile with its inner part flattened by feedback, while for DC14 the entire shape of the profile (i.e. γ, α, and β) changes as a function of the stellar-to-halo mass ratio.

Given these rather different approaches, it is interesting to investigate how the dark matter haloes inferred by the coreNFW and DC14 halo models compare across our sample. Since the physical mechanism that drives core formation is essentially the same for both models (supernova feedback after bursts of star formation), they should in principle give similar results in the overlapping halo mass range where they are both appropriate.

5.3.1. Fit quality

From Table 2 we see that the DC14 and coreNFW models generally give a similar fit quality. The DC14 model typically has a slightly lower , with only WLM and UGCA 442 breaking this trend, but in most cases both models represent the data well enough and it is not meaningful to classify one as better than the other. The latter also holds for DDO 168, which is poorly fitted by both models. As discussed in Sect. 5.1 this is probably caused by a problem with the data.

The only cases in which DC14 performs better than coreNFW are NGC 3741, DDO 87, and DDO 154. For NGC 3741 and DDO 87 the rotation curves are fitted significantly better by the DC14 model than by the coreNFW model, although the latter is also still acceptable. For DDO 154 the fit quality is good in both cases, but the coreNFW halo needs an unphysically high concentration to achieve this. coreNFW models with a lower concentration do not fit the data well. In each of these cases the problem seems to be related to the connection between the coreNFW core radius and the stellar half-mass radius. For NGC 3741 the rotation curve suggests an extended dark matter core, while the coreNFW model has a cuspy NFW shape in all but the most central region because of the tiny half-mass radius of the stars. For the other two galaxies the rotation curves require a smaller core and a less “linear” dark matter contribution than derived from R_1/2.

These issues can be resolved by making η a free parameter in the coreNFW fits, but this would break with the prescription of Read et al. (2016a) and in a way make the comparison with the DC14 model unfair. For NGC 3741 this might be justifiable. Indeed, while the stellar distribution of NGC 3741 is unusually compact, the DC14 model finds sufficient stellar mass (i.e. supernova feedback energy) to form a substantial core that is in agreement with the data. Furthermore, the distribution of the gas is much more extended than that of the stars. Since the gravity of the outflowing gas after a supernova explosion drives the formation of a dark matter core, we could interpret the extended gas distribution in NGC 3741 as a sign that the dark matter core radius is actually larger than that inferred from the half-mass radius of the stars. Therefore we performed a second coreNFW fit for NGC 3741 with η as an additional free parameter (using a flat prior such that 0 <r_c(kpc) < 7). The decomposed rotation curve is given in Fig. 9. With a larger core radius of 4.06 kpc (η = 11.2) the coreNFW model is now in excellent agreement with the data and with the DC14 model. This value of η is significantly above the upper limit η = 2.75 derived by Read et al. (2017). However, this upper limit was derived under the assumption that R_1/2 ~ 0.015r₂₀₀, which is also significantly larger than the value R_1/2 = 0.341 kpc that is measured from the stellar distribution.

For DDO 154 and DDO 87 the stellar distribution does not seem particularly unusual, although the inferred stellar half-mass radii are significantly larger by a factor of 3.9 and 1.8, respectively than those reported by Read et al. (2017). Nevertheless, this is also the case to some degree for CVnIdwA and DDO 52 (factors of 1.6 and 1.4), for which DC14 and coreNFW give similarly good fits. It is beyond the scope of this work to investigate whether this points to a problem with the Oh et al. (2015) data, and we limit ourselves to the conclusion that, based on the surface brightness profiles from Oh et al. (2015), the DC14 model performs better than the coreNFW model for the rotation curves of DDO 154 and DDO 87.

5.3.2. Best-fit parameters

In Fig. 10 we compare the main parameters of the best-fit DC14 and coreNFW models for the galaxies in our sample with M_halo ≲ 7 × 10¹⁰M_⊙. With the exception of the unphysically high concentration of the coreNFW model of DDO 154, we find a fairly good agreement between the best-fit parameter values from both models. The stellar and virial masses follow the identity line, albeit with some scatter and the two models are generally consistent within the error bars. The DC14 model typically finds somewhat less extended cores (as can be seen by comparing the individual rotation curve decompositions) and more concentrated haloes, although except for DDO 154 the concentrations are consistent within the errors. For the galaxies where the concentration difference is the highest, the DC14 model typically also finds somewhat lower stellar and halo masses than the coreNFW model, which is consistent with the cosmological scaling relations. The bottom right panel in Fig. 10 shows the log slope of the dark matter density profiles measured at the innermost point of each rotation curve. This parameter strongly depends on the local shape of the density profile and should not be over-interpreted for the comparison between the two haloes. The main conclusion to be drawn is that both models find relatively cored dark matter haloes for all the galaxies with inner log slopes that are generally not too far apart and often consistent within the uncertainties. Unlike what is claimed by Read et al. (2017), we therefore find that the DC14 model can still show significant cusp-core transformations for galaxies with M_halo ≲ 10¹⁰M_⊙.

5.4. Problematic galaxies

For completeness we show the individual fits to the problematic rotation curves in Appendix A. Using the flexible (α, β, γ) profile, Hague & Wilkinson (2014) find good fits to the rotation curves of NGC 925 and NGC 7793 and also report tight constraints on the inner log slope of NGC 3521 despite the poor fit of their model to the rising part of the rotation curve. In contrast, we find that the two physically motivated halo models investigated here cannot reproduce any of these rotation curves. For NGC 3521 we even perform additional fits at the much smaller distance of 7.7 Mpc reported by SPARC, but with the same result. Finally the cuspy rotation curve of UGC 8490 is well reproduced by the DC14 model, but coreNFW strongly underestimates the rising part. However, for both these models the parameters lie well outside the cosmological scaling relations.

6. Conclusions

For a compact sample of 13 galaxies, spanning the mass range M_halo ~ 4 × 10⁹−7 × 10¹⁰M_⊙, we have used MCMC to construct dynamical models of the rotation curves based on two recently proposed dark matter density profiles: the DC14 halo and coreNFW halo. We further applied the DC14 halo to an additional set of higher mass galaxies with M_halo ~ 10¹¹−9 × 10¹¹M_⊙. The coreNFW and DC14 halo models both use supernova feedback after bursts of star formation to transform primordial dark matter cusps into flatter cores. Although they were derived from simulations in complementary mass ranges (M_halo ~ 10⁷−10⁹M_⊙ for coreNFW versus 10¹⁰−8 × 10¹¹M_⊙ for DC14), both models should be valid for halo masses of 5 × 10⁹M_⊙ ≲ M_halo ≲ 5 × 10¹⁰M_⊙. With this analysis we investigated whether the coreNFW and DC14 halo models converge to the same solutions in this overlapping mass range and whether their predictions agree with cosmological scaling relations.

We found that both models are generally able to reproduce the rotation curves in our sample and find dark matter haloes that are in good agreement with the cosmological M_halo−c and M_∗−M_halo relations, alleviating the cusp-core controversy. This confirms the results from Read et al. (2017) for coreNFW and from Katz et al. (2017) for DC14. On the other hand, we find no evidence of the huge scatter in concentrations or the disagreement of the DC14 predictions with the M_∗−M_halo relation that were recently claimed by Pace (2016), even if a similar modelling strategy is used.

The two models generally give similarly good fits to the rotation curves, although the DC14 model does perform better in three cases. For NGC 3741 and DDO 87 the rotation curves are fitted significantly better by the DC14 model than by the coreNFW model, although the latter is also still acceptable. For DDO 154 the fit quality is good in both cases, but the coreNFW halo needs an unphysically high concentration to achieve this. In each of these cases the problem for coreNFW is related to the connection between the core size and the stellar half-mass radius.

For NGC 3741, a galaxy with a very compact stellar distribution but a remarkably extended atomic gas disk, we argue that it is justifiable to relax this connection and use η as a free parameter in the coreNFW fit. With a larger core radius the coreNFW model is in excellent agreement with the data.

The DC14 and coreNFW haloes generally converge to (approximately) the same solution, as they should. Both models find cored dark matter haloes, and while DC14 tends to predict somewhat less extended cores and more concentrated haloes, the stellar masses, halo masses, and concentrations from both models are generally comparable and agree within the errors.

---

Acknowledgments

We thank the anonymous referee for the useful comments and suggestions that significantly improved this paper. We also wish to thank E. de Blok for sharing his rotation curves and surface density profiles of the THINGS galaxies and D. Hunter and S.H. Oh for sharing the stellar and gas surface brightness profiles of the Little THINGS galaxies. F.A., M.B., and G.G. acknowledge the support of the Flemish Fund for Scientific Research (FWO-Vlaanderen). This research has made use of NASA’s Astrophysics Data System and the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

References

Appendix A: Fits to the problematic rotation curves

All Tables

All Figures