# The diversity of rotation curves of galaxies in the NewHorizon cosmological simulation

R A Jackson, J F Navarro, I M E Santos-Santos, S Kaviraj, S K Yi, S Peirani, Y Dubois, G Martin, J E G Devriendt, A Slyz, et al.

► To cite this version:

R A Jackson, J F Navarro, I M E Santos-Santos, S Kaviraj, S K Yi, et al. The diversity of rota- tion curves of galaxies in the NewHorizon cosmological simulation. Monthly Notices of the Royal Astronomical Society, 2025, 539, pp.3797-3807. 10.1093/mnras/staf667. insu-05075594

HAL Id: insu-05075594

https://insu.hal.science/insu-05075594v1

Submitted on 20 May 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scie- ntific 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The diversity of rotation curves of galaxies in the NEWHORIZON cosmological simulation

R. A. Jackson<sup>1,★</sup>, J. F. Navarro,<sup>1</sup> I. M. E. Santos-Santos<sup>1,2</sup>, S. Kaviraj<sup>1,3</sup>, S. K. Yi,<sup>4</sup> S. Peirani,<sup>5,6</sup>

Y. Dubois,<sup>7</sup> G. Martin<sup>1,8,9,10</sup> J. E. G. Devriendt,<sup>11</sup> A. Slyz,<sup>11</sup> C. Pichon,<sup>7,12</sup> M. Volonteri,<sup>7</sup> T. Kimm<sup>13</sup>
and K. Kraljic<sup>13</sup>

Accepted 2025 April 13. Received 2025 April 13; in original form 2024 July 22

## ABSTRACT

We use the cosmological hydrodynamical simulation NEWHORIZON to study the effects of the baryonic component on the inner mass profile of dark matter haloes of isolated galaxies ($10^6 < M_*/M_\odot < 10^{11.5}$). Dark matter deficits ('cores') develop only in galaxies in a narrow range of stellar mass, $5 \times 10^7 < M_*/M_\odot < 3 \times 10^9$. The lower stellar mass limit arises because a minimum amount of star formation is required to drive the baryonic outflows that redistribute dark matter and create a core. The upper limit roughly coincides with the total amount of dark matter initially contained within the innermost 2 kpc ($1 \pm 0.5 \times 10^9 M_\odot$), which roughly coincides with the stellar half-mass radius of these dwarfs. This enclosed mass is quite insensitive to the total virial mass of the system. The same upper limit applies to other simulations, like NIHAO and EAGLE-CHT10, despite their rather different galaxy formation efficiencies. This suggests that it is the galaxy total stellar mass that determines when a core is formed, and not the galaxy-to-dark halo mass ratio, as argued in earlier work. This is consistent with a back-of-the-envelope estimate for a SN-induced rate of orbital diffusion. Although NEWHORIZON dwarfs reproduce the observed diversity of rotation curves better than other simulations, there are significant differences in the gravitational importance of baryons in the inner regions of dwarfs compared to observations. These differences prevent us from concluding that cosmological simulations are currently fully able to account for the observed diversity of rotation curve shapes.

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

## 1 INTRODUCTION

A current key tension in astronomy is the ‘core-cusp’ problem (Flores & Primack 1994; Moore 1994) and its derivative issues, which describes the observed discrepancy between the inner dark matter (DM) content of low mass haloes in cosmological simulations and that inferred from the rotation curves of dwarf galaxies (de Blok 2010). LCDM simulations predict a cuspy DM density profile, where density increases towards the centre of the halo ($\rho \propto r^{-1}$) (Springel et al. 2008). This is commonly described by a Navarro–Frenk–White

(NFW) profile (Navarro, Frenk & White 1996b, 1997), although other fits have also been suggested (e.g. Einasto 1965; Navarro et al. 2010). In contrast, observations of dwarf galaxies often find slowly rising rotation curves which suggest that the density profile has a much shallower (or even flat) density profile (i.e. a constant density ‘core’) in their inner regions (e.g. Burkert 1995; de Blok et al. 2001, 2008; Gilmore et al. 2007; Kuzio de Naray, McGaugh & de Blok 2008; Kormendy et al. 2009; Oh et al. 2011, 2015; Lelli, McGaugh & Schombert 2016).

Given that dwarf galaxies remain difficult to study, both observationally and theoretically, due to the shallow detection limits of recent wide-area surveys such as the SDSS (Alam et al. 2015) and the relatively low resolution of large volume simulations; e.g.

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

<footnote n="2"><sup>2</sup>Institute for Computational Cosmology, Durham University, South Road, Durham DH1 3LE, UK</footnote>

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

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

<footnote n="5"><sup>5</sup>ILANCE, CNRS – University of Tokyo International Research Laboratory, Kashiwa, Chiba 277-8582, Japan</footnote>

<footnote n="6"><sup>6</sup>Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan</footnote>

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

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

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

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

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

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

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

<footnote>* E-mail: ryanjackson@uvic.ca</footnote>
Horizon-AGN (Dubois et al. 2014; Kaviraj et al. 2017), EAGLE (Schaye et al. 2015), Illustris (Vogelsberger et al. 2014), APOSTLE (Sawala et al. 2016), Simba (Davé et al. 2019), or IllustrisTNG (Nelson et al. 2019). It is probably unsurprising that such tensions between observations and LCDM predictions have been reported in the low-mass galaxy regime (e.g. Bullock & Boylan-Kolchin 2017; Sales, Wetzel & Fattahi 2022).

To try to resolve this, many studies have shown that it is possible to turn the DM cusps into cores via baryonic processes (Navarro, Eke & Frenk 1996a; Gnedin & Zhao 2002; Read & Gilmore 2005; Mashchenko, Wadsley & Couchman 2008; Peirani & de Freitas Pacheco 2008; Pontzen & Governato 2012; Garrison-Kimmel et al. 2013; Teyssier et al. 2013; Di Cintio et al. 2014; Chan et al. 2015; Tollet et al. 2016; Fitts et al. 2017; Peirani et al. 2017; Lazar et al. 2020; Burger & Zavala 2021; Jackson et al. 2024). Feedback processes driven by active galactic nuclei (AGNs) and/or supernovae (SNe) are capable of redistributing mass in DM haloes in different ways. In low mass galaxies SN feedback has been shown to 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. 2021). This gas removal has been invoked to explain the formation of cores in several simulations, via non-adiabatic fluctuations of the potential (Navarro et al. 1996a; Pontzen & Governato 2012). Such fluctuations lead to energy gains for DM particles, which then migrate to orbits with larger apocentres, lowering the central DM densities.

However, even models that successfully produce ‘baryon-induced cores’ (hereafter BICs) have not been able to fully reproduce the observed diversity of dwarf galaxy rotation curves (Santos-Santos et al. 2020). This raises the question of whether baryonic outflows driven by SN feedback are truly responsible for the observed diversity.

Alternative solutions include possible systematic errors (or underestimated uncertainties) in the measurements of dwarf galaxy rotation curves. There are several known sources of error (such as the failure to account for non-circular motions) which have been discussed and studied using simulations and mock observations (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). For example Roper et al. (2023) show, using variations of the EAGLE model, that simulated dwarfs with similar circular velocity curves can actually show a diversity of ‘observed’ rotation curves because of the presence of non-circular motions.

Another potential solution involves modifications to the LCDM model, such as invoking self-interacting DM (SIDM; Spergel & Steinhardt 2000; Yoshida et al. 2000; Vogelsberger, Zavala & Loeb 2012; Rocha et al. 2013; Tulin & Yu 2018) or ‘fuzzy’ DM models. SIDM, in particular, has been repeatedly championed as a viable solution of the core-cusp problem (Burger & Zavala 2019; Burger et al. 2022).

The ability of either BICs or SIDM to explain the diversity of rotation curves has been studied in detail by Santos-Santos et al. (2020). These authors show that existing simulations, such as the NIHAO series (Dutton et al. 2016, 2019; Tollet et al. 2016) and the EAGLE-CHT10 simulation (Benítez-Llambay et al. 2019), can yield large BICs, but are still unable to fully reproduce the rotation curve diversity observed in dwarf galaxies (Schaye et al. 2015; Sawala et al. 2016; Benítez-Llambay et al. 2019; Nelson et al. 2019). Models like SIDM can also produce a wide variety of rotation curve shapes, but the diversity is driven by the wide range of assumed halo concentrations rather than by the properties of SIDM. The key problem is that simulations that produce BICs are unable to reproduce

a well-defined observational trend, where dwarf galaxies with slowly rising rotation curves (normally indicative of a ‘core’) are highly baryon-dominated in the inner regions, as opposed to rapidly rising curves (indicative of a ‘cuspy’ inner profile), where baryons play a less prominent role.

In this study, we examine the formation of BICs and the diversity of galaxy rotation curves in the NEWHORIZON (NH) hydrodynamical cosmological simulation (Dubois et al. 2021) and compare these results to observational data. NH is a zoom-in of an average density region taken from the Horizon-AGN simulation (Dubois et al. 2014; Kaviraj et al. 2017), which has a volume of (142 cMpc)$^3$. NH has a maximum spatial resolution of 34 pc, with mass resolutions of $\sim 10^4$ M$_\odot$ and $\sim 10^6$ M$_\odot$ for the stars and DM, respectively.

These make NH an ideal tool to study the inner regions of galaxies over a statistically significant sample of galaxies due to the volume of the simulation. Recently, Jackson et al. (2024) have shown that the NH simulation is able to produce DM cores over a range of stellar/halo masses, but did not verify that those cores lead to a diversity of rotation curves in agreement with observed galaxies. The latter issue is the main goal of this study.

This paper is structured as follows. In Section 2, we briefly describe the NH simulation. In Section 3, we present our main findings. We summarize our conclusions in Section 4.

## 2 THE NEWHORIZON SIMULATION

The NH cosmological, hydrodynamical simulation (Dubois et al. 2021) is a high-resolution simulation of a zoomed-in region of the Horizon-AGN simulation (H-AGN hereafter; Dubois et al. 2014; Kaviraj et al. 2017). H-AGN employs the adaptive mesh refinement code RAMSES (Teyssier 2002) and utilizes a grid that evolves a 142 comoving Mpc$^3$ volume (with spatial resolution of 1 kpc). The initial conditions are generated using using MPGRAFIC (Prunet et al. 2008) using 1024$^3$ initially uniformly distributed cubic cells.

For NH, this grid is resampled at higher resolution (using 4096$^3$ initially uniformly distributed cubic cells), using the same cosmology as that used in H-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). This higher-resolution volume is a sphere which has a radius of 10 comoving Mpc which is centred on a region of average density within H-AGN. NH has a stellar mass resolution of $10^4$ M$_\odot$, DM mass resolution of $10^6$ M$_\odot$, and a maximum spatial resolution of 34 pc (compared to 1 kpc in H-AGN).$^1$ Though the simulation has been evolved down to $z = 0.18$, this study uses the $z = 0.25$ output due to the availability of data at the commencement of the project.

### 2.1 Star formation and stellar feedback

NH includes gas cooling via primordial Hydrogen and Helium, which is gradually enriched by metals produced by stellar evolution (Sutherland & Dopita 1993; Rosen & Bregman 1995). An ambient ultraviolet (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{ cm}^{-3}$ and a temperature lower than $2 \times 10^4$ K, following a Schmidt relation (Schmidt 1959; Kennicutt 1998). The star-formation efficiency is dependent on the local turbulent Mach number and the virial parameter $\alpha = 2E_k/|E_g|$, where $E_k$ is the kinetic energy of the gas and $E_g$ is the gas gravitational binding

<footnote>$^1$The gravitational force softening is equal to the local grid size.</footnote>
energy (Kimm et al. 2017). The probability of forming a star particle is then drawn at each time-step using a Poissonian sampling method, as described in Rasera & Teysier (2006).

A star particle represents a coeval stellar population, with 31 per cent of the stellar mass assumed to explode as Type II SNe 5 Myr after its birth its formation. This fraction is calculated using a Chabrier initial mass function, with upper and lower mass limits of 150 M⊙ and 0.1 M⊙ (Chabrier 2005). SN feedback takes the form of both energy and momentum, with the final radial momentum capturing the snowplough phase of the expansion (Kimm & Cen 2014). Initially the energy of each SN is $10^{51}$ erg, with a progenitor mass of 10 M⊙. Pre-heating of the ambient gas by UV radiation from young O and B stars is included by augmenting the final radial momentum from SNe following Geen et al. (2015).

## 2.2 Identification of galaxies and merger trees

For this study, we use the ADAPTAHOP algorithm to identify DM haloes (Aubert, Pichon & Colombi 2004; Tweed et al. 2009). ADAPTAHOP is efficient at removing subhaloes from host 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 on star particles (Eisenstein & Hut 1998). As HOP does not remove substructures from the main structure, which would result in star-forming clumps being removed from galaxies. Merger trees are then created for each galaxy at the final snapshot at $z = 0.25$, with an average time-step of ~15 Myr. This enables us to track the main progenitor of every galaxy with high temporal resolution.

As NH is a high resolution zoom of Horizon-AGN, higher-mass DM particles may enter the high resolution region of NH from the surrounding lower-resolution regions. With the mass differences being large, these DM particles could interact with galaxies they are passing through in unusual ways, producing spurious effects. Even though the vast majority of galaxies affected by low DM purity exist at the outer edge of the NH sphere, we exclude any galaxy/halo that contains any low-resolution particles for this study.

We also match haloes between the NH simulation and a dark-matter-only (DMO) run of the same volume; this gives us a sample of DM haloes both 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 and that the total mass of the haloes in the two runs is similar (in more than 90 per cent of cases the two agree to within a factor of 2, see Peirani et al. 2017; Jackson et al. 2024 for more details). For the purposes of this study, we define a halo virial mass as that contained within a radius where the overdensity is $200 \times$ the critical density for closure. Throughout this paper virial quantities are labelled with a ‘200’ subscript.

# 3 RESULTS

## 3.1 The stellar mass-halo mass relation

Fig. 1 shows the stellar mass, $M_*$, as a function of halo mass, measured by $V_{max}$ (the maximum circular velocity) for all central galaxies (i.e. isolated, not satellites) in NH (blue circles). We include only galaxies whose DM haloes have $M_{200} > 10^9$ M⊙ to ensure that all haloes have at least 1000 DM particles. For comparison we show the median $M_*-V_{max}$ relations from the APOSTLE (Fattahi et al. 2018) (grey line), EAGLE-CHT10 (Benítez-Llambay et al. 2019) (cyan line), and NIHAO simulations (Wang et al. 2015) (purple line).

<figure><img src="image_2.png" /><figcaption>Figure 1. The $V_{max}-M_*$ relation for central galaxies in the NH simulation (blue points). We also show the same relationship from the abundance-matching model of Moster et al. (2013) (red dashed line), as well as from the APOSTLE (grey), EAGLE-CHT10 (cyan), and NIHAO (purple) simulations. NH is fairly consistent with APOSTLE at low $V_{max}$ though for higher $V_{max}$ it tends to yield higher stellar masses at given halo mass than the other simulations.</figcaption></figure>

We also include, for reference, the results of the abundance matching study of Moster, Naab & White (2013). For further comparisons of the NH galaxies to observational results, in particular the baryonic Tully-Fisher relation see Dubois et al. (2021) fig. 20.

We find that the results for APOSTLE and NH are similar in the dwarf galaxy regime (we define the dwarf galaxy regime as roughly $M_* < 10^9$ M⊙). However, dwarfs in NIHAO and EAGLE-CHT10 are substantially less massive, at given halo mass. At higher stellar masses, NH galaxies are more massive than in the other simulations, and agree better with the abundance-matching results.

## 3.2 Rotation curve shapes and baryonic effects

A simple characterization of the shape of the rotation curve of a galaxy is provided by the ratio between the maximum circular velocity of a system and the circular velocity at a suitably defined inner ‘fiducial’ radius, $V_{fid} = V_c(r_{fid})$. For this study we use the fiducial radius defined by Santos-Santos et al. (2020):

$$r_{fid} = 2(V_{max}/70 \text{ km s}^{-1}) \text{ kpc.} \qquad (1)$$

For comparison, a DMO LCDM halo with $70 \text{ km s}^{-1}$ virial velocity reaches its peak circular velocity at $r_{max} \approx 12 \text{ kpc}$, so the ratio, $\eta_{rot} = V_{fid}/V_{max}$, provides a good measure of how rapidly the circular velocity rises in the inner regions to reach its maximum. Galaxies with $\eta_{rot} \approx 1$ have rapidly rising rotation curves, and those with ratios well below unity have slowly rising rotation curves.

More importantly, the self-similar nature of LCDM haloes implies that, in the absence of baryons, the $\eta_{rot}$ ratio is expected to be roughly constant for all haloes, regardless of mass ($\eta_{rot} \sim 0.65$). This implies that values of $\eta_{rot}$ significantly different from 0.65 signal galaxies where baryons have had significant influence shaping the rotation curve.
<figure><img src="image_3.png" /></figure>

**Figure 2.** Circular velocity profiles for four example galaxies with log M<sub>☉</sub> = 7.04, 8.46, 9.37, and 10.74 M<sub>☉</sub>. In each example we show the velocity contributions corresponding to the total mass (thick black dot-dashed), DM (thick blue), DMO run (thick orange), gas (purple), and stars (red). We also include a best-fitting NFW profile for the DMO run (thin black dashed line). The black arrow indicates the fiducial radius, r<sub>fid</sub>, for each halo. We also indicate where the profiles go inside the convergence radius (Power et al. 2003; Schaller et al. 2015) by showing the profile as a dotted line. From these examples DM ‘cores’ (i.e. reductions in DM relative to the DMO case) are clearly visible in two galaxies (2 and 3). For these cases, we fit the cored-NFW parametrization from Zavala, Vogelsberger & Walker (2013), equation (2), and show the result as the brown dashed line. Additionally, we can see that in the lowest mass halo the DM component is higher than the DMO run due to the presence of significant amounts of gas, whereas in the most massive example this occurs due to the presence of stars.

Values of $\eta_{\text{rot}}$ well below 0.65 indicate a deficit of mass in the inner regions relative to the DMO run (a clear indication of a ‘core’ in the DM inner mass profile), whereas values higher than 0.65 indicate cases where the total mass inside the fiducial radius has increased, either as a result of the assembly of the baryonic component of the galaxy and/or of the consequent halo contraction. Although a DM core may still be present in the latter case, it is not straightforward to detect because the inner DM deficit in these cases could have been compensated by the luminous component.

A better measure of the effect of baryons on the halo is provided by how much the DM mass inside $r_{\text{fid}}$ has been changed by the assembly of the galaxy, compared with a DMO run. This can be measured by the ratio, $v_{\text{fid}}$, between $V_{\text{fid}}$ in the DMO case and the DM contribution to $V_{\text{fid}}$ in the hydrodynamical simulation; i.e. $v_{\text{fid}} = V_{\text{fid,DM}}/V_{\text{fid,DMO}}$. Values of $v_{\text{fid}}$ substantially lower than unity indicate the formation of a core in the DM; values above unity signal halo contraction due to

the gravitational effect of the baryonic component. We shall hereafter identify galaxies with a ‘core’ as those with $v_{\text{fid}} < 0.75$.

### 3.3 Rotation curves of galaxies in NEWHORIZON

We begin by examining the rotation curves of four NH galaxies in Fig. 2. These examples correspond to systems spanning a wide range of stellar mass (log M<sub>☉</sub>/M<sub>☉</sub> = 7.04, 8.46, 9.37, and 10.74). In each panel we show the DMO circular velocity curve (orange), fitted with an NFW profile (black dashed line), as well as the contributions to the circular velocity of the DM (blue), gas (purple), and stellar (red) components in the simulation including baryons.

The thick black dot-dashed line shows the total circular velocity in the simulation with baryons. We indicate $r_{\text{fid}}$ in each case using a black arrow. The direction of the arrow (upward/downward) indicates whether the DM enclosed within $r_{\text{fid}}$ has increased/decreased, relative
<figure><img src="image_4.png" /><figcaption>**Figure 3.** The reduction in DM mass inside the fiducial radius (as measured by $v_{\text{fid}} = V_{\text{fid,DM}}/V_{\text{fid,DMO}}$) as a function of galaxy stellar mass, $M_*$ (left), and of galaxy formation 'efficiency', $M_*/M_{\text{halo}}$ (right). We indicate our definition of galaxies with 'cores' as those below the black dashed line indicating $v_{\text{fid}} = 0.75$. The results for APOSTLE (grey), NIHAO (purple), and EAGLE-CHT10 (cyan) galaxies are also shown. In all simulations there is a clear trend with increasing stellar mass, with low-mass galaxies having values of $v_{\text{fid}} \approx 1$, then gradually decreasing as cores for, peaking around $M_* \sim 10^9 M_{\odot}$ before increasing again. The same trend is seen in all simulations except APOSTLE, where cores do not form. Although similar trends are seen for $v_{\text{fid}}$ as a function of galaxy formation 'efficiency', $M_*/M_{\text{halo}}$, there is no agreement between simulations using this parameter. This suggests that $M_*/M_{\text{halo}}$ is not the best metric for identifying over which stellar mass range galaxies form cores. Red circles highlight the galaxies whose rotation curves are shown in Fig. 2.</figcaption></figure>

to the DMO run. We indicate the Power et al. (2003) convergence radius for each halo, where the curves transition from dotted to solid lines.

These example circular velocity profiles show how diverse rotation curves are when compared to those in the DMO simulation, which are all reasonably well approximated by NFW profiles. Indeed, some systems have slowly rising circular velocity profiles consistent with cores, whereas others have steeply rising curves indicative of a highly concentrated baryonic component and/or a DM cusp.

In the lowest-mass system (top-left panel in Fig. 2), gas contributes far more than stars to the central mass budget, and contracts the halo, yielding a fairly flat circular velocity curve with $\eta_{\text{rot}} \approx 0.9$ and $v_{\text{fid}} \sim 1.4$. In the highest mass galaxy (bottom-right), stars dominate the inner regions, also contracting the halo in the process ($\eta_{\text{rot}} \approx 0.95$, $v_{\text{fid}} \sim 1.4$).

On the other hand, the two intermediate-mass systems show evidence of an inner mass deficit compared with the DMO run, indicating the presence of a core. The brown dashed lines in Fig. 2 show fits to the inner DM profile using the density profile proposed by Zavala et al. (2013):

$$ \rho(r) = \frac{\rho_0 r_s^3}{(r + r_c)(r^2 + r_s^2)}. \qquad (2) $$

The four examples in Fig. 2 show clearly that the assembly of the galaxy can alter substantially the circular velocity profiles of NH galaxies, yielding a wide diversity of rotation curve shapes.

## 3.4 Core formation and galaxy stellar mass

Which galaxies develop baryon-induced DM cores? The left-hand panel of Fig. 3 shows, as a function of galaxy stellar mass, the change in the DM content within $r_{\text{fid}}$, as measured by the parameter $v_{\text{fid}}$

introduced above, for all NH galaxies. The right-hand panel shows $v_{\text{fid}}$ as a function of the galaxy formation 'efficiency', defined as the ratio of galaxy stellar mass to halo virial mass.

The results for NH galaxies are shown by the navy blue circles (the example galaxies shown in Fig. 2 are highlighted with red circles). The 'V'-shaped distribution of blue symbols highlight the fact that cores do not form in NH at either the low-mass or high-mass regime. Indeed, cores are only seen in a well-defined range of stellar mass, $10^{7.7} < M_*/M_{\odot} < 10^{9.5}$, with the largest cores in galaxies with $M_* \sim 10^9 M_{\odot}$.

A similar trend has been previously reported in other simulations, like NIHAO (shown in purple) and EAGLE-CHT10 (cyan). Note that APOSTLE galaxies (in grey) show no evidence of a core (Bose et al. 2019), a result that has been traced to the low gas density threshold adopted in those simulations, which prevents the gas from becoming gravitationally dominant over the DM in the dwarf galaxy regime (Benítez-Llambay et al. 2019).

The agreement between the $M_*$ versus $v_{\text{fid}}$ trends in NH, NIHAO, and EAGLE-CHT10 is at face value somewhat surprising, because it had been argued that baryon-induced cores should depend mainly on galaxy formation efficiency (i.e. $M_*/M_{\text{halo}}$), and not on $M_*$ alone (Di Cintio et al. 2014; Tollet et al. 2016; Dutton et al. 2019). The right-hand panel of Fig. 3 shows, however, that this is not the case. The $v_{\text{fid}}$ trend with $M_*/M_{\text{halo}}$ varies significantly from simulation to simulation, but the correlation with $M_*$ is quite similar for all three sets of runs.

What causes this? One clue is provided by Fig. 4, where we show the enclosed DM profiles of all NH galaxies. The three colours distinguish between systems with BICs (i.e. $v_{\text{fid}} < 0.75$, shown in blue) from those without an inner DM deficit ($v_{\text{fid}} > 0.75$, shown in green for low-mass systems and orange for high-mass ones).

These systems span almost four decades in virial mass, from roughly $10^9$ to $10^{13} M_{\odot}$. However, the range in enclosed dark mass in
<figure><img src="image_5.png" /></figure>

**Figure 4.** Lines show cumulative DM profiles for all haloes in the NH simulation. We colour-code haloes by $v_{\text{fid}}$; blue indicates galaxies with a core ($v_{\text{fid}} < 0.75$), orange and green denote galaxies with high mass (log $M_*$ > 9.5) or low mass (log $M_*$ < 8) that do not form a core (i.e. $v_{\text{fid}} > 0.75$). The coloured symbols show the effective radius versus galaxy stellar mass (coloured crosses) for each galaxy. Note that although the haloes span a wide range of virial mass (nearly four decades) the dark mass in the regions inhabited by stars (i.e. the inner 2–3 kpc) varies much less. In particular, within 2 kpc, the dark mass of systems with a core peaks at log M = 9.0 with an rms of only 0.5 dex.

the region where most of the stars reside (the ‘effective’ or half-mass radius of the stellar component, indicated by the coloured symbols in Fig. 4) is much narrower, less than 1 dex at 2–4 kpc, the typical effective radius of NH galaxies.

At such radii most haloes have enclosed dark masses of order ~$10^9 M_{\odot}$, in particular those of systems with cores (see blue lines in Fig. 4). At $r = 2$ kpc, for example the distribution of total enclosed mass in all haloes with baryon-induced cores is quite narrow: $1 \pm 0.5 \times 10^9 M_{\odot}$. In addition, the gas-to-stellar mass ratio of galaxies (which decreases steadily with increasing $M_*$) approaches unity at $M_* \sim 10^9 M_{\odot}$ (see Fig. 5).

This suggests a relatively straightforward interpretation for why cores do not form in galaxies with $M_* > 10^9 M_{\odot}$: in such systems the total mass in stars within $r_{\text{fid}}$ begins to exceed that of DM, and also that of the gaseous component. This means that the baryonic component becomes self-gravitating, and stars begin to dominate the potential. Since core formation requires fluctuations in the potential driven by the inflow and outflow of gas, such fluctuations become less and less important in more massive galaxies. Core formation thus ceases in galaxies with $M_* > 10^9 M_{\odot}$, and in more massive systems the DM inside $r_{\text{fid}}$ increases gradually as the halo contracts due to the deepening of the potential.

To show that this explanation does not only hold for NH, we show in Fig. 6 (left panel) the ratio between the stellar mass inside a fixed radius of 2 kpc (which corresponds roughly to the half mass radius of the stellar component) and the total mass inside the same radius in the DMO simulation. Blue symbols identify systems with cores in NH and in NIHAO.

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

**Figure 5.** Stellar mass versus gas mass within 15 per cent of the virial radius for haloes in NH. Haloes are colour-coded by their $v_{\text{fid}}$ value. As in Fig. 4, blue indicates systems with a core ($v_{\text{fid}} < 0.75$), orange and green denote galaxies with high mass (log $M_* > 9.5$) or low mass (log $M_* < 8$) that do not form a core (i.e. $v_{\text{fid}} > 0.75$). We include a 1:1 line to guide the eye. The gas-to-stellar mass ratio decreases with increasing stellar mass and begins to approach unity at $M_* \sim 10^9 M_{\odot}$. At larger stellar masses the gas becomes subdominant, and stars begin to dominate over the gas and the DM in the inner regions, preventing the formation of cores in massive systems.

Note that in neither simulation cores form when the total mass in stars exceeds the initial DM mass inside that radius. This holds despite the fact that the galaxy formation efficiency differs markedly in both simulations (see Fig. 1). Indeed, the right-hand panel of that figure shows the same ratio as a function of halo virial mass, and shows that, at fixed stellar mass, NIHAO galaxies inhabit substantially more massive haloes than NH.

The nice overlap between NIHAO and NH results in the left-hand panel stems from the fact that even haloes with substantially different virial mass will have similar DM mass within 2 kpc. We conclude that galaxy stellar mass is the key indicator of whether a core can form or not, and argue that cores (meant as a reduction of the inner DM content, rather than as a change in the slope of the density profile) should not, in general, be present in galaxies with $M_* > 10^9 M_{\odot}$.

## 3.5 Core sizes

We examine next the size of the DM cores, which we show, as a function of galaxy stellar mass, in Fig. 7. We begin by noting that core sizes depend sensitively on how a core is defined. For the definition we adopted above (see equation 2), we see that core sizes scale approximately with the half-mass radius of the stars, but are a factor of 2 or so smaller. This is not unexpected, since it is the formation of the stellar component that leads to core formation. (Note that we show results *only* for systems deemed to have a core; i.e. $V_{\text{fid,DM}}/V_{\text{fid,DMO}} < 0.75$).

A more robust measure of a core could be defined by locating the radius where the DM circular velocity profile starts to drop below the DMO profile by more than 25 per cent. These radii are marked in Fig. 2 for the two intermediate galaxies by a thin vertical segment,
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Left: galaxy stellar mass versus stellar-to-DM mass within 2 kpc for galaxies in both the NH and NIHAO simulations. Galaxies with a core are shown by blue symbols, dots for NH and crosses for NIHAO. High mass cusps (log M* > 9.5) are shown as orange dots for NH and crosses for NIHAO and low mass cusps (log M* < 8) are green. We see a clear trend in both simulations that as the galaxy stellar mass increases the stellar-to-halo mass ratio within 2 kpc also increases, approaching unity at M* ~ 10⁹ M⊙. At larger masses the stellar mass component becomes comparable or larger than the DM and we no longer see cores. Right: same as left panel, but as a function of halo virial mass. Colour-coding and symbols are the same as in the left panel. Different simulations form cores in haloes of different virial mass. This supports the interpretation that it is the galaxy stellar mass that dominates the core formation process, and not the galaxy-to-halo mass ratio.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Blue circles show the effective radius versus galaxy stellar mass for all 'cored' galaxies (i.e. those with v<sub>fid</sub> < 0.75) in NH. Red circles show core radii estimated by fitting equation (3) in Zavala et al. (2013) to the DM profile. Note the large scatter when using this core radius estimator. A more robust measure is provided by the radius where V<sub>DM</sub>/V<sub>DMO</sub> = 0.75, which is shown by the open red circles. This estimate of the core radii correlates well with the galaxy effective radius, suggesting that the effects of baryons on the profile are largely limited to the region inhabited by the galaxy stellar component.</figcaption></figure>

and are shown in Fig. 7 by open circles. This shows again that core
sizes scale tightly with stellar half-mass radius, a relationship that
could, in principle, be tested by observations.

3.6 The diversity of rotation curves in NEWHORIZON

The detailed information about the DM profiles before and after
galaxy formation used in earlier sections to evaluate the formation of
baryon-induced cores cannot be compared directly with observations.
For the latter, the main observable is the total circular velocity profile
of baryons plus DM. Although in principle it is possible to derive the
contribution of the DM component by subtracting the contribution
of the baryons, this procedure requires extra assumptions (such as a
mass-to-light ratio for the stars), which introduce further uncertainty.
That is why Oman et al. (2015) and Santos-Santos et al. (2020) argued
that one should adopt a simple characterization of the shape of the
circular velocity curve (like the η<sub>rot</sub> parameter introduced above),
and to compare that directly with observations to gauge the presence
of cores.

As described in that work, two observations are particularly
intriguing. One is that at fixed stellar mass/maximum rotation
velocity there is a wide spread of rotation curve shapes (i.e. a wide
range of η<sub>rot</sub>). Another is that, for dwarf galaxies (more precisely,
M* < 10<sup>9.5</sup> M<sub>☉</sub>), rotation curve shape correlates strongly with the
gravitational importance of baryons inside the fiducial radius, r<sub>fid</sub>.

The latter can be measured using the parameter η<sub>bar</sub> = (V<sub>bar,fid</sub>/V<sub>fid</sub>)<sup>2</sup>, where V<sub>bar,fid</sub> measures the baryon contribution to the circular velocity at r<sub>fid</sub>. As discussed by Santos-Santos et al. (2020), rapidly rising rotation curves occur mainly in galaxies where baryons are unimportant within r<sub>fid</sub>, whereas the most prominent ‘cores’ (i.e. low values of η<sub>rot</sub>) are seen in galaxies where baryons are quite important in the inner regions. This η<sub>rot</sub>–η<sub>bar</sub> trend has proved quite
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Rotation curves of dwarf galaxies with $75 < V_{\text{max}}/\text{km s}^{-1} < 95$, chosen to illustrate the diversity of rotation curve shapes seen (grey curves). Coloured lines indicate observed dwarf irregular galaxies and grey lines denote galaxies from NH with comparable $V_{\text{max}}$. We include an representative example DMO halo with $V_{\text{max}} = 80 \text{ km s}^{-1}$ (thick black curve) for reference. It is clear that NH is able to produce a variety of rotation curve shapes for galaxies in a narrow range of halo masses, similar to what is seen in observations.</figcaption></figure>

difficult to reproduce in simulations where the diversity of dwarf galaxy rotation curves is due to the presence of BICs (see Santos-Santos et al. 2020, for a more detailed discussion).

We begin by analysing the diversity of rotation curves over a narrow range of circular velocities in Fig. 8, where we compare four observed galaxies with NH systems with $75 < V_{\text{max}}/\text{km s}^{-1} < 95$. The thick solid black line shows the median rotation curve expected for DMO haloes in this mass range. Clearly, observed galaxies span a wide range of rotation curve shapes, some rising more slowly, and others more rapidly than expected from DM alone. This is likely driven by the influence of baryons; indeed, the NH galaxies (grey curves) show a similar diversity of rotation curve shapes as observed, with even one simulated galaxy close to matching IC2574, the dwarf galaxy with the most prominent core known to date.

## 3.7 Rotation curve shape versus central baryon dominance

The comparison with the observed correlation between rotation curve shape and inner baryon dominance is more nuanced. The circles in the right-hand panel of Fig. 9 shows the $\eta_{\text{bar}}$ versus $\eta_{\text{rot}}$ correlation for observed galaxies, coloured by $M_*$. The galaxy sample shown here is the same as that compiled by Santos-Santos et al. (2020), and includes rotation curves from the SPARC survey (Lelli et al. 2016), THINGS HI Nearby Galaxy Survey (de Blok et al. 2008), the Local Irregulars That Trace Luminosity Extremes, The HI Nearby Galaxy Survey (Oh et al. 2015), and other work from Adams et al. (2014) and Relatores et al. (2019).

The correlation alluded to above is apparent for (dwarf) galaxies with $M_* < 10^{9.5} M_{\odot}$. These galaxies (shown in orange and red) trace a wide swath, from the top left to the bottom right of the panel. Massive galaxies (shown in yellow), on the other hand, sit at the top

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Filled circles show the rotation curve shape parameter ($\eta_{\text{rot}} - V_{\text{fid}}/V_{\text{max}}$) versus the gravitational importance of the baryonic component at the fiducial radius ($\eta_{\text{bar}}$) for galaxies in NH (left) and observations (right), colour-coded by galaxy stellar mass. Galaxies in NH spread over the same range of parameter space seen in observations. Low $\eta_{\text{bar}}$ galaxies typically have fast-rising rotation curves, consistent with a cuspy NFW profile (or steeper). As $\eta_{\text{bar}}$ increases the diversity in rotation curve shapes also increases. We see that the most massive galaxies occupy the top right of the parameter space in both samples: these are baryon-dominated and rapidly rising rotation curve systems. For dwarfs (i.e. $M_* < 10^{9.5} M_{\odot}$), observations suggest a strong trend from top left to bottom right in this plot. This trend is not reproduced in NIHAO galaxies (blue crosses, left) or EAGLE-CHT10 (green crosses, right). NH fare much better, in comparison. See the text for further details.</figcaption></figure>
right of the panel; these are galaxies where baryons clearly dominate at the centre (or at least inside $r_{fid}$) and whose rotation curves are all sharply rising in the inner regions.

Surprisingly, dwarfs whose baryons dominate inside $r_{fid}$ have slowly rising curves, and show clear evidence of a core (bottom right of the panel) whereas baryons are unimportant in dwarfs with less compelling evidence of a core (i.e. those with $\eta_{rot} \sim 0.65$ or greater). As discussed by Santos-Santos et al. (2020), accounting for this correlation is an important challenge to any proposed explanation of the origin of the diversity of dwarf galaxy rotation curves. For example, the green crosses in the right-hand panel of Fig. 9 show the results for EAGLE-CHT10 galaxies. Although BICs do form, the simulated galaxies do not show the same $\eta_{bar}$ versus $\eta_{rot}$ correlation as observed. The same is true for NIHAO galaxies, which are shown by the blue crosses in the left-hand panel of Fig. 9. Neither simulation has galaxies with sharply rising rotation curves where baryons are unimportant, and few or no galaxies with obvious cores where baryons dominate the inner regions.

The left panel of Fig. 9 also shows the results for NH galaxies (circles coloured by stellar mass). NH galaxies fare better, with a clear top left to bottom right correlation for dwarfs which traces closely the observed one. Upon closer inspection, however, some differences become clear. The most apparent is that the correlation in NH is due to the fact that many dwarfs around $10^9 M_{\odot}$ have cores (i.e. low $\eta_{rot}$) whereas most galaxies at low mass do not. In other words, the trend observed between $\eta_{rot}$ and $\eta_{bar}$ reflects a mass trend where, among simulated dwarfs, those at the faint end have rising curves and are not baryon dominated, and those at the bright end have cores and are baryon-dominated within $r_{fid}$. Such mass-driven trend does not seem to be present in the observational sample, where dwarfs across the whole $10^5-10^9 M_{\odot}$ stellar mass range appear to trace the same $\eta_{rot}-\eta_{bar}$ correlation.

These differences may also be seen in Fig. 10, where we show separately $\eta_{rot}$ and $\eta_{bar}$ as a function of galaxy stellar mass. Blue circles show NH galaxies, grey ones indicate the observations.

This figure makes clear that the range of stellar masses where cores are apparent (i.e. systems below the horizontal shadowed area) is similar in observations and in NH. However, it also shows that NH shows a strong correlation between inner baryon dominance and stellar mass. The trend among observed dwarfs is not as strong, if there is a trend at all. The lack of NH galaxies in the range $10^{8.5} < M_{*}/M_{\odot} < 10^{10.5}$ where baryons are unimportant is one of the most obvious differences between simulated and observed galaxies in this sample.

Given this difference, we are hesitant to conclude that the observed diversity of rotation curves has been successfully reproduced by BICs in cosmological hydrodynamical simulations like the ones we consider here.

## 3.8 Interpretation: SNe diffusion and adiabatic contraction

Consider the following simple toy model which should capture qualitatively the impact of SNe feedback on the orbital diffusion of DM cusps (see also Kipper et al. 2024). Given the following approximations (i) the total potential is assumed to be NFW-like $-V_0^2 \log(1+r/a)a/r \approx V_0^2(r/2a - 1)$ where a is the scale length of the halo; (ii) the DM particles are assumed to move on radial orbits; (iii) the gas density distribution times the star formation and SN efficiency converted into number of bubbles is assumed to follow a given cored profile $\Sigma/(1+r/b)$, with b the scale length of this distribution; (iv) the bubble power spectrum in k, $\omega$ is ap-

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Galaxy stellar mass versus the gravitational importance of the baryonic component at the fiducial radius ($\eta_{bar}$, top panel) and versus the rotation curve shape parameter ($\eta_{rot}$) for galaxies in NH (blue points) and observations (grey points). We see a clear trend with $\eta_{bar}$ and stellar mass for NH with increasing stellar mass resulting in increasing $\eta_{bar}$. There is some evidence of this in the observations, but it is systemically at higher masses than in NH.</figcaption></figure>

proximated by a shifted Lorenzian, $1/(1+(k R_{\text{bub}} - 2)^2 + (\omega T_{\text{bub}})^2)$, where $R_{\text{bub}}$ and $T_{\text{bub}}$ are respectively the typical bubble size and its lifespan. These assumptions are sufficient to compute analytically the (orbit averaged) orbital diffusion coefficient, $D_E$, as a function of $R_{\text{bub}}$, $T_{\text{bub}}$, a, b, $\Sigma$, $V_0$ which reads at maximum

$$D_{\text{max}} = \frac{10\pi}{9} \Sigma \sqrt{\frac{b}{a} \frac{R_{\text{bub}} V_0^2 T_{\text{bub}}^2}{a^2}} \quad (3)$$

From this expression, we conclude that at low mass, for gas dominated galaxies ($\Sigma \sim 1, b \sim a$), diffusion can be significant, so orbital diffusion may proceed to replace the cusp by a core while increasing the eccentricity of DM particles within b, while at higher mass, the galaxy becomes relatively gas poorer and the gas is distributed more compactly ($\Sigma \ll 1, b \ll a$), so diffusion is less efficient and adiabatic contraction proceeds (Young 1980). The former process is most efficient if the orbital period of the DM particle, $a/V_0$ match at apocentre the bubble lifespan $T_{\text{bub}}$ and the
bubble size, $R_{\text{hub}}$ and bubble distribution scale length $a$ match the NFW scale length $a$.

## 4 SUMMARY AND CONCLUSIONS

We have used the NH cosmological simulation, along with its DMO counterpart, to examine how the assembly of the baryonic component of a galaxy affects the inner DM distribution, and its effects on the shape of the galaxy’s rotation curve. In particular, we have analysed how these effects depend on galaxy stellar mass, and whether NH galaxies can successfully reproduce the observed diversity of dwarf galaxy rotation curves. Our main conclusions may be summarized as follows.

(i) In contrast with the DMO case, the rotation curves of NH galaxies show a wide diversity of rotation curves. This diversity is driven by the assembly of the baryonic component of the galaxy, and shows a strong dependence on galaxy stellar mass.

(ii) Measuring the effect of baryons by the change in the inner DM content of a galaxy (compared with the DMO case, within a fiducial radius, $r_{\text{fid}}$), we find, as in earlier work, that DM ‘cores’ only form over a limited range of galaxy stellar mass ($10^{7.5} < M_{\star}/M_{\odot} < 10^{9.5}$).

(iii) The $M_{\star}$ range over which cores are produced by baryons is similar to that of other cosmological hydrodynamical simulations, such as NIHAO and EAGLE-CHT10, despite the fact that the galaxy formation efficiency varies widely from simulation to simulation.

(iv) Below the lower limit of the range, galaxies with $M_{\star} < 10^{7.5} M_{\odot}$ form too few stars to drive the strong baryonic inflows/outflows needed to create a core.

(v) The upper limit of the range seems to be set by the total amount of DM in the region (i.e. inner 2–3 kpc) where stars form. Indeed, haloes with baryon-induced cores all have $1 \pm 0.5 \times 10^9 M_{\odot}$ inside their inner 2 kpc in the DMO case. Galaxies more massive than this form enough stars for them to become self-gravitating, deepening substantially the gravitational potential, and leading to a contraction of the halo rather than the formation of a core.

(vi) We find that galaxy stellar mass is the key indicator of whether a core can form or not, rather than the ratio between stellar and virial mass. This stems from the fact that even haloes with substantially different virial mass will have similar DM mass within 2 kpc. It is the ratio between masses in this central region which determines whether a galaxy will have a core or a cusp.

(vii) The size of the core (when present) is well correlated with the size of the stellar component. We emphasize that ‘core radius’ estimates depend sensitively on the definition of a core, so caution should be exercised when comparing different simulations.

(viii) NH galaxies come closer than other cosmological hydrodynamical simulations to matching the observed relation between inner baryon dominance ($\eta_{\text{bar}}$) and rotation curve shape ($\eta_{\text{rot}}$). However, the correlation in NH between $\eta_{\text{bar}}$ and galaxy stellar mass is not seen in the observations.

The caveat in the latter item prevents us from concluding that the observed diversity of dwarf galaxy rotation curves can be fully explained by baryon-induced cores in simulations like NH. In particular, observed galaxies differ from NH in the galaxy stellar mass dependence of the inner baryon dominance. The discrepancy, however, is not large, and it is possible that further iterations of this (and other) simulation series may lead to galaxies as diverse as observed dwarfs in terms not only of their rotation curves, but also in the role of baryons in the inner regions.

Finally, we provided a back of the envelope toy model which qualitatively motivates the observed transition.

## 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). 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 was 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 No. 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 NEWHORIZON 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.

## DATA AVAILABILITY

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

## REFERENCES

Adams J. J. et al., 2014, ApJ, 789, 63

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

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

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

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

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, 126

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, in Corbelli E., Palle F., eds, The Initial Mass Function: From Salpeter 1955 to 2005. Springer, Dordrecht, 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, 478, 1094

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

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. J., 2008, AJ, 136, 2648

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, L1

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

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

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

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

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

Dutton A. A., Macciò A. V., Buck T., Dixon K. L., Blank M., Obreja A., 2019, MNRAS, 486, 655

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

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

Fattahi A., Navarro J. F., Frenk C. S., Oman K. A., Sawala T., Schaller M., 2018, MNRAS, 476, 3816

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

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

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, 948

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

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

Harvey D., Revaz Y., Robertson A., Hausammann L., 2018, MNRAS, 481, L89

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

Jackson R. A. et al., 2024, MNRAS, 528, 1655

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

Kennicutt R. C. Jr, 1998, ApJ, 498, 541

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

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

Kipper R., Vurm I., Niovi Triantafyllaki A., Tenjes P., Tempel E., 2024, A&A, 693, A206

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, 920

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

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

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

Moore B., 1994, Nature, 370, 629

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

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, 563

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

Navarro J. F. et al., 2010, MNRAS, 402, 21

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

Oh S.-H., de Blok W. J. G., Brinks E., Walter F., Kennicutt Robert C. J., 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

Peirani S., de Freitas Pacheco J. A., 2008, Phys. Rev. D, 77, 064023

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

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

Power C., Navarro J. F., Jenkins A., Frenk C. S., White S. D. M., Springel V., Stadel J., Quinn T., 2003, MNRAS, 338, 14

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

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

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

Relatores N. C. et al., 2019, ApJ, 887, 94

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, 634

Sales L. V., Wetzel A., Fattahi A., 2022, Nat. Astron, 6, p. 897–910

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, 243

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, 123

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

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

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

Young P., 1980, ApJ, 242, 1232

Zavala J., Vogelsberger M., Walker M. G., 2013, MNRAS, 431, L20

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