The diversity of rotation curves of simulated galaxies with cusps and cores

We use $\Lambda$CDM cosmological hydrodynamical simulations to explore the kinematics of gaseous discs in late-type dwarf galaxies. We create high-resolution 21-cm 'observations' of simulated dwarfs produced in two variations of the EAGLE galaxy formation model: one where supernova-driven gas flows redistribute dark matter and form constant-density central 'cores', and another where the central 'cusps' survive intact. We 'observe' each galaxy along multiple sight lines and derive a rotation curve for each observation using a conventional tilted-ring approach to model the gas kinematics. We find that the modelling process introduces systematic discrepancies between the recovered rotation curve and the actual circular velocity curve driven primarily by (i) non-circular gas orbits within the discs; (ii) the finite thickness of gaseous discs, which leads to overlap of different radii in projection; and (iii) departures from dynamical equilibrium. Dwarfs with dark matter cusps often appear to have a core, whilst the inverse error is less common. These effects naturally reproduce an observed trend which other models struggle to explain: late-type dwarfs with more steeply-rising rotation curves appear to be dark matter-dominated in the inner regions, whereas the opposite seems to hold in galaxies with core-like rotation curves. We conclude that if similar effects affect the rotation curves of observed dwarfs, a late-type dwarf population in which all galaxies have sizeable dark matter cores is most likely incompatible with current measurements.


INTRODUCTION
The cusp-core problem has been an enduring one in cosmology, and represents an important challenge to our current understanding of dark matter (DM) and galaxy formation in a ΛCDM universe.The structure of cold DM halos has be studied extensively using cosmological N-body simulations, and is now well understood (Navarro, Frenk & White 1996b, 1997;Wang et al. 2020).Such simulations typically produce haloes with density profiles rising to an asymptotic 'cusp' in the centres of all haloes expected to host galaxies.In contrast, many observed dwarf galaxies have rotation curves implying a constant central 'core' in DM density (Flores & Primack 1994;Moore 1994, and see de Blok 2010 for a review).A resolution of this apparent overprediction of the central DM density in dwarfs remains elusive.
Many possible explanations for this conflict have been proposed.Oman et al. (2015) argued that in the context of a DM cosmology, the possibilities fall into three categories: (i) the DM differs from the fiducial hypothesis of cold, collisionless particles; (ii) DM cusps are transformed into cores, e.g. through gravitational coupling to ★ E-mail: finn.roper@ed.ac.uk, kyle.a.oman@durham.ac.uk violent gas motions driven by supernova (SN) explosions; (iii) the DM distributions within galaxies have been incorrectly inferred from kinematic data.Of course, combinations of effects in more than one of these categories are also possible.
It has been shown that allowing for departures from the fundamental assumption in ΛCDM cosmology that DM is a cold, collisionless fluid can affect the internal kinematics of dwarf galaxies.Selfinteracting dark matter (SIDM) models have been proposed (Spergel & Steinhardt 2000), where scattering interactions between the cold DM particles are introduced.Interaction cross-sections per unit mass of / ∼ 1 cm 2 g −1 lead to the thermalisation of the centres of DM haloes needed to produce cores in dwarfs, whilst leaving the successes of CDM on larger scales essentially intact (see e.g.Rocha et al. 2013;Elbert et al. 2015;Ren et al. 2019).It has been argued that SIDM models may explain the broad diversity in dwarf rotation curve shapes, although additional effects due to galaxy formation physics complicate this picture (Creasey et al. 2017;Kaplinghat et al. 2020).
DM cores may also be created through 'baryonic' processes.Navarro, Eke & Frenk (1996a) showed that a sudden loss of mass from the inner halo, which may occur when gas is expelled by supernova explosions, can redistribute DM away from the halo centre, creating a core.This process is reversible as baryons may reaccrete over time and recreate a cusp.However, if these blowouts are sufficiently frequent and numerous, the end result may be a shallower central density profile (Pontzen & Governato 2012).For significant baryon-induced core creation (BICC) to occur, the baryons must be sufficiently gravitationally dominant to cause the halo to contract before the subsequent expulsion of gas from the centre (Benítez-Llambay et al. 2019).These results suggest that DM cores may form in dwarf galaxies with sufficiently variable and energetic central star formation (SF) feedback (Faucher-Giguère 2018), and this has been seen in multiple cosmological galaxy formation simulations (e.g.Di Cintio et al. 2014;Chan et al. 2015;Tollet et al. 2016;Verbeke et al. 2017;Hopkins et al. 2018;Benítez-Llambay et al. 2019;Jahn et al. 2023).
In principle, 'beam smearing' (see Swaters et al. 2009, and references therein) or ambiguity in the mass-to-light ratios of galaxies can lead to an incorrect inference of the DM distribution, giving the appearance of a core where none exists.However, cores are still inferred to exist from high-resolution 21-cm radio observations of dwarf galaxies, such as those taken for the THINGS (Walter et al. 2008) and LITTLE THINGS (Hunter et al. 2012) surveys, where concerns around beam smearing are greatly alleviated.The mass-to-light ratios of galaxies in these surveys were derived from an empirical relation based on the optical colours of galaxies (Bell & de Jong 2001;Bruzual & Charlot 2003), which has been robustly validated using multiple techniques (see e.g.Oh et al. 2011Oh et al. , 2015)).
Nevertheless, the possibility that observed cores are a mere symptom of systematic issues in observation or modelling has not been definitively ruled out.The shapes of DM haloes have been predicted to be triaxial in shape since some of the earliest CDM simulations (Davis et al. 1985;Frenk et al. 1988).This generates an aspherical gravitational potential, with the central regions having the greatest asphericity (Hayashi, Navarro & Springel 2007).A disc galaxy at the centre of such a halo, typically aligned on a plane close to that defined by the intermediate and major axes of the halo (Hayashi & Navarro 2006), therefore inhabits a non-axisymmetric potential.Non-circular motions (NCMs) in the disc are induced as gas orbits are elongated in the potential, manifesting as mainly bisymmetric ( = 2 harmonic) fluctuations in azimuthal velocity as a function of projection angle (Oman et al. 2019, hereafter O19; see also Marasco et al. 2018).Even very slight asphericity in the potential can result in significant NCMs in the disc (e.g.Hayashi & Navarro 2006).Although harmonic decomposition of the velocity fields of observed dwarfs suggest that the amplitudes of NCMs are relatively small (Trachternach et al. 2008), these could easily have been underestimated during model fitting.Often, variations in velocity associated with NCMs are absorbed by other model parameters that are degenerate with NCMs (see Sec. 3.2).This may lead to surprisingly large errors in the measurement of rotation curves, or, in other words, cause large discrepancies between the inferred and true circular velocity at the same radius (O19).
The weak, bisymmetric, bar-like velocity patterns discussed here are distinct from the far stronger patterns characteristic of barred spiral galaxies (e.g.Machado & Athanassoula 2010).Those patterns form through an instability occurring in massive, self-gravitating discs (Ostriker & Peebles 1973;Algorry et al. 2017), but the weaker bar-like perturbations we are concerned with are instead suppressed in such systems (see Marasco et al. 2018) because massive discs 'sphericalise' their initially triaxial host haloes as they form (see Abadi et al. 2010, and references therein).In DM-dominated dwarf galaxies, this process is less effective and so significant velocity perturbations (and therefore NCMs) due to halo triaxiality can persist.
O19 showed that NCMs may affect the rotation speed inferred at any given radius, depending on the orientation of the projection axis relative to the main kinematic axes of a gaseous disc.They also showed that conventional kinematic models systematically underestimate the circular velocity in the inner regions of galaxies, due to their inability to account for the geometric thickness of the gas discs.These effects can, and likely do, lead to many measured rotation curves rising more slowly than they should, leading to systematic underestimation of the central DM density.When such systematic effects are accounted for, the diversity of inferred rotation curve shapes of simulated ΛCDM galaxies may be reconciled with that observed (although this depends on the galaxy formation model assumed).
Santos-Santos et al. (2020, hereafter SS20) investigated multiple explanations for the observed diversity of rotation curve shapes using cosmological simulations: BICC, SIDM, and NCMs, as well as its connection with the mass discrepancy-acceleration relation (MDAR;McGaugh 2004;McGaugh, Lelli & Schombert 2016).Although the presence of cores created by BICC or SIDM generated some diversity in the shapes of dwarf rotation curves, models in both categories failed to span the full range observed.They found that these explanations for the diversity in rotation curve shapes could not account for the observation that dwarf galaxies with more steeply-rising rotation curves typically have higher central DM fractions, and vice versa.Systematic errors caused by NCMs, however, were able to account for this, and apparently also for much of the observed diversity in rotation curve shapes, albeit based on an analysis of only two simulated galaxies with DM cusps.
In this work, we investigate whether DM haloes with cusps or cores (or a mix of both) are compatible with the observed diversity both in the shapes of dwarf galaxy rotation curves and the inferred baryonic matter fraction in galaxies' inner regions, once plausible systematic effects in the kinematic modelling process have been accounted for.We use selections of simulated galaxies realised with two galaxy formation models (Secs.2.1-2.2) -one whose SF prescription causes the creation of cores in DM haloes hosting dwarfs, the other whose prescription preserves their initial cusps.We observe these galaxies multiple times each along independent sight lines (Secs.2.3-2.4).These galaxy selections are both shown to produce realistic dwarf galaxy populations as judged from well-known scaling relations (Sec.3.1) and being kinematically and geometrically similar to observed dwarfs ).We quantify the discrepancy between their circular velocity curves and their 'observed' rotation curves, and compare the characteristics of the two 'observed' populations (Sec.4).We discuss the physical interpretation, caveats, and implications of our results in Sec. 5, and summarise in Sec. 6.Throughout this paper, we assume a cosmology with H 0 = 70.4km s −1 Mpc −1 and all other cosmological parameters consistent with WMAP-7 values (Jarosik et al. 2011).

Simulations
The simulations in this work were carried out using a modified version of the N-Body, 'pressure-entropy' (Hopkins 2013) smoothed particle hydrodynamics code Gadget-3, itself a modified version of the Gadget-2 code (Springel 2005).The galaxy formation model is that of the Evolution and Assembly of GaLaxies and their Environments (EAGLE) project (Schaye et al. 2015;Crain et al. 2015).The EAGLE model includes 'subgrid' prescriptions for radiative cooling, star formation, stellar mass loss, feedback and metal enrichment The spherically averaged DM density profiles of the 21 dwarf galaxies selected for kinematics analysis from the LT simulation are shown with blue solid lines.The lines are thinner inside the 'convergence radius' (Power et al. 2003, eq. 20, with  relax ( ) ≳ 0.6 0 ), outside which the circular velocity curve converges to better than ∼ 10 per cent.Two NFW density profiles, with  c,max = 50 and 120 km s −1 , and the median concentration given their mass (Ludlow et al. 2016) are shown with dashed black lines.Lower panel: As upper panel but for the 11 dwarfs selected from the HT simulation, shown in red.LT dwarf haloes retain their  ∼  −1 cusps and therefore closely follow the NFW profile shape at all resolved radii.HT dwarf haloes, however, develop  ∼  0 cores with typical sizes of ∼ 2-3 kpc. of surrounding gas, supermassive black hole (SMBH) gas accretion and mergers, and active galactic nuclei (AGN) feedback.The model's parameters were calibrated against present-day observations of the galaxy stellar mass function, the galaxy size distribution, and the amplitude of the SMBH-stellar mass relation.The simulations' initial conditions are set at  = 127 (see Schaye et al. 2015, appendix B, for details)) and are analysed at  = 0.
We use two variations based on the 'Recal' calibration of the EA-GLE model, originally presented by Benítez-Llambay et al. (2019).The main difference between the two is the value chosen for the parameter controlling the hydrogen number density threshold above which SF is allowed to occur,  thr .The 'low threshold' variant, which we label LT, fixes  thr = 0.1 cm −3 .This value is close to typical values of the metallicity-dependent threshold proposed by (Schaye 2004) used in the fiducial Recal model (see Schaye et al. 2010, eq. 4).The second, 'high threshold' (HT) model increases the threshold to a constant  thr = 10 cm −3 .The increased threshold allows gas to accumulate in the centres of dwarf galaxies, often becoming the locally dominant contributor to the gravitational force, before SF begins and subsequent SN feedback then expels substantial amounts of gas from their central regions.Gas can then reaccrete and the cycle re-peats, with the bursty SF driving much greater fluctuations in the central gravitational force than occur in the LT variant where more continuous SF does not allow such significant gas accumulation and subsequent violent blowouts.In dwarf galaxies in the HT variant, the mechanism identified by Navarro et al. (1996a, see also Read & Gilmore 2005;Pontzen & Governato 2012) operates efficiently, triggering the formation of DM cores (Benítez-Llambay et al. 2019).In the LT variant, all dwarf galaxies retain their DM cusps.
We show the DM density profiles of all dwarf galaxies selected for kinematic analysis (see Sec. 2.2 below) from the LT and HT simulations in the upper and lower panels of Fig. 1, respectively.The DM density profiles of LT dwarfs have shapes closely resembling the Navarro, Frenk & White (1996b, hereafter NFW) profile (assuming a typical concentration; Ludlow et al. 2016) at all well-resolved (Power et al. 2003) radii: the dashed black lines show NFW profiles with  c,max = 50 and 120 km s −1 .The density profiles of HT dwarfs, on the other hand, diverge from the NFW profile shape, becoming nearly flat in the inner few kiloparsecs.We note: (i) that the density profiles of HT galaxies begin to flatten well outside their Power et al. (2003) 'convergence radii', inside which the lines are drawn thinner; (ii) that this convergence criterion is rather conservative as it is defined in terms of a cumulative quantity ( circ ), but density is a differential quantity; and (iii) that the difference with respect to the LT dwarfs, which form ∼ 200 pc cores due to numerical effects, is very clear.We are therefore confident that the cores in the HT dwarfs are not numerical artefacts.
Note that both  thr = 0.1 and 10 cm −3 are far below the densities at which stars realistically form ( thr ∼ 10 4 cm −3 ; Dutton et al. 2019).Increasing  thr is used in this work as a straightforward means of ensuring that gas accumulated in haloes becomes self-gravitating before SF occurs and so induces core formation whilst minimally altering other observables.In other simulations where gas physics is modelled differently (e.g.Jahn et al. 2023), it has been shown that a high  thr is not a necessary condition for baryon-induced cores to form, and others suggest that it may not be a sufficient condition either (e.g.Benítez-Llambay et al. 2019;Dutton et al. 2020).As such, it is primarily the resulting BICC, and not the value of  thr itself, that we are concerned with in this work.
The only other difference between the two model variants is that in the HT model AGN feedback is disabled, while it is enabled in the LT model.However, AGN feedback has negligible influence on EAGLE galaxies in the range of masses on which we focus in this work ( vir ≲ 10 11.5 M ⊙ ; Crain et al. 2015;Bower et al. 2017) -in fact, no SMBHs are seeded in haloes of  vir < 10 10 M ⊙ ℎ −1 (Crain et al. 2015).
Both simulations (LT and HT) use identical initial conditions for a periodic cube with a side length of 12 Mpc (comoving).The DM resolution elements have masses of  DM = 3.2 × 10 6 M ⊙ , and the gas particles have initial masses of  gas = 5.3 × 10 5 M ⊙ , similar to the values used in the fiducial EAGLE-Recal simulations (Schaye et al. 2015).Simulated galaxies are identified using a two-step process.First, a friends-of-friends (FoF) algorithm (Davis et al. 1985) identifies spatially coherent structures, with baryonic particles assigned to the FoF group of their nearest DM particle.Each FoF group is then analysed independently to search for gravitationally self-bound objects within them, using the subfind algorithm (Springel et al. 2001;Dolag et al. 2009).These self-bound regions are referred to as 'subhaloes' and typically each contains a single galaxy.The subhalo containing the particle at minimum gravitational potential within each FoF group is defined as being the 'central' subhalo, with the remainder being labelled 'satellite' subhaloes.

Observational comparison sample
For comparison with our sample of simulated galaxies, we select a subset of the observed galaxies compiled in SS20, table A1, originally from the THINGS (Walter et al. 2008;de Blok et al. 2008) and LITTLE THINGS (Hunter et al. 2012;Oh et al. 2015) surveys, and the SPARC compilation (Lelli, McGaugh & Schombert 2016a), retaining only those with  > 30 • , and  last > 2 fid (see Eq. 2 for definition of  fid ;  last is defined in Sec.2.3).We refer to the subset of this sample with 50 <  ,max /km s −1 < 120 as the 'selected' observational sample throughout this work, whereas those outside this interval in  ,max are labelled 'unselected'.In some cases, our analysis requires the observed data cubes, or at least the derived velocity maps.Such data are not available for all of the galaxies in the full comparison sample above.In such cases, we use a subset further restricted to the THINGS and LITTLE THINGS galaxies for which they are available (see also O19, table A2).
We stress that these comparison samples are only loosely matched to our sample of simulated galaxies; in particular, the selections of observed galaxies are highly heterogeneous in the sense that they are assembled from multiple surveys and observations of individual objects, each with their own criteria for which galaxies should be observed, whereas the simulated sample is drawn from a well-defined volume.

Selection of simulated galaxy samples
We select galaxies from the LT and HT simulations to obtain a sample analogous to observed late-type dwarf galaxies.We also exclude some systems not well suited to kinematic modelling, e.g.ongoing or recent mergers (see Appendix B for details).
Our initial selection is comprised of central (not satellite; see Sec. 2.1) galaxies with maximum circular velocities in the range 50 <  c,max /km s −1 < 120 (halo masses 2 × 10 10 ≲  200 /M ⊙ ≲ 10 11.5 ), corresponding to dwarfs but avoiding the lowest-mass systems, which are not adequately resolved in these simulations to carry out our analysis below.We then impose a minimum H i gas mass,  HI > 10 8 M ⊙ , to restrict the sample to late-types, loosely analogous to those observationally accessible to high-resolution H i imaging surveys.We find that the remaining 46 galaxies from the LT simulation and 43 from the HT simulation are isolated by at least 100 kpc from any companion satisfying the same restrictions, minimising the influence of recent/ongoing tidal interactions.
After the creation of synthetic, or 'mock' observations (Sec.2.3) of this initial sample, we proceed to a visual inspection of the systems' H i surface density and intensity-weighted mean velocity maps as seen from multiple angles, and remove those obviously unsuited to analysis with a tilted-ring model (see Sec. 2.4).Excluded here are galaxies where the gas morphology is very irregular (e.g.obvious ongoing merger, or no clear rotating disc is visible), those with very strong warps, those with strong and obvious radial gas flows within the disc, or those with a companion system that is clearly gravitationally interacting.Full details of the 25 LT and 32 HT galaxies excluded in this step are included in Appendix B.
Finally, after carrying out the kinematic modelling step in our analysis (Sec.2.4), the models for a few galaxies (5 LT and 3 HT) were found to completely fail to describe their corresponding mock observations; we also omit these from further consideration.The disparity in the final sample sizes (21 from the LT and 11 from the HT simulation) is due to the lower number of HT systems with quiescently rotating discs.This in turn is due to the repeated, violent outflows of gas driven by SN feedback (and the subsequent re-assembly of the gas discs) leaving the gas frequently very obviously out of kinematic equilibrium.

Mock observations
The purpose of mock observations is to provide an equivalent to telescope observations for simulated galaxies, such that the same analysis routines can be employed.Here, we mimic the characteristics of the THINGS (Walter et al. 2008) and LITTLE THINGS (Hunter et al. 2012) 21-cm radio surveys.
For each of our selected galaxies, we produce a set of 24 mock observations, each corresponding to a different viewing angle.We begin by measuring the angular momentum vector of the inner 30 per cent (by mass) of the H i disc to define a fiducial disc plane while avoiding the influence of relatively common warps in the outer discs.The 24 sight lines are then distributed at a constant inclination of  = 60 • and spaced by Φ = 15 • in azimuth 1 .The inclination angle is within the optimal range for tilted-ring models, which struggle to capture accurately the kinematics of nearly face-on or edge-on discs.All of our mock observations have a nominal position angle (the angle clockwise from North to the approaching half of the disc) of 270 • .This is enforced only by aligning the component of the angular momentum vector mentioned above in the plane of the sky to point North; individual mock observations of warped or otherwise asymmetric discs may have apparent kinematic position angles differing from this value.We note that pairs of mock observations separated by 180 • in Φ are not equivalent: the discs are geometrically thick and not vertically and azimuthally symmetric, which breaks this rotational symmetry.We refer to (Di Teodoro & Fraternali 2015, fig.2) for a schematic of the geometry of an observed disc.
The 24 orientations for each galaxy are then processed with martini 2 (Oman 2019), a modular Python package for the creation of spatially resolved spectral line observations from smoothed-particle hydrodynamics simulation input (for a detailed description, see O19).This produces mock H i spectroscopic 'data cubes' with two spatial and one spectral axis.The mock observations are constructed assuming a distance of 12.9 Mpc such that a 2 arcsec pixel corresponds to a physical size of 125 pc, and a peculiar velocity of zero.The extent of the spatial axes is the least of 128, 256, or 512 pixels such that the extent is larger than  last , defined as the radius of the sphere enclosing 90 per cent of the subhalo's H i mass.( last for HT-17 was set to 20 kpc since gas from a companion within the same FoF group caused the radius enclosing 90 per cent of the H i mass to be very large.)We use velocity channels with a width of 4 km s −1 , and use enough channels to contain the full H i velocity spectrum of each system.Each data cube is convolved along the spatial axes with a 12 arcsec (FWHM) circular Gaussian beam, giving an effective physical resolution (FWHM) of approximately 750 pc.The 21-cm emission of each gas particle is calculated according to its H i mass fraction.These 1 The reference azimuthal direction corresponding to Φ = 315 • is defined as as a unit vector, where ì   is the angular momentum vector defining the disc orientation.Since the galaxies are effectively randomly oriented within the simulation volume, and the vector (1, 1, 1) (in the simulation coordinates) is arbitrarily chosen, the direction Φ315 , and therefore Φ0 , is likewise arbitrary.mass fractions are corrected for self-shielding from the metagalactic ionising background following Rahmati et al. (2013), and incorporate a pressure-dependent correction for the molecular gas fraction from Blitz & Rosolowsky (2006).We label particles with H i mass fractions above one per cent 'H i-bearing'.The H i gas is assumed to be optically thin -optically thick H i typically occurs in compact (∼ 100 pc) clouds with column densities exceeding 10 22 atoms cm −2 (Braun 2012, see also Allen et al. 2012) which the LT and HT simulations do not resolve.Fig. 2 shows example visualisations of mock observations of LT-26 and HT-28 at an inclination of  = 60 • , and a projection angle of Φ = 0 • .LT-26 and HT-28 are typical of the selected LT and HT galaxy samples, respectively, and of the effect of baryon feedback on dwarfs in the samples.They are also analogous systems, with approximately equal  c,max and disc mass and size (see Appendix B) and so are chosen as illustrative examples throughout this work.The left panels show H i surface density calculated on a grid by a direct summation over simulation particle masses weighted by H i fraction.The centre panels show the flux density (0th moment) of the corresponding data cubes after mock observation.The right panels show the flux-weighted mean line-of-sight velocity with respect to the systemic velocity of the system (1st moment), masked to show only pixels where the flux density exceeds 10 −3 Jy beam −1 , or approximately Σ HI = 1 M ⊙ pc −2 ≈ 10 20 atoms cm −2 , corresponding to the typical limiting depth of observations in the THINGS and LITTLE THINGS surveys.LT-26 appears as though its inclination is less than 60 deg.When orientated at Φ = 90 deg however (see Appendix E) it is clear that this is due to a warp in the outer disc.The warp is also visible in Fig. 2 as a twist in the outer iso-velocity contours.The central region is inclined by 60 deg with a position angle of 270 deg, as intended.

Model fitting
To extract physical parameters such as rotation velocities from the data cubes, a model fitting algorithm is needed.For this, we use the 3D tilted-ring modelling tool 3D barolo 3 (Di Teodoro & Fraternali 2015), a publicly available tool designed to derive rotation curves of galaxies using emission-line observations, such as H i emission.It models the disc as a number of concentric rings, each having its own set of physical parameters (inclination, rotation speed, etc.).The algorithm is '3D' in the sense that it evaluates a figure of merit for a model by comparing the actual and predicted emission in the full data cube, rather than in integrated maps of velocity, dispersion, etc. (the classical '2D' approach).
To analyse our mock observations, we use rings of 750 pc radial width (6 pixels, equivalent to 1 beam width, i.e. consecutive rings contain nearly independent emission), and a number of rings sufficient to enclose  last .Each ring is centered within the data cube, which by construction is the location of the potential minimum of the target galaxy.This potential minimum is generally well traced (within about 150 pc ≈ 3 arcsec) by the peak of the stellar light distribution (see e.g.O19, fig.2), implying that the centres that we assume for our simulated galaxies are equivalent to those which can be measured for observed galaxies.The systemic velocity is fixed to the Hubble Table 1.Ring parameters in 3D barolo.Columns from left to right: the parameters describing each ring; the symbols used for each parameter; whether the parameter is fixed or is free to be adjusted during fitting; the value at which the parameter is fixed or the bounds (if applicable) if the parameter is free; and lastly, the units used for the values given.

Ring parameter
Symbol Fixed?Value Units Centre coordinates .The H i column density is fixed such that the integrated flux along the spectral axis in each pixel in the model is equal to the integrated flux in the corresponding pixel in the mock observations, i.e. 3D barolo's 'local normalisation' option.This prevents strongly non-uniform regions of the disc, e.g.gas overdensities or holes, from disproportionately biasing the model (Lelli et al. 2012).We have tested that adopting 3D barolo's 'azimuthal normalisation' option instead produces qualitatively equivalent results (see Di Teodoro & Fraternali 2015, for details of normalisation options).
Each ring is fully specified by the eight parameters listed in Table 1, four of which are free to vary during the fitting process: the inclination angle, position angle, rotation velocity and velocity dispersion.The inclination and position angle are constrained to remain within ±15 • and ±20 • , respectively, of their 'true' values (60 • and 270 • , respectively).We provide initial guesses of  = 60 • , PA = 270 • ,   = 30 kms −1 , and  = 8 km s −1 .The output of 3D barolo is minimally sensitive to the initial guesses, but is sensitive to the allowed ranges for position angle and, especially, inclination (Di Teodoro & Fraternali 2015), motivating us to enforce these relatively narrow allowed ranges.The model optimisation proceeds in two stages.In the first stage, all four of these parameters are free to vary.Once each ring has been independently optimised in this way, the geometric parameters (inclination and position angle) are 'regularised': the values as a function of radius are fit with Bézier functions.In the second stage, the geometric parameters are kept fixed to the value of their Bézier function at the corresponding ring radius and the rotation speed and velocity dispersion are re-optimised to construct the final model.Full details of the 3D barolo configuration used are listed in Appendix C.
We tested some alternative choices for 3D barolo's configuration parameters.We have checked that replacing the Bézier regularisation functions with polynomials of 0 th , 1 st , 2 nd , or 3 rd order produces qualitatively similar results.We have also checked that allowing 3D barolo to fit a radial velocity for each ring, and removing the S/N mask also give qualitatively equivalent results.We are therefore confident that our conclusions are not sensitively dependent on our chosen configuration.

REALISM OF SIMULATED GALAXIES
Throughout this work, simulated galaxies are mock observed and analysed to determine how increased baryon feedback and 'observation' affect their inferred rotation curves.Whether our conclusions below are applicable to real galaxies depends on whether our simulated galaxies are sufficiently similar to these, such that the observation and modelling process may plausibly lead to similar consequences.We therefore compare the overall properties and structure of our sample of simulated galaxies to those observed in this section.

Scaling relations
The simulated galaxies lie on, or close to, key scaling relations, as shown in Fig. 3. Galaxies in our final selected sample (squares) are shown in context with all central galaxies in the simulations (points), both those lying within the ranges 50 <  c,max /km s −1 < 120 and  HI > 10 8 M ⊙ but removed from our sample after visual inspection (large points; see Supplementary Tables B2 and B3) and those outside these ranges (small points).These are compared to the observational comparison sample (Sec.2.2.1), including both galaxies within the same  ,max range (large crosses) and outside (small crosses).
The left panel of Fig. 3 shows the baryonic Tully-Fisher relation (BTFR), where our selection of galaxies from both the LT (blue symbols) and HT (red symbols) simulations align very well with the observational data.Here, the definition of  max varies by dataset: for simulated systems,  max =  c,max is the maximum of the circular velocity curve4 ,  circ () = √︁ G (< )/, while for observed systems  max =  ,max is equal to the maximum velocity in their rotation curve 5 .In this figure,  bar =  ★ + 1.4 HI for all systems (see e.g.

McGaugh 2012).
The centre panel shows the H i mass-stellar mass relation.The stellar masses for observed galaxies assume 'diet Salpeter' and Chabrier IMFs for systems taken from the THINGS (Walter et al. 2008) and LITTLE THINGS (Hunter et al. 2012) surveys, respectively.For systems in the SPARC database (Lelli et al. 2016a), mass-to-light ratios at 3.6 µm of 0.7 and 0.5 M ⊙ L ⊙ −1 are assumed for bulge and disc components, respectively (Lelli et al. 2016b).For the simulated galaxies,  ★ is a direct summation of the masses of gravitationally bound 'star' particles.The simulated systems in the LT sample (blue squares) approximately follow the observed relation, whereas our selection of galaxies in the HT model (red squares) are systematically offset from both the observed relation (crosses) and the ensemble of central galaxies in the same simulation (red squares and points): they have H i masses approximately 0.5 dex higher than typical systems of the same stellar mass.The HT galaxies excluded from our selection have in many cases recently lost large amount of H i, much more than their increase in stellar mass over the same period -the same events (e.g.bursts of star formation) that triggered this loss of gas are likely responsible for their more disturbed morphologies and kinematics.
Those included in our selection, on the other hand, are typically undergoing a period of rapid accretion and slow star formation.These will likely eventually undergo a burst of star formation, increasing their stellar mass and depleting their H i, to bring them back towards the median relation.
The right panel shows the H i mass-size relation, where size is defined as the radius at which the H i surface density falls to Σ HI = 1 M ⊙ pc −2 ≈ 10 20 atoms cm −2 .The majority of LT systems, at fixed  HI , have approximately 0.2 dex larger sizes than observed, while HT systems fall within the observed scatter.Fig. 4 shows the gas mass-specific angular momentum relation, often referred to as a Fall (1983) relation.The empirical  gas - gas relation of Mancera Piña et al. (2021) based on a sample of observed dwarf galaxies from SPARC (Lelli et al. 2016a) and LIT-TLE THINGS (Hunter et al. 2012, see also Sec. 2.2.1) as well as further late-type dwarfs from LVHIS (Koribalski et al. 2018), VLA-ANGST (Ott et al. 2012), and WHISP (van der Hulst, van Albada & Sancisi 2001) observations, is shown with the black dashed line.These observations sample the mass range 10 7.5 <  gas < 10 10.5 .The selected LT sample follows the observed relation reasonably well whereas the HT sample, and particularly the subsample selected for kinematic modelling, is systematically offset to higher  gas and/or lower  gas values when compared to the LT sample.The offset between the LT and HT galaxies in this space is consistent with the typical differences between their H i masses and sizes (see Fig. 3).The somewhat larger H i disc sizes (at fixed H i mass) of LT galaxies relative to observed galaxies are compensated by their slightly lower H i masses (at fixed  c,max ) such that their specific angular momenta are similar to those of observed galaxies (  gas ≈  HI  c,max ).Although HT galaxies have H i sizes very similar to observed galaxies (at fixed H i mass), their substantially higher H i masses (at fixed  c,max ) push them well off the observed  gas - gas relation.As discussed by Benítez-Llambay et al. (2019), the large H i mass of HT galaxies is due to the EAGLE feedback implementation, which becomes inefficient at pushing gas out for the HT galaxies.Furthermore, HT galaxies selected for kinematic modelling (red squares) have comparatively high average  gas when compared to the full HT sample.This appears to be a selection bias, as galaxies with recent accretions of gas (following previous violent gas ejections) are more likely to be more suitable for modelling than those with recent turbulent ejections (see above discussion of the centre panel of Fig. 3).
Both simulated populations follow the stellar  ★ - ★ relation (not shown) from Mancera Piña et al. (2021), though the selected LT and HT samples have approximately 3 and 5 times greater scatter than observed, respectively.

Non-circular motions
There is no doubt that the presence of NCMs in the gas discs of observed dwarf galaxies influences the measurement of their rotation curves (see e.g.Valenzuela et al. 2007;Oman et al. 2015;Marasco et al. 2018; O19; SS20).It is therefore important that our simulated galaxies replicate the observed amplitudes of NCMs.Following Schoenmakers, Franx & de Zeeuw (1997), the amplitudes of NCMs are often quantified in terms of the amplitudes, either individually or collectively, of the terms   ′ >0 in the Fourier expansion of the line-of-sight velocity,  LoS , as a function of azimuthal angle6 , , around a ring of gas: with the overall amplitude of a given order defined as We recall that a term of order  in the Fourier expansion of the de-projected velocity field contributes to the  ′ =  ± 1 terms in Eq. 1 -for instance, a bisymmetric pattern such as might be induced by a bar is reflected in the amplitudes  1 ′ and  3 ′ (e.g.Schoenmakers, Franx & de Zeeuw 1997).
Multiple attempts to quantify the amplitude of NCMs from observational data have been made.Trachternach et al. (2008), for instance, found a median radially-averaged amplitude of NCMs in THINGS galaxies of 6.7 km s −1 , with approximately 90 per cent of galaxies with measurements of less than ∼ 9 km s −1 (see also Oh et al. 2011, for further discussion of a subset of the same sample; they find amplitudes ≲ 10 km s −1 using both a similar and a complementary methodology).Oh et al. (2015) find similarly low amplitudes (≲ 10 km s −1 in most cases, up to 20 km s −1 in a few cases) for the first three harmonics in LITTLE THINGS galaxies, but also note a crucial point: galaxy-scale NCMs are strongly degenerate with many other parameters such as the centroid, systemic velocity, position angle, and inclination (see also Schoenmakers et al. 1997;Wong, Blitz & Bosma 2004;Spekkens & Sellwood 2007).It is thus all too easy to minimise residuals without need for significant NCMs and so underestimate the true amplitudes of NCMs in observed galaxies.Effects associated with NCMs are then not included in the resulting model and subsequent analysis, meaning inferred rotation velocities can be profoundly in error.O19 illustrated this further: they found that parameter degeneracies left them entirely unable to measure the amplitude of the dominant (∼ 15 km s −1 )  = 2 harmonic in two simulated dwarf galaxies from mock observations, even though these caused errors of up to 10 km s −1 (≥ 20 per cent) in their rotation curves.
Even though the true amplitudes of a harmonic expansion of the gas velocities in a galaxy are seemingly observationally inaccessible, we may still ask: are the apparent amplitudes of NCMs seen in our simulations comparable to those in observed dwarfs?This must, however, be based on the application of an identical measurement to both observed and mock observed simulated data.This was the approach adopted by Marasco et al. (2018): in their fig.8, they show the relative amplitudes7  3 ′ / 1 ′ of the  ′ = 3 and  ′ = 1 harmonics as a function of  1 ′ (corrected for inclination) at a given radius.The ratio of these two harmonics is indicative of the strength of the bisymmetric ( = 2) component of the de-projected velocity field relative to the rotation velocity at the same radius.We reproduce in Fig. 5 their measurements at 2 kpc for THINGS and LITTLE THINGS galaxies (we omit the SPARC sample as the required observational data cubes are typically not available) with 50 <  ,max /km s −1 < 120.We have repeated the same measurement for our mock observations of LT and HT galaxies, shown with blue and red symbols in Fig. 5, for the Φ = 0 • and 90 • orientations.Overall, the amplitude of the  ′ = 3 harmonic,  3 ′ , in the simulated galaxies from both models is very similar to what is measured in observed galaxies.In the simulated galaxies, this is typically the harmonic with the largest amplitude (e.g.Marasco et al. 2018;O19).The simulated LT galaxies seem to have slightly higher  3 ′ on average than both the HT and observed galaxies8 , perhaps by a factor of 2, although the scatter in the distributions, and some cases the measurement uncertainties, are large.

Disc thickness
Conventional kinematic models also struggle to account for the thickness of observed discs as this implies that emission from multiple rings contributes to emission along any single sight line through an inclined disc, thus introducing many degeneracies into the model (see Sec. 5.1.2for further discussion).To avoid this, models typically assume thicknesses far smaller than those of observed discs which galaxy formation model, which leads to BICC, induces much stronger NCMs than a reference model where BICC does not occur.This may, however, be attributable to the presence of highly-disturbed galaxies in their sample, which would be obviously unsuitable for mass modelling, and which we have excluded from Fig. 5 (the HT model has a larger fraction of such galaxies).2021) for an expanded sample of observed galaxies relative to that used in this work: log(  gas ) = 1.02(log(  gas ) − 10) + 3.64; with an orthogonal intrinsic scatter of 0.15 dex shown by the grey band.The selected LT sample follow the observed  gas - gas relation well, with slightly larger scatter, however, the selected HT sample lies at systematically higher  gas than the observed relation.
induces errors in the measured rotation curves -i.e. if disc thicknesses were accurately known and incorporated into these models, the resulting rotation curves would be better representations of the discs' true rotation.Correspondingly, it is important that our simulated galaxies have realistic H i disc thicknesses so that we obtain results comparable to analyses of observed galaxies.Our sample of LT galaxies have very similar vertical gas structure to the APOSTLE galaxies studied by O19; we show their height profiles in Fig. 6 (c.f.O19, fig.12).We also show the height profiles of our sample of HT galaxies in the same figure.The HT galaxies are on average a factor of ∼ 2 thinner at all radii than their LT counterparts, likely due to their higher mid-plane gas densities (see e.g.Benítez-Llambay et al. 2018), but still are considerably thicker than the 100 pc assumed in our 3D barolo configuration.Observed late-type dwarfs have estimated half-mass heights of about 0.2-1.0kpc at a radius of 5 kpc (O'Brien, Freeman & van der Kruit 2010, see their fig.24; Banerjee et al. 2011; see also Peters et al. 2017); this range is broadly representative of further estimations of H i disc heights of nearby dwarfs of similar masses (see e.g.Patra 2020; Bacchini et al. 2022;Mancera Piña et al. 2022;Li et al. 2022, noting that differing definitions of H i height are used).We stress that these measurements are uncertain and can require many assumptions of, e.g., vertical geometry, gravitational contributions of gas, stars, and DM, and/or hydrostatic equilibrium (see O19, sec.6.2, for further discussion).O19 established that APOSTLE dwarf galaxies (and therefore by extension the similar LT dwarfs) are likely somewhat thicker than these observations imply, by a factor of ≲ 2. Our sample of HT dwarf galaxies therefore have gas discs with thicknesses fully consistent with these observational constraints.2011) are shown as the black solid and dashed bars, respectively.The atomic gas discs of galaxies selected from the HT simulation are typically approximately half as thick as those of those selected from the LT simulation over the radial range shown, making HT dwarfs fully consistent with observation.

Velocity dispersion
The thickness of an atomic gas disc is closely related to its velocity dispersion: all else being equal, a more turbulent disc is also thicker.O19 (see their fig.3) found that APOSTLE dwarf galaxies typically have somewhat larger velocity dispersions than observed galaxies at a given H i mass.In Fig. 7, we show that this is also true of our sample of LT galaxies selected for kinematic analysis (blue squares).The black crosses mark the measurements for the observed sample of THINGS and LITTLE THINGS galaxies shown in Fig. 5. HT galaxies selected for kinematic analysis (red squares) have somewhat lower velocity dispersions than their LT counterparts, and align very well with the observed sample.When we consider instead velocity dispersion profiles (along with H i surface density profiles, see Appendix E), this is found to be true at all radii.
This likely reflects their more realistic H i disc heights.We note that this observational sample contains only galaxies of sufficient kinematic regularity to be selected for mass modelling by the survey teams.Therefore, the LT and HT galaxies that we exclude from our Figure 7. Median line-of-sight velocity dispersion (calculated across pixels) as a function of H i mass.Measurements for mock observed galaxies with 50 <  c,max /km s −1 < 120 from the LT and HT simulations are shown in blue and red, respectively, with those selected for kinematic modelling shown with squares and those excluded shown with points.Measurements of selected observed galaxies from the THINGS and LITTLE THINGS surveys are shown as black crosses.For mock observed galaxies, we show measurements for two projection angles, Φ = 0 • and 90 • , and the median is calculated across all pixels with 21-cm flux density above 10 −3 Jy beam −1 (corresponding approximately to H i surface density above 1 M ⊙ pc −2 ).For observed galaxies, it is calculated across all pixels in the S/N-masked, natural-weighted 2 nd moment maps provided in the survey data releases.Mock observed, simulated galaxies have median velocity dispersions similar to those of the observed sample, although the LT sample has, on average, slightly higher velocity dispersions at a given H i mass.
sample for kinematic analysis (blue and red points), which often have much higher velocity dispersions, should not be compared directly with this sample (see also discussion of the centre panel of Fig. 3 in Sec.3.1).Observed rotation curves are often corrected for support by a radial pressure gradient9 .This correction is estimated directly by 3D barolo.In most cases the correction increases the rotation curve by only 5 to 10 per cent, essentially independently of radius, though in some rare cases it can be as large as ≳ 40 per cent.These corrections thus do not play a substantive role in our results.

Circular velocity and rotation curves of simulated galaxies
We show the pressure support-corrected rotation curves,   , from 3D barolo model fitting for each of the 24 mock observation orientations for LT-26 and HT-28 in Fig. 8 with the pale blue lines -one curve for each galaxy is highlighted in solid magenta, corresponding to the Φ = 0 • orientation visualised in Fig. 2. For this orientation, we also show the rotation curve before pressure-support correction (dotted magenta line; see Sec. 3.4).The uncertainty associated with   is taken to be the uncertainty on the rotation velocity calculated by 3D barolo (see Di Teodoro & Fraternali 2015, sec.2.5), with no contribution from the uncertainty on the pressure-support correction, as this is typically negligible.This uncertainty is shown for the Φ = 0 • orientations with the shaded magenta region 10 .
These pressure support-corrected   curves are compared with the circular velocity curve (determined from the gravitational acceleration field in the disc midplane; black line) and the median azimuthal velocity as a function of radius of H i-bearing simulation particles ( az ; green line) of each system.Using the magenta curves, we illustrate our definition of an outer, 'maximum' rotation velocity,  ,max , and an inner, 'fiducial' rotation velocity,  ,fid (dashed magenta lines, as labelled).
For simulated galaxies,  ,max is determined as the asymptotic flat rotation velocity, or if this is ill-defined,  ,max ≡   (10 kpc) (or   ( last ), if  last < 10 kpc; see Appendix D for details of the definition of  ,max ).For circular velocity curves, the asymptotic flat rotation velocity is generally equivalent to the true maximum velocity of the curve and can be used interchangeably, but for rotation curves this is often not the case: they are more irregular, and a local upward fluctuation can often cause the maximum of the curve to overestimate the asymptotically flat velocity.Our adopted definition minimises the influence of such local fluctuations in the rotation curves.
For  fid , we follow SS20 and define: (2) This definition adapts to the typical size of galaxies of given  max to yield an inner rotation velocity tracing their central matter content. fid is also typically a radius that can be resolved by most highresolution H i observations of local dwarfs (see e.g.Hunter et al. 2012), enabling comparison with observations.

Rotation curve shapes
While a qualitative impression of the structure of the galaxies can be obtained by inspection of the curves in Figs. 1 and 8, in order to quantify the connection between the dark and baryonic matter distributions we again follow SS20.We show in the left panel of Fig. 9 (similar to their fig.2)  fid plotted against  max , which quantifies how steeply the velocity curve rises.A very steeply-rising curve could already reach an amplitude of  max at the (relatively small) radius  fid -such a rotation curve would correspond to a point lying on the dashed line along  fid =  max in the figure .A matter distribution with a cusp corresponds to a relatively steeply-rising circular velocity curve, such that by  fid it is already approaching  max .For an NFW halo model appropriate for a dwarf galaxy with typical concentration (Ludlow et al. 2016),  c,fid is typically about 0.65 c,max ; the precise relation is shown with a solid black line in the figure .A constant-density core in the matter distribution instead leads to a shallower, linear rise in rotation velocity within the core radius.Cores with sizes ≳  fid therefore have lower  fid values.We note that  max and  fid are set by the shape of the total (dark plus baryonic) matter distribution.However, since the DM is often the dominant contributor to the gravitational force at all radii in dwarfs, 10 3D barolo uses an iterative algorithm to estimate errors which very occasionally fails to converge.This is the case for the outermost ring in the model for HT-28 at Φ = 0 • .

HT-28
Figure 8. Orbital velocities as a function of radius for the gas discs of simulated galaxies LT-26 and HT-28, the same example galaxies shown in Fig. 2. The circular velocity,  circ () = sgn(  () ) √︁  |   () |, where   is the radial component of the gravitational acceleration field in the disc midplane, is shown in black.The median azimuthal velocity of H i-bearing particles,  az , reflecting the actual motion of the gas, is shown in green.The 24 solid blue and magenta curves show the rotation velocity,   , obtained by modelling mock observations of the galaxies seen along 24 different sight lines at fixed inclination (60 • ).These have been corrected for pressure support.The pressure support-corrected   from mock observations for the projection angle Φ = 0 • , corresponding to the images shown in Fig. 2, is shown in solid magenta, with uncertainty shown by the magenta shading.The corresponding  ,max ,  ,fid , and  ,fid are shown with dashed magenta lines (see Eq. 2 and Appendix D for definitions).The rotation velocity without pressure support-correction for Φ = 0 • mock observations is shown with dotted magenta lines. in many cases these two parameters can in principle provide a direct measure of the DM density profile shape.
The galaxy LT-26 has a DM cusp; its steeply-rising circular velocity curve (Fig. 8, black curve in upper panel) reflects this.The values of  c,fid and  c,max measured from this curve are plotted in the left panel of Fig. 9 with a blue square symbol lying nearly directly on the line corresponding to an NFW halo model.The corresponding galaxy in the high SF density threshold simulation, HT-28, has had its DM cusp transformed into a core after repeated gas accretion and expulsion cycles.This is reflected in the more slowly-rising shape of its circular velocity curve (Fig. 8, black curve in lower panel).Since  c,max is left essentially unchanged by the formation of the DM core, the slower rise of the rotation curve leads to a lower  c,fid at approximately fixed  c,max , and the red square corresponding to the circular velocity curve of this galaxy lies below the black solid line in the left panel of Fig. 9.
When these two galaxies are mock observed and modelled, the resulting rotation curves do not exactly recover the corresponding circular velocity curves -in some cases, the discrepancy can be quite severe, as is visible in Fig. 8.In the left panel of Fig. 9, each of the 24 red and blue points corresponds to one of the 24 rotation curves measured for LT-26 and HT-28, respectively.Consecutive observations (offset by 15 • in Φ) are linked by lines.The orientation labelled Φ = 0 • , which we recall has no particular significance but is highlighted in Fig. 8, is larger than the other points and shown with error bars.The fractional uncertainty in both  ,max and  ,fid has no significant dependence on Φ and so is not shown for all points.
For LT-26, the blue points lie systematically below the blue square symbol, i.e.  c,fid is systematically underestimated by  ,fid when modelling the mock observations, and there is considerable scatter across viewing angles in both  ,fid and  ,max -a significantly greater scatter than the uncertainty associated with individual points.O19 (see also Marasco et al. 2018) attributed the scatter to NCMs and the systematic underestimate of the thickness of the H i discs of dwarf galaxies in the APOSTLE simulations.These conclusions also apply to the LT simulations, which are very similar to APOSTLE.The left panel of Fig. 9 illustrates that similar scattering effects are at play in our kinematic analyses of galaxies from the HT simulation however the underestimate of  c,fid is not present, likely due to the thinner gas discs of HT systems (see Sec. 3.3).

Rotation curve shape and central baryon dominance
We also adopt two informative velocity ratios used by SS20: the rotation curve shape parameter,  rot , and the central baryonic matter fraction parameter,  bar . rot condenses the information contained in the left panel of Fig. 9 into a single parameter, and therefore is a measure of how steeply a rotation curve rises.It is defined as: The circular velocity curve of an NFW DM halo has  rot ∼ 0.65 (i.e. they lie near the solid black lines in the left and centre panels of Fig. 9), while haloes hosting cores have correspondingly lower values.
The second velocity ratio we use is  bar , defined as: where  b,fid =  bar ( fid ) is the contribution to the circular velocity due to baryonic matter (stars and gas) at  fid .Under the assumption of spherical symmetry,  bar is equal to the baryon mass fraction within the fiducial radius.This assumption never holds exactly for late-type dwarfs, however, DM-dominated systems are likely reasonably close to spherical, such that  bar serves as an approximate tracer of the importance of the baryons in setting the kinematics near the galactic centres.For simulated galaxies, we derive  b,fid from the median gravitational acceleration due to baryon (star and gas) particles at a sample of points evenly distributed in azimuth in the disc midplane at  fid .For observed galaxies, it is derived from stellar (photometric) and H i surface density profiles under the assumption of either razorthin disc or vertically extended disc geometry (for details, see  8) as a function of outer 'maximum' velocity (see Appendix D).The square markers show the values from the circular velocity curves (derived from the gravitational acceleration in the disc midplane) of galaxies LT-24 (blue) and HT-24 (red).The blue and red points show the values from each of the 24 rotation curves obtained from mock observations for the same two galaxies.The measurement for projection angle Φ = 0 • , shown in Fig. 2 and highlighted in Fig. 8, is shown as an enlarged point with error bars.This is connected to all measurements for values of Φ in 15 We plot these two velocity ratios against each other in the centre panel of Fig. 9.As outlined above, points corresponding to steeplyrising rotation curves lie above those corresponding to slowly rising ones in this figure, and systems with a larger (apparent) contribution by baryons to the central mass content lie to the right of those centrally dominated by DM.
The values as determined from the circular velocity curves (and baryonic mass profiles) of LT-26 and HT-28 are shown with the blue and red squares, respectively.In addition to having a more slowlyrising circular velocity curve than LT-26 (as already seen in the left panel), HT-26 is also less centrally DM-dominated, with the baryons accounting for ∼ 40 per cent of the mass within  fid , compared to only ∼ 12 per cent in LT-26.As will be confirmed below, this is a feature typical of late-type dwarfs in the HT model: the higher gas density threshold for SF allows gas to condense in the galactic centres and become more tightly gravitationally bound, allowing a dense, massive gas disc to assemble11 (see also Fig. 3, where it can be seen that the HT dwarfs that we select for kinematic analysis are typically somewhat more gas-rich, and yet have smaller  HI , than their LT counterparts).
When  rot and  bar are plotted against each other for the ensemble of mock observations of a single galaxy, an anti-correlation emerges: rotation curves that rise more rapidly also give the appearance of a more DM-dominated central region.The reason for this is straightforward: if  c,fid is underestimated by  ,fid , for example, the rotation curve rises more slowly, but the absolute contribution of the baryons to the circular velocity at the fiducial radius,  b,fid , remains fixed, so the galaxy appears more centrally baryon-dominated, and vice-versa for an overestimate of  c,fid . b,fid will remain fixed whether it is determined directly from simulation, as in this work, or from the H i flux density and stellar photometry, as neither method is dynamical:  b,fid is independent of  ,fid .Some scatter is introduced in the anti-correlation due to the fact that  c,max may also be over-or underestimated by  ,max during kinematic modelling,12 represented by the error bars shown for the enlarged Φ = 0 • point.
The right panels of Fig. 9 illustrate that the values of  rot and  bar do not vary randomly as a function of the viewing angle Φ: both oscillate in an approximately sinusoidal pattern with a period of 180 • .This demonstrates that the measurements of the rotation curves near the centres of these galaxies are affected by the presence of predominantly bisymmetric NCMs.Similar patterns are seen in all of the galaxies in our samples of LT and HT galaxies (see additional figures in Appendix E).This is reminiscent of the analysis of APOSTLE galaxies by O19 (see for instance their fig.6) -it is unsurprising that the very similar LT simulation which we use here exhibits the same behaviour.Apparently, galaxies in the HT model also have significant bisymmetric distortions of the gas flows near their centres.

Rotation curve shape and baryon surface density
Rotation curve shape is known to correlate with the surface brightness of galaxies (e.g. de Blok, McGaugh & van der Hulst 1996;Swaters et al. 2009Swaters et al. , 2012;;Lelli, Fraternali & Verheĳen 2013), indicating that baryons may play a significant role in shaping their mass distributions.In Fig. 10 we show this correlation for both observed (black crosses) and simulated (blue and red markers) dwarf Figure 10.Rotation curve shape parameter,  rot , as a function of effective baryon surface density, Σ bar .Measurements for selected galaxies from the LT and HT simulations are shown in blue and red, respectively.Measurements for selected observed galaxies from the THINGS and LITTLE THINGS surveys are shown with black crosses.Σ bar for simulated galaxies is calculated directly from the simulation (i.e.not from mock observations), and  rot for the square markers is obtained from the circular velocity curves (derived from the gravitational acceleration field in the disc midplane).However, points show  rot from each of the 24 rotation curves obtained from mock observations of each galaxy.The solid black line shows the average  rot expected for an NFW halo model and the grey band shows the 10 th -90 th percentile scatter around this value (Ludlow et al. 2016).With mock observation, LT and HT values align reasonably well with the weak positive correlation seen in the observed dwarfs.
galaxies from our selected samples (see Sec. 2.2).We use  rot to capture the shape of the circular velocity curve derived from the gravitational acceleration field in the disc midplane (squares) and from rotation curves from mock observations (points) as described in Sec.4.2.The 'effective' baryon surface density, Σ bar , of each galaxy is Σ bar =  bar /2 2 b,half , where  bar is the total baryonic mass and  b,half is the radius of a sphere enclosing 0.5 bar -both calculated directly from the distribution of particles in the simulations 13 .
In the observed population the correlation between  rot and Σ bar is primarily due to the fact that  rot has a strong dependence on  max : in general, massive galaxies (i.e.high  max ) are observed to have steeply-rising rotation curves (large  rot ; see e.g.SS20, fig.3).Observed dwarf galaxies ( ,max ≲ 120/ km s −1 ), however, have a wide range of  rot at fixed Σ bar .Baryon surface density therefore cannot be used to predict the presence of a core or cusp reliably in any particular dwarf galaxy.We consider now whether this picture changes when the simulated galaxies are 'observed'.
HT galaxies within our interval of interest in  c,max typically have intrinsically higher Σ bar than LT galaxies (although there are two outliers), as expected due to their denser H i discs.Both LT 13 Calculation of Σ bar from mock 21-cm observations alone is not possible.However, we find that Σ HI derived from a mock observation is typically very close to the value obtained directly from the simulation particle distribution.We therefore have no reason to suppose that mock observation would significantly affect Σ ★ , and therefore the Σ bar values shown in Fig. 10. and HT galaxies, as measured using their 'true' circular velocity curves, replicate the observed positive correlation between  rot and Σ bar .However, neither replicate the observed scatter.The points in Fig. 10 show the approximate distribution of  rot obtained for each galaxy when their rotation curves are instead measured from mock observations.Broadly speaking, these points better cover the observed range in  rot and Σ bar , although it seems that the HT model struggles to reproduce galaxies observed to have high central baryon densities and very steeply rising rotation curves ( rot ∼ 1).

Which simulation better reproduces the observed galaxy population?
In the upper panels of Fig. 11 we show measurements from the full observational comparison sample in both the  ,fid - ,max (upperleft panel) and  bar - rot planes.The observed 'selected' galaxies have the additional constraint 50 <  ,max /km s −1 < 120 (the same as for the simulated galaxies) and are shown with crosses coloured by  ,max .These are repeated in the lower panels with small cross symbols.Galaxies with  ,max above and below this range are shown with upward-and downward-pointed triangles, respectively, and are not repeated in the other panels.
Observed galaxies scatter widely around the solid black line marking the expectation for an NFW halo model (corresponding to  rot ≈ 0.65).The scatter is significantly greater than that expected from the expected scatter in halo concentration (gray shaded band) -this is the diversity of rotation curve shapes highlighted by Oman et al. (2015).
Selected galaxies are also observed to have a diversity of central baryon mass fractions between  bar ≈ 0.05 and 0.7 (upper-right panel).Observed dwarfs (darker symbols) exhibit a (weak) anticorrelation reminiscent of that seen in Fig. 9, such that more centrally DM-dominated galaxies (i.e.lower  bar ) have more steeply-rising rotation curves (i.e. higher  rot ), and vice versa.SS20 note that this trend seems inconsistent, a priori, with BICC and SIDM models where almost all galaxies are predicted to have cores and cusps only form in baryon-dominated systems, and showed that mock observation alone could replicate this trend.
For any given observed galaxy, we are of course limited to a single sight line, such that an observational measurement leading to an equivalent of the points linked by lines in Fig. 9 is impossible.However, it is plausible that an ensemble of measurements for a collection of at least loosely similar galaxies (e.g.dwarfs) could preserve the trends outlined in Sec.4.2, independently of whether the galaxies have DM cusps, or cores produced by BICC.
In Fig. 11 we show the  c,max and  c,fid values of the circular velocity curves (derived from the gravitational acceleration field in the disc midplane) of the 21 LT and 11 HT galaxies in our mock observed sample with the square symbols in the centre-left and lowerleft panels, respectively.The measurements for the LT galaxies lie close to the solid black line marking the relation for an NFW halo of typical concentration: they all retain their central DM cusps (in the few galaxies with  c,max ≳ 90 km s −1 , the haloes are contracted, driving up  c,fid ).The measurements for the HT galaxies, on the other hand, lie systematically below the same line, reflecting the central DM deficit resulting from BICC.
If we assume that observational measurements (crosses) accurately trace the circular velocity curves of the observed galaxies, clearly neither the LT model nor the HT model is able to explain the observed diversity in terms of the circular velocity curves of the galaxies (squares).LT galaxies (blue squares) show much smaller scatter, and do not extend to values of  fid as low as those of observed  systems.HT galaxies (red squares), on the other hand, are all below the solid black line because of their pronounced cores.When the scatter induced by NCMs and the biases induced by disc thickness are taken into account, however, a different picture emerges.We turn our attention to the distribution in  ,max and  ,fid of the rotation curves of the simulated galaxies (i.e.mock observed and modelled).The measurement corresponding to the orientation Φ = 0 • for each galaxy is shown with a point with black border in the centre-left and lower-left panels of Fig. 11.We recall that the direction defined by Φ = 0 • has no particular significance, so this is equivalent to choosing a random orientation (at fixed  = 60 • ) for each galaxy.There are 24 ( ,max ,  ,fid ) points for every square symbol, each corresponding to a different viewing angle, Φ.This provides a reasonable estimate of the full range in the parameter space which can be reached by galaxies in the LT and HT models when they are 'observed'.The scatter around the values derived from the circular velocity curves (squares) is substantial for both models.While the galaxies from the LT model essentially fully cover (and even somewhat exceed) the extent of the observed distribution, the HT galaxies only rarely reach the region at the highest  ,max and  ,fid sampled by the observations.The LT sample therefore seems to better align with observations in this space, suggesting that a model where every galaxy has a cusp, but appears at times to have a core, may be more plausible than a model where every galaxy has a core as large as those created in the HT model.
LT galaxies also reproduce the anti-correlation between  rot and  bar seen in the observed population of dwarfs quite naturally, spanning the full range covered by observations.The HT model, on the other hand, seems to struggle to reproduce observed dwarfs with steeply rising rotation curves ( rot ∼ 1) where the inner regions are DM-dominated ( bar < 0.3).Furthermore, the HT galaxies with the highest  ,max values in our sample preferentially lie toward lower  rot and  bar , a feature absent for the LT galaxies, and running contrary to the observations: more massive galaxies typically have more steeply-rising rotation curves and are more centrally baryondominated (e.g.Broeils 1992;de Blok et al. 2008, see also the upperright panel of Fig. 11).The discrepancy between the 'observed' rotation curves of simulated galaxies and their circular velocity curves driven primarily by NCMs and thick H i discs seems better able to reproduce the broad trends of the observed distribution of dwarfs in  rot and  bar than any of the other possible explanations considered by SS20 (i.e.SIDM, BICC without modelling errors, or the MDAR), particularly in the LT model where all galaxies have DM cusps.
In light of the results presented in this section, we argue that if the rotation curves of observed late-type dwarfs are significantly affected by errors due to NCMs and the geometric thickness of their H i discs, the data favour a population of late-type dwarfs predominantly having DM cusps.In this interpretation, the observed DM 'cores' in such objects are merely symptomatic of the failure of a rotation curve to trace the circular velocity curve of a galaxy.We note that we have not here explored the cores created by mechanisms other than BICC, such as SIDM.

DISCUSSION
In Sec.4.4, we argued that over-and underestimates of the central portions of rotation curves are common in our analysis of mock observations, and are the dominant driver of the scatter in the mock observed values of  rot and  bar along the direction of anti-correlation between  rot and  bar .Discrepancies extending over the entire radial extent of the rotation curve are also relatively common, and lead to scatter in  bar while leaving  rot approximately constant.Overestimates of  bar and/or underestimates of  rot seem somewhat more common than the converse in our analysis; these can arise if the rotation curve is underestimated either in the centre, or at all radii.We also noted a preference for HT galaxies with higher  ,max to lie at lower  rot and lower  bar , a feature absent for LT galaxies.All of these trends and their connection to rotation curve shapes are summarised in Fig. 12.We next turn our attention to the possible physical origins of such effects, both for our simulated galaxies and their observed counterparts.

Non-circular motions
If a model assuming pure circular rotation (as in our 3D barolo analyses) is fit to a ring of gas including NCMs, each term in a Fourier series expansion of the NCMs will cause either an over-or underestimate of the rotation velocity, depending on the alignment between the phase angle and the line of sight.This implies that, provided the sight line is randomly chosen, there is no net preference for an error in the positive or negative direction, regardless of the values of the amplitudes   ′ and   ′ (see Eq. 1).
As noted in our discussion of the right panels of Fig. 9 above, bisymmetric NCMs near the centres of our simulated galaxies are the driver of the bulk of the scatter in  rot and  bar in both the LT and HT models.These cannot, however, account for the fact that the values of these parameters derived from modelling mock observations lie preferentially below (for  rot ) and above (for  bar ) the corresponding 'true' values obtained from the circular velocity curves of the galaxies.

Disc thickness
The approach of most tilted-ring modelling algorithms, including 3D barolo, is to constrain the parameters of each ring one at a time14 .When the model emission for a given ring is compared to the measured emission coming from the corresponding region on the sky, the assumption is that the emission comes from a fixed radial interval within the disc.For a razor-thin disc, the overlap between consecutive rings is minimal, even if a warp is present.For a geometrically thick disc, however, the emission in the region corresponding to a model ring includes mid-plane emission at the corresponding radius, but also emission coming from above and below the mid-plane, and therefore from other radii.
Attempting to make the model rings geometrically thick can further exacerbate the problem: adjacent rings will overlap in the plane of the sky and emission from a given location will be used to constrain multiple rings, even though the model implicitly assumes that the emission is to be decomposed into independent regions corresponding to each ring (see e.g.Iorio et al. 2017, secs. 3.2, 7.1, for further discussion).Ideally, all rings would be optimised simultaneously such that the overlap between geometrically thick rings could be accounted for, but the difficulty of parameter searches in highdimensional spaces with multiple modes and strong degeneracies between parameters have so far precluded this.11: in the LT model, more massive dwarfs (increasing  c,max , paler blue colours) are increasingly centrally baryon-dominated (increasing  bar ) and have increasingly steeply-rising rotation curves (increasing  rot ).In the HT model the trend runs in approximately the opposite direction: more massive dwarfs are increasingly centrally DM-dominated and have more slowly-rising rotation curves.Deviations in rotation curve shape similar to those illustrated in the left panel drive the bulk of the scatter along the diagonal (magenta and yellow arrows, see also Fig. 9).Deviations similar to those illustrated in the centre panel displace the corresponding point horizontally.For LT galaxies, this approximately opposes the intrinsic scalings with  c,max , e.g. an overestimate in  max ( ,max >  c,max ; orange arrow) causes a shift towards a region occupied by galaxies with typically lower values of  c,max (darker blue).For HT galaxies, however, an error in  ,max instead exaggerates the intrinsic scalings with  c,max , explaining the origin of high- ,max mock observations lying at preferentially low- rot and low- bar in the HT model -a feature absent in both the LT model and observations.Iorio et al. (2017, see their sec. 7.1) find that the net effects of 3D barolo's inability to correctly account for disc thickness during modelling are an underestimation of the rotation velocities at inner radii and an overestimation further out.This inner underestimation is clearly seen in the  ,fid values of mock observed points being, on average, lower than the 'true' squares in the left panels of Fig. 11 especially in the LT sample, which has greater disc thicknesses (see Fig. 6).Examples of the systematic overestimation of outer rotation velocities (i.e. ,max >  c,max ) can be seen in Fig. 8 (blue and magenta lines above green and black at high ); this effect is not as strong as the former, however, it is still apparent in the left panels of Fig. 11, where points are, on average, shifted to the left of squares.The net underestimation of  c,fid by  ,fid and overestimation of  c,max by  ,max indicate that disc thickness effects are the main driver of the net underestimation of  rot (and so overestimation of  bar , see Fig. 12) which, in turn, preferentially gives the appearance of DM cores.

Inclination
An error in the overall inclination of a galaxy has a straightforward effect on the rotation curve, whose amplitude is inversely proportional to sin .An overestimate of  therefore leads to both  ,max and  ,fid underestimating  c,max and  c,fid . rot is nearly unchanged, except that the decrease in  ,max reduces  ,fid , such that  c,fid is underestimated by a slightly larger factor than  c,max (see Fig. 12, centre panel), while  bar is overestimated.Global inclination errors are therefore not the primary driver of the scatter of the measurements based on mock observations in Fig. 11.
The influence of radially-localised errors in inclination on correlations between  ,fid ,  ,max ,  rot , and  bar are less easily predicted.However, 'direct' misestimates of the local inclination in tilted ring modelling are relatively rare: tracing warps is, after all, precisely the task these models were designed to fulfill.Instead, local inclination errors most commonly arise in these models as a consequence of the models' inability to capture azimuthal asymmetries in the kinematics or, in other words, NCMs.Harmonic perturbations of order  = 2, in particular, are strongly degenerate with the inclination angle (e.g.Schoenmakers et al. 1997).An azimuthally-symmetric tiltedring model therefore responds to an  = 2 harmonic perturbation by adjusting the inclination angle to minimise the model residuals.
If the perturbation has an amplitude and/or phase angle that vary with radius, the inclination error caused will also vary with radius.We prefer to think of such errors as being due to the non-circular nature of the gas orbits -exactly of the sort described in Sec.3.2 -of which locally-incorrect inclination measurements are merely a symptom.A similar argument applies to degeneracies between the other geometric parameters (centre, systemic velocity, and position angle) and harmonic perturbations of other orders (for an extensive discussion of such degeneracies, see Schoenmakers et al. 1997).

Out-of-equilibrium kinematics
That a rotation curve can be used to draw conclusions about the DM distribution within a galaxy is fundamentally predicated on the assumption that the motions of the kinematic tracers observed (in this context, H i gas) are actually in dynamical equilibrium with the underlying gravitational potential.Simply put, if the gas is not orbiting at the local circular speed, a mass model based on a gas rotation curve loses its meaning.It is therefore important to know whether the H i gas in the simulated galaxies actually orbits at the local circular speed.
We showed two examples in Fig. 8.We recall that the green curves show the median azimuthal velocity of H i-bearing gas particles 15 as a function of radius for galaxies LT-26 and HT-28.In LT-26 this reaches, and follows for several kiloparsecs, the flat portion of the circular velocity curve, while in HT-28, it flattens near the edge of the H i disc, which occurs before the circular velocity fully flattens.In both cases the media azimuthal velocity of the gas is close to the local circular velocity at  fid , but this is not the case for all galaxies in our sample.Similar figures for the full sample of galaxies selected for kinematic modelling are included in Appendix E. Inspecting these by eye, we find that only ∼ 7/21 of the sample drawn from the LT simulation have azimuthal velocity curves following the circular velocity curve at least as closely as the examples shown in Fig. 8, and we would describe only 2/21 of them as following the circular velocity curve closely at all radii within the H i disc.In the sample drawn from the HT simulation, the same fractions are somewhat higher: 10/11 and 5/11, respectively.However, we note that in both cases we have already removed from our sample those galaxies obviously unsuited to kinematic modelling following a visual inspection: 25 from the LT simulation and 32 from the HT simulation.These removed galaxies are preferentially ones where  az is a poor tracer of  circ , and more galaxies are removed for the HT than for the LT model.The fraction of the total population of gas-rich, reasonably isolated dwarfs in these simulations where the H i rotation curve is a good tracer of the circular velocity curve is therefore only about 1 5 (in both simulations), and then only if the tilted ring model accurately recovers the azimuthal velocity of the H i gas.In addition to the geometric thickness of the gas discs (see Sec. 5.1.2),the tendency for  az to underestimate  circ , especially in the inner regions, is a contributor to the systematic underestimates of  c,fid reflected in Fig. 11.
Forming an impression of how many observed galaxies may have H i azimuthal velocities differing significantly from their circular velocities is challenging, chiefly for two reasons.First, our impression is that the origins of such differences are diverse: mergers, interactions with companions, bulk flows driven by supernovae, gas accretion, an elongated potential due to a triaxial DM halo, disc instabilities, and interaction with the intergalactic medium can all disturb the ideal circular flow pattern.Ruling out each of these as a concern is labour intensive and realistically needs to be done on an object-by-object basis.Second, samples of observed galaxies selected for mass modelling are preferentially chosen to at least appear to be undisturbed and close to equilibrium.As outlined above, this mitigates how large a fraction of galaxies may be affected, but seems unlikely to reduce it to a negligible level.These issues are clearly of crucial importance in the interpretation of H i rotation curve measurements; we plan to explore them further in future work (see also Jahn et al. 2023, who find that BICC in their SMUGGLE galaxy formation model leads to  az systematically underestimating  circ ).

Imperfect model optimisation
A source of error essentially unique to our analysis of mock observations, that is, not affecting the observed galaxies to which we compare them, is that we have largely allowed 3D barolo to proceed 15 The green curves are not corrected for pressure support, but the expected corrections are small, typically increasing  az by only 5 to 10 per cent, nearly independently of radius (e.g.Oman et al. 2019, fig. 4).We reiterate that all mock observed rotation curves in this work, on which all our main conclusions rely, are corrected using the 'asymmetric drift correction' computed by 3D barolo.
automatically.That is, we do not revise and refine individual kinematic models to achieve the best possible description of the data in each case, as is customary in analyses of real galaxies.Our decision to largely automate the modelling process is driven primarily by necessity.We wish to consider of order 1000 mock analyses to help bring out more subtle trends.For comparison, the few hundred rotation curves of highly-resolved galaxies in the literature have taken the field decades to produce.We note, however, that as the survey speeds of new telescopes continue to increase, so will the volume of available data to be modelled.Indeed, new surveys rely increasingly on automated kinematic modelling pipelines (e.g.Kamphuis et al. 2015;Ponomareva et al. 2021).
We do not attempt to estimate directly how much such 'automation errors' contribute to the overall scatter and/or systematic trends shown in Fig. 11 (but see Kamphuis et al. 2015).We focus instead on the physical drivers of the scatter and trends.As discussed above, NCM-and disc thickness-driven errors seem to be the dominant effects.

Implications for other galaxy formation models in simulations
We next consider whether the trends for HT and LT galaxies illustrated in Fig. 11 may be generic features of galaxy formation models where BICC does or does not occur in dwarfs, respectively, by comparing against two other models.For this, we use selections of galaxies from the IllustrisTNG and NIHAO simulation suites, without mock observation.The TNG50-1 simulation (Nelson et al. 2019;Pillepich et al. 2019) from the IllustrisTNG suite (Nelson et al. 2018;Pillepich et al. 2018) is a high-resolution cosmological hydrodynamical simulation of a (50 cMpc) 3 volume in a flat ΛCDM (Planck Collaboration et al. 2014) cosmology, performed with the moving mesh code arepo (Springel 2010).It has a baryon mass resolution of 8.5 × 10 4 M ⊙ .The implementations of hydrodynamics, star formation, and supernova feedback are significantly different from those in the EAGLE model on which the LT and HT simulations are based.However, since SF in the TNG model occurs stochastically in regions exceeding a threshold density of  thr = 0.1 cm −3 , BICC does not operate efficiently, and in this sense the model is similar to the LT simulation.There are approximately 4500 dwarf centrals with 50 <  c,max /km s −1 < 120 and  HI > 10 8 M ⊙ in the TNG50 volume; we select 100 of these at random as a representative sample.
The NIHAO 16 project (Wang et al. 2015) consists of a suite of ∼ 100 high-resolution cosmological zoom-in simulations in a flat ΛCDM Planck Collaboration et al. ( 2014) cosmology.The targets of the zoom regions are relatively isolated galaxies.The regions are evolved with a version of the ESF-gasoline2 code (Wang et al. 2015;Wadsley, Keller & Quinn 2017).NIHAO implements both supernova and early stellar feedback (see Stinson et al. 2006Stinson et al. , 2013, respectively), respectively).Star formation follows a Kennicutt-Schmidt SF law for  < 15000 K gas, with a density threshold of  thr = 10.3 cm −3 .This is very close to the value in the HT simulation, and results in BICC in dwarf galaxies (Dutton et al. 2016;Santos-Santos et al. 2018;Dutton et al. 2019;SS20).We select the 35 galaxies from the NIHAO suite with 50 <  c,max /km s −1 < 120.
In the left panel of Fig. 13 we show  rot against  bar for LT, HT, TNG50, and NIHAO dwarfs, in all cases computed from their (spherically-averaged) circular velocity curves.The TNG50 galaxies invariably have steeply rising circular velocity curves ( rot ≳ 0.7), consistent with the presence of DM cusps.In the middle and right panels of the figure, we show  rot and  bar , respectively, as a function of  c,max .This reveals that TNG50 dwarfs have more steeply rising rotation curves across the entire range in  max considered.This may be due to their also having significantly higher central baryon fractions than LT galaxies across the same range in  c,max , which could contract their haloes.Given these qualitative similarities to the LT galaxies, we expect that if the TNG50 galaxies were mock observed, the  rot and  bar values obtained for their rotation curves would likely span a similar range as those in our sample of LT galaxies, offset to preferentially slightly higher  rot and  bar .Our assessment is that this would likely result in covering the range spanned by the observed galaxies similarly to the galaxies from the LT model.It is therefore reasonable to speculate that our conclusions based on the LT model may be generalised to other qualitatively similar (in the sense that BICC does not occur) galaxy formation models.
Whereas the TNG50 and LT galaxies are broadly similar in the context of Fig. 13, comparing the NIHAO and HT galaxies in the same space reveals more pronounced qualitative differences between these two models.There are some similarities: in both models, BICC results in slowly rising rotation curves ( rot ≲ 0.65), except in a few of the most massive NIHAO galaxies ( c,max > 100 km s −1 )this is likely related to their high central baryon fractions.However, whereas HT galaxies have slightly decreasing  rot with increasing  c,max and slightly decreasing  bar with increasing  c,max , both of these trends are increasing for NIHAO galaxies.A consequence is that our finding that the most massive dwarfs in the HT model are preferentially found towards lower  rot and  bar , which is difficult to reconcile with observations, does not apply to NIHAO galaxies.In addition, some NIHAO dwarfs occupy a region around ( bar ,  rot ) ∼ (0.15, 0.6), where no HT galaxies are found.We therefore speculate that if the NIHAO galaxies were observed, the resulting distribution would likely extend further to low  bar and high  rot than that for HT galaxies, and thus better span the range covered by observed galaxies 17 .In summary, it seems as though the conclusions which we draw based on the HT model cannot be generalised to all galaxy formation models where BICC occurs.

Correlation between core size and stellar properties of galaxies
Other authors have highlighted that the size of cores created through BICC in simulations correlate with the stellar mass (Governato et al. 2012), stellar-to-halo mass ratio (Di Cintio et al. 2014) or star formation history (Read et al. 2016), with the largest cores found in simulated galaxies with  ★ ∼ 10 8 − 10 9 M ⊙ ,  ★ / 200 ∼ 10 −2 , and/or the most temporally extended star formation histories.If such correlations are clearly measured in observed galaxies, this will lend weight to interpretations involving BICC.While models parametrising these effects provide good descriptions of simulated galaxies with baryon-induced cores (references above, and see also Tollet et al. 2016;Benítez-Llambay et al. 2019;Lazar et al. 2020), and it is clear that they achieve a good description of some observed galaxies (e.g.Oh et al. 2015;Read et al. 2016Read et al. , 2019)), they so far fail to provide a comprehensive explanation for their overall scatter in central DM density slope, central rotation speed, or other measures of central DM content for galaxies in the regime where core formation peaks in these models (Santos-Santos et al. 2018;Relatores et al. 2019;Santos-Santos et al. 2020;Frosst et al. 2022, andsee Sales et al. 2022, for a review).This leaves space for additional scatter due to effects such as those discussed in Sec.5.1.
velocity curves of NIHAO galaxies do not span the ranges observed, for instance near  rot ∼ 1.0 and  bar ∼ 0.15 (see also SS20, sec.4.2.2).Some additional source of scatter is therefore needed if this model is to explain the observations.

SUMMARY AND CONCLUSIONS
The diverse kinematics of observed local dwarf irregular galaxies is difficult to replicate directly in ΛCDM hydrodynamical simulations (Oman et al. 2015).SIDM has been proposed as a possible explanation (see e.g.Kaplinghat et al. 2020).More commonly, gravitational coupling of violent gas motions to the DM (BICC; Navarro et al. 1996a;Read & Gilmore 2005;Pontzen & Governato 2012) have been invoked to explain the observed diversity (e.g.Di Cintio et al. 2014;Brook 2015;Chan et al. 2015;Tollet et al. 2016;Read et al. 2016;Pacheco-Arias et al. 2021;Jahn et al. 2023), with some concluding that this may not be sufficient (e.g.Bose et al. 2019;Benítez-Llambay et al. 2019).However, to our knowledge, no ΛCDM or SIDM galaxy formation model has yet been identified that can reproduce a realistic galaxy population including late-type dwarf galaxies with both the most steeply-and the most slowlyrising rotation curves observed.The observation that more centrally baryon-dominated late-type dwarfs seem to have more slowly-rising rotation curves (as quantified by the parameters  rot and  bar ), highlighted by SS20, proves to be especially challenging for BICC models to reproduce.This discrepancy has been framed as a major challenge to the fiducial ΛCDM cosmology (e.g.Flores & Primack 1994;Moore 1994;Oman et al. 2015;Bullock & Boylan-Kolchin 2017;Sales et al. 2022).However, recent studies have shown that systematic errors in the kinematic modelling of late-type dwarfs may be key to resolving this discrepancy (Read et al. 2016;Iorio et al. 2017;Pineda et al. 2017;O19;Genina et al. 2020).O19 and SS20 showed that discrepancies between the rotation curves and circular velocity curves of simulated late-type dwarfs induced primarily by the presence of non-circular flows in their H i discs can plausibly reproduce the observed kinematic diversity, as well as the observed trend in  rot and  bar , for galaxies with DM cusps drawn from the APOSTLE simulations.In the present study we build on this work by extending a similar analysis to many more simulated galaxies created using a galaxy formation model similar to APOSTLE (with low SF density threshold; LT), and a second model where BICC creates DM cores (with high SF density threshold; HT).
The two galaxy formation models are based on the EAGLE project (Crain et al. 2015;Schaye et al. 2015), with the gas density threshold for star formation modified to control the creation of DM cores (Benítez-Llambay et al. 2019).We select galaxies with maximum circular velocities between 50 and 120 km s −1 , H i masses exceeding 10 8 M ⊙ , and which are isolated from massive companions.By visual inspection, we remove those clearly disturbed which are unsuited to kinematic analysis, leaving a sample of galaxies with properties similar to those of observed late-type dwarfs.The simulated galaxies are mock observed from multiple viewing angles to produce H i spectroscopic data cubes which are then analysed using a conventional 3D tilted-ring modelling tool in order to extract rotation curves.
We quantify rotation curve shapes by their amplitudes at outer radii ( ,max ) and at an inner radius ( ,fid ), and the ratio of the two ( rot ).In addition, we quantify the dynamical importance of baryons near the galactic centre by the baryonic-to-total mass ratio in this region ( bar ).Our main conclusions are summarised as: (i) The diversity in the rotation curve shapes as measured from the 'true' circular velocity curves of galaxies (i.e.directly from the simulation outputs) from both the LT and the HT models falls far short of the observed diversity.
(ii) Assuming pure circular motion in galactic H i discs in the presence of non-circular motions with realistic amplitudes introduces significant scatter in the recovered rotation curves for different viewing angles of the same galaxy.This scatter is similar in amplitude to the observed diversity in  ,fid at fixed  ,max , and naturally leads to a correlation of the form  rot ∝  −0.5 bar .This scatter due to non-circular motions is symmetric around the true values of  rot and  bar .
(iii) Assuming too thin an H i disc typically causes an underestimation of the inner circular velocity curve ( ,fid <  c,fid ) and an overestimation of the outer circular velocity curve ( ,max >  c,max ), creating a bias towards slowly rising rotation curves (low  rot ) and centrally high baryon fractions (high  bar ).This bias is stronger in the LT galaxies, which have thicker H i discs compared to the HT galaxies in our sample.
(iv) The combination of these two effects, which are dominant in our analysis of simulated galaxies, means that it is much more common for a galaxy with a DM cusp to appear to have a core than vice versa.It is likewise unusual, albeit not impossible, for simulated galaxies with cores to appear to have high central DM fractions.If the rotation curves of observed late-type dwarfs are similarly impacted -which seems likely -a late-type dwarf population where most galaxies have DM cusps is easier to reconcile with the available data than one where most galaxies have sizeable DM cores.
(v) The mock observations of galaxies from the HT model with higher  ,max lie preferentially at lower  rot and/or lower  bar .No such trend is apparent in the LT model, or in observed late-type dwarfs.This lends some additional weight to a scenario where most late-type dwarfs host DM cusps.
(vi) An important caveat is that the HT simulation has significant relevant qualitative differences with respect to another galaxy formation model where BICC occurs (NIHAO) -our conclusions based on this model are therefore unlikely to hold for all such models.The LT simulation, however, is broadly similar to another simulation where BICC does not occur (TNG50), suggesting that the conclusions based on this model may apply more generally to such models.
Our results demonstrate that the possibility that the measured rotation curves of late-type dwarfs differ significantly from their circular velocity curves due to a combination of non-circular motions and geometrically thick H i discs must be taken seriously.Indeed, these effects may plausibly be the dominant contributors to the scatter observed in the rotation curve shapes of these objects regardless of whether their DM haloes have central DM cusps, or cores.Our analysis shows that a scenario where all dwarf galaxies have central DM cusps remains viable; it may even be preferred.

Figure 1 .
Figure 1.Upper panel:The spherically averaged DM density profiles of the 21 dwarf galaxies selected for kinematics analysis from the LT simulation are shown with blue solid lines.The lines are thinner inside the 'convergence radius'(Power et al. 2003, eq.20, with  relax ( ) ≳ 0.6 0 ), outside which the circular velocity curve converges to better than ∼ 10 per cent.Two NFW density profiles, with  c,max = 50 and 120 km s −1 , and the median concentration given their mass(Ludlow et al. 2016) are shown with dashed black lines.Lower panel: As upper panel but for the 11 dwarfs selected from the HT simulation, shown in red.LT dwarf haloes retain their  ∼  −1 cusps and therefore closely follow the NFW profile shape at all resolved radii.HT dwarf haloes, however, develop  ∼  0 cores with typical sizes of ∼ 2-3 kpc.

Figure 2 .
Figure 2. Visualisations of the H i gas distribution in simulated galaxies LT-26 (upper panels) and HT-28 (lower panels).Left panels: H i surface density, Σ HI , calculated directly from the simulation particle distribution.The spatial resolution is chosen for visual clarity.Centre panels: 21-cm flux density  HI , from mock observed data cubes.Right panels: Intensity-weighted mean velocity, from mock observed data cube.The velocity map is masked to show only pixels where the flux density exceeds 10 −3 Jy beam −1 (indicated by the magenta contour in  HI ; approximately Σ HI = 1 M ⊙ pc −2 ≈ 10 20 atoms cm −2 ).Velocity contours correspond to tick locations on the colour bar.Both systems are mock observed at an inclination of  = 60 • (LT-26 appears nearly face-on due to a warped outer disc, see Sec. 2.3), projection angle of Φ = 0 • , position angle of 270 • , and assuming a distance of 12.9 Mpc.The mock observations are convolved with a 12 arcsec (750 pc) FWHM circular Gaussian beam.

Figure 3 .
Figure 3. Left panel:The baryonic Tully-Fisher relation (BTFR) for all simulated galaxies in the low-SF density threshold (LT) model (blue symbols), and the high threshold (HT) model (red symbols).Galaxies with both maximum circular velocity 50 <  c,max /km s −1 < 120 and total H i mass  HI > 10 8 M ⊙ are indicated with larger symbols (points and squares); galaxies selected for mock observation and kinematic analysis are indicated with square symbols.The region outside of this velocity range is shaded grey.The baryonic mass is defined as  bar =  ★ + 1.4 HI .A selection of observed galaxies from the compilation of SS20 are shown with crosses, with larger symbols used for those within the same range in  ,max .Centre panel: H i mass-stellar mass relation; symbols are as in the left panel.The region corresponding to  HI < 10 8 M ⊙ is shaded grey.Right panel: H i mass-size relation; symbols and region shading are as in the centre panel. HI is defined as the galactic radius at which the H i surface density falls to 1 M ⊙ pc −2 .

Figure 5 .Figure 6 .
Figure5.The strength of the bisymmetric NCMs (traced by the  ′ = 3 harmonic, see Eq. 1) relative to the rotation speed (traced by the  ′ = 1 harmonic), both measured at a radius of 2 kpc, as a function of the inclination-corrected amplitude of the rotation speed at the same radius.Measurements for mock observed galaxies from the LT and HT simulations selected for kinematic modelling are shown with blue and red points, respectively.Measurements for selected observed galaxies from the THINGS and LITTLE THINGS surveys, also at a radius of 2 kpc, are shown with black crosses.We show two measurements for each simulated galaxy, for projection angles of Φ = 0 • and 90 • .Error bars show formal errors on  3 ′ / 1 ′ (seeMarasco et al. 2018, for details), and arrows indicate upper limits.

V
max [km s 1 ]

Figure 11 .
Figure11.Upper-left: Inner 'fiducial' velocity,  ,fid , (see Eq. 2 and Fig.8) as a function of maximum velocity,  ,max , for a selection of galaxies observed in H i with 50 <  ,max /km s −1 < 120 shown as crosses, and unselected galaxies with  ,max above and below this range shown as upward-and downward-pointed triangles, respectively (see Sec. 2.2.1).Cross markers are coloured by  ,max (redundant with the horizontal axis, but included for consistency with the upperright panel).The dashed black line shows  fid =  max , while the solid black line shows the relation for an average NFW profile assuming the mass-concentration relation ofLudlow et al. (2016), with 10 th -90 th percentile scatter shown as a grey band.Upper-right: Rotation curve shape parameter,  rot , as a function of central baryon-to-total mass fraction  bar (see Sec. 4.2).Symbols and colours are as in the upper-left panel.The solid black line shows the average  rot expected for NFW cusps and the grey band shows corresponding 10 th -90 th percentile scatter.Centre-left: Similar to upper-left panel; lines and cross symbols (without colouration) are repeated.The square symbols are measured from the circular velocity curves ( c,max ,  c,fid , derived from the gravitational acceleration field in the disc midplane) of the 21 simulated LT galaxies selected for kinematic modelling.For each galaxy, there are 24 points for the 24 rotation curves ( ,max ,  ,fid ), one for each mock observation orientation, with points corresponding to Φ = 0 • given black borders.These give an impression of the range in rotation curve shapes which might be measured given the intrinsic shapes reflected by the square markers.Darker and lighter error bars show median fractional uncertainties (for Φ = 0 • mock observations) for galaxy subsamples with  c,max < 75 km s −1 and  c,max ≥ 75 km s −1 , respectively.Centre-right: Similar to upper-right; lines and cross symbols (without colouration) are repeated.Blue squares and points are as in centre-left panel.Lower-left: As centre-left, but for the 11 selected, simulated HT galaxies.Lower-right: As centre-right, but for selected HT galaxies.

Figure 12 .
Figure12.A schematic representation of the connection between rotation curve shape and the parameters  rot and  bar .Left: An overestimate (magenta) or underestimate (yellow) of the inner rotation curve, such as might be induced by NCMs near the galactic centre.Since  ,max is unaffected (and  b,fid is assumed held fixed), the corresponding point in the  rot - bar plane is deflected along a line  rot ∝  −0.5 bar (inset panel).Centre: A global overestimate (orange) or underestimate (cyan) of the circular velocity curve, such as might be induced by a constant inclination error across the entire galactic disc.Since  ,fid and  ,max change approximately proportionally to one another,  rot remains approximately constant, while  bar is inversely proportional to  2 ,fid (inset panel).Right: The coloured patches are schematic representations of the distributions of square symbols in the centre-right and lower-right panels of Fig.11: in the LT model, more massive dwarfs (increasing  c,max , paler blue colours) are increasingly centrally baryon-dominated (increasing  bar ) and have increasingly steeply-rising rotation curves (increasing  rot ).In the HT model the trend runs in approximately the opposite direction: more massive dwarfs are increasingly centrally DM-dominated and have more slowly-rising rotation curves.Deviations in rotation curve shape similar to those illustrated in the left panel drive the bulk of the scatter along the diagonal (magenta and yellow arrows, see also Fig.9).Deviations similar to those illustrated in the centre panel displace the corresponding point horizontally.For LT galaxies, this approximately opposes the intrinsic scalings with  c,max , e.g. an overestimate in  max ( ,max >  c,max ; orange arrow) causes a shift towards a region occupied by galaxies with typically lower values of  c,max (darker blue).For HT galaxies, however, an error in  ,max instead exaggerates the intrinsic scalings with  c,max , explaining the origin of high- ,max mock observations lying at preferentially low- rot and low- bar in the HT model -a feature absent in both the LT model and observations.
Gas specific angular momentum,  gas =  HI , as a function of gas mass,  gas = 1.4HI , for galaxies from the LT and HT simulations selected for kinematic analysis, shown with blue and red square markers, respectively, and excluded galaxies, shown with points of the same colours.The dashed black line shows a best-fitting  gas - gas relation found byMancera Piña  et al. ( Ludlow et al. (2016)d lines.The solid black line shows the relation for an average NFW halo model assuming the mass-concentration relation ofLudlow et al. (2016), and the shaded band illustrates the 10 th -90 th percentile scatter around this relation.The dashed black line shows  fid =  max .Centre panel: Rotation curve shape parameter,  rot , as a function of central baryon-to-total mass fraction,  bar (see Sec. 4.2.2).The solid black line shows the average  rot expected of NFW cusps and the grey band shows the approximate 10 th -90 th percentile scatter around this value.Symbol shapes and colours, and connecting lines, are as in the left panel.Right panels:  bar and  rot as functions of viewing angle, Φ. Symbols and the solid lines have the same meaning as in the other panels.Dashed lines show best fitting sinusoidal functions with a 180 • period, revealing the systematic effect bisymmetric NCMs have on  bar and  rot in both models.
Figure13.Left panel: Rotation curve shape parameter,  rot , as a function of central baryonic matter fraction,  bar , both determined from the (sphericallyaveraged) circular velocity curves.Galaxies from the LT and HT simulations selected for kinematic analysis are shown with blue and red square markers, respectively, while galaxies not selected are shown with points of the same colours.A selection of galaxies from the NIHAO simulation suite is shown with orange upward-pointed triangles, and a selection from the TNG50-1 simulation with cyan downward-pointed triangles.The LT and TNG50 populations occupy broadly similar regions of this parameter space, as do the HT and NIHAO populations.Centre panel:  rot as a function of maximum circular velocity,  c,max .Symbols are as in the left panel.LT and TNG50-1 galaxies have similarly low SF density thresholds, and all have steeply-rising rotation curves suggestive of DM cusps ( rot ≳ 0.65).HT and NIHAO galaxies have similar, higher SF density thresholds, and often have slowly-rising rotation curves, suggestive of DM cores.Right panel:  bar against  c,max , both from circular velocity curves.Symbols are as in the other panels.The galaxy population in the LT simulation has a strong positive gradient in this space.The populations from both TNG and NIHAO have weaker positive gradients, while the population from the HT simulation has approximately constant  bar across the range in  c,max shown.