© ESO 2020

1. Introduction

Dark matter (DM), which makes up 85% of all the matter in the Universe, is still one of the most elusive mysteries in present-day physics. The existence of nonbaryonic DM is considered a fact that has been proved by Planck measurements of the cosmic microwave background (CMB), for instance. This component is thought to be made of nonrelativistic particles that can be described by a collisionless fluid that interacts with ordinary particles (baryons) through gravity alone because the cross section for nucleons is thought to be only ∼10⁻²⁶ cm² (e.g., Jungman et al. 1996).

At present, Λ cold dark matter (ΛCDM) is known as a very successful cosmological scenario in describing the formation and evolution of the large-scale structures of the Universe (Frenk et al. 1985; Bode et al. 2001). The N-body simulations in this scenario have been able to resolve the current structure of virialized objects from clusters to dwarf galaxies. The simulations have shown that the ΛCDM power spectrum of the density perturbations combined with the collisionless nature of the CDM particles lead in any virialized object from clusters to dwarf galaxies to a very peculiar cuspy dark halo density profile: ρ_DM ∝ r⁻¹. In detail, dark halos around galaxies are well described by the Navarro-Frenk-White (NFW) profile (Navarro et al. 1996, 1997).

It is also known that this sharp central cusp of the DM density distribution strong contradicts the analysis of the kinematics of disk galaxies. For the first time, Gentile et al. (2004) showed in mass models obtained from six high-quality extended RCs that combined optical and 2D HI data overwhelming evidence for cored dark halo distributions beyond any statistical or bias uncertainty. More references on this topic is given in Salucci (2019).

About 100 RCs of the above quality have been thoroughly investigated so far, and the results have led to a similar decisive support in favor of the DM halo cored distribution (e.g., Adams et al. 2012; Donato et al. 2009; Weinberg et al. 2015; Simon et al. 2005; de Blok & Bosma 2002; Spekkens et al. 2005). The current status of the analysis of individual RCs is described in Korsaga et al. (2018, 2019a), who obtained mass models of 31 spiral and irregular galaxies, derived using hybrid rotation curves that combined high-resolution optical GHASP Fabry-Perot Hα RCs with extended radio WHISP HI RCs. Moreover, for each galaxy, the analysis took advantage of high-quality 3.4 μm WISE photometry, which is a fair indicator of the stellar disk mass. Korsaga and collaborators found that independent of the value of the galaxy optical velocity (V_opt), the performance of the baryonic matter plus NFW halo in reproducing the RCs was poor. Cored DM halo models are also required in galaxies that are dominated by luminous matter inside R_opt.

This topic has been considered serious enough for alternatives to the collisionless dark particle scenario to have been proposed, including warm dark matter (WDM; e.g., Di Paolo et al. 2018), self-interacting DM (Spergel & Steinhardt 2000), and ultra-light axions (ULAs; e.g., Tulin et al. 2013); moreover, the possibility that the dark particles gain energy from the standard model particles by means of feedbacks, that is, by sudden gas outflows into the galaxy halo driven by supernova explosions, raises doubts on the collisionless status of dark particles in the ΛCDM scenario (e.g., Governato et al. 2012).

Although the cored density distribution is often considered as preponderant in disk systems, a complete and comprehensive comparison between the galaxy kinematical data and the NFW halo profile has not been carried out so far. The notable questions even in disk systems are (i) the degree of the disagreement of data with simulations as a function of the galaxy reference velocity (V_opt) and of its Hubble type, and (ii) the way in which this discrepancy can be qualitatively characterized in general and in the various different objects.

The results found in this work are a necessary first step to answer these questions. In this regard, we stress that we tested the NFW profile not only by determining whether it reproduces the kinematics of disk systems, but also by verifying whether the resulting best-fit models (i) reproduce the relationship between two profile parameters, that is, concentration (c) and virial mass (M_vir), that emerges from N-body simulations; and if they (ii) lead to fractions of dark matter in systems of different virial masses that agree with independent measurements.

Only a complete assessment of the actual inability of the virialized noncollisional matter to describe the actual density halos around galaxies can shed light on how the collisionless scenario may be changed properly. The processes within a collisionless particle scenario include adiabatic contraction, feedback from supernovae, and clump migration (see, e.g., Blumenthal et al. 1986; El-Zant et al. 2016; Macciò et al. 2012; Di Cintio et al. 2014). Moreover, newly proposed scenarios for the DM particles do not always produce cores. For instance, the ULAs scenario seems to form core radii in the DM halos with M <  10¹¹ M_⊙ (e.g., de Martino et al. 2018); in the baryonic feedback scenario, the DM halos retain their cuspy profiles if M_stellar/M_halo <  0.001 or M_stellar/M_halo >  0.005; and in WDM scenario, it is difficult to explain cores in the halo densities of dwarf galaxies (e.g., Salucci 2019).

Furthermore, we also exploit the fact that the luminous and dark masses of galaxies can be determined by methods alternative to the RC analysis and then compared with those obtained from the RC analysis. In detail, disk masses can be derived by multiplying the galaxy luminosities with their mass-to-light ratios estimated from their colors, while the virial masses are obtained by means of weak-lensing measurements (e.g., Reyes et al. 2012; Mandelbaum et al. 2016) or from the abundance-matching method (e.g., Moster et al. 2013; Shankar et al. 2006; Rodriguez-Puebla et al. 2012).

We assume Ω_total = 1; Ω_matter = 0.3; Ω_baryonic = 0.04; Ω_DM = 0.26; Λ = 0.7 and H₀ = 72 km s⁻¹ Mpc⁻¹.

The structure of this paper is as follows: in Sects. 2 and 2.1, we describe the samples we used. In Sect. 2.2 we review the process of coaddition of individual RCs, and in Sect. 2.3 we discuss the mass modeling, which includes the disk and halo components. In Sect. 3 we obtain the mass model of 26 coadded RCs and discuss their implications, in particular, the implications from the resulting best-fitting values of the model-free parameters, M_D, M_vir, and c. Finally, we draw conclusions in Sect. 4 based on our results.

2. Samples and mass model

2.1. Sample

The 26 coadded RCs that we obtained from more than 3200 individual RCs of disk systems and that we discuss in Sect. 2.2 belong to four samples: Persic et al. (1996), Catinella et al. (2006), Di Paolo et al. (2019), and Karukes & Salucci (2017). These coadded RCs have never been used for the goal of testing NFW profile. 900 RCs in Persic et al. (1996) were also independently analyzed by Catinella et al. (2006). (i) In the first sample, optical RCs out to ∼1.2 R_opt are combined with HI RCs out to 2 R_opt, where R_opt is the optical radius and defined as the radius enclosing 83% of the total stellar mass of a galaxy. (ii) The second sample is composed of optical RCs alone. (iii) The third and the fourth samples consist of a combination of about 90% optical RCs and 10% HI RCs. Then, the innermost region (R <  2/3 R_opt) of almost all the 3200 galaxies are traced by high-resolution optical RCs. In short, each coadded RC is obtained from a large number of measurements, all of which belong to RCs of galaxies of same Hubble type and with similar values of V_opt, which is the velocity at the optical radius (R_opt). Moreover, they have a very high spatial resolution, so that they are not affected by the beam-smearing effects that may distort the inner profile of low-resolution RCs. Finally, all the coadded RCs we used are extended out to at least, R ≃ 3 R_D, where R_D is the disk scale length. This is discussed in Sect. 2.2.

The first sample has been published in Persic et al. (1996, hereafter PSS96). It consists of ∼900 individual RCs, arranged into 11 M_I (the I-band absolute magnitude) bins, whose central values span from −17.5 to −23 (which corresponds to 75 km s⁻¹ to 279 km s⁻¹), and which yield 11 coadded rotation curves built by arranging the RCs in their corresponding M_I bins. In radial normalized coordinates x ≡ R/R_opt, the radial bins have a size of ∼0.1.

The second sample comes from Catinella et al. (2006, hereafter Ca06; see also Lapi et al. 2018), and it consists of ∼2200 individual RCs that also include the PSS96 sample. Individual RCs are arranged into nine optical velocity bins, whose central values span from 104 km s⁻¹ to 330 km s⁻¹, and the normalized radial bin has a size of 0.06.

The third sample, from Di Paolo et al. (2019, hereafter LSB), consists of 72 low surface brightness galaxies with velocities at the optical radius V_opt spanning from ∼24 km s⁻¹ to ∼300 km s⁻¹. They are arranged into five optical velocity bins whose central values are given in Di Paolo et al. (2019) and whose normalized radial size is ∼0.1.

The fourth sample, from Karukes & Salucci (2017, hereafter KS17), consists of 36 RCs of dwarf irregular galaxies arranged into one velocity bin whose optical velocities V_opt span from 17 km s⁻¹ to 61 km s⁻¹. Its central value is ⟨V_opt⟩ = 40 km s⁻¹, and the normalized radial average size of each bin is ∼0.1.

Samples 1 and 2 therefore include similar objects (normal spirals), and the second sample acts as an independent verification of the first. Samples 3 and 4 include peculiar disk systems (low surface brightness and dwarf disks). The LSB sample is an environment extremely useful for investigating dark matter. The disk components of LSBs have a much lower surface density than normal spirals. The dwarf disks are objects with the smallest disks.

The PSS96 and Ca06 samples include about 90% of the objects belonging to Sb-Im Hubble types. Both samples have galaxies with V_opt in the range of normal spirals and are much larger than those in the KS17 sample. The LSB includes only low surface brightness objects, which in turn are almost entirely missing in the previous three samples.

2.2. Coadded rotation curves

It is well known that each disk galaxy has two tags (properties) that specify it completely: (i) R_D the size of the stellar disk that is derived by its surface luminosity profile. In these objects it is reasonable to assume that (Freeman 1970)

(1)

I₀ is the central surface luminosity, ∼100 L_⊙ pc², and R_D = 1/1.67 R_1/2 is the scale length that determines the distribution of stars in the galaxy disk, while R_1/2 is the half-light radius. The photometric quantity R_opt, analogous to the well-known half-light radius R_1/2 in ellipticals (e.g., Salucci 2019), is related to the scale length of the Freeman disk by

(2)

Thus, for exponential thin disks: R_D = 3.2/1.67 R_1/2; moreover, R_opt and R_D are physically identical. We stress that in the case of stellar disks with (usually mildly) nonexponential profiles, R_opt also still marks the dynamical edge of the disk as the quantity 3.2 R_D in the case of Freeman disks.

(ii) Disk galaxies include irregular spirals, dwarf irregular disks and low surface brightness disks in addition to normal spirals. These types of objects are all rotationally supported, but have very different photometric properties (see Salucci 2019) and therefore are different probes for the distribution of dark matter. We specify that for a given R_D, the disk mass can vary largely, depending on the morphological type.

We can represent all the rotation curves of disk systems by means of a universal profile (e.g., Salucci et al. 2007). In detail, the universal RC (URC) is obtained from thousands of individual RCs when they are arranged in a number of coadded rotation curves. The uncertainties of coadded (stacked) rotation curves are significantly smaller than those of the individual RCs. Moreover, these coadded curves are as extended as the most extended individual RCs: r_outermost >  5R_D, where r_outermost is the center of the outermost radial bin in each coadded RC and has a very high spatial resolution (< 1/4 R_D). The process of building such curves is explained in detail in the original works, but we briefly recall it here.

The basic assumption underlying the coaddition process and verified in disk systems is that the galaxies with similar V_opt ≡ V(R_opt) (that also have similar R_opt) have similar rotation curves, especially when the latter are expressed in the normalized radial coordinate x ≡ R/R_opt (Yegorova & Salucci 2007; Gammaldi et al. 2018; López Fune 2018; Rhee 1997). This leads to similar mass distributions, with the DM density always taken as a cored distribution.

The concept of coadded RCs, implicit in Rubin et al. (1985), which was pioneered by Persic & Salucci (1991), set by Persic et al. (1996), and extended to large galactocentric radii and different Hubble types by Salucci et al. (2007), Di Paolo et al. (2019), and Karukes & Salucci (2017) allows us to derive V_coadd(R/R_opt, M_I) or V_coadd(R/R_opt, V_opt), that is, the average (coadded) RCs of galaxies with the magnitude lying in a bin centered at the I-band absolute magnitude M_I and of size δM_I (or centered at V_opt and of size δV_opt). The number of objects and the position and size of each magnitude (or velocity) bin are indicated in detail in the original works.

The coadded RCs were built in three steps in a process that is similar to the process used in cosmological numerical simulations to work out the averaged DM density distribution of halos of similar virial mass M_vir (see: Navarro et al. 1996). We describe here the details for one of the four samples we used (PSS96); the procedure is only marginally different from those adopted for the other three samples.

(i) The whole I-magnitude range is divided into 11 successive bins, each of which is centered at M_j, I (j = 1…,11) and has a size δM_j, I. The 11 values of the latter quantities, alongside with the 11 values of V_j, opt, R_j, opt, are given in Table 1 (also in Persic et al. 1996). The rotation curve V_jk(R) of each galaxy of the sample assigned to its corresponding j-magnitude bin is normalized by its V_jk, opt value. Next, by using the known galaxy value of R_jk, opt, each RC is expressed in terms of its normalized radial coordinate x_jk ≡ R/R_jk, opt. Here, j is the index of the magnitude (velocity) bins, jk is the index of the various galaxies in the above bins before coaddition, ji is the index of the radii at which jcoadded RC have measurements. We directly determined in the first sample (PSS96) the various R_jk, opt from the definition, for LSB and KS17, the Freeman disk is an excellent distribution for their stellar content, and this also holds for the spirals in Ca06, if slightly less so. In the above works, the R_jk, D values were used to derive R_jk, opt from Eq. (2). In this way, every RC of our sample was then double-normalized and tagged as

(3)

Let us notice that V is velocity in physical units while v is double-normalized velocity. (ii) Each of the N_j double-normalized individual rotation curves v_jk in each j bin was binned into 20 radial bins of length 0.1 and index ji, and the number of data N_ji was then averaged to obtain the double-normalized coadded RCs and their uncertainties computed from their rms δv_j, coadd(x_ji, M_j, I) and N_ji: δv_j, coadd(x_ji, M_j, I)/(N_ji − 2)^1/2.

(iii) The final step is that the double-normalized curves are denormalized using the averaged values of V_jk, opt and R_jk, opt.

In general, the extension of the 26 coadded RCs is x_{j, outermost} = 1.75 ± 0.25 and the positions of the various x_ji are given in the table available at the CDS for j = 1…26. For x >  1.2, the concurring RCs are mostly HI RCs. We recall that Gentile et al. (2004) have shown that HI and optical RCs describe the same gravitational potential. Noticeably, these coadded curves vary very little: almost always, δV_coadd ≤ 0.03V_coadd (e.g., Salucci 2019) because for most of the radial bins in each coadded RC, we have many measurements, from 20 to 200. When these are averaged, the statistical component of the variance in the velocity measurements is much reduced. Furthermore, v_jk(x_jk, M_j, I), the individual double-normalized RCs of objects with the same magnitude M_I are all very similar.

It is also important to stress in reference to all samples (with j = 1...26), that in each of the j luminosity (velocity) bins, the values of the normalization velocities V_jk, opt of the individual RCs assigned to the j bin are very similar. They each differ from the averaged bin value V_j, opt by less than 10%. For simplicity of notation for the quantities V_j, opt and R_j, opt we do not indicate the sign < >, but these quantities must be averaged over the index jk with fixed j. In the first luminosity (velocity) bin of each sample, the values of the self-normalization velocities V_jk, opt may vary significantly among the galaxies of the bin; however, in these cases, corresponding to the lowest luminosity (velocity) bins of the four samples, the profiles of the individual RCs, that is, d log V/d log R, which are crucial for the mass modeling, are all very similar and independent of V_jk, opt, for example. We stress that the averaged value ⟨d log V/d log R⟩ between 1 R_D and 3 R_D is lower than 0.02 for all the 26 coadded RCs, which is irrelevant in the mass modeling process.

Similarly, in each of the 26 magnitude (velocity) j bins the radial normalization quantities R_jk, opt (fixed j) are very similar, and the rms of the averaged values R_j, opt is < 25%. Moreover, because the individual RCs locally have a mostly linear profile, V(R + δR)≃V(R)+const. × δR, the variance of the various R_jk, opt(fixed j) with respect to the average R_j, opt affects the coaddition process only very mildly. Therefore, the finiteness of the 26 luminosity (velocity) bins does not affect the structure of the double-normalized coadded RCs and their subsequent denormalizations.

The URC of disk systems is instead the analytical function devised to fit the coadded RCs V_j, coadd(r, M_j, I) in Persic et al. (1996) and in the others in the subsequent works. It was chosen as the sum in quadrature of the components to the circular velocity coming from a Freeman stellar disk and a Burkert dark matter halo (Burkert 1995; Persic et al. 1996). In this work, the URC provides a crucial model comparison. In previous works, we have shown that all the 26 coadded RCs that we used here, are perfectly fit by this function, that is, in more detail, by a velocity model including a cored DM halo with and a Freeman disk. This model has the same number of free parameters as the model we analyze here: the core radius r₀, the central DM density ρ₀, and the disk mass M_D. In all 26 cases the URC fitting uncertainties were found to be lower than the rms of the V_j, coadd(x_ji, M_j, I or V_j, opt) measurements (Persic et al. 1996; Catinella et al. 2006; Di Paolo et al. 2019; Karukes & Salucci 2017; Salucci et al. 2007).

We briefly recall that at a physical radius R, the normalized radius x and the double-normalized circular velocity v(x) are defined as

(4)

For a double-normalized coadded RC j with velocity data at x_ji, with j = 1, …26, we have

(5)

which gives R_ji, the positions in physical units of the coadded velocity data of v_j(x).

Figure 1 shows the 26 coadded RCs in double-normalized coordinates. To express them in physical units, it is necessary to rescale them according to the above relations and the values in Tables 1–3. We also plot the coadded RCs in a 3D form to show that in the (v, x, V_opt) coordinates system, they clearly define a 3D surface that leads to the concept of the URC (see Fig. 2).

Fig. 1. Twenty-six double-normalized coadded RCs. The data and their uncertainties are listed in the CDS table.

Fig. 2. R/Ropt − log Vopt − V/Vopt for the 26 coadded RCs. The legend is the same as in Fig. 1. For clarity we do not plot the error related to each measurement; this can be found in the CDS table.

2.3. Mass model

In the ΛCDM scenario, we have from the N-body simulations that the virialized dark halos show a universal spherically averaged density profile,

(6)

where R is the radial spherical coordinate, and ρ_s and r_s are the characteristic density and the scale radius of the dark halo. It is useful to define the concentration c ≡ R_vir/r_s, where R_vir is the virial radius given by . The critical density is with H₀ the current value of Hubble’s parameter. We assume H₀ = 72 km s⁻¹ Mpc⁻¹. No result of this paper changes by assuming the lower value of H₀ = 67 km s⁻¹ Mpc⁻¹ favored by Plank CMB measurements (Planck Collaboration VI 2020).

In the circular velocity model, the contribution of the DM halo to the circular velocity adds in quadrature to that of the stellar disk as

(7)

where V_D and V_NFW are the contribution of the stellar Freeman disk and the dark halo, respectively. In the disk systems we consider in this work, the stellar component is described by an exponential Freeman disk (see Freeman 1970) with a surface density given by

Its contribution to the circular velocity is

(8)

where it is worth recalling that M_D is the disk mass, x ≡ R/R_opt and I_n, K_n are the modified Bessel functions of nth order.

The dark halo contribution to the circular velocity that we investigate here is the contribution generated by the NFW dark halo density profile (ρ_NFW) and is given by

(9)

where M_vir is the halo virial mass, y = R/R_vir is the normalized radial coordinate, and c = R_vir/r_s is the concentration parameter.

The velocity model that we tested with 26 coadded rotation curves obtained from more than 3200 individual RCs that consist of more than 10⁵ independent kinematical measurements is

(10)

the notation (c, M_D, M_vir) has to be considered as a 26 × 3 matrix.

Finally, we stress that we did not test an empirical dark matter density profile (as in Burkert 1995) for which we would need data out to the galaxy virial radius to prove its success, but a theoretical profile obtained from N-body simulations in the ΛCDM scenario, whose functional form is believed to be perfectly known from r = 0 to r = R_vir. In detail, the two free parameters ρ_s and r_s can be derived by the inner regions of the RCs, that is, for R <  r_s (see Eq. (6)). We do not need very extended kinematics to disprove the presence of an NFW halo. It is clear that with more external data the estimate of these parameters may change, but our aim is not to show a global failure of the NFW + disk model but only its extremely serious troubles in the region in which the model-independent analysis of RCs has revealed DM density cores, that is, for r <  r₀ ∼ r_s. We did not investigate whether for r >  R_opt the DM halo density converges to an NFW, as it appears to do in spirals (Salucci et al. 2007).

3. Results

The data fitting was performed by means of the nonlinear least-square method (χ²-Levenberg-Marquardt method). We fit each coadded rotation curve with the above model. The χ² statistic is defined as

(11)

where V_ji, coadd is the value of coadded RCs at x_ji tagged as j; same definition stands for σ_ji, coadd and V_ji, mod that are the (observational) uncertainties of the j coadded RCs value at x_jiR_j, opt and the corresponding model value, respectively. N_ji are the number of kinematical measurements of each coadded RC (see Tables 1–3), and X is the set of three free parameters of the model, X = (M_D, M_vir, c), which are to be determined through the fitting procedure. The values of R_ji = x_jiR_j, opt, V_ji, coadd, and σ_ji are given in the table available at the CDS.

These best-fit values and those of the reduced (), where and the denominator is the number of degrees of freedom of the model, are shown in Tables 1–3 for the 25 coadded RCs investigated here. In Fig. A.1 we show the triangular plots for 25 mass models.

The coadded RC of KS17 has been modeled in the original paper, where the authors found

(12)

with .

At this point, we state that in order to reclaim success, we expect from the model under investigation that it fits the data satisfactorily () and that the values of free parameters are at least in a fair agreement with the following relations holding for disk galaxies: with 0.2 <  α <  0.8 from the colors and the luminosities of spirals (e.g., Salucci et al. 2008), c ∼ (7 − 17) from N-body simulations (Bullock et al. 2001) according to the halo masses (see Eq. (13)), and M_vir ∼ 10¹² M_⊙(V_opt/(200 km s⁻¹))³ from weak-lensing measurements and abundance-matching technique (e.g., Salucci 2019).

Overall, the model performs very poorly. The most relevant problems are that there is a wide range of values for the concentration parameter c, which ranges from 1.5 to 28.3, and as indicated in boldface in Tables 1–3, in 50% of the cases it is well outside the simulation range.

There are also 12 entirely implausibly high values for the halo virial masses in the 26 best-fitting models. We stress that these high values do not originate from the lack of kinematical measurements in the outer parts of galaxy halos (i.e., for R/R_opt >  1), but from the values of the (coadded) RCs in the innermost regions of the galaxies (R <  R_opt): the RC profiles are so steep (i.e., V_coadd ∝ R) there that in order to reproduce this RC characteristic with an NFW halo profile, the best-fit model is forced to have very high values of the parameters r_s and ρ_s, which in turn yields implausibly high values for the virial halo mass, for which .

Moreover, according to Tables 1–3, in order to minimize the value, the best fit (NFW halo + Freeman disk) model forces the parameter M_D in 5 cases out of the total 26 to assume values that are far lower than those obtained from the L_B galaxy luminosities and average spiral mass-to-light ratios M_D/L_B = 2. The uncertainties on M_D propagated from those in R_opt are smaller than 25% and therefore negligible for the results of this work. From its maximum value, obtained with the assumption that at R = 2.2 R_D the disk component entirely dominates the circular velocity V(2.2 R_D), we have M_D, max ∼ G⁻¹V² (2.2 R_D) 2.2 R_D. The actual disk mass can be a fraction of it.

Furthermore, the values are not satisfactory. In detail, we find that this quantity exceeds the value of 1.4 in 16 cases. We recall that the URC velocity model (Burkert halo + Freeman disk) fits all the 26 coadded RCs we used here extremely well (see Persic et al. 1996; Lapi et al. 2018; Di Paolo et al. 2019; Karukes & Salucci 2017). However, poor fits are only a part of the failure of the model we tested.

The performance of the model under investigation was also tested by correlating the best-fit values of the two DM structural parameters c and M_vir and then comparing them with the corresponding scaling relation that emerges in the N-body simulations within the ΛCDM scenario (Bullock et al. 2001):

(13)

The cosmic variance in log c is 0.1 dex, which is negligible for the aim of our work. Figure 3 shows the c − log(M_vir) relation of our 26 best-fit models with their corresponding uncertainties. We can easily realize that most of the best-fit values lie very far from the Eq. (13) relation, of which they are unable to trace even the gross trend.

Fig. 3. Relation between c and log Mvir from the 26 mass models. Blue points indicate the results of this paper (red points indicate a very large uncertainty on their estimates). The ΛCDM outcome from simulations is also shown (green line).

Another scaling relationship of great importance is the disk mass M_D versus the virial halo mass M_vir; noticeably, this relation can be independently derived by the abundance-matching method (e.g., Moster et al. 2010). We find

(14)

where M₁, A, β, and γ are

(15)

A completely independent analysis has yielded the disk mass M_D versus virial halo mass M_vir (see Shankar et al. 2006). This is just mildly different from the mass in Eq. (14), also considering that both have an uncertainty of 0.2 dex in log (M_D),

(16)

which we also included in the comparison with the best-fit values of the model under investigation. The stellar-halo mass relationship so obtained is shown in Fig. 4, along with the above independent relationships (14) and (16). It is evident that most of the values of the (NFW + stellar disk) model lie very far from either relation. The results of KS17 are entirely inconsistent with the observational relations.

Fig. 4. MD − Mvir relationship in logarithmic scale obtained from the 25 coadded RCs (blue points). The red circles are the same as in Fig. 3. The orange and gray lines correspond to relations (14) and (16). The value of log MD for KS17 has the extremely discrepant coordinates (see the text).

The poorest performance of the model under study is probably related to the fact that this never shows a halo mass with M_vir <  10¹² M_⊙. Halo masses like this are known to exist in a great number in the ΛCDM scenario.

The failure of the model under analysis to reproduce the kinematics of spirals is summarized in Tables 1–3 and Figs. 3 and 4. In order to have a more complete view of the problem, we plot the values of three important structural parameters derived from the best fits of the 26 coadded RCs: in Fig. 5 we show c as a function of log (M_D/M_vir) and . The two surfaces restrict the volume in which we expect the values of the DM parameters to lie in, which means that the ΛCDM NFW model would be a valid representation of the 26 coadded RCs. Each sphere tags the best-fit model of a coadded RC. We realize a large failure of the NFW density profile. Nearly all of the coadded RCs lie far outside of the goal volume.

Fig. 5. Relation between c, ratio and chired. ratio ≡ log(MD/Mvir) and chired ≡ . The values are obtained from the 26 mass models (3D blue spheres) and the ΛCDM expectation (the volume between two red surfaces).

We now discuss the other two minor, but cosmologically important components of baryonic matter in disk systems. A bulge contribution to the circular velocity is present in (only) disk galaxies with the highest stellar mass. In Appendix B we include this component in the circular velocity model of the coadded RCs relative to these objects. We find that the resulting best-fit model agrees well with that obtained from Eq. (10). The other component is a gaseous HI disk in the outermost parts of the least luminous objects. In Appendix C we present the model for the coadded RCs of the latter. We add an HI gaseous disk contribution and compare the resulting best-fit model with that of Eq. (10). We find a good agreement in this case as well. This shows that including these two minor baryonic components does not change the results of this work.

It is important to compare the results of our work with a previous study of individual RCs (Korsaga et al. 2019b). It is worth noticing that the coadded RCs with respect to the individual RCs show on average higher values of c and a stronger correlation with the halo mass as likely expected by the coaddition process, which has eliminated a part of the random errors of the individual RCs. In Fig. 6 we plot our 26 c versus V_opt along with those obtained from individual RCs in Korsaga et al. (2019b). The agreement further reinforces the statement of the gross inability of NFW halos to reproduce the observed disk kinematics.

Fig. 6. Relation between concentration c and optical velocity Vopt. The quantity Vopt is derived using the relation MB = −7log Vopt − 4.9 from Ponomareva et al. (2017). The blue points are the results of this work, and the orange points are the results of Korsaga et al. (2019b) using the fixed M/L fitting technique.

4. Conclusion

More than 3200 RCs of disk systems with different V_opt, ranging from 20 km s⁻¹ to 330 km s⁻¹, were stacked into 26 coadded RCs that extended out to R_opt and beyond. These data represent the whole kinematics of disk systems well and were used to investigate a particular mass distribution. This is characterized by a dark halo with an NFW profile and a stellar Freeman exponential thin disk. The model has three free parameters: the disk mass M_D, the halo virial mass M_vir, and the halo concentration c. These were derived by means of best-fitting the coadded RCs. Because the NFW velocity model is an analytical model that is valid from r = 0 and r = R_vir, we can derive all the related ρ_s, r_s, and R_vir from the two halo parameters above.

The aim of our investigation was to determine whether dark halos around spirals possess cuspy inner profiles. The outcomes crucially add to the results that were obtained in the past 20 years by means of the analysis of several dozen individual RCs. Although these have given an unequivocal answer against this possibility, they cannot be considered either entirely unbiased or complete.

By means of the above 26 coadded RCs we have completed the investigation about the core-cusp question in disk systems of any morphology and luminosity (mass). We found that generally, the model under investigation fails to reproduce the observational data. Not only does it fit the coadded RCs pooly, implying high values for , but in many cases, the best-fit values of the structural parameters c, M_vir and M_D are also very different from the expectations from simulations and from measurements independent of RCs data. Noticeably, the model fails at any reference velocity from (V_opt = 20 km s⁻¹) to (V_opt >  300 km s⁻¹). One exception may be for objects with 190 km s⁻¹ ≤ V_opt ≤ 210 km s⁻¹, for which the model appears to perform as well as the Burkert cored halo + Freeman disk model; in this case, only individual HI RCs with the optical velocity in this range and with very much extended kinematics can solve this riddle (Karukes & Salucci 2017). We stress that according to our results, we conclude that within the ΛCDM scenario, the corification of initial cusps must occur in the smallest and largest disk systems.

A detailed investigation shows that the failure of the model has different aspects: the concentration c takes a wide range of values from 1.5 to 28.3 that is well outside the range found in simulations. The dark matter halo masses in the 50% of the cases reach implausibly high values. In some cases, the values of the disk masses are ridiculously low. The value, on the other hand, can be acceptable for only 38% of the cases.

From the best-fit parameters we reconstructed the model scaling relations c − M_vir and M_D − M_vir that emerge to clearly contradict those obtained from observations and from N-body simulations (Salucci et al. 2007; Shankar et al. 2006; Bullock et al. 2001). The 3D relation of in Fig. 5 clearly shows that most of the best-fit values lie outside the volume that indicates a successful representation of the kinematics of disk galaxies.

This result must be gauged by the fact that, in previous works, we found that the cored halo+ Freeman disk is able to reproduce all these 26 coadded RCs with very low values (≃1) and yields meaningful relations among the model structural parameters (see Karukes & Salucci 2017; Lapi et al. 2018; Di Paolo et al. 2019).

In this work, however, we make no claim about (i) the ΛCDM scenario, in that cores in the DM density might be produced during the cosmological evolution of galaxies, (ii) the DM distribution at large radii r >  100 kpc, in galaxies of other Hubble types, in spirals at high redshifts z >  0.5, and in disk systems just after their formation for the evident lack of available measurements. Only a large number of high-quality individual RCs, that is, 100 per magnitude, per Hubble type, per 0.5 redshift can resolve all these questions.

Acknowledgments

We thank the referee for suggestions that have improved the quality of the paper and Chiara Di Paolo for useful discussions.

References

Appendix A: Best-fit rotation curves and triangle plots

In this section we plot the 25 velocity best-fitting models to their corresponding coadded RCs (left panels of Fig. A.1). The triangular plots for the 25 mass models are also shown in this figure (right panels of Fig. A.1).

Fig. A.1. RC modeling based on the NFW halo profile and the stellar Freeman exponential disk (left panels) and their corresponding triangle plots of 2D invariance between parameters with 1σ uncertainty (right panels). Left panel: the red circles are observed data with error bars, the blue line is the best-fitting model, the dashed green line shows the NFW halo component, and the dot-dashed red line is the disk component.

Fig. A.1. continued.

Fig. A.1. continued.

Fig. A.1. continued.

Fig. A.1. continued.

Fig. A.1. continued.

Fig. A.1. continued.

Fig. A.1. continued.

Fig. A.1. continued.

Appendix B: Bulge contribution to the rotation curves of disk systems

The largest disk galaxies with V_opt >  220 km s⁻¹ have an additional inner stellar component, a bulge. The rotation curve of these systems is then composed of a spherical bulge, a disk, and a dark halo,

where V_D and V_NFW are the disk and dark halo components of the rotation curve, respectively, and are defined in Eqs. (8) and (9). V_B is the bulge contribution. This component is smaller than the stellar disk by at least a factor of 3. We assumed that we can represent this component by the following velocity profile:

where M_B is the bulge mass. In the above, the bulge behaves like a point mass at large distances, for R ≪ 1/3 R_D, it increases like R².

The velocity model therefore has four free parameters. Table B.1 shows the best-fitting results for the most massive bins of the three samples Ca06, PSS96, and LSB where this component is (mildly) relevant. The corresponding RCs are shown in Fig. B.1. The bold and red values are the same as in Tables 1–3. We realize that including the very central bulge does not change the results of this work. It is worth noting that we did not investigate the bulge-dominated Sa spirals here.

Fig. B.1. Observed and fit RC modeling: The components of the NFW dark halo (green dashed), disk (dot-dashed red), and bulge (dashed black) that contribute to the total rotation curve and the observed data (red circles) with error bars are shown in the left panels. Their corresponding triangle plots of the 2D invariance between parameters with 1σ uncertainty are shown in the right panels.

Appendix C: Gas contribution to the rotation curves

Disk galaxies also receive an HI contribution to their circular velocities. We did not consider this here because (i) this component is relevant only in small objects and outside R_opt, (ii) in any case, the effects of incorporating this gaseous component results in NFW cusps that are even more unlikely (e.g., Di Paolo et al. 2019; Karukes & Salucci 2017).

In this section we add these component to the stellar disk + NFW halo model to the circular velocity of galaxies for the lowest velocity bins of our sample. The HI surface density can approximately be reproduced by a Freeman distribution with a scale length R_HI that is about three times larger than that of the stellar disk (see Evoli et al. 2011; Wang et al. 2014). We then take for the contribution of the HI disk to the circular velocity

where and M_HI is the mass of HI disk, considered as a free parameter as M_D and I_n and K_n are the modified Bessel functions calculated at 0.53x. The expected range of the free parameter values of the stellar+ gaseous disk and the dark halo are reported in Table C.1, and the mass model best-fitting result is shown in Fig. C.1. Compared to the results we obtained when we neglected the gas contribution (see Fig. A.1), it produces only a mild difference in the mass model when we include the gas.

Fig. C.1. Best-fit RC of the first bin of the LSB sample by considering the gaseous component. The dashed green, dot-dashed red, dotted black, and solid blue lines stand for the DM halo, the stellar disk, the gaseous disk, and the total contributions to the circular velocity.

All Tables

All Figures

	Fig. 1. Twenty-six double-normalized coadded RCs. The data and their uncertainties are listed in the CDS table.
In the text

	Fig. 2. R/Ropt − log Vopt − V/Vopt for the 26 coadded RCs. The legend is the same as in Fig. 1. For clarity we do not plot the error related to each measurement; this can be found in the CDS table.
In the text

	Fig. 3. Relation between c and log Mvir from the 26 mass models. Blue points indicate the results of this paper (red points indicate a very large uncertainty on their estimates). The ΛCDM outcome from simulations is also shown (green line).
In the text

	Fig. 4. MD − Mvir relationship in logarithmic scale obtained from the 25 coadded RCs (blue points). The red circles are the same as in Fig. 3. The orange and gray lines correspond to relations (14) and (16). The value of log MD for KS17 has the extremely discrepant coordinates (see the text).
In the text

	Fig. 5. Relation between c, ratio and chired. ratio ≡ log(MD/Mvir) and chired ≡ . The values are obtained from the 26 mass models (3D blue spheres) and the ΛCDM expectation (the volume between two red surfaces).
In the text

	Fig. 6. Relation between concentration c and optical velocity Vopt. The quantity Vopt is derived using the relation MB = −7log Vopt − 4.9 from Ponomareva et al. (2017). The blue points are the results of this work, and the orange points are the results of Korsaga et al. (2019b) using the fixed M/L fitting technique.
In the text

Fig. A.1. RC modeling based on the NFW halo profile and the stellar Freeman exponential disk (left panels) and their corresponding triangle plots of 2D invariance between parameters with 1σ uncertainty (right panels). Left panel: the red circles are observed data with error bars, the blue line is the best-fitting model, the dashed green line shows the NFW halo component, and the dot-dashed red line is the disk component.

In the text

	Fig. A.1. continued.
In the text

	Fig. A.1. continued.
In the text

	Fig. A.1. continued.
In the text

	Fig. A.1. continued.
In the text

	Fig. A.1. continued.
In the text

	Fig. A.1. continued.
In the text

	Fig. A.1. continued.
In the text

	Fig. A.1. continued.
In the text

	Fig. B.1. Observed and fit RC modeling: The components of the NFW dark halo (green dashed), disk (dot-dashed red), and bulge (dashed black) that contribute to the total rotation curve and the observed data (red circles) with error bars are shown in the left panels. Their corresponding triangle plots of the 2D invariance between parameters with 1σ uncertainty are shown in the right panels.
In the text

	Fig. C.1. Best-fit RC of the first bin of the LSB sample by considering the gaseous component. The dashed green, dot-dashed red, dotted black, and solid blue lines stand for the DM halo, the stellar disk, the gaseous disk, and the total contributions to the circular velocity.
In the text