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

R Jackson, S Kaviraj, S Yi, S Peirani, Y Dubois, G Martin, J Devriendt, A
Slyz, C Pichon, M Volonteri, et al.

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

R Jackson, S Kaviraj, S Yi, S Peirani, Y Dubois, et al.. The formation of cores in galaxies across cosmic time – the existence of cores is not in tension with the ΛCDM paradigm. Monthly Notices of the Royal Astronomical Society, 2024, 528 (2), pp.1655-1667. 10.1093/mnras/stae056. insu-04423473

HAL Id: insu-04423473

https://insu.hal.science/insu-04423473

Submitted on 29 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# The formation of cores in galaxies across cosmic time – the existence of cores is not in tension with the ΛCDM paradigm

R. A. Jackson<sup>1,2,3★</sup> S. Kaviraj<sup>1</sup>,<sup>2</sup> S. K. Yi,<sup>3</sup> S. Peirani,<sup>4,5</sup> Y. Dubois,<sup>5</sup> G. Martin<sup>6,7,8</sup> J. E. G. Devriendt,<sup>9</sup> A. Slyz,<sup>9</sup> C. Pichon,<sup>5,10</sup> M. Volonteri,<sup>5</sup> T. Kimm<sup>10,3</sup> and K. Kraljic<sup>11</sup>

Accepted 2024 January 5. Received 2024 January 2; in original form 2023 September 29

## 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_{\star} \ge 10^{10.2} M_{\odot}$) and lower ($M_{\star} \le 10^{8} M_{\odot}$) ends of the stellar mass distribution contain cusps. However, Haloes containing galaxies with intermediate ($10^{8} M_{\odot} \le M_{\star} \le 10^{10.2} M_{\odot}$) stellar masses are generally cored, with typical halo masses between $10^{10.2} M_{\odot}$ and $10^{11.5} M_{\odot}$. 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_{\star} \ge 10^{9.5} M_{\odot}$) 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.

**Key words:** galaxies: formation – galaxies: evolution – galaxies: dwarf – galaxies: haloes – methods: numerical.

## 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

<footnote n="1"><sup>1</sup>Department of Physics and Astronomy, University of Victoria, BC V8P 5C2, Canada</footnote>

<footnote n="2"><sup>2</sup>Centre for Astrophysics Research, School of Physics, Astronomy and Mathematics, University of Hertfordshire, Hatfield, AL10 9AB, UK</footnote>

<footnote n="3"><sup>3</sup>Department of Astronomy and Yonsei University Observatory, Yonsei University, Seoul 03722, South Korea</footnote>

<footnote n="4"><sup>4</sup>Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Bd de l'Observatoire, Université Côte d'Azur, CS 34229, F-06304 Nice Cedex 4, France</footnote>

<footnote n="5"><sup>5</sup>Institut d'Astrophysique de Paris, Sorbonne Universités, UMPC Univ Paris 06 et CNRS, UMP 7095, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="6"><sup>6</sup>School of Physics and Astronomy, University of Nottingham, University Park, Nottingham NG7 2RD, UK</footnote>

<footnote n="7"><sup>7</sup>Korea Astronomy and Space Science Institute, 776 Daedeokdae-ro, Yuseong-gu, Daejeon 34055, South Korea</footnote>

<footnote n="8"><sup>8</sup>Steward Observatory, University of Arizona, 933 N. Cherry Avenue, Tucson, AZ 85719, USA</footnote>

<footnote n="9"><sup>9</sup>Department of Physics, University of Oxford, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="10"><sup>10</sup>School of Physics, Korea Institute for Advanced Study (KIAS), 85 Hoegiro, Dongdaemun-gu, Seoul 02455, South Korea</footnote>

<footnote n="11"><sup>11</sup>Observatoire Astronomique de Strasbourg, Université de Strasbourg, CNRS, UMR 7550, F-67000 Strasbourg, France</footnote>

<footnote>*E-mail: ryanjackson@uvic.ca</footnote>
# The formation of cores in galaxies across cosmic time – the existence of cores is not in tension with the ΛCDM paradigm

R. A. Jackson<sup>1,2,3★</sup> S. Kaviraj<sup>1</sup>,<sup>2</sup> S. K. Yi,<sup>3</sup> S. Peirani,<sup>4,5</sup> Y. Dubois,<sup>5</sup> G. Martin<sup>6,7,8</sup> J. E. G. Devriendt,<sup>9</sup> A. Slyz,<sup>9</sup> C. Pichon,<sup>5,10</sup> M. Volonteri,<sup>5</sup> T. Kimm<sup>10,3</sup> and K. Kraljic<sup>11</sup>

Accepted 2024 January 5. Received 2024 January 2; in original form 2023 September 29

## 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_{\star} \ge 10^{10.2} M_{\odot}$) and lower ($M_{\star} \le 10^{8} M_{\odot}$) ends of the stellar mass distribution contain cusps. However, Haloes containing galaxies with intermediate ($10^{8} M_{\odot} \le M_{\star} \le 10^{10.2} M_{\odot}$) stellar masses are generally cored, with typical halo masses between $10^{10.2} M_{\odot}$ and $10^{11.5} M_{\odot}$. 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_{\star} \ge 10^{9.5} M_{\odot}$) 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.

**Key words:** galaxies: formation – galaxies: evolution – galaxies: dwarf – galaxies: haloes – methods: numerical.

## 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

<footnote n="1"><sup>1</sup>Department of Physics and Astronomy, University of Victoria, BC V8P 5C2, Canada</footnote>

<footnote n="2"><sup>2</sup>Centre for Astrophysics Research, School of Physics, Astronomy and Mathematics, University of Hertfordshire, Hatfield, AL10 9AB, UK</footnote>

<footnote n="3"><sup>3</sup>Department of Astronomy and Yonsei University Observatory, Yonsei University, Seoul 03722, South Korea</footnote>

<footnote n="4"><sup>4</sup>Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Bd de l'Observatoire, Université Côte d'Azur, CS 34229, F-06304 Nice Cedex 4, France</footnote>

<footnote n="5"><sup>5</sup>Institut d'Astrophysique de Paris, Sorbonne Universités, UMPC Univ Paris 06 et CNRS, UMP 7095, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="6"><sup>6</sup>School of Physics and Astronomy, University of Nottingham, University Park, Nottingham NG7 2RD, UK</footnote>

<footnote n="7"><sup>7</sup>Korea Astronomy and Space Science Institute, 776 Daedeokdae-ro, Yuseong-gu, Daejeon 34055, South Korea</footnote>

<footnote n="8"><sup>8</sup>Steward Observatory, University of Arizona, 933 N. Cherry Avenue, Tucson, AZ 85719, USA</footnote>

<footnote n="9"><sup>9</sup>Department of Physics, University of Oxford, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="10"><sup>10</sup>School of Physics, Korea Institute for Advanced Study (KIAS), 85 Hoegiro, Dongdaemun-gu, Seoul 02455, South Korea</footnote>

<footnote n="11"><sup>11</sup>Observatoire Astronomique de Strasbourg, Université de Strasbourg, CNRS, UMR 7550, F-67000 Strasbourg, France</footnote>

<footnote>*E-mail: ryanjackson@uvic.ca</footnote>
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 $\Lambda$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 ($\rho \propto r^{-1}$, 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 $\rho \propto 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 $\Lambda$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
$\Lambda$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 \sim 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
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 $\Lambda$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 ($\rho \propto r^{-1}$, 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 $\rho \propto 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 $\Lambda$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
$\Lambda$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 \sim 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)$^3$. NH offers a maximum spatial resolution of 34 pc and mass resolutions of ~$10^4$ M$_\odot$ and ~$10^6$ M$_\odot$ 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)$^3$ volume (achieving a spatial resolution of 1 kpc), using 1024$^3$ 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$^3$ uniformly distributed cubic cells), with the same cosmology as that used in Horizon-AGN ($\Omega_m = 0.272$, $\Omega_b = 0.0455$, $\Omega_\Lambda = 0.728$, $H_0 = 70.4 \text{ km s}^{-1} \text{ Mpc}^{-1}$, 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^6$ M$_\odot$ (compared to $8 \times 10^7$ M$_\odot$ in Horizon-AGN), stellar mass resolution of $10^4$ M$_\odot$ (compared to $2 \times 10^6$ M$_\odot$ in Horizon-AGN) and a maximum spatial resolution of 34 pc (compared to 1 kpc in Horizon-AGN).<fnref n="1" /> The simulation has been performed down to $z \sim 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 \text{ H cm}^{-3}$ and a temperature lower than $2 \times 10^4$ K, following a Schmidt relation (Schmidt 1959). The star formation efficiency is dependent on the local turbulent Mach number and the virial parameter $\alpha_{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^4$ M$_\odot$ 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$_\odot$) 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$_\odot$ and $0.1$ M$_\odot$. 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^{51}$ erg and the supernova has a progenitor mass of $10$ M$_\odot$. Preheating 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^4$ M$_\odot$. 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

<footnote n="1">¹The gravitational force softening is equal to the local grid size.</footnote>
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)$^3$. NH offers a maximum spatial resolution of 34 pc and mass resolutions of ~$10^4$ M$_\odot$ and ~$10^6$ M$_\odot$ 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)$^3$ volume (achieving a spatial resolution of 1 kpc), using 1024$^3$ 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$^3$ uniformly distributed cubic cells), with the same cosmology as that used in Horizon-AGN ($\Omega_m = 0.272$, $\Omega_b = 0.0455$, $\Omega_\Lambda = 0.728$, $H_0 = 70.4$ km s$^{-1}$ Mpc$^{-1}$, 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^6$ M$_\odot$ (compared to $8 \times 10^7$ M$_\odot$ in Horizon-AGN), stellar mass resolution of $10^4$ M$_\odot$ (compared to $2 \times 10^6$ M$_\odot$ in Horizon-AGN) and a maximum spatial resolution of 34 pc (compared to 1 kpc in Horizon-AGN).<fnref n="1" /> The simulation has been performed down to $z \sim 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$^{-3}$ and a temperature lower than $2 \times 10^4$ K, following a Schmidt relation (Schmidt 1959). The star formation efficiency is dependent on the local turbulent Mach number and the virial parameter $\alpha_{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^4$ M$_\odot$ 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$_\odot$) 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$_\odot$ and 0.1 M$_\odot$. 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^{51}$ erg and the supernova has a progenitor mass of 10 M$_\odot$. Preheating 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^4$ M$_\odot$. 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

<footnote n="1">¹The gravitational force softening is equal to the local grid size.</footnote>
galaxy in the final snapshot at $z \sim 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*<sub>NH</sub>/*M*<sub>NHDM-only</sub> < 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<sup>7</sup> M<sub>☉</sub> *M*<sub>⋆</sub> < 10<sup>11.5</sup> M<sub>☉</sub>, thus
ensuring that we have a statistically significant sample.

We start by showing, in Fig. 1, how the *M*<sub>Halo</sub> versus *M*<sub>⋆</sub> and *M*<sub>⋆</sub> versus *R*<sub>eff</sub> 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*<sub>⋆</sub> versus *R*<sub>eff</sub> 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.2} M_{\odot}$) which, in the NH full run, contains a galaxy with $M_{\star} = 10^{10.9} M_{\odot}$. 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.31} M_{\odot}$), which contains a dwarf galaxy with $M_{\star} = 10^{8.15} M_{\odot}$ 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 (C/r<sup>α</sup>, 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 ($\alpha$), 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 \sim 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 $\alpha < -0.65$), whereas in intermediate mass Haloes (with DM masses between $10^{10.2} M_{\odot}$ and $10^{11.5} M_{\odot}$) the inner DM profiles are flatter and core-like ($\alpha > -0.65$). We pick this cut-off as no DM-only Haloes have profiles $\alpha > -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 $\alpha$ as a function of stellar and halo mass

Fig. 4 shows how the measured $\alpha$ for each halo varies with stellar mass (left panel), halo mass (centre panel), and $M_{\star}/M_{\text{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^9 M_{\odot}$), 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.
galaxy in the final snapshot at $z \sim 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_{\text{NH}}/M_{\text{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^7 M_{\odot} < M_{\star} < 10^{11.5} M_{\odot}$, thus ensuring that we have a statistically significant sample.

We start by showing, in Fig. 1, how the $M_{\text{Halo}}$ versus $M_{\star}$ and $M_{\star}$ versus $R_{\text{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_{\star}$ versus $R_{\text{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_{\text{halo}} = 10^{12.2} M_{\odot}$) which, in the NH full run, contains a galaxy with $M_{\star} = 10^{10.9} M_{\odot}$. 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_{\text{halo}} = 10^{10.31} M_{\odot}$), which contains a dwarf galaxy with $M_{\star} = 10^{8.15} M_{\odot}$ 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 ($C/r^{\alpha}$, 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 ($\alpha$), between 1–2 per cent of $R_{\text{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 \sim 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 $\alpha < -0.65$), whereas in intermediate mass Haloes (with DM masses between $10^{10.2} M_{\odot}$ and $10^{11.5} M_{\odot}$) the inner DM profiles are flatter and core-like ($\alpha > -0.65$). We pick this cut-off as no DM-only Haloes have profiles $\alpha > -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 $\alpha$ as a function of stellar and halo mass

Fig. 4 shows how the measured $\alpha$ for each halo varies with stellar mass (left panel), halo mass (centre panel), and $M_{\star}/M_{\text{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^9 M_{\odot}$), 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.
<figure><img src="image_1.png" /><figcaption>**Figure 1.** Left: $M_{\text{Halo}}$ versus $M_{\star}$ for the NH galaxy sample used in this study (blue circles). We compare to the abundance matching studies from Behroozi et al. (2013) and Moster et al. (2013) and other simulations: APOSTLE (Sawala et al. 2016), EAGLE CHT10 (Benfitez-Llambay et al. 2019), and NIHAO (Wang et al. 2015). On average, NH Haloes contain more stellar mass than other simulations, with the observational data lying between NH and the other simulations. Right: $M_{\star}$ versus $R_{\text{eff}}$ for the NH galaxy sample (blue circles). We also show the best fit observational results from Mowla et al. (2019) ($z = 0.37$). The sizes of NH galaxies are well-matched to the observations, particularly in the stellar mass range of interest in this study. For more details of how NH matches to observational results, see Dubois et al. (2021).</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>**Figure 2.** The DM density profiles of two galaxies in the NH simulation (blue lines) and their counterparts in the NH DM-only run (orange lines) at $z = 0.18$. The left panel shows an example of an NFW-like halo, where the density continues to increase towards the centre, with a slope of 1.68. In the DM-only run, the central regions have lower densities compared to the NH run, in which the presence of baryons in these central regions further deepens the potential. The right panel shows an example of a core-like halo with a slope of 0.43, in which the DM has been redistributed from the centre to larger radii in the NH run compared to the DM-only run. The NH simulation is therefore able to produce both cuspy and core-like profiles in the presence of baryons. In both examples, we indicate the region within which the slope is calculated using the dashed red line.</figcaption></figure>

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_{\odot}$ we find Haloes with both cores and cusps in NH unlike in Tollet et al. (2016). In the right panel we consider $M_{\star}/M_{\text{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
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Left: $M_{\text{Halo}}$ versus $M_{\star}$ for the NH galaxy sample used in this study (blue circles). We compare to the abundance matching studies from Behroozi et al. (2013) and Moster et al. (2013) and other simulations: APOSTLE (Sawala et al. 2016), EAGLE CHT10 (Benfitez-Llambay et al. 2019), and NIHAO (Wang et al. 2015). On average, NH Haloes contain more stellar mass than other simulations, with the observational data lying between NH and the other simulations. Right: $M_{\star}$ versus $R_{\text{eff}}$ for the NH galaxy sample (blue circles). We also show the best fit observational results from Mowla et al. (2019) ($z = 0.37$). The sizes of NH galaxies are well-matched to the observations, particularly in the stellar mass range of interest in this study. For more details of how NH matches to observational results, see Dubois et al. (2021).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** The DM density profiles of two galaxies in the NH simulation (blue lines) and their counterparts in the NH DM-only run (orange lines) at $z = 0.18$. The left panel shows an example of an NFW-like halo, where the density continues to increase towards the centre, with a slope of 1.68. In the DM-only run, the central regions have lower densities compared to the NH run, in which the presence of baryons in these central regions further deepens the potential. The right panel shows an example of a core-like halo with a slope of 0.43, in which the DM has been redistributed from the centre to larger radii in the NH run compared to the DM-only run. The NH simulation is therefore able to produce both cuspy and core-like profiles in the presence of baryons. In both examples, we indicate the region within which the slope is calculated using the dashed red line.</figcaption></figure>

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 $\sim 10^{10.2} M_{\odot}$ we find Haloes with both cores and cusps in NH unlike in Tollet et al. (2016). In the right panel we consider $M_{\star}/M_{\text{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
<figure><img src="image_3.png" /><figcaption>**Figure 3.** The slope of the density profile for galaxies in the full NH (i.e. DM + baryons) simulation, and its DM-only counterpart, as a function of halo mass. Galaxies from the full run are colour-coded by their stellar mass. We show our definition of a cored halo ($\alpha = -0.65$) using a dashed red line. Galaxies are generally cusps at the lowest halo masses, whilst at intermediate masses cores form. As halo masses continue to grow, profiles return to cusps at the highest masses. There is a transition halo mass ($M_{\text{halo}} \sim 10^{10.2} M_{\odot}$) where the profile slope generally flattens, compared to the DM-only run suggesting that baryonic processes are key for the creation of cores.</figcaption></figure>

(see Oh et al. (2011) for more details). Whilst these three studies find a similar relationship between $\alpha$ and $M_*/M_{\text{halo}}$, we find that NH galaxies show a different trend. In NH, cores are found across a range of $M_*/M_{\text{halo}}$ values. Unlike many other studies, which find a dependence on $M_*/M_{\text{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_{\text{halo}}$ relation where galaxies in NH only contain a core. Instead, it appears that once the halo crosses a threshold mass ($\sim 10^{10.2} M_{\odot}$) and enough stars have formed ($\sim 10^8 M_{\odot}$), 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_{\text{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 $\alpha-M_*/M_{\text{halo}}$ relationship. This gives credence to the idea that core formation is not solely dictated by $M_*/M_{\text{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 $\alpha > -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 ~70 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 $\sim 10^{10.3} M_{\odot}$. 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_{\odot}) < 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 $\alpha$. 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
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The slope of the density profile for galaxies in the full NH (i.e. DM + baryons) simulation, and its DM-only counterpart, as a function of halo mass. Galaxies from the full run are colour-coded by their stellar mass. We show our definition of a cored halo ($\alpha = -0.65$) using a dashed red line. Galaxies are generally cusps at the lowest halo masses, whilst at intermediate masses cores form. As halo masses continue to grow, profiles return to cusps at the highest masses. There is a transition halo mass ($M_{\text{halo}} \sim 10^{10.2} M_{\odot}$) where the profile slope generally flattens, compared to the DM-only run suggesting that baryonic processes are key for the creation of cores.</figcaption></figure>

(see Oh et al. (2011) for more details). Whilst these three studies find a similar relationship between $\alpha$ and $M_*/M_{\text{halo}}$, we find that NH galaxies show a different trend. In NH, cores are found across a range of $M_*/M_{\text{halo}}$ values. Unlike many other studies, which find a dependence on $M_*/M_{\text{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_{\text{halo}}$ relation where galaxies in NH only contain a core. Instead, it appears that once the halo crosses a threshold mass ($\sim 10^{10.2} M_{\odot}$) and enough stars have formed ($\sim 10^8 M_{\odot}$), 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_{\text{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 $\alpha-M_*/M_{\text{halo}}$ relationship. This gives credence to the idea that core formation is not solely dictated by $M_*/M_{\text{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 $\alpha > -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 ~70 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 $\sim 10^{10.3} M_{\odot}$. 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_{\odot}) < 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 $\alpha$. 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
<figure><img src="image_4.png" /></figure>

**Figure 4.** DM density profile slope between 1–2 per cent of the virial radius as a function of the stellar mass (right), halo mass (centre), and stellar/halo mass (right) of galaxies in the full NH (i.e. DM + baryons) run. Galaxies from the NH simulation are shown in red in all the plots. In all panels we compare to the simulation results of Tollet et al. (2016) (shown as blue diamonds) and find reasonable agreement in the stellar and halo mass ranges where cores are found. We compare to observations from Oh et al. (2015) in the left panel (purple stars), finding that, whilst observed cores have similar values for $\alpha$, they are generally found at slightly lower stellar masses than in NH. In the right panel (stellar/halo mass) we also compare to the simulation results of Di Cintio et al. (2014, green squares) and Lazar et al. (2020, black crosses) who find a clear trend between $M_*/M_{\text{Halo}}$ and $\alpha$. We do not see a similar relationship in NH, where $M_*/M_{\text{Halo}}$ increases with $\alpha$.

<figure><img src="image_5.png" /></figure>

**Figure 5.** The percentage of galaxies that contain a cored DM density profile (defined as $\alpha > -0.65$) as a function of look-back time in the NH simulation. Initially there are few to no cores in Haloes but, as star formation starts to increase, this value rises rapidly before remaining roughly constant until the end of the simulation. This suggests that for many galaxies cores are not destroyed and/or cores are constantly being created in galaxies.

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 \text{ H cm}^{-3}$) 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

<figure><img src="image_6.png" /></figure>

**Figure 6.** The halo mass of galaxies with a cored DM density profile (defined as $\alpha > -0.65$) as a function of look-back time. The darker line shows the evolution of the median halo mass and the edges of the shaded region represents the 16th and 84th percentiles. The median halo mass remains roughly constant throughout time, indicating that there may be (a) a threshold halo mass that allows for the creation and maintenance of a core and (b) a secondary threshold above which Haloes can no longer support a core and so the central regions returns to a cusp.

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,
<figure><img src="image_6.png" /></figure>

**Figure 4.** DM density profile slope between 1–2 per cent of the virial radius as a function of the stellar mass (right), halo mass (centre), and stellar/halo mass (right) of galaxies in the full NH (i.e. DM + baryons) run. Galaxies from the NH simulation are shown in red in all the plots. In all panels we compare to the simulation results of Tollet et al. (2016) (shown as blue diamonds) and find reasonable agreement in the stellar and halo mass ranges where cores are found. We compare to observations from Oh et al. (2015) in the left panel (purple stars), finding that, whilst observed cores have similar values for $\alpha$, they are generally found at slightly lower stellar masses than in NH. In the right panel (stellar/halo mass) we also compare to the simulation results of Di Cintio et al. (2014, green squares) and Lazar et al. (2020, black crosses) who find a clear trend between $M_*/M_{\text{Halo}}$ and $\alpha$. We do not see a similar relationship in NH, where $M_*/M_{\text{Halo}}$ increases with $\alpha$.

<figure><img src="image_7.png" /></figure>

**Figure 5.** The percentage of galaxies that contain a cored DM density profile (defined as $\alpha > -0.65$) as a function of look-back time in the NH simulation. Initially there are few to no cores in Haloes but, as star formation starts to increase, this value rises rapidly before remaining roughly constant until the end of the simulation. This suggests that for many galaxies cores are not destroyed and/or cores are constantly being created in galaxies.

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_{\text{H}} > 10 \text{ H cm}^{-3}$) 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

<figure><img src="image_8.png" /></figure>

**Figure 6.** The halo mass of galaxies with a cored DM density profile (defined as $\alpha > -0.65$) as a function of look-back time. The darker line shows the evolution of the median halo mass and the edges of the shaded region represents the 16th and 84th percentiles. The median halo mass remains roughly constant throughout time, indicating that there may be (a) a threshold halo mass that allows for the creation and maintenance of a core and (b) a secondary threshold above which Haloes can no longer support a core and so the central regions returns to a cusp.

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,
<figure><img src="image_7.png" /><figcaption>**Figure 7.** Halo mass versus stellar mass for centrals (left) and satellites (right), colour coded by the slope of the DM density profile. Haloes with cores are found in both populations which indicates that tidal interactions are not solely responsible for the creation of cores. Cores are found at lower halo and stellar masses in satellites than centrals, likely due to the effect of tidal stripping.</figcaption></figure>

<table><caption>**Table 1.** Median values of the masses of different components within galaxies in the mass range $8 < \log(M_s/M_\odot) < 9$, measured within both the central 2 kpc and the virial radius, for galaxies with cuspy (top row) and cored (bottom row) profiles, respectively. The stellar, DM and gas masses are consistent between both profile types, both within 2 kpc and the whole halo. However, the amount of star-forming gas is elevated by an order of magnitude in cored galaxies, indicating that star formation and the abundance of dense, cold gas plays a key role in the formation of cores.</caption><thead><tr><th rowspan="2">Profile</th><th>Log stellar mass [M<sub>☉</sub>]</th><th>Log DM mass [M<sub>☉</sub>]</th><th>Log gas mass [M<sub>☉</sub>]</th><th>Log SF gas mass [M<sub>☉</sub>]</th></tr><tr><th>2 kpc | Total</th><th>2 kpc | Total</th><th>2 kpc | Total</th><th>2 kpc | Total</th></tr></thead><tbody><tr><td>Cusp</td><td>7.64 | 8.27</td><td>8.58 | 10.31</td><td>8.44 | 9.29</td><td>6.98 | 6.98</td></tr><tr><td>Core</td><td>7.80 | 8.56</td><td>8.57 | 10.37</td><td>8.47 | 9.41</td><td>7.83 | 7.97</td></tr></tbody></table>

are hosted in a similarly sized Haloes (log $M_{\text{Halo}} \sim 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
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Halo mass versus stellar mass for centrals (left) and satellites (right), colour coded by the slope of the DM density profile. Haloes with cores are found in both populations which indicates that tidal interactions are not solely responsible for the creation of cores. Cores are found at lower halo and stellar masses in satellites than centrals, likely due to the effect of tidal stripping.</figcaption></figure>

<table><caption>**Table 1.** Median values of the masses of different components within galaxies in the mass range $8 < \log(M_s/M_\odot) < 9$, measured within both the central 2 kpc and the virial radius, for galaxies with cuspy (top row) and cored (bottom row) profiles, respectively. The stellar, DM and gas masses are consistent between both profile types, both within 2 kpc and the whole halo. However, the amount of star-forming gas is elevated by an order of magnitude in cored galaxies, indicating that star formation and the abundance of dense, cold gas plays a key role in the formation of cores.</caption><thead><tr><th rowspan="2">Profile</th><th>Log stellar mass [M<sub>☉</sub>]</th><th>Log DM mass [M<sub>☉</sub>]</th><th>Log gas mass [M<sub>☉</sub>]</th><th>Log SF gas mass [M<sub>☉</sub>]</th></tr><tr><th>2 kpc | Total</th><th>2 kpc | Total</th><th>2 kpc | Total</th><th>2 kpc | Total</th></tr></thead><tbody><tr><td>Cusp</td><td>7.64 | 8.27</td><td>8.58 | 10.31</td><td>8.44 | 9.29</td><td>6.98 | 6.98</td></tr><tr><td>Core</td><td>7.80 | 8.56</td><td>8.57 | 10.37</td><td>8.47 | 9.41</td><td>7.83 | 7.97</td></tr></tbody></table>

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
<figure><img src="image_8.png" /><figcaption>**Figure 8.** Top row: Enclosed mass (gas, DM, stellar) within 1 and 5 physical kpc versus look-back time for a cuspy halo (left) and a cored halo (right). The mass within 1 kpc is shown using the solid lines, while the mass within 5 kpc is shown using the dashed lines. Bottom row: sSFR of the galaxy (red) and inner DM density profile slope (black dashed line) versus look-back time (note that the slope is only shown when the inner DM density profile is above the resolution criteria set in subsection 3.1). Cored Haloes experience bursty star formation over an extended period of time which results in the removal of gas from within the central 1 kpc. When this occurs the DM responds to this removal of mass and also moves out of the central regions. The halo experiences multiple episodes of this process and eventually forms a core. Conversely, low mass cuspy Haloes have much briefer periods of star formation at high redshift before becoming quenched. This stops the removal of gas from within 1 kpc and means that the gravitational potential is not significantly altered and a core is never formed.</figcaption></figure>

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 ($\Phi_0$) 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
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Top row: Enclosed mass (gas, DM, stellar) within 1 and 5 physical kpc versus look-back time for a cuspy halo (left) and a cored halo (right). The mass within 1 kpc is shown using the solid lines, while the mass within 5 kpc is shown using the dashed lines. Bottom row: sSFR of the galaxy (red) and inner DM density profile slope (black dashed line) versus look-back time (note that the slope is only shown when the inner DM density profile is above the resolution criteria set in subsection 3.1). Cored Haloes experience bursty star formation over an extended period of time which results in the removal of gas from within the central 1 kpc. When this occurs the DM responds to this removal of mass and also moves out of the central regions. The halo experiences multiple episodes of this process and eventually forms a core. Conversely, low mass cuspy Haloes have much briefer periods of star formation at high redshift before becoming quenched. This stops the removal of gas from within 1 kpc and means that the gravitational potential is not significantly altered and a core is never formed.</figcaption></figure>

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 ($\Phi_0$) 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
<figure><img src="image_9.png" /></figure>

**Figure 9.** The gravitational potential versus log radius, at different look-back times (11, 10.87, 10.69, 10.57, and 3.01 Gyrs), for the halo of an example galaxy that forms a core. Dashed lines indicate the contribution to the potential from the DM whereas solid lines show the gas potential. These look-back times correspond to starbursts that contribute to the creation of the core. All values are normalized by their maximum, in order to compare the changing shape of the potential across cosmic time. We see that even on these short timescales, the bursts of star formation and subsequent movement of gas are able to alter the potential in the central regions of the halo, enabling the formation of a core.

galaxy, containing a cusp at the final timestep in NH, with a stellar and halo mass of $10^{11.3} M_{\odot}$ and $10^{12.9} M_{\odot}$, 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

<figure><img src="image_10.png" /></figure>

**Figure 10.** Top panel: The enclosed mass (within 1 kpc) for $M_{\star}$ (blue), $M_{\text{gas}}$ (green), and $M_{\text{DM}}$ (orange) versus look-back time (Gyr) for a massive halo ($M_{\text{Halo}} = 10^{12.9} M_{\odot}$) with a cusp at the final timestep of the full NH (i.e. DM + baryons) run. We also show the enclosed mass within 5 kpc for $M_{\star}$ (blue dashed line) to explore whether stars are being formed *in situ* or not. Bottom panel: The evolution of the sSFR (red) with look-back time and $\alpha$ (black dashed line) for this galaxy. At large look-back times, this halo had a cored profile when the gas and DM masses were roughly equal in the centre. However, ongoing star formation with a high sSFR changes the makeup of the central masses. As stars begin to dominate the mass in the centre of the halo, the core disappears as a cusp is reformed.

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 ~30 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.
<figure><img src="image_11.png" /></figure>

**Figure 9.** The gravitational potential versus log radius, at different look-back times (11, 10.87, 10.69, 10.57, and 3.01 Gyrs), for the halo of an example galaxy that forms a core. Dashed lines indicate the contribution to the potential from the DM whereas solid lines show the gas potential. These look-back times correspond to starbursts that contribute to the creation of the core. All values are normalized by their maximum, in order to compare the changing shape of the potential across cosmic time. We see that even on these short timescales, the bursts of star formation and subsequent movement of gas are able to alter the potential in the central regions of the halo, enabling the formation of a core.

galaxy, containing a cusp at the final timestep in NH, with a stellar and halo mass of $10^{11.3} M_{\odot}$ and $10^{12.9} M_{\odot}$, 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

<figure><img src="image_12.png" /></figure>

**Figure 10.** Top panel: The enclosed mass (within 1 kpc) for $M_{\star}$ (blue), $M_{\text{gas}}$ (green), and $M_{\text{DM}}$ (orange) versus look-back time (Gyr) for a massive halo ($M_{\text{Halo}} = 10^{12.9} M_{\odot}$) with a cusp at the final timestep of the full NH (i.e. DM + baryons) run. We also show the enclosed mass within 5 kpc for $M_{\star}$ (blue dashed line) to explore whether stars are being formed *in situ* or not. Bottom panel: The evolution of the sSFR (red) with look-back time and $\alpha$ (black dashed line) for this galaxy. At large look-back times, this halo had a cored profile when the gas and DM masses were roughly equal in the centre. However, ongoing star formation with a high sSFR changes the makeup of the central masses. As stars begin to dominate the mass in the centre of the halo, the core disappears as a cusp is reformed.

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 ~30 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 $\Lambda$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:

(i) The central DM density slope ($\alpha$) 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 $\alpha > -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.

(ii) 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.

(iii) 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.

(iv) Haloes that contain galaxies in the upper ($M_{*} \ge 10^{10.2} M_{\odot}$) and lower ($M_{*} \le 10^{8} M_{\odot}$) ends of the stellar mass distribution contain cusps. However, in Haloes that contain galaxies that have intermediate ($10^{8} M_{\odot} \le M_{*} \le 10^{10.2} M_{\odot}$) stellar masses, $\alpha$ 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_{\odot}$ and $10^{11.5} M_{\odot}$.

(v) The evolution of $\alpha$ 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).

(vi) The fraction of galaxies that contain a core becomes roughly constant (~15 per cent) after a redshift of $z \sim 1.5$, with cores not forming until the cosmic star formation rate increases around a redshift of $z \sim 2$. The mass of Haloes that contain a core is roughly constant in the mass range $10^{10.2} M_{\odot} < M_{*} < 10^{11} M_{\odot}$, with a median value of $10^{10.3} M_{\odot}$, 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 $\Lambda$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 $\Lambda$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.
# 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 $\Lambda$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:

(i) The central DM density slope ($\alpha$) 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 $\alpha > -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.

(ii) 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.

(iii) 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.

(iv) Haloes that contain galaxies in the upper ($M_{*} \ge 10^{10.2} M_{\odot}$) and lower ($M_{*} \le 10^{8} M_{\odot}$) ends of the stellar mass distribution contain cusps. However, in Haloes that contain galaxies that have intermediate ($10^{8} M_{\odot} \le M_{*} \le 10^{10.2} M_{\odot}$) stellar masses, $\alpha$ 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_{\odot}$ and $10^{11.5} M_{\odot}$.

(v) The evolution of $\alpha$ 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).

(vi) The fraction of galaxies that contain a core becomes roughly constant (~15 per cent) after a redshift of $z \sim 1.5$, with cores not forming until the cosmic star formation rate increases around a redshift of $z \sim 2$. The mass of Haloes that contain a core is roughly constant in the mass range $10^{10.2} M_{\odot} < M_{*} < 10^{11} M_{\odot}$, with a median value of $10^{10.3} M_{\odot}$, 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 $\Lambda$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 $\Lambda$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.

# REFERENCES

Alam S. et al., 2015, ApJS, 219, 12

Aubert D., Pichon C., Colombi S., 2004, MNRAS, 352, 376

Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, L57

Benítez-Llambay A., Frenk C. S., Ludlow A. D., Navarro J. F., 2019, MNRAS, 488, 2387

Benson A. J., Bower R. G., Frenk C. S., Lacey C. G., Baugh C. M., Cole S., 2003, ApJ, 599, L38

Bode P., Ostriker J. P., Turok N., 2001, ApJ, 556, L93

Bondi H., Hoyle F., 1944, MNRAS, 104, 273

Bose S. et al., 2019, MNRAS, 486, 4790

Bower R. G., Benson A. J., Malbon R., Helly J. C., Frenk C. S., Baugh C. M., Cole S., Lacey C. G., 2006, MNRAS, 370, 645

Bullock J. S., Boylan-Kolchin M., 2017, ARA&A, 55, 343

Burger J. D., Zavala J., 2019, MNRAS, 485, 1008

Burger J. D., Zavala J., 2021, ApJ, 921, L126

Burger J. D., Zavala J., Sales L. V., Vogelsberger M., Marinacci F., Torrey P., 2022, MNRAS, 513, 3458

Burkert A., 1995, ApJ, 447, L25

Chabrier G., 2005, The Initial Mass Function: From Salpeter 1955 to 2005. p. 41

Chan T. K., Kereš D., Oñorbe J., Hopkins P. F., Muratov A. L., Faucher-Giguère C. A., Quataert E., 2015, MNRAS, 454, 2981

Chan T. K., Kereš D., Wetzel A., Hopkins P. F., Faucher-Giguère C. A., El-Badry K., Garrison-Kimmel S., Boylan-Kolchin M., 2018, MNRAS, 1094

Chua K. T. E., Pillepich A., Vogelsberger M., Hernquist L., 2019, MNRAS, 484, 476

Cole S., Lacey C. G., Baugh C. M., Frenk C. S., 2000, MNRAS, 319, 168

Croton D. J. et al., 2006, MNRAS, 365, 11

Davé R., Anglés-Alcázar D., Narayanan D., Li Q., Rafieferantsoa M. H., Appleby S., 2019, MNRAS, 486, 2827

Davis F. et al., 2022, MNRAS, 511, 4109

de Blok W. J. G., 2010, Adv. Astron., 2010, 789293

de Blok W. J. G., McGaugh S. S., Bosma A., Rubin V. C., 2001, ApJ, 552, L23

de Blok W. J. G., Walter F., Brinks E., Trachternach C., Oh S. H., Kennicutt R. C., Jr, 2008, AJ, 136, 2648

Dekel A., Silk J., 1986, ApJ, 303, L39

Dekel A. et al., 2021, MNRAS, 508, 999

Di Cintio A., Brook C. B., Macciò A. V., Stinson G. S., Knebe A., Dutton A. A., Wadsley J., 2014, MNRAS, 437, 415

Di Cintio A., Brook C. B., Dutton A. A., Macciò A. V., Obreja A., Dekel A., 2017, MNRAS, 466, 1

Dodelson S., Widrow L. M., 1994, Phys. Rev. Lett., 72, 17

Downing E. R., Oman K. A., 2023, MNRAS, 522, 3318

Dubois Y., Teyssier R., 2008, A&A, 482, 13

Dubois Y., Devriendt J., Slyz A., Teyssier R., 2012, MNRAS, 420, 2662

Dubois Y., Volonteri M., Silk J., 2014a, MNRAS, 440, 1590

Dubois Y. et al., 2014b, MNRAS, 444, 1453

Dubois Y. et al., 2021, A&A, 651, 109

Dutton A. A. et al., 2016, MNRAS, 461, 2658

Dutton A. A., Buck T., Macciò A. V., Dixon K. L., Blank M., Obreja A., 2020, MNRAS, 499, 2648

Einasto J., 1965, Trudy Astrofizicheskogo Instituta Alma-Ata, 5, 87

Eisenstein D. J., Hut P., 1998, ApJ, 498, L137

Fitts A. et al., 2017, MNRAS, 471, 3547

Flores R. A., Primack J. R., 1994, ApJ, 427, L1

Freundlich J., Dekel A., Jiang F., Ishai G., Cornuault N., Lapiner S., Dutton A. A., Macciò A. V., 2020, MNRAS, 491, 4523

Garrison-Kimmel S., Rocha M., Boylan-Kolchin M., Bullock J. S., Lally J., 2013, MNRAS, 433, 3539

Geen S., Rosdahl J., Blaizot J., Devriendt J., Slyz A., 2015, MNRAS, 448, 3248

Genina A. et al., 2018, MNRAS, 474, 1398

Gilmore G., Wilkinson M. I., Wyse R. F. G., Kleyna J. T., Koch A., Evans N. W., Grebel E. K., 2007, ApJ, 663, L948

Gnedin O. Y., Zhao H., 2002, MNRAS, 333, 299

Governato F. et al., 2010, Nature, 463, 203

Haardt F., Madau P., 1996, ApJ, 461, L20

Harvey D., Revaz Y., Robertson A., Hausmann L., 2018, MNRAS, 481, 89

Hoyle F., Lyttleton R. A., 1939, Mathematical Proceedings of the Cambridge Philosophical Society, 35, 592

Jackson R. A. et al., 2021a, MNRAS, 502, 1785

Jackson R. A. et al., 2021b, MNRAS, 502, 4262

Jahn E. D. et al., 2023, MNRAS, 520, 461

Kaviraj S. et al., 2017, MNRAS, 467, 4739

Kimm T., Cen R., 2014, ApJ, 788, L121

Kimm T., Katz H., Haehnelt M., Rosdahl J., Devriendt J., Slyz A., 2017, MNRAS, 466, 4826

Komatsu E. et al., 2011, ApJS, 192, 18

Kormendy J., Fisher D. B., Cornell M. E., Bender R., 2009, ApJS, 182, 216

Kuzio de Naray R., McGaugh S. S., de Blok W. J. G., 2008, ApJ, 676, L920

Lazar A. et al., 2020, MNRAS, 497, 2393

Lelli F., McGaugh S. S., Schombert J. M., 2016, AJ, 152, 157

Li Z., Dekel A., Mandelker N., Freundlich J., François T. L., 2023, MNRAS, 518, 5356

Madau P., Dickinson M., 2014, ARA&A, 52, 415

Martin G. et al., 2019, MNRAS, 485, 796

Martin G. et al., 2021, MNRAS, 500, 4937

Martizzi D., Teyssier R., Moore B., 2013, MNRAS, 432, 1947

Mashchenko S., Couchman H. M. P., Wadsley J., 2006, Nature, 442, 539

Mashchenko S., Wadsley J., Couchman H. M. P., 2008, Science, 319, 174

McKinney J. C., Tchekhovskoy A., Blandford R. D., 2012, MNRAS, 423, 3083

Moore B., 1994, Nature, 370, 629

Moster B. P., Naab T., White S. D. M., 2013, MNRAS, 428, 3121

Mowla L., van der Wel A., van Dokkum P., Miller T. B., 2019, ApJ, 872, L13

Navarro J. F., Eke V. R., Frenk C. S., 1996a, MNRAS, 283, L72

Navarro J. F., Frenk C. S., White S. D. M., 1996b, ApJ, 462, L563

Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, L493

Nelson D. et al., 2019, Comput. Astrophys. Cosmol., 6, 2

Oh S.-H., de Blok W. J. G., Brinks E., Walter F., Kennicutt R. C., Jr, 2011, AJ, 141, 193

Oh S.-H. et al., 2015, AJ, 149, 180

Oman K. A. et al., 2015, MNRAS, 452, 3650

Oman K. A., Marasco A., Navarro J. F., Frenk C. S., Schaye J., Benítez-Llambay A., 2019, MNRAS, 482, 821

Ostriker E. C., 1999, ApJ, 513, L252

Peirani S. et al., 2017, MNRAS, 472, 2153

Pontzen A., Governato F., 2012, MNRAS, 421, 3464

Pontzen A., Governato F., 2014, Nature, 506, 171

Prunet S., Pichon C., Aubert D., Pogosyan D., Teyssier R., Gottloeber S., 2008, ApJS, 178, 179

Rasera Y., Teyssier R., 2006, A&A, 445, 1

Read J. I., Gilmore G., 2005, MNRAS, 356, 107

Rocha M., Peter A. H. G., Bullock J. S., Kaplinghat M., Garrison-Kimmel S., Oñorbe J., Moustakas L. A., 2013, MNRAS, 430, 81

Roper F. A., Oman K. A., Frenk C. S., Benítez-Llambay A., Navarro J. F., Santos-Santos I. M. E., 2023, MNRAS, 521, 1316

Rosen A., Bregman J. N., 1995, ApJ, 440, L634

Sales L. V., Wetzel A., Fattahi A., 2022, Nature Astron.

Santos-Santos I. M. E. et al., 2020, MNRAS, 495, 58

Sawala T. et al., 2016, MNRAS, 457, 1931

Schaller M. et al., 2015, MNRAS, 451, 1247

Schaye J. et al., 2015, MNRAS, 446, 521

Schmidt M., 1959, ApJ, 129, L243

Shakura N. I., Sunyaev R. A., 1973, A&A, 500, 33
# 6 DATA AVAILABILITY

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

# REFERENCES

Alam S. et al., 2015, ApJS, 219, 12

Aubert D., Pichon C., Colombi S., 2004, MNRAS, 352, 376

Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, L57

Benítez-Llambay A., Frenk C. S., Ludlow A. D., Navarro J. F., 2019, MNRAS, 488, 2387

Benson A. J., Bower R. G., Frenk C. S., Lacey C. G., Baugh C. M., Cole S., 2003, ApJ, 599, L38

Bode P., Ostriker J. P., Turok N., 2001, ApJ, 556, L93

Bondi H., Hoyle F., 1944, MNRAS, 104, 273

Bose S. et al., 2019, MNRAS, 486, 4790

Bower R. G., Benson A. J., Malbon R., Helly J. C., Frenk C. S., Baugh C. M., Cole S., Lacey C. G., 2006, MNRAS, 370, 645

Bullock J. S., Boylan-Kolchin M., 2017, ARA&A, 55, 343

Burger J. D., Zavala J., 2019, MNRAS, 485, 1008

Burger J. D., Zavala J., 2021, ApJ, 921, L126

Burger J. D., Zavala J., Sales L. V., Vogelsberger M., Marinacci F., Torrey P., 2022, MNRAS, 513, 3458

Burkert A., 1995, ApJ, 447, L25

Chabrier G., 2005, The Initial Mass Function: From Salpeter 1955 to 2005. p. 41

Chan T. K., Kereš D., Oñorbe J., Hopkins P. F., Muratov A. L., Faucher-Giguère C. A., Quataert E., 2015, MNRAS, 454, 2981

Chan T. K., Kereš D., Wetzel A., Hopkins P. F., Faucher-Giguère C. A., El-Badry K., Garrison-Kimmel S., Boylan-Kolchin M., 2018, MNRAS, 1094

Chua K. T. E., Pillepich A., Vogelsberger M., Hernquist L., 2019, MNRAS, 484, 476

Cole S., Lacey C. G., Baugh C. M., Frenk C. S., 2000, MNRAS, 319, 168

Croton D. J. et al., 2006, MNRAS, 365, 11

Davé R., Anglés-Alcázar D., Narayanan D., Li Q., Rafieferantsoa M. H., Appleby S., 2019, MNRAS, 486, 2827

Davis F. et al., 2022, MNRAS, 511, 4109

de Blok W. J. G., 2010, Adv. Astron., 2010, 789293

de Blok W. J. G., McGaugh S. S., Bosma A., Rubin V. C., 2001, ApJ, 552, L23

de Blok W. J. G., Walter F., Brinks E., Trachternach C., Oh S. H., Kennicutt R. C., Jr, 2008, AJ, 136, 2648

Dekel A., Silk J., 1986, ApJ, 303, L39

Dekel A. et al., 2021, MNRAS, 508, 999

Di Cintio A., Brook C. B., Macciò A. V., Stinson G. S., Knebe A., Dutton A. A., Wadsley J., 2014, MNRAS, 437, 415

Di Cintio A., Brook C. B., Dutton A. A., Macciò A. V., Obreja A., Dekel A., 2017, MNRAS, 466, 1

Dodelson S., Widrow L. M., 1994, Phys. Rev. Lett., 72, 17

Downing E. R., Oman K. A., 2023, MNRAS, 522, 3318

Dubois Y., Teyssier R., 2008, A&A, 482, 13

Dubois Y., Devriendt J., Slyz A., Teyssier R., 2012, MNRAS, 420, 2662

Dubois Y., Volonteri M., Silk J., 2014a, MNRAS, 440, 1590

Dubois Y. et al., 2014b, MNRAS, 444, 1453

Dubois Y. et al., 2021, A&A, 651, 109

Dutton A. A. et al., 2016, MNRAS, 461, 2658

Dutton A. A., Buck T., Macciò A. V., Dixon K. L., Blank M., Obreja A., 2020, MNRAS, 499, 2648

Einasto J., 1965, Trudy Astrofizicheskogo Instituta Alma-Ata, 5, 87

Eisenstein D. J., Hut P., 1998, ApJ, 498, L137

Fitts A. et al., 2017, MNRAS, 471, 3547

Flores R. A., Primack J. R., 1994, ApJ, 427, L1

Freundlich J., Dekel A., Jiang F., Ishai G., Cornuault N., Lapiner S., Dutton A. A., Macciò A. V., 2020, MNRAS, 491, 4523

Garrison-Kimmel S., Rocha M., Boylan-Kolchin M., Bullock J. S., Lally J., 2013, MNRAS, 433, 3539

Geen S., Rosdahl J., Blaizot J., Devriendt J., Slyz A., 2015, MNRAS, 448, 3248

Genina A. et al., 2018, MNRAS, 474, 1398

Gilmore G., Wilkinson M. I., Wyse R. F. G., Kleyna J. T., Koch A., Evans N. W., Grebel E. K., 2007, ApJ, 663, L948

Gnedin O. Y., Zhao H., 2002, MNRAS, 333, 299

Governato F. et al., 2010, Nature, 463, 203

Haardt F., Madau P., 1996, ApJ, 461, L20

Harvey D., Revaz Y., Robertson A., Hausmann L., 2018, MNRAS, 481, 89

Hoyle F., Lyttleton R. A., 1939, Mathematical Proceedings of the Cambridge Philosophical Society, 35, 592

Jackson R. A. et al., 2021a, MNRAS, 502, 1785

Jackson R. A. et al., 2021b, MNRAS, 502, 4262

Jahn E. D. et al., 2023, MNRAS, 520, 461

Kaviraj S. et al., 2017, MNRAS, 467, 4739

Kimm T., Cen R., 2014, ApJ, 788, L121

Kimm T., Katz H., Haehnelt M., Rosdahl J., Devriendt J., Slyz A., 2017, MNRAS, 466, 4826

Komatsu E. et al., 2011, ApJS, 192, 18

Kormendy J., Fisher D. B., Cornell M. E., Bender R., 2009, ApJS, 182, 216

Kuzio de Naray R., McGaugh S. S., de Blok W. J. G., 2008, ApJ, 676, L920

Lazar A. et al., 2020, MNRAS, 497, 2393

Lelli F., McGaugh S. S., Schombert J. M., 2016, AJ, 152, 157

Li Z., Dekel A., Mandelker N., Freundlich J., François T. L., 2023, MNRAS, 518, 5356

Madau P., Dickinson M., 2014, ARA&A, 52, 415

Martin G. et al., 2019, MNRAS, 485, 796

Martin G. et al., 2021, MNRAS, 500, 4937

Martizzi D., Teyssier R., Moore B., 2013, MNRAS, 432, 1947

Mashchenko S., Couchman H. M. P., Wadsley J., 2006, Nature, 442, 539

Mashchenko S., Wadsley J., Couchman H. M. P., 2008, Science, 319, 174

McKinney J. C., Tchekhovskoy A., Blandford R. D., 2012, MNRAS, 423, 3083

Moore B., 1994, Nature, 370, 629

Moster B. P., Naab T., White S. D. M., 2013, MNRAS, 428, 3121

Mowla L., van der Wel A., van Dokkum P., Miller T. B., 2019, ApJ, 872, L13

Navarro J. F., Eke V. R., Frenk C. S., 1996a, MNRAS, 283, L72

Navarro J. F., Frenk C. S., White S. D. M., 1996b, ApJ, 462, L563

Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, L493

Nelson D. et al., 2019, Comput. Astrophys. Cosmol., 6, 2

Oh S.-H., de Blok W. J. G., Brinks E., Walter F., Kennicutt R. C., Jr, 2011, AJ, 141, 193

Oh S.-H. et al., 2015, AJ, 149, 180

Oman K. A. et al., 2015, MNRAS, 452, 3650

Oman K. A., Marasco A., Navarro J. F., Frenk C. S., Schaye J., Benítez-Llambay A., 2019, MNRAS, 482, 821

Ostriker E. C., 1999, ApJ, 513, L252

Peirani S. et al., 2017, MNRAS, 472, 2153

Pontzen A., Governato F., 2012, MNRAS, 421, 3464

Pontzen A., Governato F., 2014, Nature, 506, 171

Prunet S., Pichon C., Aubert D., Pogosyan D., Teyssier R., Gottloeber S., 2008, ApJS, 178, 179

Rasera Y., Teyssier R., 2006, A&A, 445, 1

Read J. I., Gilmore G., 2005, MNRAS, 356, 107

Rocha M., Peter A. H. G., Bullock J. S., Kaplinghat M., Garrison-Kimmel S., Oñorbe J., Moustakas L. A., 2013, MNRAS, 430, 81

Roper F. A., Oman K. A., Frenk C. S., Benítez-Llambay A., Navarro J. F., Santos-Santos I. M. E., 2023, MNRAS, 521, 1316

Rosen A., Bregman J. N., 1995, ApJ, 440, L634

Sales L. V., Wetzel A., Fattahi A., 2022, Nature Astron.

Santos-Santos I. M. E. et al., 2020, MNRAS, 495, 58

Sawala T. et al., 2016, MNRAS, 457, 1931

Schaller M. et al., 2015, MNRAS, 451, 1247

Schaye J. et al., 2015, MNRAS, 446, 521

Schmidt M., 1959, ApJ, 129, L243

Shakura N. I., Sunyaev R. A., 1973, A&A, 500, 33
Somerville R. S., Primack J. R., 1999, MNRAS, 310, 1087

Spergel D. N., Steinhardt P. J., 2000, Phys. Rev. Lett., 84, 3760

Springel V. et al., 2008, MNRAS, 391, 1685

Strigari L. E., Frenk C. S., White S. D. M., 2017, ApJ, 838, L123

Sutherland R. S., Dopita M. A., 1993, ApJS, 88, 253

Teyssier R., 2002, A&A, 385, 337

Teyssier R., Pontzen A., Dubois Y., Read J. I., 2013, MNRAS, 429, 3068

Tollet E. et al., 2016, MNRAS, 456, 3542

Tulin S., Yu H.-B., 2018, Phys. Rep., 730, 1

Tweed D., Devriendt J., Blaizot J., Colombi S., Slyz A., 2009, A&A, 506, 647

Vogelsberger M., Zavala J., Loeb A., 2012, MNRAS, 423, 3740

Vogelsberger M. et al., 2014, Nature, 509, 177

Walker M. G., Peñarrubia J., 2011, ApJ, 742, L20

Wang L., Dutton A. A., Stinson G. S., Macciò A. V., Penzo C., Kang X., Keller B. W., Wadsley J., 2015, MNRAS, 454, 83

Weinberg M. D., 2001, MNRAS, 328, 311

White S. D. M., Rees M. J., 1978, MNRAS, 183, 341

Wright A. H. et al., 2017, MNRAS, 470, 283

Yoshida N., Springel V., White S. D. M., Tormen G., 2000, ApJ, 544, L87

## 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^8 M_{\odot} < M_{\star} < 10^{10.2} M_{\odot}$ 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.

<figure><img src="image_11.png" /></figure>

**Figure A1.** Inner DM density profile slope ($\alpha$) versus stellar mass. We calculate the slope of the inner DM density profile using three different radial ranges $0.01 < r/R_{vir} < 0.02$ (blue), $1 < r/kpc < 2$ (orange), and $1 < r/r_{soft} < 5$ (purple) where $r_{soft}$ is ~500 pc. All three radial ranges produce consistent values for $\alpha$, suggesting that the results are relatively independent of which range is selected.

This paper has been typeset from a TEX/LATEX file prepared by the author.
Somerville R. S., Primack J. R., 1999, MNRAS, 310, 1087

Spergel D. N., Steinhardt P. J., 2000, Phys. Rev. Lett., 84, 3760

Springel V. et al., 2008, MNRAS, 391, 1685

Strigari L. E., Frenk C. S., White S. D. M., 2017, ApJ, 838, L123

Sutherland R. S., Dopita M. A., 1993, ApJS, 88, 253

Teyssier R., 2002, A&A, 385, 337

Teyssier R., Pontzen A., Dubois Y., Read J. I., 2013, MNRAS, 429, 3068

Tollet E. et al., 2016, MNRAS, 456, 3542

Tulin S., Yu H.-B., 2018, Phys. Rep., 730, 1

Tweed D., Devriendt J., Blaizot J., Colombi S., Slyz A., 2009, A&A, 506, 647

Vogelsberger M., Zavala J., Loeb A., 2012, MNRAS, 423, 3740

Vogelsberger M. et al., 2014, Nature, 509, 177

Walker M. G., Peñarrubia J., 2011, ApJ, 742, L20

Wang L., Dutton A. A., Stinson G. S., Macciò A. V., Penzo C., Kang X., Keller B. W., Wadsley J., 2015, MNRAS, 454, 83

Weinberg M. D., 2001, MNRAS, 328, 311

White S. D. M., Rees M. J., 1978, MNRAS, 183, 341

Wright A. H. et al., 2017, MNRAS, 470, 283

Yoshida N., Springel V., White S. D. M., Tormen G., 2000, ApJ, 544, L87

## 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^8 M_{\odot} < M_{\star} < 10^{10.2} M_{\odot}$ 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.

<figure><img src="image_13.png" /></figure>

**Figure A1.** Inner DM density profile slope ($\alpha$) versus stellar mass. We calculate the slope of the inner DM density profile using three different radial ranges $0.01 < r/R_{vir} < 0.02$ (blue), $1 < r/kpc < 2$ (orange), and $1 < r/r_{soft} < 5$ (purple) where $r_{soft}$ is ~500 pc. All three radial ranges produce consistent values for $\alpha$, suggesting that the results are relatively independent of which range is selected.

This paper has been typeset from a TeX/LATEX file prepared by the author.