Review of Solutions to the Cusp-Core Problem of the CDM Model

This review aims at proposing to the field an overview of the Cusp-core problem, including a discussion of its advocated solutions, assessing how each can satisfactorily provide a description of central densities. Whether the Cusp-core problem reflects our insufficient grasp on the nature of dark matter, of gravity, on the impact of baryonic interactions with dark matter at those scales, as included in semi-analytical models or fully numerical codes, the solutions to it can point either to the need for a paradigm change in cosmology, or to to our lack of success in ironing out the finer details of the


Introduction
Notwithstanding the many successes of the concordance cosmology model, also known as the ΛCDM model, fundamental issues remain open. Most of the cosmology community considers the ΛCDM paradigm to accurately describe the Universe's large scale structure formation and evolution, its early state and the proportions of its content in matter and energy [1][2][3][4][5]. This view has been reinforced with each new discovery, such as the kinetic Sunyaev-Zel'Dovich (SZ) effect, the Cosmic Microwave Background (CMB) B-mode polarisation [6], or its lensing [7,8]... Among its most famous resisting issues, the questions of the nature, "fine tuning" and "cosmic coincidence" problems of the cosmological constant Λ remain open [9,10].
The "cosmological constant fine tuning problem" emerges from the large discrepancy between the observed value of Λ and the huge values predicted by quantum field theories for the present quantum vacuum energy, the latter exceeding the former by more than 100 orders of magnitude [9][10][11]: the cosmological observation upper bounds are limited to ρ Λ 10 −47 GeV 4 , while the naive theoretical expectations obtain ρ Λ 10 71 GeV 4 . This is the most extreme fine tuning problem known to physics. As Λ is constant, the concordance model also presents a strong fine tuning of the dark energy (DE) initial conditions. The "cosmic coincidence problem" questions the relative coincidence of the orders of magnitude at present of dark matter (DM) and DE energy densities [12].
Further issues of the ΛCDM paradigm revolve around observations of the CMB or of smaller scales structures. Unexplained statistical anomalies in the CMB large-angle fluctuations contradict the ΛCDM CMB multipoles statistical independence, which reflects the assumption that the content of our Universe should realise a Gaussian and statistically isotropic random field. These anomalies comprise detection in the CMB of a power hemispherical asymmetry [13][14][15][16][17][18], of a cold spot, marking an underdense region [19][20][21], and of a quadrupole-octupole alignment [22][23][24][25][26]. Such anomalies could either mark new physics, be caused by statistical effects [27], or reflect unknown systematics.
Furthermore, puzzling tensions have been found between the cosmological fluctuations σ 8 parameters measured by cluster number counts and weak lensing, compared with the 2013 Planck results [28], as well as between the 2013 Planck H 0 Hubble parameter, compared with SNIa data. Those tensions have persisted between the Planck 2015 data and the σ 8 growth rate [29], as well as with the CFHTLenS weak lensing [30].
In addition to those CMB issues, the description of small scale structures in the ΛCDM model also raise problems e.g., [31][32][33][34][35][36]. They count the "missing satellite problem" MSP [32,37], the "Too Big To Fail" problem TBTF [34,35], the angular momentum catastrophe AMC [38,39], the satellites planes problem [40], and the baryonic Tully-Fisher relation reproduction problem [41]. Among those small-scale problems, the so called Cusp-core problem CC [31,42] retains a central role. The CC coins the gap between the flat central slopes of density profiles in dwarf galaxies (i.e., dwarfs), Low Surface Brightness (LSBs), and Irregulars galaxies on one hand, while they all are DM dominated, and on the other hand, the cuspy central slopes of dissipationless N-body simulations density profile predictions [43][44][45].
Besides the aforementioned issues, further problems of the ΛCDM paradigm [46][47][48][49][50] are left out of this review, as we focus here on the CC problem. Note however that the formation of cusps and issues of satellites in ΛCDM are tightly connected, as shown through the mechanism in some unified models [51][52][53] that can turn cusps into cores as well as alleviate the MSP, since the shape of the parent halo determines strongly its tidal effect on satellites e.g., [54][55][56], with cored profiles even capable of stripping the satellite out of existence [56].
We organise the discussion as follows: starting with a short exposition of the initial discussions on the CC problem, Section 2 will summarise its issues, while early solutions involving baryons, DM matter type changes or gravity modifications will be presented in Section 3. A review of the baryonic solutions to SPH simulations collapse models through gas heating feedback effects from dynamical friction of baryonic clumps will be discussed in Section 4, including a discussion in Section 4.1, where the supernovae explosions [57] solutions of core heating will be presented, a Section 4.2 that will discuss both baryonic approaches and a unified baryonic solution [53] and a presentation of the resulting mass dependent density profile in Section 4.3. Finally, Section 5 will review and discuss solutions modifying the DM nature, before to conclude in Section 6.

Exposition of the Cusp/Core Problem
The CC problem emerged after Flores & Primack [42] and Moore [31] showed DDO galaxies' rotation curves could be better fitted with cored (or pseudo-)isothermal density profiles, ruling out cuspy profiles. As dissipationless CDM model simulations persisted to produce cuspy profiles (see Figure 1), this was coined as the CC problem.
Understanding of the construction of RCs was refined with their decomposition for some spiral galaxies into their components in terms of stars, gas and DM, fitted with various models in Gentile et al. [98][99][100]. The results favoured constant density core profiles over cuspy models, in the same way as the THINGS (The HI Nearby Galaxy Survey) 7 dwarf galaxies processed by Oh et al. [36]. Figure 3 compares the RCs of the latter, among others, with those of NFW and pseudo-isothermal (ISO) profiles, clearly showing the ISO fits them best. The LITTLE THINGS survey also produced similar results [101].
However, as the previous results emanate mainly from low brightness objects such as dwarfs or LSBs, while large galaxies or high-surface brightness objects inner structure determination is more complex, general statements on inner density profiles of all galaxies, in particular on their cored or cuspy nature, are far from settled. Indeed, while the high-surface brightness galaxies studied by [102] present cored centres, their analysis by e.g., [68,78,[103][104][105] disagree. Other surveys contain a mix of cored and cuspy objects. In particular, in the THING sample, galaxies with low luminosity, M B > −19, tend to follow a cored ISO profile, while luminous galaxies, M B < −19, are equally well fitted by cored or cuspy profiles.
Such confused picture is reinforced by the lack of systematically flat slopes in dwarf galaxies: in Simon et al. [103], a large scatter in the inner slopes of the low mass spirals NGC2976, NGC4605, NGC5949, NGC5963 and NGC6689, was obtained, with a cored profile α −0.01 for NGC2976, a cuspy slope α −1.28 for NGC5963, while α −0.88 for NGC4605, α −0.88 for NGC5949 and α −0.80 for NGC6689. We illustrate this in Figure 4: its top left panel confronts NGC5963's DM halo RC (black dots with error-bars) to RCs computed from (a) the [74] model that accounts for baryonic physics (yellow dashed line), (b) a fitted NFW profile (cyan line), (c) a fitted pseudo-isothermal profile (ISO, shortdashed magenta line), while its top right panel confronts the same kinds of fits to the RC from NGC5949, and the bottom panel only displays NGC2976's RC with a flat power law (black line) and model (a) (dashed line). The NGC5963 RC is well fitted by the cuspy NFW profile's. The NGC5949 RC is equally well approximated by RCs from a cuspy NFW or ISO cored profiles. Finally, the NGC2976 RC reflects a very flat inner profile (α 0.01). The three RCs displayed in Figure 4 all agree with the [76] model.  [106]. Note that the Einasto profile fits also well those simulations [44]. The same comparison for the Moore profile is displayed in the top right panel, and each panel represents simulations of dwarfs galaxies (red line), medium size galaxies (green line), and galaxy clusters (blue line). The bottom left and right panels focus on the residuals between the corresponding theoretical profiles and the compared simulations this figure is reproduced from Figure 1 in [106] .
It appears then that if a cored profile agrees with RCs of many dwarfs, some remain clearly cuspy.
In agreement with the preceding considerations, galaxies' observed RCs were recognised by [109] to exhibit much more diversity than what simulations produced.
The conclusion this brings emphasizes the difficulty of galaxies' inner slope assessment, including for dwarfs. This moreover calls for a redefinition of the CC problem in terms of inner mass, rather than inner slope of galaxies density or RCs. Indeed, as shown by several studies, including those mentioned here, and despite kinematic maps recent improvements, there exists no agreement on DM slope distributions correlated with morphologies [36,103,108], but rather a range of profiles.
This lack of agreement appears even more obvious for the extreme ends of the mass spectrum of galaxies: for larger masses, dominated by stars, as in spiral galaxies, and for lighter galaxies, such as dwarf spheroidals (dSphs), which models require the account of biases, that result in opposite consequences.
To interpret and assess these results on dSphs, a range of approaches have been tried. As stellar orbit anisotropy and mass are degenerate in the spherical Jeans equation [110], its results highly vary with assumptions. Similar degeneracies also hinder the approach applying parameter space maximum likelihood to Jeans modelling [111][112][113]. With Schwarzschild modelling, cored profiles are found [114][115][116][117], for, e.g., Sculptor or Fornax dSphs. More recently a study by [118] has found that the majority of dSPhs favor cuspy profiles. Multiple stellar populations methods also obtained cored profile [119][120][121][122] for Fornax (measuring its slope at 1 kpc) and Sculptor (for its slope at 500 pc). Nevertheless, the lack of agreement on a single inner distribution for dSphs still remains, with contradictory results such as the cusp obtained from Schwarzschild modelling for Draco [117].
Such scatter preventing a unified understanding occurs as well for galaxy clusters scales. In clusters MACS1206, MS2137-23, RXJ1133, A383, A1201, A963 and RXJ1133, studies with a combination of BCG's (Brightest Central Galaxy) stars velocity dispersion, weak and strong lensing, Sand et al. [123] found cores, and only the RXJ1133 profile could be compatible with the NFW model. Such flatter profiles were also obtained, in similar studies (see also [124]), for A611 by Newman et al. [125], 2009, for A383 by Newman et al. [126], 2011 and for MS2137, A963, A383, A611, A2537, A2667, A2390 by Newman et al. [127], 2013. This contradicts other studies such as the cuspy profile in the strong lensing and X-ray observations of A611, among seven cuspy or flat profiles of relaxed massive clusters with average slope α = 0.50 ± 0.1, by Donnaruma et al. [128].
Despite the agreement of the IC2574 galaxy's RC with the simulations from Ref. [109] (presented in Figure 5, shaded green band) for radii > 6 kpc, the galaxy's inner region behaviour differs completely from their simulations. Such divergence emphasizes that correct RC slope predictions for some inner regime doesn't entail that the successful simulation remains accurate in all regimes of the RC. As the main characteristic of the CC problem resides in the inner galaxy's mass deficit [109], the focus should lie in the DM excess in the central distribution produced by CDM rather than on the density profile shape. The observation vs prediction discrepancy already appears where the asymptotic value of the circular velocity is reached [138].
It is therefore unfortunate that practically no observational papers in the field depart from the inner slope estimation method using α. A notable exception can be found in the inner slope estimation via the integrated indicator Γ ≡ ∆ log M/∆ log r < 3 − α performed by [120] in the cases of the Fornax and Sculptor galaxy clusters. As Γ is integrated, contrary to the local α, its evaluation is easier, but provides a lesser constraint on the inner slopes.
Finally, an alternate probe of the DM inner profile flattening to rotation curves consists in the dynamical friction timing (DFT) technique, that is assessing the time an object like a globular clusters (GC) requires to move inside the distribution of mass. It clearly represents an indirect method to obtain the DM inner slope, however it presents the advantage of directly probing the density, rather that the integrated mass furnished by rotation curves. Thanks to the presence of GCs in dwarf galaxies, like Sagittarius or Fornax, Ref. [86] was able to show analytically that the inner density profile must be cored. This behaviour focused on GCs was confirmed by several simulations [139][140][141][142]. The key interests of the DFT technique are twofold: (1) it concentrates on smallest scale galactic systems like dwarf galaxies, on which the proposed solutions to the CC problem show their main drawbacks, and (2) its method is independent from that of rotation curves, giving information on the structure of the mass distribution in the studied systems. Recently, DFT studies have shifted their focus from GCs to DF on stars. Using analytical approaches, as in Hernandez [143], or numerical techniques, in Inoue [144], on ultra-faint satellite galaxies of our galaxy, was shown that DM fractions imply DFT for the stars themselves and are incompatible with cusped dark matter haloes, at these tiny scales, out of reach of SNF models, and even maybe of DFBC solutions.
Since smaller mass objects are more likely to harbour inner profiles in agreement with predictions of dissipationless, N-body, codes (cuspy), a definite verdict on the nature of the inner slope, cored or cuspy, of the smaller end of the mass spectrum, i.e., dSphs, is critical.  [72]. The green shaded region and line, as well as the black line, respectively stands for the scatter of one and both median RCs for simulated galaxies from which this figure, Figure 1 is extracted from Ref. [109].

Initially Proposed CC Problem Solutions
In the wealth of propositions to solve the CC problem, the initial attempts were determined by the limited state of knowledge and observations. Systematic effect and poor resolution were first invoked as the origin for the CC inconsistencies e.g., [89,145]: they blamed beam smearing, slit-misplacement, non-circular motions or beam off-centering, as they are sources of systematically lowering slopes.
For instance, finite HI emission beam size tends to smear the HI observations, inducing a larger measured disk, depending on the HI distribution, beam size, intrinsic velocity gradients and galaxy inclination angle. With modern high spatial resolution (<1 kpc, see below), the smear can be resolved.
The slit misplacement for an Hα emitting galaxy can entail missing the galaxy's dynamic centre, inducing a flatter measured profile. One of the solutions to avoid this slope underestimate have used 3d spectroscopy [94,102]. As Hα observations assume circular motions for the gas, any non-circular orbits underestimate the slope. Some have found [97] that those departures from circular velocity remain within the order of a few km/s. Current high resolution rotation curve data are able to discriminate cuspy from cored haloes with asymptotic inner slopes [84].
Projection effects were put forward by several papers [92,106,146] to explain the discrepancy of their simulations with observations, despite their claims of good agreement of their profiles, becoming shallower in the range the virial radius towards the centre. Their argument lies on the practice of representing the 3D DM haloes motion with a spherically averaged circular velocity, traced by the gas disk rotation speeds. Their claim is to blame the disagreement with data on the correspondence fitting formulae rather than their simulations [92]. Current DM-only high resolution simulations clearly reject this explanation, as their inner profiles cannot reach shallower slope than −0.8 [70], compared with high resolution galaxy observations revealing much smaller values.
More generally, simulations limitations, such as over-merging or lack of relaxation or resolution, have been blamed to source the CDM model's failure [82,147,148]. This was already questioned then by Ref. [149]'s convergence tests proving the dissipationless N-body simulations of the time were adequate to properly obtain the CDM density profiles, despite their neglecting the galaxies inner kpc regions effects of baryons. Such effects were shown to dominate the DM generated profiles in clusters' central 10 kpc regions [123,[125][126][127]. Modern simulations no longer suffer from the blamed limitations and hydrodynamic cosmological high resolution models now account for baryons, as will be discussed in the next section on baryonic resolutions of the CC Problem.
Questioning of the very paradigm of CDM opened the door to speculations on the nature of DM, such as warm WDM [150], fluid [151], repulsive [152], fuzzy [153], decaying [154], self-annihilating [155], or self-interacting [156] DM. More fundamental questioning led to altered small scale DM power spectrum e.g., [157] or modifications of gravity itself, such as the so-called f (R) [158,159] or f (T) [160][161][162][163] theories, or the so-called Modified Newton Dynamics MOND, [164,165] 1 . Such paradigmatic modifications, whether touching DM nature, initial conditions or the nature of gravity, will hereafter be denoted, in their impact the CC problem, as "cosmological solutions", and developed in Section 5.

Baryonic Resolutions of the CC Problem
Before resorting to the above mentioned "cosmological solutions" to deal with the CC problem, a thorough check, whether neglected or poorly understood local physics could explain it, would be wise. Indeed, discarding the concordance, i.e., ΛCDM, model would risk to significantly alter the otherwise very successful predictions ΛCDM provided, thus loosing the explanations for many observations of our Universe. This is the aim of the current section.
Keeping the framework of the ΛCDM model, "astrophysical solutions" alter the results of pure CDM structure formation with some "heating" mechanism applied to inner DM distribution, originated in the effects of baryons. The inner density profile gets then flatter than for pure CDM evolution. Such "heating" mechanism could be provided by 1. accounting for the action, in structure formation, of angular momentum; 2. the induced dynamics from a central black hole, claimed to shallow the DM cusp [170][171][172]; 3. the modification of equilibrium from the presence of a rotating bar; 4. the mechanism of dynamical friction that induces baryons-DM transfer of angular momentum [74,173,174]; 5. the supernova (SN) explosions-generated bulk motions in gas, dubbed AGN feedback [54,55,72].
Dynamical friction was shown to flatten, and even erase, the DM formed cusps in both galaxy clusters and dwarfs by transferring energy from baryon clumps to DM [173,174]. Conversely, DM adiabatic contraction (AC) in haloes steepens their density profiles [193,194]. 4 The above effects see also [196] were all combined in the model from Del Popolo [74], including: -tidal torques caused ordered angular momentum acquired during the structure formation; -its obtained random angular momentum; -its dynamical friction induced baryons-DM energy and angular momentum exchanges; -and the effects of its DM adiabatic contraction (hereafter, AC).
In addition to inducing flatter density profiles in the center of galaxies and clusters, dynamical friction and angular momentum also tend to modify their overall internal structure [197][198][199][200]. The simple confrontation of real structures, including baryons, with the outputs of dissipationless simulations is demonstrably incorrect [74]: the non-negligible action of baryons in the inner regions of structures explains the failures of N-body simulations to reproduce observed density profiles. This is illustrated in Figure 6, showing halo evolution for two masses, 10 9 and 10 14 M , modeled using the "Dynamical Friction from Baryonic Clumps" scenario (hereafter DFBC scenario) from Ref. [74] and discussed below.
The impact of baryons content, environment and formation history on dwarf galaxies density profiles was demonstrated with the Del Popolo DFBC model [76].
The "heating" mechanism 5 was first put forward by Navarro et al. [201], flattening the inner profile from the feedback of SN explosion into the medium, as discussed below.
They will be discussed in more details in Sections 4.1 and 4.2, and are put in context below. In general, baryons and DM usually "interact" in smooth distributions via AC, that can be counterbalanced by transferring energy back from baryons to DM. AC leads DM to collapse towards the centre of the halo, steepens its profile and increase its central density. The counteracting transfer has the opposite effect. It can occur via 1. Dynamical friction transfer of incoming clumps orbital energy to DM.
Such source expelling from, or relocating baryons in, the halo (for instance, bulk motions induced from SN explosions, [54,55]) can temporarily flatten the gravitational potential, shift DM particles outwards, and flatten the cusp.
We will now discuss the two favoured "heating" mechanisms in the CC problem solution, and how those solutions are synthesised into the proposal of a mass dependent density profile for DM haloes with baryons.
Profile flattening by a single sudden baryons expulsion into the halo was first shown by Ref. [201], and found to be most effective in dwarfs and other shallow potential structures. However, energy from a single explosion was shown [214] to remain insufficient to drive the profile to a core, while repeated moderately violent explosive events would be required see [215] for a different point of view.
A more detailed study of the [201] approach, focused on evolving from NFW profiles to DDO154's RC via a gas outflow event modeled from an abrupt change in the disk potential, by Gelato and Sommer-Larsen [202], established that at least 75% of the disk mass should be expelled that way to generate this RC.
In the case of larger galaxies and deeper potential structures, a core was shown by Read and Gilmore [203] to result from repeated and alternate outflows and re-accretions of gas.
In primordial galaxies, a core was determined to follow from SN-explosions-driven gas random bulk motions [54,55], as confirmed in [72] simulations. The average slope of the latter simulations were compared with observations of THINGS dwarfs in Refs. [36,85]. Producing simulations for larger galaxies than [72], and comparing them with observations [73], correlation between the stellar mass M * and the inner slope were obtained for M * > 10 7 M galaxies were obtained. These results were obtained from simulated galaxies thanks to the N-Body+SPH code GASOLINE [216]. With an 86 kpc softening, gas particles and DM particles resolution were reduced to M p,gas = 3 × 10 3 M and M p,DM = 1.6 × 10 4 M , reps., thanks to the "zoom" technique of Ref. [217]. Two realisations were produced, distinguished by their star formation thresholds: a High Threshold run (HT), with star forming hydrogen cloud densities >100/cm 3 , and a Low Threshold run (LT), with hydrogen densities >0.1/cm 3 . Both SNF blast wave mechanism [218] and/or early stellar feedback [219] were implemented in those simulations, as later by [69] in the same way. The Governato simulations resulting energy feedback to the interstellar medium (ISM) amounted to 10 51 ergs from >8M stars, while the ISM gas ejecting energy was coupled to the SN expelled energy via the coefficient esf . The fiducial esf = 0.1 was selected for the MaGICC simulations [219].
Comparably, a combination of bursty star formation and SNF, inducing inner (1 kpc) halo potential fast oscillations and expanding gas bubbles when the central cold gas density reaches >100/cm 3 and stars start to form 5 , was shown to flatten the cusp [221]. However in this model, no noticeable modifications of the DM inner density profiles are induced for smaller densities (e.g., 0.1/cm 3 ). The case of larger mass galaxies was treated by Governato et al. [73] with the same method (also see Figure 7 top panels). Agreement with the results of [221] (Figure 7 bottom left panel) was obtained from introducing a new SNF scheme in the RAMSES adaptive mesh refinement code by Teyssier et al. [212], showing flat galaxies inner profile for M * > 10 7 M . Similar outcomes with a lower threshold at M * > 10 6 M were produced by Onorbe et al. [222] (Figure 7 bottom right panel) 6 .  Since their simulations found little effect of baryons on rotation curves of galaxies with V max < 60 km/s, even in their inner regions, this perspective was denounced by [109]. They argued that the cores obtained by [221], since no evidence of core formation was detected in their [109] work, were inherently reflecting an ad hoc choice of parameters.
Such understanding agrees with the findings of Ref. [223], which 7 high-resolution dwarfs simulations, included in halos with 1-2 ×10 10 M and following different assembly history, could not reveal any inner core flattening. A very extended assembly history dwarf, thus also having an extended star formation rate (SFR) history, exhibits their lowest detected inner slope: they measured it to reach, at 0.01-0.02 R vir , up to −0.8. It also remains consistent with the core-less simulated realistic galaxies of Ref. [224].
In addition to those contradictions in core formation between various high resolution hydrodynamical simulations, several criticisms have arisen on the methods and assumptions used in simulation models of the types of [72,225], as they obtain faulty 1. core formation energetics: the number of stars, obtained in their M * < 10 7 M galaxies, generates an insufficient energy to flatten the galaxies cusp [226]. Moreover, their required core-forming baryonic mass marginally exceeds the dSphs observed baryon content [215]. That problem in Penarrubia [226] is illustrated in Figure 8 , presents: (a) the range of minimum SN explosions energy required to obtain a core of size r c given within 0.1 < r c /r s < c, with c and r s , respectively being the concentration parameter and the scale radius of the NFW profile at given virial mass M vir , as the red shaded area; (b) for a fixed core size r c = 1 kpc, the dotted black line is obtained; (c) the SN explosion energy outputs, indicated from two different studies [227,228] as dot-dashed green and dashed blue lines, denote those compatible with the "missing satellite" [32,37] problem and reveal, for haloes with M vir < 10M , the tension with the "core/cusp" problem. The right panel, extracted from Figure 3 in [229], shows: (a) the left panel red shaded area of Penarrubia et al. [226]'s energy estimates, as shaded grey area; (b) solid lines with symbols for the conversion energy of cuspy to pseudo-isothermal density profiles at fixed core sizes, indicated in the legend by the symbol types; (c) the solid black line for the ∆E scaling with halo mass M h , while the ratio of the cusp mass redistribution limit radius, r m , over the halo radius, r h , is fixed; (d) the dotted, dashed, and dot-dashed lines, respectively, for the Behroozi et al. [230], Kravtsov [228], and Moster et al. [231] M − M h relations.
2. energy coupling: their SN 0.4 exceeds the 0.05 coupling that studies such as [232] deduced; 3. star formation threshold: the [72] results have been found in [109,233] to require a very high star formation threshold; 4. tensions with TBTF solution: those models simply encounter adverse conditions [215,234,235] when trying to solve the TBTF problem; 5. feedback resolution: accurately describing the cusp transformation into a core via their feedback processes would require better resolution than they present [236][237][238][239].
Despite the SNF limitations, 8 for M * > 10 6 M , correlations between M * and galaxies inner slope have been found, comparing THINGS dwarves average slopes with [72]'s simulations in Refs. [36,85], and similarly in [73] for larger objects.
Inversely, hydrodynamic simulations for M * < 10 6 M produce cuspy profiles, in contradiction with the Fornax and Sculptor cored inner structure found by Ref. [120] thanks to their more reliable DM haloes' slopes from mass profiles Γ ≡ d log M d log r < 3 − α. Only with 50% more SN energy injections could Ref. [240] claim agreement with [120]'s cores.
Assuming SN feedback and other specific assumptions, a Semi-Analytic model from Gnedin and Zhao [214] claimed impossibility to obtain cores. As, among details of stellar feedback, the core formation impact of baryonic clumps, such as in Ref. [74], is ignored, such objection is weakened.
However, there exists one core-forming SNF dwarf galaxies simulation for masses < 10 6 M : with the P-SPH mode in the GIZMO code, Ref. [222] obtained cores in a much less natural way than in the DFBC scenario, discussed in the next section.

Discussion of Gas Clumps Dynamical Friction
As discussed in the introduction of the present Section 4, Refs. [173,174] proposed an alternate path from SNF for the turning of cusps into cores: baryon clumps with mass ratios to the system above 0.01% can transfer, through dynamical friction (DF), clumps orbital energy to DM. Such heating induces similar effects to SNF without the need to wait for a full massive stellar cycle, moving to outer orbits the central DM particles and thus flattening the inner DM density profile. It is more efficient on earliest, smallest haloes of the hierarchical structure formation.
DFBC was first shown to function in galaxies [173], then in clusters [174]. Cluster C0337-2522 evolution with such model [205] showed the BCG generation precedes DFBC inner slope formation, resulting in flatter than NFW inner profiles with 0.49 < α < 0.90. Implemented in a hybrid N-body/SPH simulation, DFBC galaxy evolution showed heated up cusps by baryons subhalos into a core of 1 kpc, and was compared with DM-only systems and mixed DM/baryons systems [206].
Refs. [208][209][210]'s simulations confirmed such outcomes. The main steps of the DFBC can be outlined in the following: • the linear phase develops from initial DM and diffuse gas proto-structures.
• As in local star-generating systems and following the Kennicutt-Schmidt law, a few percents of those clumps' gas converts to stars [251], while DF migrates them to the DM halo (galactic) center. • The first collapse phase's AC compressed baryons [193,194] and e.g., 10 9 M galaxy at z 5 in our Figure 6, reenforcing the DM cusp. • The previous effect is countered by DF-induced migration of clumps to galaxy center, as DF transfers energy and angular momentum from baryons to DM, heating up the cusp into a core. • SNF only intervenes later (e.g., around z = 2), when SN explosions expel gas, decreasing surrounding stellar density. Moreover, such feedback destroys the smallest gas clumps as soon as a small fraction of their mass turns to stars. 9 Although sharing some common features (e.g., clumps to DM energy transfer via gravitational interaction), leading some to consider SNFF and DFBC as the same idea implemented differently, they fundamentally differ on two crucial points: 1. their respective profile flattening epochs differ markedly: while pure SNFF DM profile does not depart noticeably from NFW at z 3 [222,225,252], DFBC flattening initiate at higher redshifts (z < 5, see Figure 6); 2. their clumps moving energy sources diverge in nature: in Refs. [54,55]'s definition, while clumps motion is driven, in the SNFF model, by SN explosions energy, the DFBC views them to just "passively" infall to the halo centre.
The two competing scenarios predictions for the central density slope vs. stellar mass and vs. circular velocity relationships were confronted in Refs. [253,254] with a series of high resolution data predictions from the Milky Way, Sculptor, Fornax, LITTLE THINGS [101], THINGS dwarfs [36,85] and Refs. [103,108]. Such confrontation showed slightly better performance from the DFBC model compared with the SNFF scenario, as the former predicts core emergence at smaller stellar masses than 10 6 M , limiting the latter. Although very small dwarfs (M * ≤ 10 4 M ) cannot produce cores in DFBC, this nevertheless agrees with results from Ref. [255].
Both competing model face the same caveats in their resolution of the CC problem: galaxies containing a bulge have been demonstrably able to reconstitute their cusp [79]. This effect can even been found in dwarf galaxies [237].

Mass Dependent DM Density Profiles
Assuming that baryons are indeed responsible for some cored DM haloes, while allowing some others to be cuspy, and using the results of SPH simulations, the SNFF or DFBC models, some authors synthesised the diversity of DM density profiles and generalized the usual density profiles, giving the mass dependent empirical law that would solve the CC problem. For example, Ref. [256] generalized the Einasto profile into a mass dependent parameterised profile, encompassing both cuspy and cored cases, while [257] obtained such mass dependent profile by modifying the Zhao [258] profile. In the following we will give the details concerning the works of [257,259].
In previous discussions, we saw that the DM density profile is well approximated by the NFW law, and even better fitted by the Einasto profile. Those two profiles are obtained by DM-only simulations. There, the hydrodynamical processes, fundamental in the determination of the structure of the density profile, are neglected. For example AC [193,194] produced by gas cooling, which strengthens cusps, is not taken into account by N-body, DM-only simulations. At the same time, the presence of baryons can expand haloes. As discussed in Sections 4.1 and 4.2, stellar feedback and DFBC can produce a halo expansion and a flattening of the inner density profile. As discussed there, Ref. [69] showed how the inner slope is modified by SNF. The results of this paper were used to obtain a DM halo whose slope depends from the ratio between stellar and DM masses, The generalised NFW profile (gNFW) is a specific form of a law depending on three parameters (α, β, γ): being r s the scale radius and ρ s the scale density. −γ and −β, are the logarithmic slopes of the inner and outer regions respectively, while the sharpness of the transition from inner to outer region is specified by α. The peculiar case (α, β, γ) = (1, 3, 1) corresponds to the NFW profile. The goal of [259] was that of making the profile given in Equation (1) mass dependent, in the sense that its three parameters would depend on the ratio M * M halo , using haloes obtained, in [69] from SPH simulations. The density profile of the haloes were fitted, using Equation (1). The behavior of the three parameters was captured using fitting laws as a function of M * M halo . The outer slope β was modeled with a parabola. In the cases of the inner slope, γ, and of the transition parameter, α, double power law model were used.
The best fits for these parameters to the values for SPH galaxies are shown as dotted lines in Figure 9. The functional forms of the fits are:  Given the stellar-to-halo mass ratio of a galaxy, by means of Equations (1) and (2), it is possible to compute the entire DM profiles.
Ref. [259] also showed that the concentration parameter is dependent on M *  Figure 10. Each simulation is represented by its symbol and size as described in Table 1 of [259]. The dependence of the ratio between concentration parameters obtained from SPH and pure DM N-body simulations, C SPH /C DM , on M * M halo is nearly exponential. Its best fit is: where X + = log 10 ( M * M halo ) + 4.5. C DM has been derived fitting a NFW profile to the DM-only version of each galaxy, while C SPH has been obtained from hydrodynamical runs where their galaxies DM halo model profile had their scale radius correspondence converted between r s into r −2 .
A comparison of the model with simulated data of the density profiles, and rotation velocities gives a good agreement.
The mass dependent density profile idea from Ref. [259] was applied by several other authors. As mentioned above, it was applied by Lazar et al. [256] to the Einasto profile and by Freundlich et al. [257] to the so-called Dekel-Zhao density profile, a peculiar case of Equation (1). We now focus on another popular profile, similar in performance to the Einasto profile, the modified NFW model, adopted by Freundlich et al. [257] in the form from Zhao [258] where the radial dependence is encoded in x = r/r c , with respect to the characteristic radius r c , and the density scale is set by the typical density ρ c . It is equivalent to the generalised NFW profile, noting the parameters α = 1/b, β = g and γ = a. Note that, for any natural numbers n and k, setting g = 3 + k/n and b = n, it is possible to integrate the mass and obtain the corresponding equilibrium velocity dispersion and gravitational potential analytically. Some choices of integers actually yield very good DM simulations fits, and are even able to accommodate baryons' presence, as well as cusps or core. The profile referred to as the Dekel-Zhao (DZ) profile, found for k = 1 and n = 2, which corresponds in Equation (4) to b = 2 and g = 3.5, and leaves two free parameters, a and a concentration parameter c, have been found, using the NIHAO suite of simulations [260], particularly efficient at representing cored profiles [261], but also cuspy profiles.
Following Equation (4), it takes the form using x = r/r c and defining the concentration c = R vir /r c from the virial equilibirum radius of the spherical mass distribution R vir , where the average density at that radius computes as ρ vir = 3M vir /4πR 3 vir , while the shape factor µ = c a−3 (1 + c 1/2 ) 2(3−a) is involved in the profile mass integration, yielding the average critical density ρ c = c 3 µρ vir , and thus relating the density scale to the shape parameters, c and the central slope a, as ρ c = (1 − a/3)ρ c .
As measurements are made at fixed given resolution, the central slope is actually an integrated slope with respect to the actual profile. Freundlich et al. [257] explicits, for a fixed resolution r 1 , an inner logarithmic slope and the corresponding averaged concentration parameter The physically motivated restrictions to positive density and inner slope result in the conditions a ≤ 3 and a + 3.5c 1/2 (r 1 /R vir ) 1/2 ≥ 0, respectively. One can reobtain the intrinsic parameters (a, c), employed in analytic expressions, from the observed (s 1 , c 2 ), usually bestowed to numerical tests, by inverting Equations (6) and (7), as and Equations (18) and (21) and Section 2.3 in Freundlich et al. [257] respectively provided analytic expressions for the DZ profile's gravitational potential, velocity dispersion and lensing properties, as well as, using some NiHAO SPH simulations suite's profiles (see their Figure 3), compared the DZ, gNFW and Einasto profiles through fits, assuming their parameters free, and finally finding the DZ profile performs significantly better than Einasto's and marginally better than the gNFW on NiHAO.
In this context, the mass dependent, DZ-based profile from Freundlich et al. [257], was constructed, as in Di Cintio et al. [259], from the NIHAO simulations suite [260], to set the DZ parameters functional forms. SPH simulations were fitted by Di Cintio et al. [69] with a similar profile to Equation (4), to set the free parameters forms as a function of M star /M vir , as did Lazar et al. [256] with the Einasto profile and Freundlich et al. [257] with the DZ profile. However, the mass dependence introduced, similarly to Equations (2) by Di Cintio et al. [259], to improve galaxies and clusters density profiles fits, such as in the modified Einasto profile from Lazar et al. [256] or in the modified NFW from [262], does not allow to keep any analytic form for the mass, as detected through lensing, velocity dispersion or gravitational potential. On the contrary, the DZ profile keeps that possibility intact. This is the advantage of the DZ profile fit by Freundlich et al. [257], applied to simulated NIHAO haloes' density profile logarithms, optimising the parameters, spacing the N 100 grid profile radii r in the range 0.01 R vir -R vir , through a least-square minimization. Setting the virial mass M vir , the profile can set the scale by extracting the corresponding R vir , leaving the remaining free parameters a and c. The star distribution can be samely fitted to get the star mass at virial radius M star and obtain the corresponding M star /M vir . The theoretical slope and concentration are related to the integrated (observed) ones through Equations (6)- (9), respectively, as shown in Figure 11.  (47) and (48) of [257] and their best fits, Equations (45) and (49) on their Figure 9, here marked as plain black lines. The best fit residual rms σ is shown in gray. The dashed black line in the c 2 /c DMO panel represents the best-fit function from Di Cintio et al. [259]. Color coding, indicated on the right side in each panel, for each halo compounds the information of the other panel.
This fitting procedure allowed to obtain, for each NIHAO haloes, their s 1 and c 2 as a function of M star /M vir , as seen Figure 8, black dots in Freundlich et al. [257].
The obtained s 1 and c 2 behaviour were synthesised in two parameterised functions of , where x 0 , s , s , and ν are ajustable and which used values were given in Table 1 of Freundlich et al. [257] s 1 (x) = s while the parameter c is introduced for The resulting mass-dependent profile presents the advantages, over the Di Cintio et al. [259] approach, to require less free parameters as well as to allow integration into analytic expressions for lensing properties, dispersion velocity and gravitational potential.

Cosmological Solutions to the CC Problem
The CC problem, together with the other small scale problems encountered by the ΛCDM model, open the question of a more general failure of the CDM's paradigm. This justifies models in which the nature of DM or gravity itself are modified. Such models encountered various degrees of success in their already widely checked instances. We will here focus on alternate models of DM. Several approaches have been used and we will detail the two main avenues, before discussing other possibilities.

Warm DM
Warm DM (WDM) presents the simplest departure from CDM, endowing DM with a small velocity dispersion (σ 100 m/s nowadays), also linked, for thermally produced DM, to a small DM relic mass [150,264]. Since velocity dispersion is expected to decrease with time, it should be more important in the past and thus smear small scale structure formation. Such smear induces similar DM "heating" as the baryonic solutions of Section 4: in this case the WDM particles higher velocity than for CDM smears the haloes to flatter profiles as well as produce fewer low mass haloes. These effects of WDM on structure formation have been copiously simulated e.g., [252,[265][266][267].
Despite solving the CC problem for haloes with their corresponding scales on a case by case basis by tuning its particle mass, WDM does not manage to solve it with one mass for all galaxies, nor in the entire CC problem mass range [268]. As shown in Ref. [269], neither WDM or pure collisionless CDM models can match observed disk galaxies rotation curves, while they can agree with hydrodynamical simulations including baryons (see Section 2's discussion).
In addition, Ref. [265], for a m = 2 keV thermal relic, and several other authors e.g., [270][271][272], shown that, for instance, the strong-lensing subhalo fraction is too high compared with WDM produced subhaloes. Furthermore, WDM thermal mass displays tensions between expected core size (1 kpc core corresponding to 0.1 keV) and large scale structure m 1-2 keV leads to 10-20 pc cores, see Figure 12 and Ref. [273]. In conclusion, because its power spectrum falls off too sharply, WDM does not improve on CDM [274]. Finally, since WDM structure formation is also shown to modify the Lyman-α forest [275], the situation is even worse, as WDM solution to the CC problem in our vicinity cannot, at the same time, be consistent with observed high redshift Lyman-α forest.  [273]. The WDM core formation is illustrated against the CDM density profile (solid black line) by showing profiles from a range of five WDM particle masses, from m = 2 keV (WDM1, as per the legend) to 0.05 keV (WDM5), in the left panel. The right panel synthesises the WDM core radius obtained for its given particle mass with other constraints: while the solid black line marks the core size dependence at fixed halo mass, the shaded area limits the plane to the cosmological constraints for 0.15 < Ω m < 0.6 and the vertical dashed line sets the upper limit from large scale observations.

Self-Interacting DM
The second simplest variation away from CDM, coined Self-interacting DM hereafter SIDM [156], consists in granting to DM a self-interaction cross-section, within the same magnitude of nucleon-nucleon interactions ( (m/g) −1 cm 2 ). 10 The resulting elastic scattering in regions of highest density (inner galactic regions) redistributes energy and angular momentum so that tri-axiality is reduces and a Burkert profile core forms [276].
The CC problem has been claimed to be solved in dwarfs, MW-sized galaxies, and galaxy clusters by SIDM from some cosmological simulations [127,[277][278][279] consistent with merging galaxy clusters observations [280][281][282], 0.1-0.5 cm 2 /g cross sections. Figure 13 presents simulations of haloes from galaxies to clusters masses, with σ/m set within two values, in the case of Ref. [278]. Although SIDM appears to solve the CC problem, since its cored subhaloes are more sensitive to disruption and tidal stripping than CDM's, and to improve on WDM [278,279], as enough subhaloes survive, a different point of view still remains [268].
The SIDM's attraction does not entirely resides in its structure formation effects. As "hidden sector" particle models naturally produce them e.g., [156,283,284], it also retains appeal from the particle theory point of view.  Figure 4 in Ref. [278]. Two SIDM models, designated in the legend as SIDM 1 and SIDM 0.1 , are shown, using DM particle masses σ/m = 1, marked with blue stars, and 0.1, characterised with green triangles, respectively. They are ovelayed with the NFW, plotted as a black line, and the Burkert, traced in blue, density profiles. The core radii of Burkert's profiles are located with arrows.

Other DM Models
Other proposals to change from the CDM model are more or less loose declinations altering around the SIDM concept. Here is a list of the main variations: RDM switching interaction to negative scattering results in repulsive DM RDM, Ref. [152].
The next two variations from SIDM are very commonly related to indirect DM detection: SADM if self-interaction produces a DM particles annihilation, the result is coined Self-Annihilating DM SADM, Ref. [155] proposing cross section-velocity σv 10 −29 (m/GeV) cm 2 . Such annihilation decreases the dense regions' particle numbers, in particular in the halo's centre, reducing central gravity, therefore allowing to expand central particles' orbits and thus to flatten the central profile. At the same time, annihilation results in radiation emissions that are possibly detectable.
DDM Alternately, self-interaction can result in DM decay into relativistic particles, designated as Decaying DM DDM, Ref. [154]. Their gravitational effect on structures is similar to SADM, as they equally deplete galaxies' central density, since the relativistic particles escape away, while larger scales structures behaviour remain similar to CDM. Similarly, the relativistic decay products produce radiations that can also be detected.
The last main variations are more loosely related to SIDM: BCDM Superfluid behaviour of non-relativistic, massive boson condensates in haloes centre can also result in smoothing down their profiles from cusp to cored [285]. Recent structure formation simulations [286] demonstrated that the Bose condensate DM (BCDM) small scale gravity vs. uncertainty principle opposition reduced substructures and flattened the density profile at those scales, while producing larger scale structures indistinguishable from CDM outputs.
SFDM As scalar field condensation also forms a Bose condensate, similar flat galaxies inner profiles were obtained from such implementation of the BCDM [287], that was called Scalar Field Dark Matter (SFDM) FDM Further implementation of the gravity vs. uncertainty principle exploited the waveparticle duality to obtain Fuzzy DM FDM, Ref. [153], using galactic core sized Compton wavelength, ultra-light (m 10 −22 eV) scalar particles. As they cannot be "squeezed" below their Compton wavelength, they also develop flatter profiles and form less substructures.

Conclusions
Since the first assessment by Öpik of the amount of matter surrounding our solar system through vertical star motion with respect to the plane of the ecliptic more than a century ago [288], current data indicates a dominant amount of matter should be in the form of DM. However, such dynamics and geometric data has still not been confirmed, to date, by direct or indirect DM particles detection [4,5].
Grasping the essence of DM from a gravitational point of view should privilege the study of galaxy satellites, dwarf galaxies and other small scale structures, dominated by DM, all the more as they present significant departures from the ΛCDM model predictions.
In this review, we focused on the departure constituted by the CC problem, formulated more than twenty years ago. Although we now have some competing explanations, there is no concensus on the CC problem's cause.
In this review, we have discussed the focus on baryonic effects on small scale structures. We recalled that the DFBC model can explain cores in M > 10 5 M dwarfs from baryon clumps interactions with DM, while cusps form in MW-type galaxies [253]. Although the SNFF model gets similar results, its prediction cover structures in a more restricted range M * > 10 6 M [73]. Note that both models, by depleting the DM density in central halo regions, also explain other small scale problems such as the TBTF or the subhaloes abundance problems. However, baryonic solutions to the CC and other small scale problems are still debated, as some issues remain, such as with the SNFF in classical MW dwarfs [215,226], or the possible future discovery of M * < 10 6 M dwarfs with flat inner profiles which would drive the SNFF model to conclude that DM is actually not cold 11 , or even that no baryonic physics solutions clearly appear for the TBTF problem in isolated galaxies. Escaping this debate would entail understanding how the SIDM model or one its variants, as we have discussed, could solve the CC problem. Future surveys, including measurements of star velocities, would anyhow be required to further our understanding of the CC problem and other small scale issues.
From an observational point of view, as discriminating cores from cusps in dSphs remains delicate, the discussion on their structure persists see the introduction and, e.g., [289]. Although the Subaru Hyper-Supreme-Camera [290] and GAIA [291] were foreseen as providing possible resolution for the question, their abilities would not exceed larger dwarf galaxies, of the size of Sagittarius [292]. Based on Jeans equations, their determination of density profile using stars line of sight velocity and 2D projection radius are plagued with degeneracy with the anisotropy parameter. Notwithstanding the challenges GAIA faces from it [292], improvements to measure dwarfs's stars proper motions of the [293], that yield density slope at half-light radius, were recommended, that focused on only one component out of the three velocities and on just two components out of three positions see [294]. However successful these efforts may be, understanding dwarfs inner structure will be crucial. In contrast, the difficulties of the SNFF tensions do not appear fundamental: baryons, through the more efficient DFBC model, can still turn cusps into cores sooner than gas turns into stars. Consequently, should cored dwarfs with M * < 10 6 M be discovered, despite their inconsistency with the SNFF, provided M * > 10 5 M , the ΛCDM model correct with baryonic effects of the DFBC model would remain consistent with observations.
Beyond the CC problem and the other documented issues in small scales of the ΛCDM model, further testing could involve • checking the predicted large subhaloes number within galaxies' virial radius, barely agreeing with anomalies in the gravitational lensing fluxes [295]; • measuring subhaloes perturbations of the MW's cold tidal streams see [296].
Finally, the nature and existence of DM would find a definitive solution in detection of DM particles, whether in colliders, direct or indirect means. Despite the hopes raised by a diphoton excess decay at 750 GeV at 3.9 σ significance in the LHC, as it was later ruled out as a statistical fluctuation, no evidence of super-symmetry (SUSY), that would allow WIMPS particles DM, have so far emerged and collider detection of DM is yet elusive. Although controversial, the annual modulation claimed in DAMA/LIBRA/CoGeNT remains the only possible direct DM particle detection. To date, no incontrovertible direct, indirect or collider evidence for the nature or existence of DM has been announced.
Author Contributions: All authors have contributed equally in this work. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest. 1 Generalisations of Einstein's General Relativity, f (R), and f (T) theories modify gravity in distinct ways. Replacing, in their Lagrangian, the Ricci scalar by a function f (R), the first type dates from the Buchdahl [158] introduction in 1970 but picked up interest with Starobinsky [159] a decade later. Developed from the Teleparallel Equivalent of GR, f (T) theories also modify their Lagrangian from the torsion scalar to a function. On one hand, both types were introduced to absorb the origin of cosmic acceleration into gravity and dispense from the need for DE [160]. On the other hand, MOND, introduced to fit galaxies rotation curves by Milgrom in 1983 [164,165], purports to replace DM. The DM effects of the MOND Newtonian modification also lead to other GR generalisations with the same aim [166][167][168][169] 2 A General Relativistic version of secondary infall models is at the root of works by a group including Mimoso and Le Delliou [188][189][190][191]. 3 A proto-structure's peak height measures the ratio of its central peak overdensity δ(0) to its mass variance σ see [192], as ν = δ(0)/σ. More massive objects are characterised by larger ν. 4 As computed with iterative methods e.g., [195]. 5 This process agrees with the assumptions from Ref. [72]; the bulk gas flows start beyond >10/cm 3 [220]. 6 Essentially, M * < 10 6 M galaxies are not able to turn cusps into cores under the SNF mechanism. 7 "Too Big to Fail" is used in the context of Milky Way (MW) satellite simulations producing bigger objects than observed MW satellites without any mechanism explaining why such object would fail to be detected. 8 Recall footnote 11. 9 Star formation is not an efficient process. 10 Note that 1 cm 2 /g 1 barn/GeV, so multiplying by the DM particle mass yields the cross-section. 11 Note that such conclusion is not true for the DFBC model.