Issue 
A&A
Volume 643, November 2020



Article Number  A161  
Number of page(s)  21  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201937079  
Published online  19 November 2020 
NavarroFrenkWhite dark matter profile and the dark halos around disk systems^{⋆}
^{1}
SISSA, Via Bonomea 265, 34136 Trieste, Italy
email: salucci@sissa.it
^{2}
Faculty of Physics, Semnan University, Semnan 3513119111, Iran
email: RaziehDehghani@semnan.ac.ir, hghafarnejad@semnan.ac.ir
Received:
7
November
2019
Accepted:
25
July
2020
Context. The Λ cold dark matter (ΛCDM) scenario is able to describe the Universe at large scales, but clearly shows some serious difficulties at small scales. The corecusp question is one of these difficulties: the inner dark matter (DM) density profiles of spiral galaxies generally appear to be cored, without the r^{−1} profile that is predicted by Nbody simulations in the above scenario.
Aims. It is well known that in a more physical context, the baryons in the galaxy might backreact and erase the original cusp through supernova explosions. Before the efficiency and the presence of this effect is investigated, it is important to determine how wide and frequent the discrepancy between observed and Nbodypredicted profiles is and what its features are.
Methods. We used more than 3200 quite extended rotation curves (RCs) of good quality and high resolution of disk systems that included normal and dwarf spirals as well as low surface brightness galaxies. The curves cover all magnitude ranges. All these RCs were condensed into 26 coadded RCs, each of them built with individual RCs of galaxies of similar luminosity and morphology. We performed mass models of these 26 RCs using the NavarroFrenkWhite (NFW) profile for the contribution of the DM halo to the circular velocity and the exponential Freeman disk for the contribution of the stellar disk.
Results. The fits are generally poor in all the 26 cases: in several cases, we find χ_{red}^{2} > 2. Moreover, the bestfitting values of three parameters of the model (c, M_{D}, and M_{vir}) combined with those of their 1σ uncertainty clearly contradict wellknown expectations of the ΛCDM scenario. We also tested the scaling relations that exist in spirals with the outcome of the current mass modeling: the modeling does not account for these scaling relations.
Conclusions. The results of testing the NFW profile in disk systems indicate that this DM halo density law cannot account for the kinematics of the whole family of disk galaxies. It is therefore mandatory for the success of the ΛCDM scenario in any disk galaxy of any luminosity or maximum rotational velocity to transform initial cusps into the observed cores.
Key words: galaxies: kinematics and dynamics / dark matter
The 26 stacked rotation curves are only available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/cat/J/A+A/643/A161
© 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 presentday 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^{−26} cm^{2} (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 largescale structures of the Universe (Frenk et al. 1985; Bode et al. 2001). The Nbody 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^{−1}. In detail, dark halos around galaxies are well described by the NavarroFrenkWhite (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 highquality 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 highresolution optical GHASP FabryPerot Hα RCs with extended radio WHISP HI RCs. Moreover, for each galaxy, the analysis took advantage of highquality 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), selfinteracting DM (Spergel & Steinhardt 2000), and ultralight 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 bestfit models (i) reproduce the relationship between two profile parameters, that is, concentration (c) and virial mass (M_{vir}), that emerges from Nbody 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; ElZant 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^{11} 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 masstolight ratios estimated from their colors, while the virial masses are obtained by means of weaklensing measurements (e.g., Reyes et al. 2012; Mandelbaum et al. 2016) or from the abundancematching method (e.g., Moster et al. 2013; Shankar et al. 2006; RodriguezPuebla et al. 2012).
We assume Ω_{total} = 1; Ω_{matter} = 0.3; Ω_{baryonic} = 0.04; Ω_{DM} = 0.26; Λ = 0.7 and H_{0} = 72 km s^{−1} Mpc^{−1}.
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 bestfitting values of the modelfree 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 highresolution 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 beamsmearing effects that may distort the inner profile of lowresolution 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 Iband absolute magnitude) bins, whose central values span from −17.5 to −23 (which corresponds to 75 km s^{−1} to 279 km s^{−1}), 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^{−1} to 330 km s^{−1}, 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^{−1} to ∼300 km s^{−1}. 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^{−1} to 61 km s^{−1}. Its central value is ⟨V_{opt}⟩ = 40 km s^{−1}, 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 SbIm 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)
I_{0} is the central surface luminosity, ∼100 L_{⊙} pc^{2}, 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 halflight radius. The photometric quantity R_{opt}, analogous to the wellknown halflight radius R_{1/2} in ellipticals (e.g., Salucci 2019), is related to the scale length of the Freeman disk by
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 Iband 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 Imagnitude 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 jmagnitude 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 doublenormalized and tagged as
Values of the bestfit parameters of the mass modeling with their 1σ uncertainties for the PSS96 sample.
Let us notice that V is velocity in physical units while v is doublenormalized velocity. (ii) Each of the N_{j} doublenormalized 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 doublenormalized 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 doublenormalized 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 doublenormalized 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 selfnormalization 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 doublenormalized 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_{0}, the central DM density ρ_{0}, 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 doublenormalized circular velocity v(x) are defined as
For a doublenormalized coadded RC j with velocity data at x_{ji}, with j = 1, …26, we have
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 doublenormalized 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. Twentysix doublenormalized coadded RCs. The data and their uncertainties are listed in the CDS table. 
Fig. 2. R/R_{opt} − log V_{opt} − V/V_{opt} 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 Nbody simulations that the virialized dark halos show a universal spherically averaged density profile,
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_{0} the current value of Hubble’s parameter. We assume H_{0} = 72 km s^{−1} Mpc^{−1}. No result of this paper changes by assuming the lower value of H_{0} = 67 km s^{−1} Mpc^{−1} 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
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
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
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^{5} independent kinematical measurements is
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 Nbody 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 modelindependent analysis of RCs has revealed DM density cores, that is, for r < r_{0} ∼ 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 leastsquare method (χ^{2}LevenbergMarquardt method). We fit each coadded rotation curve with the above model. The χ^{2} statistic is defined as
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_{ji}R_{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_{ji}R_{j, opt}, V_{ji, coadd}, and σ_{ji} are given in the table available at the CDS.
These bestfit 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
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 Nbody simulations (Bullock et al. 2001) according to the halo masses (see Eq. (13)), and M_{vir} ∼ 10^{12} M_{⊙}(V_{opt}/(200 km s^{−1}))^{3} from weaklensing measurements and abundancematching 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 bestfitting 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 bestfit 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 masstolight 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^{−1}V^{2} (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 bestfit values of the two DM structural parameters c and M_{vir} and then comparing them with the corresponding scaling relation that emerges in the Nbody simulations within the ΛCDM scenario (Bullock et al. 2001):
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 bestfit models with their corresponding uncertainties. We can easily realize that most of the bestfit 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 M_{vir} 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 abundancematching method (e.g., Moster et al. 2010). We find
where M_{1}, A, β, and γ are
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}),
which we also included in the comparison with the bestfit values of the model under investigation. The stellarhalo 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. M_{D} − M_{vir} 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 M_{D} 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^{12} 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 bestfit 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(M_{D}/M_{vir}) 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 bestfit 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 bestfit 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 V_{opt}. The quantity V_{opt} is derived using the relation M_{B} = −7log V_{opt} − 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^{−1} to 330 km s^{−1}, 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 bestfitting 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 corecusp 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 bestfit 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^{−1}) to (V_{opt} > 300 km s^{−1}). One exception may be for objects with 190 km s^{−1} ≤ V_{opt} ≤ 210 km s^{−1}, 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 bestfit 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 Nbody 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 bestfit 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 highquality 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
 Adams, J. J., Gebhardt, K., Blanc, G. A., et al. 2012, ApJ, 745, 92 [Google Scholar]
 Blumenthal, G., Faber, S., Flores, R., & Primack, J. 1986, ApJ, 301, 27 [Google Scholar]
 Bode, P., Ostriker, J., & Turok, N. 2001, ApJ, 556, 93 [Google Scholar]
 Bullock, J. S., Kolatt, T. S., Sigad, Y., et al. 2001, MNRAS, 321, 559 [Google Scholar]
 Burkert, A. 1995, ApJ, 447, 25 [Google Scholar]
 Catinella, B., Giovanelli, R., & Haynes, M. P. 2006, ApJ, 640, 751 [Google Scholar]
 de Blok, W. J. G., & Bosma, A. 2002, A&A, 385, 816 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 de Martino, I., Broadhurst, T., Tye, S. H. H., et al. 2018, Galaxies, 6, 10 [Google Scholar]
 Di Cintio, A., Brook, C. B., Macciò, A. V., et al. 2014, MNRAS, 437, 415 [Google Scholar]
 Di Paolo, C., Nesti, F., & Villante, F. L. 2018, MNRAS, 475, 5385 [Google Scholar]
 Di Paolo, C., Salucci, P., & Erkurt, A. 2019, MNRAS, 490, 4551 [Google Scholar]
 Donato, F., Gentile, G., Salucci, P., et al. 2009, MNRAS, 397, 1169 [Google Scholar]
 ElZant, A. A., Freundlich, J., & Combes, F. 2016, MNRAS, 461, 1745 [Google Scholar]
 Evoli, C., Salucci, P., Lapi, A., & Danese, L. 2011, ApJ, 743, 45 [Google Scholar]
 Freeman, K. C. 1970, ApJ, 160, 811 [Google Scholar]
 Frenk, C., White, S. D. M., Efstathiou, G., & Davis, M. 1985, Nature, 317, 595 [Google Scholar]
 Gammaldi, V., Karukes, E., & Salucci, P. 2018, Phys. Rev. D, 98, 3008 [Google Scholar]
 Gentile, G., Salucci, P., Klein, U., Vergani, D., & Kalberla, P. 2004, MNRAS, 351, 903 [Google Scholar]
 Governato, F., Zolotov, A., Pontzen, A., et al. 2012, MNRAS, 422, 1231 [Google Scholar]
 Jungman, G., KamionKowski, M., & Griest, K. 1996, Phys. Rev., 267, 195 [Google Scholar]
 Karukes, E. V., & Salucci, P. 2017, MNRAS, 465, 4703 [Google Scholar]
 Korsaga, M., Carignan, C., Amram, P., Epinat, B., & Jarrett, T. 2018, MNRAS, 478, 50 [Google Scholar]
 Korsaga, M., Amram, P., Carignan, P., & Epinat, B. 2019a, MNRAS, 482, 154 [Google Scholar]
 Korsaga, M., Epinat, B., Amram, P., et al. 2019b, MNRAS, 490, 2977 [Google Scholar]
 Lapi, A., Salucci, P., & Danese, L. 2018, ApJ, 859, 2 [Google Scholar]
 López Fune, E. 2018, MNRAS, 475, 2132 [Google Scholar]
 Macciò, A. V., Stinson, G., Brook, C. B., et al. 2012, ApJ, 744, 1 [Google Scholar]
 Mandelbaum, R., Wang, W., Zu, Y., et al. 2016, MNRAS, 457, 320 [Google Scholar]
 Moster, B. P., Somerville, R. S., Maulbetsch, C., et al. 2010, ApJ, 710, 903 [Google Scholar]
 Moster, B. P., Naab, T., & White, S. D. 2013, MNRAS, 428, 3121 [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563 [Google Scholar]
 Navarro, J., Frenk, C., & White, S. 1997, ApJ, 490, 493 [Google Scholar]
 Persic, M., & Salucci, P. 1991, ApJ, 368, 60 [Google Scholar]
 Persic, M., Salucci, P., & Stel, F. 1996, MNRAS, 281, 27 [Google Scholar]
 Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ponomareva, A. A., Verheijen, M. A. W., Peletier, R. F., & Bosma, A. 2017, MNRAS, 469, 2387 [NASA ADS] [CrossRef] [Google Scholar]
 Reyes, R., Mandelbaum, R., Gunn, J. E., et al. 2012, MNRAS, 425, 2610 [NASA ADS] [CrossRef] [Google Scholar]
 Rhee, M. H. 1997, ASP Conf. Ser., 117, 90 [Google Scholar]
 RodriguezPuebla, A., Drory, A., & Avila, R. 2012, ApJ, 756, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Rubin, V. C., Burstein, D., Ford, W. K., Jr., & Thonnard, N. 1985, ApJ, 289, 81 [NASA ADS] [CrossRef] [Google Scholar]
 Salucci, P. 2019, A&ARv, 27, 2 [CrossRef] [Google Scholar]
 Salucci, P., Lapi, A., Tonini, C., et al. 2007, MNRAS, 378, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Salucci, P., Yegorova, I., & Drory, N. 2008, MNRAS, 388, 159 [NASA ADS] [CrossRef] [Google Scholar]
 Shankar, F., Lapi, A., Salucci, P., De Zotti, G., & Danese, L. 2006, ApJ, 643, 14 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Simon, J. D., Bolatto, A. D., Leroy, A., Blitz, L., & Gates, E. L. 2005, ApJ, 621, 757 [NASA ADS] [CrossRef] [Google Scholar]
 Spekkens, K., Giovanelli, R., & Haynes, M. P. 2005, AJ, 129, 2119 [NASA ADS] [CrossRef] [Google Scholar]
 Spergel, D. N., & Steinhardt, P. J. 2000, Phys. Rev. Lett., 84, 3760 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Tulin, S., Yu, H., & Zurek, K. M. 2013, Phys. Rev. D, 87, 115007 [CrossRef] [Google Scholar]
 Wang, J., Fu, J., Aumer, M., et al. 2014, MNRAS, 441, 2159 [NASA ADS] [CrossRef] [Google Scholar]
 Weinberg, D. H., Bullock, J. S., Governato, F., Kuzio de Naray, R., & Peter, A. H. G. 2015, Proc. Natl. Acad. Sci., 112, 12249 [Google Scholar]
 Yegorova, I., & Salucci, P. 2007, MNRAS, 377, 507 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Bestfit rotation curves and triangle plots
In this section we plot the 25 velocity bestfitting 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 bestfitting model, the dashed green line shows the NFW halo component, and the dotdashed 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^{−1} 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^{2}.
The velocity model therefore has four free parameters. Table B.1 shows the bestfitting 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 bulgedominated Sa spirals here.
Fig. B.1. Observed and fit RC modeling: The components of the NFW dark halo (green dashed), disk (dotdashed 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. 
Values of the bestfit parameters of the mass model along with their 1σ uncertainties.
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 bestfitting 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. Bestfit RC of the first bin of the LSB sample by considering the gaseous component. The dashed green, dotdashed 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. 
Coadded RC model result with gaseous disk.
All Tables
Values of the bestfit parameters of the mass modeling with their 1σ uncertainties for the PSS96 sample.
Values of the bestfit parameters of the mass model along with their 1σ uncertainties.
All Figures
Fig. 1. Twentysix doublenormalized coadded RCs. The data and their uncertainties are listed in the CDS table. 

In the text 
Fig. 2. R/R_{opt} − log V_{opt} − V/V_{opt} 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 M_{vir} 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. M_{D} − M_{vir} 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 M_{D} for KS17 has the extremely discrepant coordinates (see the text). 

In the text 
Fig. 5. Relation between c, ratio and chired. ratio ≡ log(M_{D}/M_{vir}) 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 V_{opt}. The quantity V_{opt} is derived using the relation M_{B} = −7log V_{opt} − 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 bestfitting model, the dashed green line shows the NFW halo component, and the dotdashed 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 (dotdashed 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. Bestfit RC of the first bin of the LSB sample by considering the gaseous component. The dashed green, dotdashed 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.