The formation of cores in galaxies across cosmic time – the existence of cores is not in tension with the ΛCDM paradigm

ABSTRACT

The ‘core-cusp’ problem is considered a key challenge to the ΛCDM paradigm. Haloes in dark matter only simulations exhibit ‘cuspy’ profiles, where density continuously increases towards the centre. However, the dark matter profiles of many observed galaxies (particularly in the dwarf regime) deviate strongly from this prediction, with much flatter central regions (‘cores’). We use NewHorizon (NH), a hydrodynamical cosmological simulation, to investigate core formation, using a statistically significant number of galaxies in a cosmological volume. Haloes containing galaxies in the upper (M_⋆ ≥ 10^10.2 M_⊙) and lower (M_⋆ ≤ 10⁸ M_⊙) ends of the stellar mass distribution contain cusps. However, Haloes containing galaxies with intermediate (10⁸ M_⊙ ≤ M_⋆ ≤ 10^10.2 M_⊙) stellar masses are generally cored, with typical halo masses between 10^10.2 M_⊙ and 10^11.5 M_⊙. Cores form through supernova-driven gas removal from halo centres, which alters the central gravitational potential, inducing dark matter to migrate to larger radii. While all massive (M_⋆ ≥ 10^9.5 M_⊙) galaxies undergo a cored-phase, in some cases cores can be removed and cusps reformed. This happens if a galaxy undergoes sustained star formation at high redshift, which results in stars (which, unlike the gas, cannot be removed by baryonic feedback) dominating the central gravitational potential. After cosmic star formation peaks, the number of cores, and the mass of the Haloes they are formed in, remain constant, indicating that cores are being routinely formed over cosmic time after a threshold halo mass is reached. The existence of cores is, therefore, not in tension with the standard paradigm.

1 INTRODUCTION

The study of dwarf galaxies is a key area of research in the current astronomical landscape. Dwarfs are important because they dominate the galaxy number density, across all environments and at all epochs (Wright et al. 2017; Martin et al. 2019; Jackson et al. 2021b). Given their low stellar masses and high dark matter (DM) content (Di Cintio et al. 2017; Chan et al. 2018; Jackson et al. 2021a), dwarfs also provide a unique laboratory for testing our current understanding of DM and the ΛCDM structure formation paradigm (Davis et al. 2022). However, dwarf galaxies remain difficult to study, both observationally and theoretically, due to the shallow detection limits of recent wide-area surveys such as the SDSS (Alam et al. 2015) and the relatively low resolution of large volume simulations (e.g. Horizon-AGN Dubois et al. 2014b; Kaviraj et al. 2017, EAGLE Schaye et al. 2015, Illustris Vogelsberger et al. 2014, APOSTLE Sawala et al. 2016, Simba Davé et al. 2019, IllustrisTNG Nelson et al. 2019). It is, therefore, unsurprising that many tensions between observations and ΛCDM predictions are found in the dwarf regime (e.g. Bullock & Boylan-Kolchin 2017; Sales, Wetzel & Fattahi 2022).

One key tension is the ‘core-cusp’ problem (Flores & Primack 1994; Moore 1994), which describes the observed discrepancy between the measured dark matter (DM) density slopes in cosmological simulations and those inferred from the rotation curves of local dwarf galaxies (de Blok 2010). Simulations of cold DM predict a cuspy DM density profile, where density increases towards the centre of the halo (ρ ∝ r⁻¹, e.g. White & Rees 1978; Springel et al. 2008). This is commonly described by a Navarro-Frenk-White (NFW) profile (Navarro, Frenk & White 1996b, 1997), although other fits to this profile have also been suggested (e.g. Einasto 1965). In contrast, observations of dwarf galaxies often find slowly rising rotation curves which would indicate that the profiles instead have a much shallower (or even flat) density slope (i.e. a core) in their inner regions with ρ ∝ r^∼0 (e.g. Burkert 1995; de Blok et al. 2001; Gilmore et al. 2007; Kuzio de Naray, McGaugh & de Blok 2008; de Blok et al. 2008; Kormendy et al. 2009; Oh et al. 2011, 2015; Lelli, McGaugh & Schombert 2016).

This discrepancy has three main suggested solutions. First, it is possible that the discrepancy is driven by systematic errors (or unknown uncertainties) in the measurement of dwarf galaxy rotation curves. Given the difficulty in observing dwarf galaxies and extracting accurate rotation curves, this is a plausible solution. Many sources of error (such as the failure to account for non-circular motions) have been discussed and studied using simulations (Strigari, Frenk & White 2017; Genina et al. 2018; Harvey et al. 2018; Oman et al. 2019; Santos-Santos et al. 2020; Downing & Oman 2023; Roper et al. 2023). Another solution involves modifications to the ΛCDM model, such as invoking warm DM (Dodelson & Widrow 1994; Bode, Ostriker & Turok 2001) or self-interacting DM (SIDM, Spergel & Steinhardt 2000; Yoshida et al. 2000; Vogelsberger, Zavala & Loeb 2012; Rocha et al. 2013; Tulin & Yu 2018). SIDM has, in particular, been recently proposed as a viable alternative which solves the core-cusp problem (Burger & Zavala 2019; Burger et al. 2022). Finally, there is a potential solution within ΛCDM itself, which relies on the effect that baryons and baryonic processes may themselves have on DM density profiles.

Baryonic feedback processes that operate within galaxies, driven by active galactic nuclei (AGN) and/or supernovae (SN), are capable of redistributing mass within the DM Haloes, which can then alter the shape of the DM density profile via noise driven orbital shuffling (Weinberg 2001). Using the Horizon-AGN simulation, Peirani et al. (2017) have shown that AGN can affect the shape of a halo’s DM density profile in massive galaxies. They show that DM density profiles in a Horizon-AGN run that includes baryons are initially steeper than their DM-only counterpart. However, by z ∼ 1.5 AGN activity is able to flatten these profiles. After AGN activity ceases the profile then returns to being a cusp. Other work has also shown that AGN can create cores in massive galaxies at high redshift (Dekel et al. 2021; Li et al. 2023). More commonly in dwarf galaxies (where many cores are found in observations), supernova (SN) feedback is known to similarly be able to remove gas from the central regions of galaxies (Dubois & Teyssier 2008; Di Cintio et al. 2017; Chan et al. 2018; Jackson et al. 2021b).

Many theoretical studies have shown that SN feedback can create cores via violent or repeated bursts of star formation (Dekel & Silk 1986; Navarro, Eke & Frenk 1996a; Gnedin & Zhao 2002; Read & Gilmore 2005; Mashchenko, Wadsley & Couchman 2008; Pontzen & Governato 2012; Garrison-Kimmel et al. 2013; Teyssier et al. 2013; Di Cintio et al. 2014; Chan et al. 2015; Dutton et al. 2016; Tollet et al. 2016; Fitts et al. 2017; Freundlich et al. 2020; Lazar et al. 2020; Burger & Zavala 2021). One method by which this can occur, and which is invoked to explain cores in several simulations, is via the adiabatic change of the DM potential. Briefly, the gravitational potential of the DM halo is perturbed by SN-induced fluctuations which stirs the gas. DM particles on average gain energy and migrate to larger orbits due to these stochastic perturbations which resonate with their orbital frequencies, leading to lower central densities and a cored DM density profile.

Since the early attempts to model galaxies and DM Haloes with N-body simulations (e.g. Somerville & Primack 1999; Cole et al. 2000; Benson et al. 2003; Bower et al. 2006; Croton et al. 2006), computational power has increased exponentially. We are now able to run large, hydrodynamical simulations in cosmological volumes, which can model baryons and DM self-consistently, allowing us to test the interplay between them and quantify the influence of baryonic processes on the DM. Different cosmological simulations, that use a variety of underlying codes and subgrid models, have subsequently attempted to address the core-cusp problem. Those which adopt a zoom-in technique (resulting in generally higher resolution but lacking a statistical sample) have had success in both recreating the broad properties of galaxies across a range of stellar masses and creating DM density profiles with both cores and cusps. These simulations (e.g. Di Cintio et al. 2014; Chan et al. 2015; Tollet et al. 2016; Fitts et al. 2017; Lazar et al. 2020) consistently find that SN feedback is capable of driving gas outflows that alter the gravitational potential in such a way as to form DM cores in dwarf galaxies. They have also shown that cores form more efficiently in a relatively narrow range of stellar or halo mass and that the inner slope of dark matter Haloes may be halo mass dependent.

Conversely, large-volume simulations such as Illustris (Vogelsberger et al. 2014), EAGLE (Schaye et al. 2015) and APOSTLE (Sawala et al. 2016), which reproduce the properties and scaling relations that govern massive galaxies well, have been unable to reproduce cored DM density profiles in their fiducial runs, regardless of the stellar or halo mass of the galaxies in question (Schaller et al. 2015; Benítez-Llambay et al. 2019; Chua et al. 2019). However, re-simulations using the EAGLE model have managed to produce a diversity of DM density profiles when adjusting the gas density threshold required for star formation, albeit for smaller samples of galaxies (Benítez-Llambay et al. 2019).

This inability to form cores in large volume simulations is attributed to two causes: the star formation (SF) density threshold of the gas and the modelling of the inter-stellar medium (ISM). These parameters can often be linked, with a lack of spatial resolution in simulations leading to approximations in the modelling of the ISM, and a low SF density threshold leading to a lack of cores (Oman et al. 2015; Schaller et al. 2015; Bose et al. 2019). Indeed, many subsequent studies have shown that efficient core formation relies on the presence of dense gas (and therefore a high SF threshold), as well as bursty SF, in order to perturb the gravitational potential sufficiently to move DM particles to larger radii (Governato et al. 2010; Teyssier et al. 2013; Di Cintio et al. 2014; Benítez-Llambay et al. 2019; Jahn et al. 2021). Therefore, a simulation which is both capable of resolving the small-scale processes within galaxies and offers a cosmological volume containing a statistically large sample of dwarf galaxies is needed to truly address the core-cusp problem.

In this study we use the NewHorizon (NH) hydrodynamical cosmological simulation (Dubois et al. 2021), to study the DM density profiles of galaxies across all stellar masses and in a range of environments (from the field to large groups). NH is a zoom-in of an average density region taken from the larger Horizon-AGN simulation (Dubois et al. 2014b; Kaviraj et al. 2017), which has a volume of (142 cMpc)³. NH offers a maximum spatial resolution of 34 pc and mass resolutions of ∼10⁴ M_⊙ and ∼10⁶ M_⊙ in the stars and DM, respectively. In addition, it has a timestep resolution of 15 Myr. This makes it an ideal tool to study the inner DM density profiles of galaxies, particularly in the dwarf regime where cores are believed to be present (e.g. Martin et al. 2021; Jackson et al. 2021a, b), whilst also offering a statistically significant sample of galaxies due to its large volume.

This paper is structured as follows. In Section 2 we briefly describe the NH simulation and in Section 3 we show the broad properties of the DM density profiles in NH galaxies. In Section 4, we study how cores and cusps form in the simulation. We summarize our findings in Section 5.

2 SIMULATION

The NH cosmological, hydrodynamical simulation (Dubois et al. 2021), is a high-resolution simulation, produced using a zoom-in of a region of the Horizon-AGN simulation. Horizon-AGN employs the adaptive mesh refinement code ramses (Teyssier 2002) and utilizes a grid that simulates a (142 cMpc)³ volume (achieving a spatial resolution of 1 kpc), using 1024³ uniformly distributed cubic cells that have a constant mass resolution, via the mpgrafic (Prunet et al. 2008) code.

For NH, we resample this grid at higher resolution (using 4096³ uniformly distributed cubic cells), with the same cosmology as that used in Horizon-AGN (Ω_m = 0.272, Ω_b = 0.0455, Ω_Λ = 0.728, H₀ = 70.4 km s⁻¹ Mpc⁻¹, and n_s = 0.967; Komatsu et al. (2011)). The high-resolution volume occupied by NH is a sphere which has a radius of 10 comoving Mpc, centred on a region of average density within Horizon-AGN. NH has a DM mass resolution of 10⁶ M_⊙ (compared to 8 × 10⁷ M_⊙ in Horizon-AGN), stellar mass resolution of 10⁴ M_⊙ (compared to 2 × 10⁶ M_⊙ in Horizon-AGN) and a maximum spatial resolution of 34 pc (compared to 1 kpc in Horizon-AGN).¹ The simulation has been performed down to z ∼ 0.18.

2.1 Star formation and stellar feedback

NH employs gas cooling via primordial Hydrogen and Helium, which is gradually enriched by metals produced by stellar evolution (Sutherland & Dopita 1993; Rosen & Bregman 1995) allowing for metal cooling. An ambient UV background is switched on after the universe is re-ionized at z = 10 (Haardt & Madau 1996). Star formation is assumed to take place in gas that has a hydrogen number density greater than n_H > 10 H cm⁻³ and a temperature lower than 2 × 10⁴ K, following a Schmidt relation (Schmidt 1959). The star formation efficiency is dependent on the local turbulent Mach number and the virial parameter α_vir = 2E_k/|E_g|, where E_k is the turbulent kinetic energy of the gas and E_g is the gas gravitational binding energy (Kimm et al. 2017). The probability of forming a star particle of mass M_{*, res} = 10⁴ M_⊙ is drawn at each time step using a Poissonian sampling method, as described in Rasera & Teyssier (2006).

Each star particle represents a coeval stellar population, with 31 per cent of the stellar mass of this population (corresponding to stars more massive than 6 M_⊙) assumed to explode as Type II supernovae, 5 Myr after its birth. This fraction is calculated using a Chabrier initial mass function (Chabrier 2005), with upper and lower mass limits of 150 M_⊙ and 0.1 M_⊙. Supernova feedback takes the form of both energy and momentum, with the final radial momentum capturing the snowplough phase of the expansion (Kimm & Cen 2014). The initial energy of each supernova is 10⁵¹ erg and the supernova has a progenitor mass of 10 M_⊙. Pre-heating of the ambient gas by ultraviolet radiation from young O and B stars is included by augmenting the final radial momentum from supernovae following Geen et al. (2015).

2.2 Supermassive black holes and black hole feedback

Supermassive black holes (SMBHs) are modelled as sink particles. These accrete gas and impart feedback to their ambient medium via a fraction of the rest-mass energy of the accreted material. SMBHs are allowed to form in regions that have gas densities larger than the threshold of star formation, with a seed mass of 10⁴ M_⊙. New SMBHs do not form at distances less than 50 kpc from other existing black holes. A dynamical gas drag force is applied to the SMBHs (Ostriker 1999) and two SMBHs are allowed to merge if the distance between them is smaller than 4 times the cell size, and if the kinetic energy of the binary is less than its binding energy.

The Bondi–Hoyle–Lyttleton accretion rate determines the accretion rate on to SMBHs, with a value capped at Eddington (Hoyle & Lyttleton 1939; Bondi & Hoyle 1944). The SMBHs release energy back into the ambient gas via a thermal quasar mode and a jet radio mode, when accretion rates are above and below 1 per cent of the Eddington rate, respectively (Dubois et al. 2012). The spins of these SMBHs are evolved self-consistently through gas accretion in the quasar mode and coalescence of black hole binaries (Dubois, Volonteri & Silk 2014a), which modifies the radiative efficiencies of the accretion flow (following the models of thin Shakura & Sunyaev accretion discs) and the corresponding Eddington accretion rate, mass–energy conversion, and bolometric luminosity of the quasar mode (Shakura & Sunyaev 1973). The quasar mode imparts 15 per cent of the bolometric luminosity as thermal energy into the surrounding gas. The radio mode employs a spin-dependent variable efficiency and spin up and spin down rates that follow the simulations of magnetically choked accretion discs (see e.g. McKinney, Tchekhovskoy & Blandford 2012).

2.3 Identification of galaxies and merger trees

We use the adaptahop algorithm to identify DM Haloes (Aubert, Pichon & Colombi 2004; Tweed et al. 2009). adaptahop efficiently removes subHaloes from principal structures and keeps track of the fractional number of low-resolution DM particles within the virial radius of the halo in question. We identify galaxies in a similar fashion, using the HOP structure finder applied directly to star particles (Eisenstein & Hut 1998). The difference between HOP and adaptahop lies in the fact that HOP does not remove substructures from the main structure, since this can result in star-forming clumps being removed from galaxies. We then produce merger trees for each galaxy in the final snapshot at z ∼ 0.18, with an average timestep of ∼15 Myr, enabling us to track the main progenitor of every galaxy with high temporal resolution.

Given that NH is a high-resolution zoom of Horizon-AGN, it is worth considering the DM purity of galaxies, since higher-mass DM particles may enter the high resolution region of NH from the surrounding lower-resolution regions. Given the large mass difference, these DM particles may interact with galaxies they are passing through in unusual ways. The vast majority of galaxies affected by low DM purity exist at the outer edge of the NH sphere. The galaxies studied in this paper all have DM Haloes with a purity of 100 per cent at the final snapshot.

3 CORES AND CUSPS IN THE GALAXY POPULATION

It has been shown in previous studies that modern simulations are able to produce DM density profiles that correspond to both cores and cusps when there are baryons present (Di Cintio et al. 2014; Chan et al. 2015; Dutton et al. 2016; Tollet et al. 2016; Fitts et al. 2017; Benítez-Llambay et al. 2019; Lazar et al. 2020). In order to investigate this using NH we first match Haloes between a DM-only and full (i.e. DM + baryons) versions of the simulation. This gives us a sample of DM Haloes with and without baryons for comparison. The matching is confirmed by ensuring that it passes a threshold of sharing a significant number (75 per cent) of the same DM particles. In addition, if the total mass ratio of the Haloes in the two runs is outside of the range 0.1 < M_NH/M_NHDM-only < 10) then it is excluded from our analysis (see Peirani et al. (2017) for more details). With this sample of Haloes in hand, we begin our analysis by creating DM halo density profiles for all galaxies in both the full NH run and its DM-only counterpart. The volume of NH is large enough that we are able to investigate the DM Haloes of ∼1000 galaxies in the stellar mass range 10⁷ M_⊙ < M_⋆ < 10^11.5 M_⊙, thus ensuring that we have a statistically significant sample.

We start by showing, in fig. 1, how the M_Halo versus M_⋆ and M_⋆ versus R_eff relations in NH compare to those of comparable simulations (APOSTLE, EAGLE CHT10, and NIHAO) and observational results (Behroozi, Wechsler & Conroy 2013; Moster, Naab & White 2013). NH Haloes of a given mass contain more stellar mass than other simulations (though similar to APOSTLE at lower masses) and also observations although at higher stellar masses are consistent with observations. However the M_⋆ versus R_eff relation is well matched to the observation data of Mowla et al. (2019).

In Fig. 2, we show two examples of DM Haloes in the two versions (DM-only and full, i.e. DM + baryons) of NH, which exhibit differing shapes in their central regions. The left panel shows a massive halo (M_halo = 10^12.2M_⊙) which, in the NH full run, contains a galaxy with M_⋆ = 10^10.9 M_⊙. The DM density profile in both runs shows a clear cusp-like shape with no flattening in the central region. For Haloes (and galaxies) at such masses the presence of baryons increases the central density of the DM profile, due to the extra baryonic mass creating a deeper gravitational potential.

The right panel shows a low mass halo (M_halo = 10^10.31M_⊙), which contains a dwarf galaxy with M_⋆ = 10^8.15 M_⊙ in the NH full run. While the DM-only run shows a cuspy profile, the presence of baryons in this example flattens the central regions in the NH full run, as DM is redistributed within the halo from central regions to larger radii. It therefore seems that, in NH, baryons are able to influence the DM density profiles of Haloes in multiple ways, depending on the halo/galaxy mass in question.

3.1 Parametrizing DM density profile shapes

In order to study how baryonic processes affect the inner DM density profiles of our simulated galaxies, we first need to parametrize the profile shapes in these inner regions. To do this we construct the DM density profile by computing the mean density profiles (binned in spherical shells equally spaced in log r) and then fit a power law (Cr^α, where C is a constant) to obtain the slope within a given radial range in log space. In order to compare the results of this study to other work, and to ensure that we probe the same region of each halo regardless of mass, we measure the slope of the profile between 1 and 2 per cent of a halo’s virial radius. We ensure that this is above the DM spatial resolution limit (∼500 pc) and remove any galaxies that do not meet this requirement. It is also important to ensure that our selection of radial range does not influence the overall results. To account for this, we also measure the slope of the DM density profile between 1 and 2 kpc which is usually larger than 2 per cent of the virial radius. In Appendix  A, we show the different slopes obtained for all galaxies using these different radial ranges. As we find that the overall trends do not change with the radial ranges used (even if there is some variation on an individual galaxy basis) we proceed by using the fractional virial radius range for our analysis (i.e. slopes measured within 1 and 2 per cent of a halo’s virial radius).

In Fig. 3 we show the inner DM density profile slope (α), between 1–2 per cent of R_vir versus the halo mass for all galaxies in the full NH run (colour-coded by their stellar mass) and their corresponding halo in the DM-only run (orange points) at the last timestep of the simulation (which corresponds to z ∼ 0.18). The inner slopes of Haloes show a clear relationship with halo mass for NH galaxies. In the lowest and highest mass Haloes we find cusp-like profiles (defined in this study as α < −0.65), whereas in intermediate mass Haloes (with DM masses between 10^10.2 M_⊙ and 10^11.5 M_⊙) the inner DM profiles are flatter and core-like (α > −0.65). We pick this cut-off as no DM-only Haloes have profiles α > −0.65 as well as being roughly consistent with previous studies. At both extremes of halo mass, the Haloes show evidence of being even denser in their centres compared to their DM-only counterparts, with profiles that are steeper than those predicted by an NFW profile. Haloes in the DM-only run are all roughly consistent with the predicted values for an NFW-profile and so are cuspy in nature. Therefore, the presence of baryons is able to cause a halo’s centre to be both more or less dense than an NFW-profile, depending on the mass of the halo and the baryonic processes that operate in the central regions.

3.2 α as a function of stellar and halo mass

Fig. 4 shows how the measured α for each halo varies with stellar mass (left panel), halo mass (centre panel), and M_⋆/M_halo (right panel). The red points indicate the values for galaxies in the full NH run in each panel. In the left panel we compare to the theoretical results of Tollet et al. (2016) and observational results of Oh et al. (2015). We find reasonable agreement with Tollet et al. (2016) in higher mass galaxies (down to ∼10⁹ M_⊙), albeit with a larger scatter for NH galaxies. The lowest mass galaxies in NH show significantly cuspier profiles, although the broad transition region where cores are no longer found is consistent. Oh et al. (2015) find cores in a similar stellar mass range, although many cores are at lower stellar mass in their study. It is worth noting, however, that their sample only contains a small number of galaxies.

When we consider the central panel, we also find agreement between the NH distribution and Tollet et al. (2016) at high and intermediate halo masses. However, below ∼10^10.2 M_⊙ we find Haloes with both cores and cusps in NH unlike in Tollet et al. (2016). In the right panel we consider M_⋆/M_halo, which in many studies is believed to be the key determining factor for core formation. Here we include the simulation results from Tollet et al. (2016), Di Cintio et al. (2014), and Lazar et al. (2020) in our comparison. The simulation studies (Di Cintio et al. 2014; Tollet et al. 2016; Lazar et al. 2020) all measure the DM density slope within the same radius, although we note that in Oh et al. (2015) the slope is calculated within a break radius and may therefore be different (see Oh et al. (2011) for more details). Whilst these three studies find a similar relationship between α and M_⋆/M_halo, we find that NH galaxies show a different trend. In NH, cores are found across a range of M_⋆/M_halo values. Unlike many other studies, which find a dependence on M_⋆/M_halo and which is deemed to be a proxy of star formation efficiency (Di Cintio et al. 2014; Tollet et al. 2016; Dutton et al. 2020), there is no clear region of the M_⋆/M_halo relation where galaxies in NH only contain a core. Instead, it appears that once the halo crosses a threshold mass (∼10^10.2 M_⊙) and enough stars have formed (∼10⁸ M_⊙), there is a possibility a core will be created. The idea of a ‘threshold halo mass’ for core formation to take place has been previously suggested in Chan et al. (2015) but is demonstrated here, using NH, using a much larger galaxy sample.

As there is no clear trend with M_⋆/M_halo in NH, unlike previous studies, it would imply that this metric should not be used solely to determine whether a galaxy should contain a core. There is, however, a range of halo and stellar masses where most galaxies have a cored profile. Interestingly, low-mass galaxies in NH may overproduce stars compared to observed abundance matching relations (Dubois et al. ), explaining the discrepancy in the α–M_⋆/M_halo relationship. This gives credence to the idea that core formation is not solely dictated by M_⋆/M_halo, as NH creates cores at similar stellar and halo masses to other simulations.

4 THE FORMATION OF CORES ACROSS COSMIC TIME

A key question in the formation of cores is whether the halo initially forms as a core or whether a core develops due to some physical processes. In this section, we investigate the prevalence of cored DM profiles across cosmic time. Fig. 5 shows the percentage of galaxies with core-like profiles (recall that this is defined as α > −0.65) as a function of look-back time. We see that the frequency of cores has been relatively constant for the last ∼9 Gyrs, suggesting that either, once cores have formed, they are resilient and hard to destroy or cores are constantly being formed to replenish those that are destroyed. At the last timestep we find that |${\sim}70~{{\ \rm per\ cent}}$| of cores have existed for at least 3 Gyrs, implying that most cores are relatively long lived. Prior to this period, however, there were few to no cores. Indeed, it is during the period where cosmic star formation starts to increase and then peaks (Madau & Dickinson 2014) that the number of cores begins to increase rapidly. This not only allows time for Haloes to grow to a mass capable of maintaining a core but also suggests that star formation plays a key role in core formation.

As we now see that cores have existed across most of cosmic time, we can further assess the concept of a ‘threshold halo mass’, by looking at the mass of Haloes which contain cores throughout the simulation. In Fig. 6, we show the median halo mass of galaxies with cores throughout the simulation (with the edges of the shaded region indicating the 16th and 84th percentile values). The median halo mass of cores is relatively constant across time at ∼10^10.3 M_⊙. This implies that the process that creates the flattened density profile requires some mass limit to be reached and also suggests that there is an upper limit beyond which the halo centres will return to NFW-like profiles again, as cores are not found at all halo masses. This, combined with the results from Fig. 4, suggests that cores will form once Haloes are both above a given mass and when star formation is occurring within the galaxy. Given that this halo mass is roughly consistent with the other simulations we have compared our results to (Di Cintio et al. 2014; Tollet et al. 2016; Lazar et al. 2020), it appears that disparate baryonic prescriptions can give rise to core formation within different simulations.

4.1 How and why do cores form?

In order to study how and why cored profiles form, we consider a region of the stellar and halo mass parameter space where Haloes have both cored and cuspy inner profiles and examine the physical differences between galaxies in this regime. We proceed by selecting galaxies with stellar masses in the range 8 < log(M_⋆/M_⊙) < 9 in the final snapshot (z = 0.18) as this contains a large number of both cores and cusps, whilst also retaining enough resolution in the simulation to study the evolution of the galaxies.

We start by looking at external processes that could create cored profiles. Tidal interactions between galaxies have been shown to strongly affect the DM distribution of galaxies, albeit usually at the outskirts (Jackson et al. 2021a). Fig. 7 shows halo mass versus stellar mass for central galaxies (left) and satellites (right) defined by the adaptahop structure finder, colour-coded by their values of α. Both central and satellite galaxies contain a variety of core and cusp profiles, with cores in satellites being found at slightly lower halo masses (likely due to tidal stripping of the outskirts reducing their original halo mass). This suggests that becoming a satellite does not initially alter the central DM distributions of galaxies and therefore that tidal processes are not a likely cause of core creation, at least not directly. This also indicates that it is possible to find galaxies with cores as satellites of larger galaxies, such as the those in the local group, as has been previously suggested by Walker & Peñarrubia (2011).

As external processes are unlikely to be able to alter the central part of DM profiles (at least until the outskirts have been stripped), we proceed by investigating the internal physical properties of galaxies in this transition region. In Table 1, we show the median values of stellar, DM, total gas and star-forming gas mass, measured both within the central 2 kpc and within the virial radius, for all cored and cuspy galaxies in our transition region. The stellar, DM, and total gas mass are consistent, regardless of the shape of the DM density profile, with gas and DM being the dominant contributors to the central mass budget in both cases. However, in the cored profiles the amount of star-forming gas (n_H >10 H cm⁻³) is significantly elevated in the centre, suggesting that there is (or has recently been) increased star formation in these galaxies compared to their cuspy counterparts. Feedback processes from star formation are known to be able to drive outflows of gas which, if contributing a significant amount of mass, can alter a galaxy’s gravitational potential and induce the DM to move to larger radii (Navarro et al. 1996a; Read & Gilmore 2005; Mashchenko, Couchman & Wadsley 2006; Martizzi, Teyssier & Moore 2013; Pontzen & Governato 2014).

Therefore, the relative contributions to the total mass in the central regions of the halo is key to determining the shape of the profile. In order to understand how the mass budget impacts the DM profile, we compare the central mass evolution of two example core and cusp galaxies. Both galaxies are centrals and, by the end of the simulation, are hosted in a similarly sized Haloes (log M_Halo ∼ 10.3). Fig. 8 shows the enclosed mass within two radial distances (1 and 5 kpc) as a function of look-back time for a galaxy with a cuspy halo (left panel) and a cored halo (right panel). These galaxies represent typical examples of cuspy and cored Haloes in this transition regime. The blue curves show the stellar mass evolution, the green curves represent total gas and the orange curves show DM, with solid lines indicating the mass within 1 kpc and dashed lines showing the mass within 5 kpc. In the bottom panels, the red lines show the evolution of the specific star formation rate (sSFR) and the dashed black line shows the evolution of the inner DM density slope of the galaxies in question.

At large look-back times (i.e. in the early universe) both galaxies show highly elevated sSFRs, which coincides with a rapid increase in stellar mass in the central regions. At the same epoch, gas and DM constitute roughly equal amounts of mass within 1 kpc and are dominant over the stellar component. In the case of the galaxy that ends up with a cuspy halo, the gas component drops rapidly and star formation quenches, leaving the DM component at a constant level and the dominant constituent of the central mass budget.

However, in the case of the galaxy that ends up with a cored halo, we see repeated bursts of star formation, with the gas mass within 1 kpc repeatedly rising and falling, with the falls occurring shortly after the peak sSFR values. When looking at all cores in NH, we find that no galaxy has a core when bursty star formation lasts for less than 2 Gyrs. The gas and DM components within 1 kpc are similar in mass initially but, unlike in the cuspy case above, the gas mass remains high and even becomes the dominant central mass component at 8 < t < 10 Gyr. Despite star formation taking place, the amount of stellar mass being formed is less than the oscillations in gas mass, and thus the removal of gas from the centre is a combination of stars forming and gas being moved to larger radii.

As the gas mass is comparable to the DM mass, moving large amounts of gas to larger radii will alter the gravitational potential, allowing the DM to migrate to larger radii as well. For example, in the right panel of Fig. 8, as gas is removed from the centre (at around 11 Gyrs), the DM mass within 1 kpc also decreases. Although gas in-falls to the centre again and the mass increases the DM mass does not and remains relatively constant. Gas is then once again removed leading to a further decrease in the DM. In this way the DM mass is removed in stages by bursts of star formation and gas removal, which flattens the central DM density profile and creates a core. We can also calculate the dynamical time for the halo during this period, within 1 kpc, to understand how this compares to the oscillations of gas. In the example of the core in Fig. 8, we find that when the core is forming, the dynamical time is 0.2 Gyrs, which is roughly equal to the duration of the observed bursts.

In order for a core to be formed in this manner, the gravitational potential of the galaxy needs to be altered from a classical NFW profile. As previously shown, large amounts of gas could be moved from the very central regions to larger radii, thus altering the shape of the potential. In Fig. 9, we show the gravitational potential of the example cored galaxy as a function of its radius. In order to compare the shape of the potential at different times, as well as to its DM-only counterpart, we normalize the potential by its lowest value (Φ₀) and the halo’s radius by its maximum value (r_max). We show both the DM component of the potential (dashed lines) and the gas component (solid lines) for different look-back times that correspond to moments before and after bursts of star formation which create the core and the final timestep of the simulation (black). Between the first two timesteps shown, the potential of both the DM and gas is significantly altered due to the removal of gas (and subsequently DM) from the centre of the halo. After this period the potential starts to return to its previous shape until another burst of SF further alters it towards a core. Indeed, even by the end of the simulation the potential has not returned to its original shape, before the bursts of star formation occurred.

4.2 The evolution of massive galaxy profiles – the reformation of cusps

As seen in Fig. 4, above a certain halo mass most galaxies in the simulation have cuspy DM density profiles. However, if core formation is determined purely by a threshold halo mass and star formation, have these galaxies had cores in their past and, if so, what removed them? In Fig. 10, we show an example of a massive galaxy, containing a cusp at the final timestep in NH, with a stellar and halo mass of 10^11.3 M_⊙ and 10^12.9 M_⊙, respectively. In the top panel we show the stellar (blue), gas (orange), and DM (green) mass enclosed within the inner 1 kpc (and the stellar mass within the inner 5 kpc as a blue dashed line) as a function of look-back time. For the lower panel we show the sSFR (red line) of the galaxy as a proxy for SN feedback and the DM density slope within 1−2 per cent of R_vir (black dashed line).

This massive galaxy experiences a period when it has a cored DM density profile at high redshift, coinciding with high sSFR, and comparable gas/DM/stellar masses in the central 1 kpc. However, unlike in Fig. 8, the DM density profile steepens again as the stellar mass continues to increase in both the inner 1 and 5 kpc, despite ongoing star formation and high gas mass. The stellar mass increases in both radial bins, which suggests that it is increased central in situ star formation (rather than existing mass falling into the centre or ex-situ mass accretion) that drives the reformation of the cusp. Unlike in the cored example (Fig. 8), sustained star formation leads to the stellar component becoming dominant over both gas and DM in the central regions. This increased stellar component deepens and anchors the potential, reducing the outflow of gas and thus stopping the dynamical heating of DM particles. The DM reacts to this new, deeper potential removing the core and recreating the cusp. Therefore the halo’s inability to control star formation is at least partly responsible for the rejuvenation of the cusp.

In addition, just like in the cored example, we calculate the dynamical time for the halo within 1 kpc. We see a similar dynamical time (∼0.2 Gyrs) whilst this halo has a core. However, as the star formation continues the dynamical time rapidly decreases to 0.04 Gyrs meaning that the cusp can rapidly reform, as is observed in the slope of the DM density profile.

The potential role of BHs (and therefore AGN) is also worth considering, especially as it has been previously shown to affect DM density profiles in larger galaxies (Peirani et al. 2017). However, for galaxies in the mass ranges where cores are found in NH the BH occupation fraction is |${\sim}30~{{\ \rm per\ cent}}$|⁠. In addition, in those galaxies which contain BHs both cuspy and cored Haloes are found. Therefore it is logical to conclude that, at least in this mass range, BHs and AGN are not primarily responsible for the creation of cores or the reformation of cusps.

5 SUMMARY AND CONCLUSIONS

Dwarf galaxies dominate the galaxy number density in all environments across cosmic time. However, due to their low stellar masses, and consequently fainter surface brightnesses, they have remained difficult to study in past surveys. This, in part, has contributed to many tensions between observations and theoretical predictions within ΛCDM in the dwarf regime. One longstanding tension is the core-cusp problem. DM-only simulations suggest that DM Haloes should have a cuspy profile, with the density continuing to increase towards the centre. However, many observational studies of dwarf galaxies indicate that they actually possess flattened density profiles in their central regions, known as cores.

Here, we have used the NH cosmological simulation, in conjunction with its DM-only counterpart, to examine the evolution of the DM density profiles of galaxies across a large range of stellar masses and environments (from the field to large groups). This has enabled us to study the diversity of DM density profiles in both dwarf and massive galaxies in a statistical fashion, for the first time within a cosmological volume. Our main conclusions are as follows:

• The central DM density slope (α) between 1−2 per cent of R_vir can be used to probe whether a halo contains a core or a cusp. In our study we consider galaxies with α > −0.65 as having cores. While Haloes in NH exhibit DM density profiles with both cores and cusps, only cusps are present in Haloes in the DM-only run. Thus, it is the presence of baryons that drives the observed diversity of DM density profiles.

• Cores in NH form when there are repeated bursts of star formation, in a galaxy that has a central gas mass that is comparable to the DM mass. No galaxy (at least in NH) creates a core when this bursty star formation lasts for a period shorter than 2 Gyrs. The resultant SN feedback removes gas from the central regions, altering the shape of the (central) gravitational potential. The DM responds to this change and moves to larger radii. This process is typically repeated multiple times, with the central DM being removed in stages, until star formation ceases. If this process occurs for a sufficient length of time (∼2–3 Gyrs), the shape of the DM density profile is eventually transformed from an initial cusp into a core.

• Conversely cusps can be reformed if the central star formation continues rapidly after the core forms. In this scenario a core forms, albeit briefly, before being removed when stellar mass becomes dominant in the most central regions. This mass comes from central in situ star formation rather than stellar mass infall from larger radii into the central regions or ex-situ accretion (such as mergers). At this stage, stars dominate the central gravitational potential and, since they are not directly affected by SN feedback, anchor the new potential in place. This, in turn, causes the DM to adiabatically cool, reforming a cuspy DM density profile.

• Haloes that contain galaxies in the upper (M_⋆ ≥ 10^10.2 M_⊙) and lower (M_⋆ ≤ 10⁸ M_⊙) ends of the stellar mass distribution contain cusps. However, in Haloes that contain galaxies that have intermediate (10⁸ M_⊙ ≤ M_⋆ ≤ 10^10.2 M_⊙) stellar masses, α becomes more positive than in the DM-only run (i.e. the profiles become core-like), with this becoming most common in Haloes with DM masses between 10^10.2 M_⊙ and 10^11.5 M_⊙.

• The evolution of α with stellar and halo mass in NH is broadly consistent with what is seen in other comparable simulations (except at the lowest stellar masses) and also in limited available observational data (given the difficulties in inferring a DM density profile from observations). However, the trend of increasing core frequency with increasing M_⋆/M_Halo, that is observed in some other simulations, is not seen in NH, indicating that this may not be the best metric to determine whether a galaxy has a core or a cusp (as has been previously suggested).

• The fraction of galaxies that contain a core becomes roughly constant (∼15 per cent) after a redshift of z ∼ 1.5, with cores not forming until the cosmic star formation rate increases around a redshift of z ∼ 2. The mass of Haloes that contain a core is roughly constant in the mass range 10^10.2 M_⊙ < M_⋆ < 10¹¹ M_⊙, with a median value of 10^10.3 M_⊙, which indicates (1) a potential threshold halo mass below which cores cannot form and (2) a maximum halo mass above which cores are removed. Additionally at the final timestep 70 per cent of cores have existed for at least 3 Gyrs, suggesting that most cores are relatively long-lived.

We conclude that, in NH, a cosmological simulation that adopts a ΛCDM cosmology, the core-cusp problem is not a challenge to the standard paradigm. While this has been shown in previous theoretical work, we are able to perform this test with a statistically significant number of galaxies, over a larger range of galaxy masses, and within a cosmological volume. The limited amount of observational data makes comparisons to the real universe difficult, though this may be addressed in future deep-wide surveys, which will allow us to constrain the stellar masses where cores form in real galaxies. As NH has been shown to successfully recreate cores and cusps, it is a unique tool to investigate other ‘small-scale’ challenges to ΛCDM, such as the diversity of rotation curves in dwarf galaxies, which will be explored in a forthcoming paper (Jackson et al. in prep).

ACKNOWLEDGEMENTS

SK acknowledges support from the STFC [grant numbers ST/S00615X/1 and ST/X001318/1]. SK also acknowledges a Senior Research Fellowship from Worcester College Oxford. SKY acknowledges support from the Korean National Research Foundation (NRF-2020R1A2C3003769, NRF-2022R1A6A1A03053472). SKY acted as a corresponding author as the head of the group whilst most of this investigation was being performed at Yonsei University. This research has made use of the Horizon cluster on which the simulation was post-processed, hosted by the Institut d’Astrophysique de Paris. We warmly thank S. Rouberol for running it smoothly. This work is partially supported by the grant Segal ANR-19-CE31-0017 of the French Agence Nationale de la Recherche and by the National Science Foundation under Grant Number NSF PHY-1748958. This work was granted access to the HPC resources of CINES under the allocations c2016047637, A0020407637, and A0070402192 by Genci and as a ‘Grand Challenge’ project granted by GENCI on the AMD Rome extension of the Joliot Curie supercomputer at TGCC. Part of NH was performed using Nurion at KISTI under the allocation of KSC-2017-G2-0003, and the large data transfer was supported by KREONET which is managed and operated by KISTI. SKY acknowledges support from the Korean National Research Foundation (NRF-2020R1A2C3003769, NRF-2022R1A6A1A03053472). This study was funded in part by the NRF-2022R1A6A1A03053472 grant. TK was supported by the National Research Foundation of Korea (Nos. 2020R1C1C1007079 and 2022R1A6A1A03053472)

For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission.

6 DATA AVAILABILITY

All data used in this paper can be made available upon reasonable request to the first author.

Footnotes

References

APPENDIX A: RADIUS RANGE COMPARISON

In Fig. A1, we compare the slope of the inner DM density profiles using three different radial ranges as a function of stellar mass. In blue we show profiles calculated in the range 0.01 < R_vir < 0.02. The orange points show the same for 1 < r/kpc < 2, while the green points show profiles calculated in the range 1 < r/r_soft < 5, where r_soft is ∼500 pc. For most galaxies, the radius at which the inner DM density profile gradient is measured does not alter the results found. The overall trend remains the same, with cores formed in galaxies between 10⁸ M_⊙ < M_⋆ < 10^10.2 M_⊙ and cusps formed in the very low and very high mass galaxies. This also suggests that the cores produced are larger than these radii, indicating that many cores in NH are at least 2 kpc in size.