The many reasons that the rotation curves of low-mass galaxies can fail as tracers of their matter distributions

It is routinely assumed that galaxy rotation curves are equal to their circular velocity curves (modulo some corrections) such that they are good dynamical mass tracers. We take a visualisation-driven approach to exploring the limits of the validity of this assumption for a sample of $33$ low-mass galaxies ($60<v_\mathrm{max}/\mathrm{km}\,\mathrm{s}^{-1}<120$) from the APOSTLE suite of cosmological hydrodynamical simulations. Only $3$ of these have rotation curves nearly equal to their circular velocity curves at $z=0$, the rest are undergoing a wide variety of dynamical perturbations. We use our visualisations to guide an assessment of how many galaxies are likely to be strongly perturbed by processes in several categories: mergers/interactions (affecting $6$/$33$ galaxies), bulk radial gas inflows ($19$/$33$), vertical gas outflows ($15$/$33$), distortions driven by a non-spherical DM halo ($17$/$33$), warps ($8$/$33$), and winds due to motion through the IGM ($5$/$33$). Most galaxies fall into more than one of these categories; only $5$/$33$ are not in any of them. The sum of these effects leads to an underestimation of the low-velocity slope of the baryonic Tully-Fisher relation ($\alpha\sim 3.1$ instead of $\alpha\sim 3.9$, where $M_\mathrm{bar}\propto v^\alpha$) that is difficult to avoid, and could plausibly be the source of a significant portion of the observed diversity in low-mass galaxy rotation curve shapes.


INTRODUCTION
Since the discovery of flat rotation curves in galaxies (Rubin et al. 1980;Bosma 1981) leading to the widespread acceptance of dark matter (DM) theories, rotation curves have been used to study DM. Low-mass galaxies, with maximum circular velocities 120 km s −1 , are particularly well suited for such analysis because their high DM mass fractions reduce the relative gravitational influence of baryons, so that their circular velocity almost directly traces their DM content. The baryonic Tully-Fisher relation (BTFR;McGaugh et al. 2000) provides a concise summary of this trend: the baryonic (gas plus stellar) mass of galaxies is observed to be proportional to about the fourth power of their maximum rotation velocities, bar ∝ 4 max (but see Ponomareva et al. 2018), but a constant baryon-to-DM mass ratio would instead imply a shallower slope close to bar ∝ 3 max (e.g. Sales et al. 2017). The slope and scatter of the BTFR for the lowest mass galaxies ( bar 10 9 M ), however, remain challenging to constrain (Sorce & Guo 2016;Papastergis et al. 2016;Bradford et al. 2016;Verbeke et al. 2017;Ponomareva et al. 2018;Mancera Piña et al. 2019;Lelli et al. 2019;Wingfield McQuinn et al. 2022;Ball et al. 2022, andLelli 2022, a review), and leaves the connection between the luminous components of galaxies and the DM haloes in which they form at the low-mass edge of galaxy formation uncertain (Trujillo-Gomez et al. 2011;Desmond 2012;Papastergis et al. 2015;Brook et al. 2016;Oman et al. 2016;Brooks et al. 2017;Sales et al. 2017).
Studies of dwarf galaxies have revealed several potential problems in near-field cosmology (see Bullock & Boylan-Kolchin 2017;Sales et al. 2022, for reviews). One such problem that remains unresolved is the 'cusp-core' problem (Flores & Primack 1994;Moore 1994;de Blok 2010); the inner slopes of low-mass galaxy rotation curves are often slowly rising compared to the mass profile implied by the steep central density 'cusps' predicted by N-body simulations (Navarro et al. 1996b).
There have been many proposed resolutions of the cusp-core problem within the ΛCDM framework. One such proposal is that gas flows driven by supernova feedback couple gravitationally to the DM and re-distribute it, producing and maintaining a central density 'core' (e.g., Navarro et al. 1996a;Read & Gilmore 2005;Pontzen & Governato 2012, and see Pontzen & Governato 2014, for a review). The 'bursty' star formation histories arising in some galaxy formation simulations produce cores in a limited mass range (e.g., Di Cintio et al. 2014;Chan et al. 2015;Tollet et al. 2016;Jahn et al. 2021), and the conditions necessary for core formation via this mechanism are now well-understood (Bose et al. 2019;Benítez-Llambay et al. 2019). However, whether such effects can fully reproduce the diverse rotation curves observed for dwarf galaxies remains unclear (e.g. Oman et al. 2015;Santos-Santos et al. 2020;Roper et al. 2022).
Another proposed scenario involves allowing cold DM particles to scatter from each other, leading to heat transfer to the inner regions of DM haloes and redistributing the DM to produce a core (Spergel & Steinhardt 2000). Such 'self-interacting dark matter' (SIDM) models inherit the large scale successes of the standard ΛCDM model, and are able to produce a range of rotation curve shapes by including the gravitational influence of baryons, which can re-form a cusp (see Tulin & Yu 2018, for a recent review). This shows promise (e.g., Ren et al. 2019;Kaplinghat et al. 2020) however, again, concerns whether SIDM can account for the full observed diversity remain (Creasey et al. 2017;Santos-Santos et al. 2020).
More prosaically, the problem could be that the circular velocity curves of low-mass galaxies are not accurately measured by the methods used to extract them from observations. The inclination angle (possibly varying with radius), non-circular motions, potentially anisotropic velocity dispersion, and geometrically thick and/or flared nature of gas discs are just some of the challenging issues that models in principle need to account for to accurately measure a rotation curve. Strong degeneracies between parameters describing the geometry and kinematics of a gas disc further complicate matters. Attempts to model realistic galaxies with known rotation curves have revealed that the errors due to these issues can be quite severe (Read et al. 2016;Oman et al. 2019;Roper et al. 2022), although Frosst et al. (2022) argue that such effects may still fall well short of explaining the observed diversity in rotation curve shapes. There is, however, an even more worrying possibility: that the rotation curves of low-mass galaxies may in some cases not faithfully trace their circular velocity curves 1 . In this case even a perfectly accurate measurement of the rotation curve gives no meaningful information about the total matter distribution within a galaxy.
It is clear that some low-mass galaxies are not in dynamical equilibrium, and thus that their rotation curves are not reliable tracers of their circular velocity curves (and consequently of their DM content). Obvious perturbations, such as mergers or star formation-driven 'superbubbles', are easily identified, however low-mass galaxies' shallow gravitational potential wells make them especially susceptible to additional perturbations which may not be so obvious. How often these more subtle physical processes may cause departures from equilibrium in these objects remains almost unexplored in the literature (see Hayashi & Navarro 2006 on the effect of a triaxial DM halo; Valenzuela et al. 2007 on the influence of lopsided gas discs; Read et al. 2016 on the influence of the star formation cycle; Verbeke et al. 2017 sec. 4.1 for a brief exploration of the topic).
In this work we make an initial assessment of the relative importance of different types of perturbations using a sample of galaxies with maximum circular velocities 60 < max /km s −1 < 120 from the APOSTLE suite of cosmological hydrodynamical simulations. We create visualisations of the galaxies and compute their rotation and circular velocity curves at a range of times over the past ∼ 4 Gyr. We use these to investigate the kind of perturbations that affect lowmass galaxies, their frequencies, their effects on the galaxies' rotation curves, and what conditions are necessary for galaxies to actually rotate at their circular speeds.
We begin in Section 2 with a brief description of the APOSTLE simulations and our methods for calculating rotation curves and producing visualisations. In Section 3, we present our main results: we describe the perturbations affecting galaxies in our sample, and investigate their influence on key galaxy scaling relations. We summarize our conclusions and discuss their implications and applicability to real galaxies in Section 4.

The APOSTLE simulations
The APOSTLE 2 simulations Fattahi et al. 2016) are a suite of zoom-in cosmological hydrodynamical galaxy formation simulations. The suite is made up of 12 regions selected to resemble the Local Group of galaxies in terms of the masses, separation and kinematics of a pair of galaxies analogous to the Milky Way and Andromeda, and a lack of other massive galaxies within a few megaparsecs. A region about 2 − 3 Mpc in radius around each pair was simulated at multiple resolution levels (lowest 'L3' to highest 'L1') with the 'Reference' calibration (Crain et al. 2015) of the EA-GLE galaxy formation model (Schaye et al. 2015). The model is implemented using a smoothed-particle hydrodynamics framework in the pressure-entropy formulation (Hopkins 2013) and includes prescriptions for radiative cooling (Wiersma et al. 2009a), star formation (Schaye 2004;Schaye & Dalla Vecchia 2008), stellar and chemical enrichment (Wiersma et al. 2009b), thermal-mode stellar feedback (Dalla Vecchia & Schaye 2012) and cosmic reionisation (Haardt & Madau 2001;Wiersma et al. 2009b). The feedback from supermassive black hole accretion implemented in the EAGLE model has a negligible effect on the galaxies in the APOSTLE simulations . The simulations assume the WMAP-7 cosmological parameters (Komatsu et al. 2011).
Galaxies are identified in the simulations following a two-step process. First, particles are linked together by a friends-of-friends (FoF) algorithm (Davis et al. 1985). Each FoF group is independently analysed using the S halo finding algorithm (Springel et al. 2001;Dolag et al. 2009) which identifies gravitationally bound substructures. The subhalo with the minimum gravitational potential in each FoF group is labelled the 'central' galaxy of the group, while others are labelled 'satellites'. We label galaxies from the APOSTLE simulations following the same convention as Oman et al. (2019): for example, AP-L1-V6-5-0 refers to APOSTLE resolution level L1, region (volume) V6, FoF group 5, subhalo 0 (the 'central' subhalo). We always refer to the identifier of the galaxy in the last snapshot; its progenitor(s) may have different identifiers. We track the progenitors of galaxies in our sample back through time using the merger tree algorithm of Helly et al. (2003). When a galaxy has more than one progenitor at a previous time, we follow the progenitor that contributed the most particles to the descendant.
In this work we focus on a sample drawn exclusively from the highest-resolution (L1) simulations in the suite. Only regions V1, V4, V6, V10 and V11 have been simulated at this resolution. At L1 resolution, the gas (dark matter) particle mass is typically 7.4 × 10 3 M (3.6 × 10 4 M ), and the maximum gravitational softening length is ≈ 134 pc. According to the criterion of Power et al. (2003), the circular velocity curves of low-mass galaxies at this resolution level are numerically converged to better than 10 per cent at radii 700 pc. We focus on recent times, between 8.9 Gyr and 13.76 Gyr ( = 0). In this period there are 17 full simulation outputs (every ∼ 0.3 Gyr), or 'snapshots', and 147 partial outputs (every 34 Myr), or 'snipshots', where some detail -such as abundances of individual elements -is omitted.
Our sample of galaxies is the same as that used by Oman et al. (2019). The galaxies are selected to have maximum circular velocities 60 < max /km s −1 < 120, to be centrals (not satellites), and to be found in FoF groups which do not include any contaminating low-2 A Project Of Simulating The Local Environment resolution particles from outside the nominal zoom-in regions of the simulations. There are 33 such galaxies, with 4 found in simulation region V1, 5 in V4, 11 in V6, 10 in V10, and 3 in V11. All are at separated from the nearest of the pair of galaxies analogous to the Milky Way and M 31 by at least 450 kpc, and up to 3 4 Mpc.

Circular velocity and rotation curves
We calculated the total circular velocity curves of galaxies in our sample as circ = √︁ ( (< ))/ , where is the gravitational constant and (< ) is the mass enclosed within radius of the location of the particle with the minimum gravitational potential, including all particle types (DM, gas, stars, and black holes). The spherically symmetric approximation is reasonable for our sample of galaxies, which are invariably DM-dominated both globally and locally at all radii. Furthermore, as will be seen below, the actual rotation curves preferentially underestimate the (spherically averaged) circular velocity curves, so the reduction in circ by a few per cent (Binney & Tremaine 2008, sec. 2.6.1b) due to this approximation tends to slightly underestimate differences between the two 4 .
Before calculating rotation curves, we set the velocity zero point of each galaxy to the mean velocity of its 100 innermost 'atomic' gas particles. We define atomic gas particles as those with H mass fractions of greater than 0.5. The H mass fractions are calculated as detailed in Oman et al. (2019) -in brief, these assume the empirical prescription of Rahmati et al. (2013) to compute the neutral 3 The 'zoom-in' region has an irregular shape and can extend beyond the nominal radius of 2−3 Mpc. The condition that no low-resolution particles are present in the FoF group ensures that the galaxies in our sample are sufficiently far from the boundary of the 'zoom-in' region to avoid any spurious numerical effects. 4 In the central ∼ 1 kpc, the circular velocity curves are also underestimated due to the limited resolution of the simulations, but that effect actually modifies (rather than misestimates) the gravitational potential, and the rotation curve should be expected to agree with the gravitational potential actually realised in the simulations. fractions of particles, and the relation given in Blitz & Rosolowsky (2006) to partition atomic from molecular gas. We then calculate the angular momentum vector of the atomic gas disc by summing the angular momenta of the innermost 50 per cent of atomic gas particles (or 125,000, whichever is fewer). We rotate the coordinate frame so that the angular momentum vector points in the -direction, placing the disc in the -plane. We measure the median azimuthal velocity of atomic gas particles gravitationally bound to the galaxy within cylindrical annuli of 0.5 kpc width. This bin width offers a good compromise between limiting noise in the measurement and resolving the structure in the rotation curves. We measure the rotation curves out to the edge of the atomic gas disc disc , which we define as the radius enclosing 90 per cent of the H mass.
The rotation curves are not corrected for a possible radial pressure gradient in the gas disc (often incorrectly termed an 'asymmetric drift correction', see e.g. Valenzuela et al. 2007, appendix A). Such corrections for our sample of galaxies (at = 0) were computed by Oman et al. (2019) and are invariably small ( 10 per cent), except for during mergers. We have checked that the corrections are no larger even at the earliest times that we consider in our analysis. Since we focus below on links between visible (Sec. 2.3) gas kinematic features and rotation curve features we omit further discussion of pressure-support corrections for simplicity, but note that we do not expect that accounting for these would qualitatively change any of our conclusions as deviations of our galaxies' rotation curves from their circular velocity curves are well in excess of 10 per cent in most cases (see Sec. 3.1,below).
This process was repeated for the 17 snapshots between 8.88 Gyr and 13.76 Gyr to produce a set of circular velocity and rotation curves over time, for each of the 33 galaxies in our sample. Fig. 1 shows the resulting curves at each snapshot for the galaxy AP-L1-V6-5-0, as an illustrative example.

Images and videos
We use the -(Benitez-Llambay 2015) toolkit to create videos of galaxies in our sample over time to explore the kinds of perturbations that affect them and their effects on their rotation curves.
In -, the 'observer' is referred to as the 'camera'. The parameters specifying the camera position and orientation are 'anchored' at the times corresponding to snapshots. The camera is pointed at the centre of the galaxy of interest (defined as the location of the particle with the minimum gravitational potential), and placed at a distance such that an image with a 90 • field of view extends to about twice the radius of the atomic gas disc. We track both a 'faceon' view camera offset from the centre along the angular momentum vector of the disc (see Sec. 2.2), and an 'edge-on' view camera offset along an arbitrarily chosen orthogonal axis. Each galaxy is visualised at 383 times evenly spaced between 8.88 and 13.76 Gyr. Since the time when a visualisation is to be created does not in general correspond to the time of a snapshot or snipshot, particle positions are linearly interpolated between the two simulation outputs closest to the desired time. The parameters describing the camera position and orientation are also linearly interpolated to the desired time. The close spacing of the snipshots in time means that a higher-order interpolation scheme is not necessary. Finally, the normalisation of the colour scale of the images is linearly interpolated between the maximum pixel values in the first and last image (and likewise for the minima) in each series to prevent 'flickering' and over/undersaturation of the images. To focus attention on the object of interest, the contributions of simulation particles more than 50 kpc from the Selected frames from the face-on video of the galaxy AP-L1-V6-5-0, showing a gas-rich merger which strongly disrupts the gas disc and the rotation curve around ≈ 10.5 ± 0.5 Gyr (shown in Figure 1). The partially stripped, but still gas-rich, secondary halo has a second approach, once again disturbing the galaxy, around ≈ 11.7 ± 0.1 Gyr. The first two images show the DM density only (grey-scale) and gas density only (purple-orange colour map) respectively, for the same time as shown in the third panel. The further images are composites of DM density and gas density. On compatible pdf viewer software a video will play before the figure is displayed. It shows the evolution of the galaxy over ∼ 4 Gyr with side-by-side face-on (left) and edge-on (right) views of the galaxy showing the DM and gas density composite visualisation. The same video is available in the supplementary materials as AP-L1-V6-5-0-composite-edge-and-face.mp4 (see Appendix A).
centre of the object of interest is exponentially suppressed with a scale length of 50 kpc, such that anything beyond ∼ 300 kpc is essentially invisible.
We use this procedure to create videos visualising the galaxy faceon and edge-on for DM and gas particle types, and assemble these in a variety of combinations (e.g. face-on and edge-on with composite DM plus gas images; face-on showing DM and gas particles side by side, edge-on showing DM and gas particles side by side) to create an information-rich set of videos for each galaxy. We also produce a set of figures for each galaxy showing its circular velocity and rotation curve at the time of each snapshot side-by-side with an image of the galaxy at the same time. Further details are given in the Appendix. Fig. 2 shows a few example frames from a DM-plus-gas composite face-on view video (on compatible software the video itself will be shown before the figure is displayed) for the galaxy AP-L1-V6-5-0 (the same galaxy as in Fig. 1), showing a gas-rich merger.

RESULTS
We examined the videos and rotation curves for each galaxy in detail. We noted the characteristics of each galaxy, the types of perturbations visibly affecting each, and their effects on the galaxy and its rotation curve over time. We allowed the qualitative impressions formed during this process to guide the creation of a quantitative summary of the different types of perturbations affecting the galaxies. Our thoughts regarding the advantages of this approach are summarised in Sec. 4.2 below.

Quality of the rotation curve as a circular velocity tracer
In Fig. 3, we show the rotation curves at the times of the last 3 snapshots (13.10, 13.43, 13.76 Gyr) of each galaxy, as well as their circular velocity curves at the time of the last snapshot. The circular velocity curve (purple line) at the times of the two preceding snapshots is invariably very similar to that at the time of the last snapshot, so we omit them from the figure. Some rotation curves accurately trace the circular velocity curve, while others do not. Likewise, some galaxies have rotation curves that are highly variable over the ∼ 660 Myr spanned by the three snapshots, while others are quite stable. Guided by our visual impression of the curves in Fig. 3, we devised a summary statistic that captures these features. It is defined: where: circ,0 is the circular velocity curve at the time of the last snapshot, rot,0 , rot,1 and rot,2 are the rotation curves at times of the last, second-last and third-last snapshots, respectively, are the radii where the curves are sampled, and 0.75 (·) denotes the 75 th percentile. The radii are evenly spaced every 500 pc out to the radius enclosing 90 per cent of the H mass of the galaxy -since this varies with time, pairs of curves are compared using sampling points common to the pair. Conceptually, 0 measures how well the rotation curve traces the circular velocity curve (at the time of the final snapshot), while 0,1 and 1,2 measure the time variability of the rotation curve. In all cases, smaller values indicate better agreement. The 75 th percentile is used to enforce that 'agreement' between two curves must extend over most of the curves ( 3 4 of their extent) to Table 1. Summary of perturbations affecting the gas kinematics in our sample of galaxies. The rows are in order of increasing parameter (column 3) (Sec. 3.1) quantifying how closely the rotation curve traces the circular velocity curve; higher indicates poorer agreement. The range in is separated into 4 classes (column 1) from class 1, 'excellent agreement', to class 4, 'poor agreement'. The remaining columns provide quantitative estimates of the strength of various perturbations with entries corresponding to a nominal 'strong perturbation' regime shown in bold. Further details are given in the specified sections. Column (4): Time(s) since the big bang of the first pericentric passage of companions with DM mass ratio greater than 1:20. Currently strongly interacting companions are marked ‡, and the entire entry is shown in bold (Sec. 3.2.1). (5): Peak (most negative) bulk cylindrical radial atomic gas inflow rate during the last ∼ 600 Myr, values < −5 kpc Gyr −1 in bold (Sec. 3.2.2). (6): Peak bulk vertical (sgn( ) ) atomic gas expansion rate during the last ∼ 600 Myr, values > 1 kpc Gyr −1 in bold (Sec. 3.2.2).  AP-L1-V6-7-0 0.20 8.9, 11.5, 12.8 ‡ −12.8 −0. obtain a correspondingly small value. The three values are combined as a weighted sum to give , with slightly more weight placed on the agreement between the rotation curve and circular velocity curve than its time variability. Our particular definition of involves several subjective choices, but we have checked that the qualitative conclusions that we reach are not unduly influenced by these. For example, using 0 in place of (i.e. neglecting the time variability of the rotation curves), or using pairs of snapshots spaced by ∼ 1 Gyr instead of ∼ 300 Myr, lead to only small quantitative differences that do not change our interpretation of our analysis. We acknowledge that is not a perfect measure -for example galaxy AP-L1-V10-5-0 ( = 0.37) has a rotation curve that agrees closely with its circular velocity curve over about 70 per cent of its extent, but has large deviations in the remaining 30 per cent such that it is perhaps unduly penalised by our choice to use the 75 th percentile -but it neverthe-less seems to successfully broadly sort our sample of galaxies by the degree of agreement between their circular velocity and rotation curves.
The panels of Fig. 3 are arranged in order of increasing . It is visually clear that the rotation curves of galaxies with higher do not trace the circular velocity curve as closely as those of galaxies with lower , and are likewise more time-variable. We divide galaxies into 4 classes based on the statistic. Of the 33 galaxies in our sample, 3 are labelled 'class 1' ('excellent' agreement between circular velocity and rotation curves; < 0.125), 9 'class 2' ('good' agreement; 0.125 ≤ < 0.175), 7 'class 3' ('fair' agreement; 0.175 ≤ < 0.225) and 14 'class 4' ('poor' agreement; ≥ 0.225). The residuals between the rotation curves and their circular velocity curves at = 0 as a function of radius are shown in Fig. 4, grouped and coloured by class. Consistent with the definition of the parameter, the amplitude  of the residuals increases with . The rotation curves of class 1 and 2 galaxies underestimate their circular velocity curves within the central ∼ 1 kpc, while class 3 and 4 galaxies have rotation curves that systematically underestimate their circular velocity curves across the full radial extent of their discs, but especially near their centres. The magnitude and radial variation of the residuals, especially in the central regions, is in most cases significantly larger than the statistical uncertainties usually assigned to the rotation curves of observed galaxies with similar rotation speeds (e.g. the more massive galaxies in the LITTLE THINGS sample, see Oh et al. 2015or Iorio et al. 2017. These are usually of about 5 − 10 km s −1 , or 10 − 15 per cent, corresponding to about the edge of the green band in Fig. 4. It is also larger than the expected amplitude of pressure-support corrections (Sec. 2.2), with the exception of the lowest (class 1) galaxies -accounting for pressure support would therefore not qualitatively change our results. The values of galaxies are tabulated in Table 1, with rows ordered by increasing . The same table also provides a concise summary of various effects that can (and often do) perturb the rotation curves of the galaxies. We discuss each in turn in Sec. 3.2, but first give a brief qualitative overview.
The overall impression that emerges immediately on visual inspection of the videos of the galaxies in our sample is one of rich variety, both in galaxy properties and in the perturbations that they are undergoing. Whilst selected with a simple criterion: 60 < max /km s −1 < 120, there are large galaxies with gas discs extending nearly 30 kpc in radius, but also tiny galaxies which barely −0.5 resemble discs (radii as small as 2 kpc). Some galaxies have obvious, strong gas outflows, while others are rapidly accreting new gas. There are several instances of galaxies losing the majority of their gas and then accreting a new disc that is highly inclined relative to the previous disc, resulting in a strongly warped disc. Some galaxies are very elongated and/or have frequent lopsided (harmonic of order = 1) perturbations, while others have a long-term stable, circular disc. Mergers and interactions with companions are common, with a range of impact parameters. In many cases the gas merges quickly, while the secondary DM halo completes several orbits before fully merging, visibly disturbing the gas kinematics at each pericentric passage. A few complicated triple mergers are also present in the sample. Other common disruptions include non-merging interactions with gas-rich or gas-less haloes (for our selection in max , we do not find any star-less or 'dark' galaxies that significantly perturb the gas kinematics). In some cases a wind from motion through the intergalactic medium (IGM) seems to cause strong = 1 deformations of the disc.

Mergers and close companions
Close interactions and mergers with gas-rich companions cause the most obvious disturbances to rotation curves. Gas-less (but not necessarily dark) companions cause less disruption, but can still visibly disturb the gas kinematics. In most cases the effect of a gas-less companion on the rotation curve is minimal, and even in the most extreme cases the rotation curve is usually still a reasonably good tracer of the circular velocity curve, if no other perturbation is ongoing. Similar statements apply to the gas-less remnant of an initially gas-rich companion as it returns on subsequent orbital passages. We list two examples of interactions with companions (and all other types of perturbations discussed in subsections below) and where they can be most clearly seen in our collection of visualisations in Table 2.
Using the merger trees (Sec. 2.1), we identify all companion galax-ies that merged into each galaxy in our sample and their progenitors. In addition, we track the progenitors and descendants of all companion galaxies, defined as those found in the same FoF group as the galaxy at any time (but that did not later merge). For each companion and merged object, we find its maximum DM mass at any time and compare it to the maximum DM mass of the galaxy of interest, discarding any with a mass ratio less than 1:20. We found that interaction with smaller mass ratios caused little visible disturbance to the gas discs, and no noticeable perturbation to their rotation curves. We define the time of the first pericentric passage of an interaction as the time of the earliest simulation snapshot when both galaxies are found in the same FoF group and the sign of the radial velocity difference between the companion and host is positive.
In Table 1, we list the times of first pericentric passages for all such interactions, excluding those before 8 Gyr. There are 25 interactions in total, occurring in 15 galaxies. All of the tabulated companions/mergers are initially gas rich -their peak (over time) gas-to-stellar mass ratios are ≥ 1.7. We also note that all companions and mergers with mass ratios greater than 1:20 had starsperturbations due to 'dark' galaxies are unimportant for the galaxies in our sample. The collision of two gas discs in a 1:20 or greater merger invariably strongly and globally disturbs the gas morphology and kinematics, making any other possible perturbations moot. We therefore flag ongoing interactions (which may persist long after the first pericentric passage), defined as those where the closest approach of the companion occurs during the last 3 simulation snapshots (∼ 650 Myr) and is closer than 25 kpc. These are marked ‡ in Table 1 and shown with open symbols in later figures.
It is clear from Table 1 that galaxies with an ongoing interaction with a gas-rich companion have rotation curves that are poor tracers of the circular velocity (classes 3 & 4, with one exception in class 2 where the mass ratio was close to 1:20 to begin with and the companion has been heavily stripped by the end of the simulation). However, galaxies may recover quickly (in a little more than a dynamical time) from earlier interactions, depending on the mass ratio, impact parameter, and the relative inclinations of the gas discs. For example, galaxy AP-L1-V11-3-0 (class 1) finished merging with a massive ( tot ∼ 3.3 × 10 10 M , dark matter mass ratio ∼ 0.13) companion around 12.5 Gyr (see entry in Table 2 for details), but since the approach was nearly in the plane of the gas disc, the disruption of the disc was limited and the gas settled after the merger ended, such that the rotation curves trace the circular velocity very well by 12.7 Gyr and thereafter.
On the other hand, galaxy AP-L1-V6-15-0 (class 4) experienced an oblique collision with a galaxy ( tot ∼ 4.3 × 10 9 M , dark matter mass ratio ∼ 0.15), dispersing nearly all of the primary galaxy's gas on its first approach (∼ 11.5 Gyr), before the now partially stripped, but still gas-rich companion halo returns and the rest of the gas merges (∼ 12.2 Gyr). This dramatic event severely disrupts both the gas disc and DM halo, and the gas dynamics are entirely out of equilibrium until ∼ 13.0 Gyr where they begin to settle. By = 0 the rotation curve still underestimates the circular velocity overall and has radially localised features (e.g. 'wiggles').

Bulk non-circular gas flows
Bulk non-circular gas flows (e.g. radial or vertical flows) directly violate the assumption of rotational support implied by the expectation that the rotation curve of a galaxy should agree with its circular velocity curve. Bulk outflows in low-mass galaxies in APOSTLE are driven predominantly by the injection of thermal energy by supernovae and are preferentially ejected along the 'path of least resis- AP-L1-V11-3-0-gas-edge-and-face.mp4 Merger with the gas disc of a companion. The companion arrives on a prograde orbit nearly in the plane of the disc around time 10.2 (this is actually the second passage, the first was around 8.4 Gyr, at this time the gas discs interacted but did not collide); the disc survives and settles (by time ∼ 12.7 Gyr). Merger/companion ( §3.2.1) 11.2 -12.0 Gyr ( = 0.27) AP-L1-V6-15-0-gas-edge-and-face.mp4 Merger with the gas disc of a companion. The companion arrives on an oblique prograde orbit; the disc is almost completely destroyed and does not re-form until time ∼ 13.3 Gyr. Radial inflows ( §3.2.2) 13.0 -13.5 Gyr ( = 0.18) AP-L1-V4-14-0-gas-edge-and-face.mp4 Gas is visibly ejected from the disc, likely by a series of supernova explosions, around time 12.5 Gyr. This gas quickly begins to settle back onto the disc. While this is ongoing, the entire disc contracts radially. Vertical outflows ( §3.2.2) 13.0 -13.5 Gyr ( = 0.26) AP-L1-V4-13-0-gas-edge-and-face.mp4 Several prominent whisps of ejected gas are visible both above and below the disc (in the right panel of the video), launched over a period of a few hundred megayears.
The DM halo is visibly elongated throughout, driving transient lopsided ( = 1 harmonic) and bisymmetric ( = 2) deformations of the gas disc. For example, at time 13.0 Gyr, the disc is both elongated and lopsided. The position angle of the halo as seen in the visualisation rotates over time. Initially (8.9 Gyr) the major axis in projection runs from the lower-left to upperright of the image, by 10.3 Gyr it is horizontal, etc. When the = 2 mode of the gas disc is excited, its major axis is preferentially about 45 • clockwise from the major axis of the halo (e.g. horizontal in the image at 8.9 Gyr). Warped disc ( §3.2.4) 13.0 -13.5 Gyr ( = 0.18) AP-L1-V4-6-0-gas-edge-and-face.mp4 By 11 Gyr the gas disc is very small after being consumed by star formation and losing gas to supernova feedback. Between 12 and 13 Gyr a large amount of gas accretes onto the disc, misaligned with the existing disc. By 13.5 Gyr the edge-on planes of the inner and outer discs are clearly visible in the right panel of the visualisation. IGM wind ( §3.2.5) all times ( = 0.09) AP-L1-V1-4-0-gas-edge-and-face.mp4 Throughout the visualisation the diffuse gas surrounding the disc has a noticeable net flow from right to left in the image, in both the face-on and edge-on views.
tance': vertically from the disc (see Table 2 for an example). Bulk inflows within the disc, on the other hand, tend to be radial and are usually associated with gas accretion (see example in Table 2). We quantify bulk non-circular gas flows as follows. We focus on the atomic gas disc by first selecting only 'atomic' gas particles, which we recall that we define as those with H mass fractions of > 0.5, and then selecting only those particles within a cylindrical aperture with a radius equal to the radius enclosing 90 per cent of the H mass of the galaxy, and a half-height equal to the half-height enclosing 90 per cent of the H mass. We calculate the radial (in cylindrical coordinates) and vertical bulk flow rates of the selected particles as their massweighted average radial and vertical velocities. For the vertical flow rate, we use the speed towards or away from the disc midplane (i.e. sgn( ) ). Calculating these flow rates for a few consecutive simulation snapshots revealed that they are highly time-variable, motivating us to choose a summary statistic. The peak (most negative) radial inflow rates and peak (most positive) vertical outflow rates from the last 3 snapshots (13.10, 13.43 and 13.76 Gyr) are plotted in Fig. 5 against the parameter defined in equation (1). To emphasize that these flow rates capture a global contraction/expansion of the disc 5 rather than e.g. the speed of gas selected to be 'outflowing' or 'inflowing', we show values in units of kpc Gyr −1 (rather than e.g. km s −1 ). We also note that 'harmonic' non-circular motions, such as a bar-like distortion of the gas orbits, are not captured in this measurement because such distortions do not result in a net transport of gas. Fig. 5 shows that galaxies with stronger inflows and/or outflows 5 An equivalent interpretation is in terms of the total momentum of the H gas: the radial inflow rate is determined from the total (cylindrical) radial momentum of selected particles divided by the total mass of selected particles, which could be termed a 'specific linear momentum' and has the dimensions of a speed. Similarly, the vertical outflow rate is the total specific momentum towards/away from the disc midplane. The average radial velocity (in cylindrical coordinates) of 'atomic' gas particles in each galaxy in our sample is calculated at each of the last 3 snapshots (13.10, 13.43 and 13.76 Gyr), and the minimum value (i.e. peak inflow rate) is plotted on the vertical axis. Lower panel: The average vertical velocity away from the disc midplane (i.e. sgn( ) ) of 'atomic' gas particles in each galaxy in our sample is calculated at each of the last 3 snapshots, and the maximum value (i.e. peak outflow rate) is plotted on the vertical axis. Galaxies currently strongly interacting with a companion (marked ‡ in Table 1) are plotted with open symbols. Galaxies with stronger bulk flows have preferentially higher values. Galaxies with peak radial inflow rates stronger (more negative) than −5 kpc Gyr −1 and/or peak vertical outflow rates greater than 1 kpc Gyr −1 (red dashed lines; galaxies with flows stronger than either or both limits are shown with red markers) are not found in class 1 (green background), and only one galaxy without strong bulk flows is found in class 4 (red background), out of 14 galaxies in this class. In all panels, points are numbered in order of increasing , the same order in which galaxies appear in Table 1. tend to have rotation curves that are poorer tracers of their circular velocity curves. We illustrate this by plotting galaxies with radial inflow rates stronger (more negative) than −5 kpc Gyr −1 and/or vertical outflow rates greater than 1 kpc Gyr −1 (approximately the median flow rates for galaxies in our sample) with red markers. Entries in Table 1 exceeding these thresholds are also highlighted in bold face. By this measure, most galaxies in our sample (26/33) have strong bulk flows in at least one of these two directions, but no galaxies in our class 1 ( < 0.125) do. There is furthermore a clear correlation between the vertical outflow rate and , albeit with large scatter (see Table 3). While the strongest (most negative) radial inflow rates seem to occur in galaxies with higher , these two quantities are not significantly correlated (see Table 3).
The connection between bulk flows, in the vertical direction in particular, was one of the first that we noticed in our initial visual analysis of our collection of videos: a lack of visible outflows from a galaxy is a strong predictor that its rotation curve will be a good tracer of its circular velocity curve. However, given the diversity of perturbations which can cause rotation curves to differ from the Table 3. Statistical correlation coefficients for the distributions plotted in Figs. 5-8. Interacting or merging galaxies (open symbols in the figures) are excluded from the calculations. The -values reported are 'one-sided', i.e. they reflect the probability that the hypothesis that the data do not have a correlation of the given sign (negative for rows 1, 3 and 5; positive for rows 2 and 4) can be rejected. We caution that Pearson's statistic is indicative only, since we do not expect the correlations to be necessarily linear, nor the variables involved to be Gaussian-distributed. The assumptions underlying Spearman's statistic are satisfied. The correlations of peak bulk vertical outflow rate, DM halo / and warp angle with are likely significant, while the other two quantities are consistent with being uncorrelated with .

Pearson
-value Spearman -value Peak bulk radial inflow rate  Table 1) is a strong predictor of poor agreement between the rotation curve and the circular velocity curve, but a spherical halo does not guarantee close agreement. The coloured background marks the same intervals in as introduced in Fig. 3. Points are numbered in order of increasing , the same order in which galaxies appear in Table 1. circular velocity curve, having weak bulk flows does not guarantee this to be the case, as is evident from exceptions such as AP-L1-V10-16-0. This galaxy has amongst the weakest bulk flows in our sample, but falls in class 4 ( ≥ 2.225).

Dark matter halo shape
Elongated or triaxial DM haloes give rise to non-circular gas orbits, with gas often visibly sloshing around in the aspherical potential (see Table 2 for an example). In galaxies where this mechanism is effective, the rotation curves are highly variable as strong, transient lopsided (harmonic of order = 1) and bisymmetric ( = 2) modes are excited in the gas disc. Fig. 6 shows the anti-correlation between the intermediate-to-major axis ratio of the DM halo and the parameter defined in equation (1). We focus on the shape of the halo in the region occupied by the disc by calculating axis ratios using DM particles in a spherical aperture with a radius equal to twice the radius enclosing 90 per cent of the galaxy's H mass. The squares of the axis lengths are proportional to the eigenvalues of the reduced inertia tensor: where and are the coordinate vector and mass of particle , respectively.
Even very small departures from / = 1 seem to be sufficient to drive large changes in the rotation curves. The red dashed line in Fig. 6 marks / = 0.95 -no galaxies with / < 0.95 fall in our class 1 ( < 0.125), and all save one class 4 ( ≥ 0.225) galaxies have / < 0.95. We highlight the entries for galaxies with / < 0.95 in Table 1 -these make up 17 of the 33 galaxies in our sample. The anticorrelation in the figure has considerable scatter, reflecting the fact that a galaxy with a spherical halo can be perturbed by some other mechanism, but an aspherical halo seems to be a strong predictor of the rotation curve being a poor tracer of the circular velocity curve. Nevertheless, the anti-correlation is very clear (see Table 3).
Although not shown in Fig. 6, we also investigated trends in as function of the minor-to-major axis ratio ( / ) and the triaxiality parameter ( ≡ 2 − 2 2 − 2 ). These show somewhat weaker trends than that with / , suggesting that a prolate or triaxial halo shape (i.e. / ≠ 1) has a stronger perturbative effect than an oblate shape ( / ∼ 1). This agrees with intuition: a light, rotationally supported disc has possible stable configurations in the potential an oblate halo, but is unstable in a prolate or triaxial potential.
The mass of the gas disc also plays a role. A more massive disc may resist the perturbative effect of an aspherical halo, or even 'sphericalise' the halo if it is massive enough. We will return to the importance of the gas disc mass in Sec. 3.3 below.

Warped discs
Several galaxies in our sample have visible warps in their gas discs; one example of a prominent warp is listed in Table 2. We quantify the strength of a warp by the angle warp between the angular momentum vectors of the inner and outer gas discs, which we define as the inner 30 per cent and outer 60-90 per cent of the H gas by mass. We plot warp against our parameter defined in equation (1) in Fig. 7. Most galaxies in our sample have warp angles of less than ∼ 20 • and these span the entire range in , but a minority have large warp angles caused by rapid accretion of gas with angular momentum strongly misaligned with the existing disc, or an interaction with a companion -these galaxies have preferentially higher values (see Table 3) and fall in our classes 3 & 4 ( 0.175). We flag galaxies with warp > 30 • (dashed line in the figure) as strongly warped, and highlight their entries in Table 1 in bold face. There are 8 such galaxies in our sample. We note that the influence of a warp on the rotation curve is somewhat exaggerated in our analysis because we have measured the rotation curves as the median azimuthal velocity of particles in a fixed plane aligned with the inner disc. The rotation speed in a warped outer disc is therefore underestimated by about a factor of cos warp , which will be reflected in the measured value.  Table 1) is associated with a poor agreement between the rotation and circular velocity curves ( 0.2). The coloured background marks the same intervals in as introduced in Fig. 3. Points are numbered in order of increasing , the same order in which galaxies appear in Table 1.

IGM wind
In the visualisations of some of our galaxies, a 'wind' blowing against the gas disc due to its motion through the IGM is clearly visible, and appears to deform the disc, often resulting in a lopsided disc displaced 'downwind'. One example where this effect is especially clear is given in Table 2. The effect is visually similar to recent observations of the isolated Wolf-Lundmark-Melotte (WLM) dwarf galaxy (Yang et al. 2022), which also show gas clouds trailing 'downwind' from the galaxy's proper motion, providing evidence for ram-pressure stripping by the intergalactic medium. We note, however, that the simulated galaxy in question is in class 1 ( = 0.09) and has an IGM wind speed at = 0 (see below for details) close to the median in our sample of galaxies. This highlights both the difficulty in quantifying the strength of the wind and its potentially quantitatively subtle effect on the kinematics of the gas disc, despite the perturbative effect of the wind being visually very clear. As a consequence, our efforts to quantify such perturbations have yielded less clear-cut results than for the other types of perturbations discussed above, suggesting that perturbation due to a wind may be more nuanced. The effect on the rotation curve likewise often seems to be fairly subtle.
We estimate the speed of the IGM wind as follows. We first select gas particles in a spherical shell between one and two times the virial 6 radius around the galaxy. In order to avoid undue bias by other nearby galaxies, we further restrict our selection to include only those particles not gravitationally bound to any subhalo according to the halo finder. We take the median velocity in the rest frame of the galaxy (the same frame used when measuring the rotation curves) of the remaining selected particles to be the IGM wind velocity. We have verified that the conclusions that we reach are not very sensitive to 6 We define the virial radius as the radius of a sphere within which the mean matter density is 200 times the critical density crit = 3 2 0 /8 . . IGM wind speed wind , calculated as the median velocity of the gas particles between 1 and 2 times 200 of the galaxy and not belonging to any FoF group in a frame of reference where the gas disc is at rest (see Sec. 3.2.5 for details), plotted against the degree to which the rotation curve traces the circular velocity curve, (Sec. 3.1). No strong trend is visible. The dashed red line marks 50 km s −1 ; we consider galaxies above this line to be the strongest outliers in wind , and mark the corresponding values in bold in Table 1. The coloured background marks the same intervals in as introduced in Fig. 3. Points are numbered in order of increasing , the same order in which galaxies appear in Table 1. the precise radial range used (within a factor of about 3), or whether or not bound particles are included.

Rotation curves of low-mass galaxies
The speed of the wind wind is plotted against the parameter defined in equation (1) in Fig. 8. Any correlation is less clear than those seen in Figs. 5-7 above, but there is tentative evidence that the galaxies with the highest IGM wind speeds in our sample have higher values, or at least avoid the lowest values, but there is no statistically significant correlation between wind and (see Table 3. We draw a boundary at 50 km s −1 (dashed line in the figure) separating the galaxies with the highest wind speeds (5 of the 33 in our sample) from the others, and highlight the entries corresponding to the galaxies above this threshold in Table 1.
In addition to the caveats listed above, we are cautious in our interpretation of perturbations due to the IGM wind because we struggle to find a clear correspondence between the galaxies where we identified what appeared to be a wind in our visualisations and those with a high wind speed (or other similar quantitative measures that we explored). Furthermore, our impression from our visual inspection is that periods of strong IGM wind are often short-lived (the example in Table 2 is an exception to this), and any perturbation of the rotation curve does not seem to persist after the wind subsides.

Summary
Taking all of these various kinds of perturbations into account, it is perhaps unsurprising that so few galaxies in our sample have a rotation curve that closely traces their circular velocity curve at = 0. However, looking at Table 1, there are also a few galaxies that have avoided any obvious recent disturbance and yet have rotation curves that are poor tracers of their circular velocity curves. The galaxy AP-L1-V10-30 presents an intriguing case. The visualisations (e.g. AP-L1-V10-30-gas-edge-and-face.mp4) do not reveal any obvious perturbations in any of the categories discussed above at late times, except perhaps some vertical outflows from the disc. Inspecting its entry in Table 1, none of its properties exceed our (admittedly somewhat arbitrary) thresholds for strong perturbations. And yet, its late-time rotation curve (Fig. 3) significantly underestimates the circular velocity curve in the central ∼ 2 kpc, and is time-variable in the outer regions of the disc. The presence of such an example in our sample of galaxies emphasizes that we have only scratched the surface of a complex topic: it is clear that many types of perturbations significantly influence the gas kinematics in low-mass galaxies, but a more complete understanding of the prevalence and importance of each type will require further study, often on a galaxy-by-galaxy basis.
Ultimately, given the intrinsically limited information available from observations of real galaxies, a practical question to ask is: are galaxies where the rotation curve is a good tracer of the circular velocity curve separated from others in terms of observable properties? In our exploration of our sample of galaxies, it became clear very quickly that gas mass (e.g. at fixed stellar mass) plays an important role. Once galaxies with ongoing interactions or mergers are removed from consideration, galaxies with higher gas mass are more likely to have rotation curves that trace the circular velocity well, and vice versa (see further discussion in Sec. 3.3 below). This is tentatively consistent with the stabilising effects of a massive disc against some perturbations, including the influence of an aspherical DM halo (the type of perturbation with the clearest effect out of those that we investigated, save mergers), as mentioned above.

Trends with galaxy scaling relations
We plot each pair-wise relation between gas mass, stellar mass, and max at = 0 in Fig. 9 (tabulated values are available in Oman et al. 2019, table A1). We also plot the data from the SPARC compilation (Lelli et al. 2016a) with small grey points 7 . The simulated galaxies broadly follow observed trends in these relations; we do not discuss the comparison further (see Oman et al. 2019, for a detailed comparison). Points for simulated galaxies are coloured by their their class, from green (class 1) to red (class 4). We plot galaxies with a recent merger or interaction with a companion (see Sec. 3.2.1) with an open symbol. Considering the lower panels, we highlight that all non-interacting galaxies with gas mass greater than 1.5 × 10 9 M (9 galaxies) are in class 1 or 2, while only 2/17 non-interacting galaxies with lower gas masses are. We note that galaxies with recent/ongoing interactions or mergers have preferentially higher gas masses (at fixed max ) than non-interacting galaxies -this is unsurprising, since all companion galaxies of galaxies in our sample bring a lot of gas with them (see Sec. 3.2.1). There is no similarly clear separation between points of different colours in maximum circular velocity or stellar mass, besides some weak trends coming from the fact that gas mass correlates with both of these parameters.
It is tempting to attribute gas mass being more important than stellar mass in this context to the galaxies in our sample having gas masses exceeding their stellar masses, such that the gravitational influence of the stars on the gas kinematics is not dominant. However, the gas is typically more extended, so the stars can still dominate the gravitational potential near the centre, and can exert a strong nongravitational influence through supernova feedback. The strength showing galaxies with a recent gas-rich merger. SPARC data from Lelli et al. (2016a) are plotted with small grey points. The areas outside of our selection 60 < max /km s −1 < 120 are shaded in the right panels.
of supernova feedback would be expected to correlate instead with recent star formation, but we did not find any strong trend with recent star formation rate (not shown).
Considering the lower-right panel of Fig. 9, we were surprised not to find a stronger dependence on a combination of max and gas , as might be expected if the gas-to-total mass ratio was a primary driver of . We note, however, that a larger sample of simulated galaxies would be very helpful in exploring these issues further.
We next turn our attention to the possible biases introduced into observational scaling relations involving rotation curve measurements by the types of perturbations discussed above. We emphasize that we investigate here only the 'direct' impact due to the difference between the rotation curve and the circular velocity curve -the rotation curves that we consider are those that an observer with perfect knowledge of the gas kinematics would measure: the median azimuthal velocity as a function of radius. With real observations, this direct impact is likely to be compounded by additional errors induced by e.g. attempting to model a non-equilibrium system assuming equilibrium dynamics, assuming circular orbits when the actual orbits are noncircular, etc. (see e.g. Read

The BTFR
The upper panel of Fig. 10 shows the BTFR of the galaxies in our sample, along with galaxies from the SPARC compilation (Lelli et al. 2016a) and the BTFR of SPARC galaxies reported by Lelli et al. (2016b) to provide context. We do not undertake a comparison with the observed BTFR in this work, as this has previously been addressed by Oman et al. (2016) and Sales et al. (2017). In this panel, the horizontal axis shows max , the maximum of the circular velocity curve. This can be thought of as the 'truth' that is obtained in the ideal case where the gas rotation curve follows the circular velocity curve and the measurement of the rotation curve is without error. We also plot an indicative linear fit to the points in this panel as a black solid line, excluding interacting/merging galaxies (open symbols) from the calculation. The fit minimizes the sum of the squared offsets in bar from the BTFR. The best-fitting slope ( bar ∝ max ) is = 3.9. The ratio flat / max , plotted against max . The maximum gas rotation velocity systematically underestimates the maximum circular velocity ( flat / max < 1), and the underestimates get systematically worse at lower max ; this has the potential to bias both the normalisation and the slope of the BTFR. In all panels, points are numbered in order of increasing , the same order in which galaxies appear in Table 1.
In the centre panel of Fig. 10, the coloured points show a measurement of the maximum rotation velocity of the gas, determined from the flat portion of the rotation curve following the approach of Roper et al. (2022, Appendix C; if the rotation curve is still rising at the outermost point, the value at this point is used) -we label this flat . Briefly, flat is measured by iteratively searching for the velocity in the rotation curve that maximises the number of other velocities within an adaptive fractional difference from it. Each point is joined to its position in the upper panel by a solid line. Unsurprisingly, galaxies in our classes 3 & 4 move further (on average) from their positions in the upper panel than those in our classes 1 & 2. Furthermore, nearly all points shift to the left, as the rotation curves preferentially underestimate the circular velocity curves. This is emphasized in the lower panel of the figure, where the ratio flat / max is plotted against max -here it is clear that the underestimates get systematically worse towards lower max (or lower bar ). The max and flat values for each galaxy in our sample are tabulated in Table B1. The trends evident in the bottom panel of Fig. 10 mean that the BTFR is biased to a higher normalisation (because flat systematically underestimates max ), and shallower slope (because the underestimates get worse at lower max ). The change in slope is illustrated in the centre panel by the dashed line, which has a slope of = 3.1. This shows a linear fit to the filled points in this panel, similar 8 to the solid line (repeated in the upper and centre panels).
In the APOSTLE simulations, the BTFR has a steep cutoff around max = 50 km s −1 (see Oman et al. 2016;Sales et al. 2017). Replacing max with flat seems to soften the cutoff, potentially enough for the trend to become more reminiscent of the constant slope often claimed in observational studies (e.g. McGaugh et al. 2000;Ponomareva et al. 2018;Lelli et al. 2019), although an analysis including galaxies at lower max would be needed to confirm this.
If observed galaxies are subject to broadly similar perturbations as those that we observe in our simulations, which seems likely, then the observed BTFR is probably biased in a similar sense. The magnitude of the effect, however, depends on the details of each type of perturbation and their relative frequencies, which the simulations may not capture in full detail.
Interestingly, attempting to remove galaxies where the rotation curve does not trace the circular velocity curve from a sample to be used to measure the BTFR likely still results in a bias, because those galaxies where the rotation curve is a good tracer of the circular velocity curve are not an unbiased sub-sample: they tend to be the most gas-rich galaxies and to have higher bar at fixed max (see discussion of Fig. 9). There is observational evidence for such a bias: Papastergis et al. (2016) found that using a sample selected to be extremely gas-rich ( gas / ★ 2.7) yields a steeper slope for the BTFR bar ∝ ( /2) as a function of H line width than studies with less extreme selections (e.g. Zaritsky et al. 2014;Hall et al. 2012;McGaugh 2012). They find = 3.75±0.11 (rather than ∼ 3.3-3.4). Ball et al. (2022) similarly find that restricting their sample to gasrich ( HI / ★ > 2) galaxies significantly increases the slope of the BTFR, from about 3.3 to 3.9. They also find that dividing their galaxy into high-and low-baryonic mass sub-samples (at bar = 10 10 M gives different slopes, of 2.9 and 4.1, respectively. However, Gogate et al. (2022) instead find no strong dependence on gas fraction. All of these studies use spatially-integrated spectral line widths for the velocity axis of the BTFR. Searching for similar trends when spatially resolved rotation curves are used instead is an interesting avenue for future studies.

The gas specific angular momentum-mass ( gas − gas ) relation
We loosely follow Mancera Piña et al. (2021, their eq. (1)) to define the specific angular momentum of the gas, gas , in our galaxies. We evaluate it as: where we have discretized the disc into a series of radial bins (of 500 pc width, extending out to the radius enclosing 90 per cent of the H mass) labelled with gas masses gas, , rotation speeds and radii . We consider below both the case where is measured from the circular velocity curve (we label this gas,circ ), and that where it is measured from the gas rotation curve ( gas,rot ). Following Mancera Piña et al. (2021), we define gas = 1.33 HI . Each time that we make a measurement of gas , we also determine whether it is converged according to the criterion in Mancera Piña et al. (2021, sec. 3.2). The gas,circ and gas,rot values for each galaxy in our sample are tabulated in Table B1, along with their convergence diagnostics.
In the upper panel of Fig. 11 we show the gas,circ − gas relation for galaxies in our sample. For comparison, we also show the measurements of Mancera Piña et al. (2021), including only those observed galaxies where the gas measurement is converged according to their criterion. In general the simulated galaxies follow the observed relation remarkably well (see also Roper et al. 2022, fig. 4), especially when restricted to those where the gas measurement converges (larger coloured points) and that are not currently merging/interacting (filled points) -although this reduces our sample of 33 galaxies to only 7.
In the second panel of Fig. 11, we show the gas,rot − gas relation. Similarly to Fig. 10, we join the points to their positions in the upper panel with line segments, but in most cases these are too small to be visible, emphasizing that the gas − gas relation is less sensitive to rotation curves that fail to trace the circular velocity curve of their galaxy than the BTFR. We again mark the galaxies where the gas measurement converges with larger symbols -in this case, there are 9/33 (that are also not merging/interacting with companions). In the lower panel of the figure we show the ratios gas,rot / gas,circ , which shows that the specific angular momenta scatter by only about 15 per cent between the two sets of measurements, with a bias toward gas,rot < gas,circ that is more pronounced for galaxies with lower gas masses (which are also those preferentially found in our classes 3 & 4, as shown in Fig. 9). We caution, however, that there are only 4 galaxies in our sample (which are not merging/interacting with a companion) for which both the gas,rot and gas,circ are converged.
Even though the gas specific angular momentum-mass relation has less potential than the BTFR to be biased by the failure of rotation curves to trace the circular velocity curves of their galaxies, a possible 15 per cent contribution to the scatter is hardly negligible -it represents about a third of the estimated intrinsic scatter of the relation (0.15 dex or about 40 per cent, according to Mancera Piña et al. 2021). Conversely, the tightness of the gas − gas relation does not seem to have any strong implications regarding how closely rotation curves correspond to circular velocity curves.  Figure 11. The gas specific angular momentum-mass ( gas − gas ) relation. Upper panel: Gas specific angular momentum gas,circ (see Eq. (6)) calculated assuming the gas is rotating at the circular velocity against the gas mass ( gas = 1.4 HI ) for our sample of simulated galaxies. The points are coloured and open/filled as in Fig. 9. Galaxies whose gas has not converged by the edge of the gas disc according to the criterion of Mancera Piña et al. (2021, sec. 2.3) are shown with smaller markers. The measurements of Mancera Piña et al. (2021) are shown with the small grey markers and error bars, including only those galaxies where the measurement has converged, and their fit to these points is shown with the thin grey line. Centre panel: As upper panel, but with gas calculated using the actual rotation curve of the gas for the simulated galaxies. The observational measurements and fit are repeated from the upper panel. Lower panel: The ratio gas / gas,circ , plotted against the gas mass gas . Only galaxies where both gas and gas,circ have converged are shown with larger points. In all panels, points are numbered in order of increasing , the same order in which galaxies appear in Table 1.

The H mass-size relation
The H mass-size relation ( HI − HI ) relates the mass of an H disc and its size, usually defined as the radius at which the H surface density drops below 1 M pc 2 (e.g. Wang et al. 2016). We show this relation for our sample of simulated galaxies in Fig. 12 (see also Oman et al. 2019, fig. 1). Essentially all galaxies in our sample lie  Oman et al. (2019, table A1). We note that this definition differs from the radius disc enclosing 90 per cent of the H mass used elsewhere in this work. The points are coloured and open/filled as in Fig. 9. The relation for galaxies in the SPARC compilation (Lelli et al. 2016a) is also shown for comparison. Points are numbered in order of increasing , the same order in which galaxies appear in Table 1. on a tight relation offset to slightly larger sizes (by ∼ 0.2 dex) than the observed relation. The only strong outlier (AP-L1-V10-20-0, well below the relation) actually has an H surface density profile that is flat at just under 1 M pc −2 out to about 7 − 8 kpc, so its departure from the relation is simply a consequence of the somewhat arbitrary choice of threshold surface density. Stevens et al. (2019) argued that the H mass-size relation is an essentially inevitable consequence of the way that H gas settles into a disc and is converted to/from H 2 . Any departures of more than a factor of ∼ 2 (in radius) from the observed relation should only occur in simulations for quite extreme failures in modeling the interstellar medium. Fig. 12 demonstrates that even the most kinematically disturbed galaxies in our sample still lie on this relation, and therefore that a galaxy's position in the HI − HI plane apparently carries no information about the connection between its circular velocity curve and its rotation curve.

The fid − max relation
Santos-Santos et al. (2020) adapted a relation introduced by Oman et al. (2015) that relates the maximum rotation speed (or circular velocity) and the rotation speed (or circular velocity) at an inner radius, fid = ( fid ). The radius fid ≡ ( max /70 km s −1 ) 2 kpc is defined to adapt to the scale of each galaxy. This quantifies the shape of the rotation curve or circular velocity curve: a more slowly rising curve (i.e. a rotation curve with a shallow inner slope) has a lower fid at a given max than a more steeply rising curve. We plot this relation for our sample of simulated galaxies in the upper panel of Fig. 13, here using fid,circ and max measured from their circular velocity curves. We also plot measurements from the compilation of Santos-Santos et al. (2020) for context (see Oman et al. 2019;Santos-Santos et al. 2020;Roper et al. 2022, for further discussion of the comparison).
In the second panel, we plot the locations of our sample of galaxies in the same space, but measured from their = 0 rotation curves ( flat  Fig. 9. The areas outside of our selection 60 < max /km s −1 < 120 are shaded in light grey. Second panel: As upper panel, but the inner and outer rotation velocities fid,rot and flat are measured from the = 0 rotation curve. Points are joined by a line to their positions in the upper panel. Third panel: Ratio of the outer rotation velocity measured from the rotation curve and that measured from the circular velocity curve flat / max , plotted against max . Point colours and open/filled symbols as in upper panels. Lower panel: Ratio of the inner rotation velocity measured from the rotation curve and that measured from the circular velocity curve fid,rot / fid,circ , plotted against max . Point colours and open/filled symbols as in upper panels. In all panels, points are numbered in order of increasing , the same order in which galaxies appear in Table 1. is measured as in Fig. 10). The values corresponding to each plotted point for all galaxies in our sample are tabulated in Table B1. In addition to the tendency for rotation curves to underestimate the maximum circular velocity as discussed above, 29 of the 33 galaxies in our sample have a significantly lower rotation velocity than circular velocity at fid , with the effect being severe for most class 4 galaxies. In some cases the shape of the rotation curve is broadly preserved (displacements parallel to the solid grey line in the second panel), while in others the rotation curve rises much more slowly (vertical displacement downwards). Interestingly, the resulting scatter in the space of fid − flat is not dissimilar from that observed for the SPARC galaxies, with even some similarly extreme outliers. We caution, however, that in practice rotation curve measurements do not recover the median azimuthal velocity as a function of radius exactly but are subject to various systematic errors in modelling, especially in their central regions (Oman et al. 2019). The scatter in the lower panel of Fig. 13 is therefore likely a lower bound on what would be obtained were these simulated galaxies 'observed' and modelled analogously to real galaxies -as is confirmed by Oman et al. (2019) for a subset of the galaxies in our sample. The discrepancy between the rotation curves of low-mass galaxies and their circular velocity curves may be a significant contributor to the diversity in the shapes of observed dwarf galaxy rotation curves highlighted by Oman et al. (2015).

Summary
That the cold gas in some observed galaxies is out of equilibrium and is therefore a poor dynamical mass tracer is well known. However, just how rare it may be that an atomic gas rotation curve can reasonably be interpreted as a circular velocity curve has not been previously been systematically explored. Our visualisations of the gas kinematics of low-mass APOSTLE galaxies (60 < max /km s −1 < 120) over the past ∼ 4 Gyr emphasize the wide variety of processes perturbing them.
Only about a third (12/33) of the galaxies in our sample have rotation curves that we would describe as similar to their circular velocity curves, with examples of close matches being rarer still (3/33). These are found at preferentially higher gas masses ( gas 1.5 × 10 9 M ). Based on our visual inspection of galaxies and their recent history, the most frequent types of perturbations include: • Mergers and interactions with gas-rich companion galaxies (6/33).
The fractions in parentheses indicate the fraction of galaxies in our sample that exceed the thresholds for 'strong perturbations' of the given type outlined in Secs. 3.2.1-3.2.5 (entries in bold face in Table 1). Only 5/33 galaxies in our sample avoid 'strong' perturbations in all of these categories at = 0. Some of these types of perturbations (e.g. mergers) are readily identified observationally, such that the galaxy in question can be excluded from samples for kinematic analysis, but others (e.g. IGM wind, influence of triaxial DM halo) are much more subtle. Furthermore, because susceptibility to perturbation correlates with galaxy properties such as total cold gas mass, omitting perturbed galaxies from analyses introduces biases. In particular, we find that this has probably led to an underestimate of the low-velocity slope of the baryonic Tully-Fisher relation, offering a straightforward explanation for the steeper slope for gas-rich galaxies found by (Papastergis et al. 2016).
Whether our findings based on the APOSTLE simulations are applicable to observed galaxies depends on how faithfully the simulations capture the relevant physical processes. Of the main categories of perturbations that we see operating in the simulations, we would characterize only one (supernova-driven outflows) as sensitively dependent on modelling choices in which there is significant ambiguity. Other processes like the merger rate or the shapes of DM haloes are natural consequences of structure formation in a ΛCDM cosmology, and depend on physics that is well understood and straightforward to implement in the models. There is only a single galaxy in our sample that we have flagged as having strong vertical outflows, but not any other strong perturbations at = 0. Our main conclusion that a majority of galaxies in the max range of our sample have rotation curves that differ significantly from their circular velocity curves is therefore probably also applicable to real low-mass galaxies, but confirming this in other cosmological hydrodynamical galaxy formation models would reinforce this. Galaxies with non-equilibrium gas kinematics are therefore likely one of the main drivers of the observed kinematic diversity (as highlighted by Oman et al. 2015) in dwarfs.

Reflections on visualisation-driven analysis
Our starting point for all of the analysis presented above was our collection of galaxy evolution visualisations and their circular velocity and rotation curves at the corresponding times. This allowed us to build a strong intuition for the perturbations affecting the galaxies in our sample. The visualisations highlight the diversity and complexity of these low-mass galaxies in a way that cannot be fully captured by integrated properties (such as those in Figs. 5-8) and provided important context for our more quantitative analysis. We can identify several instances where we would probably have reached qualitatively different conclusions if the visualisations were not available to guide our intuition and analysis.
The sheer wealth of information represented by the visualisations and rotation curves eventually motivated our choice to focus our analysis on the current time ( = 0). The gas discs of essentially every galaxy in our sample have been subject to different perturbative processes at different times. The time dimension of our data set remains largely unexplored, offering an interesting avenue for future work.
The human eye is an exceptionally powerful tool for reducing complex visual information to simple patterns and trends. Visualisationdriven analysis of cosmological hydrodynamical simulations has, in our opinion, a largely untapped potential to advance our understanding of a wide variety of physical processes in galaxies.
for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K00042X/1, ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure. This work has made use of NASA's Astrophysics Data System.

DATA AVAILABILITY
The SPARC data are available at https://cdsarc.cds.unistra. fr/viz-bin/cat/J/AJ/152/157, with supplementary data tabulated in Santos-Santos et al. (2020), table A1. Access to the APOS-TLE simulation data is available on reasonable request to the corresponding author. Basic properties of galaxies in our sample are tabulated in Oman et al. (2019), table A1. Table B1.
Integrated properties of galaxies from our sample of simulated galaxies. The first three columns are repeated from Maximum circular velocity, max .

(9):
Flat value of the rotation curve, flat , see Sec. 3.3.1. Galaxies marked † have a rotation curve that is still rising at the outermost measured radius.