Non-circular motions and the diversity of dwarf galaxy rotation curves

TL;DR

This study uses simulated dwarf galaxy data to show that non-circular gas motions can cause significant errors in rotation curve estimates, affecting interpretations of dark matter distribution.

Contribution

It demonstrates how non-circular motions in gas lead to diverse and biased rotation curves in dwarf galaxies, highlighting limitations of the tilted-ring modeling approach.

Findings

Non-circular motions cause large diversity in rotation curves.

Tilted-ring models often underestimate true circular velocities.

Non-circular motions can mimic the presence of dark matter cores.

Abstract

We use mock interferometric HI measurements and a conventional tilted-ring modelling procedure to estimate circular velocity curves of dwarf galaxy discs from the APOSTLE suite of {\Lambda}CDM cosmological hydrodynamical simulations. The modelling yields a large diversity of rotation curves for an individual galaxy at fixed inclination, depending on the line-of-sight orientation. The diversity is driven by non-circular motions in the gas; in particular, by strong bisymmetric fluctuations in the azimuthal velocities that the tilted-ring model is ill-suited to account for and that are difficult to detect in model residuals. Large misestimates of the circular velocity arise when the kinematic major axis coincides with the extrema of the fluctuation pattern, in some cases mimicking the presence of kiloparsec-scale density 'cores', when none are actually present. The thickness of APOSTLE discs compounds this effect: more slowly-rotating extra-planar gas systematically reduces the average line-of-sight speeds. The recovered rotation curves thus tend to underestimate the true circular velocity of APOSTLE galaxies in the inner regions. Non-circular motions provide an appealing explanation for the large apparent cores observed in galaxies such as DDO 47 and DDO 87, where the model residuals suggest that such motions might have affected estimates of the inner circular velocities. Although residuals from tilted ring models in the simulations appear larger than in observed galaxies, our results suggest that non-circular motions should be carefully taken into account when considering the evidence for dark matter cores in individual galaxies.

Tables

Table 1. Table 1: Summary of the key parameters of the APOSTLE simulations used in this work. Particle masses vary by up to a factor of 2 2 2 between volumes at a fixed resolution ‘level’; the median values below are indicative only (see Fattahi et al., 2016 , for full details) . Details of the WMAP-7 cosmological parameters used in the simulations are available in Komatsu et al. ( 2011 ) .

	Particle masses ($M_{\odot }$)		Max softening
Simulation	DM	Gas	length ($\mathrm{pc}$)
AP-L3	$7.3\times {10}^6$	$1.5\times {10}^6$	$711$
AP-L2	$5.8\times {10}^5$	$1.2\times {10}^5$	$307$
AP-L1	$3.6\times {10}^4$	$7.4\times {10}^3$	$134$

Table 2. Table 2: Selected properties of the synthetically observed APOSTLE galaxies in our sample. Columns: (1)  Number corresponding to labels in Figs.  1 , 3 & 9 , entries in bold face are those referred to as ‘equilibrium’ systems in the text and shown in Fig.  9 ; (2)  Object label as described in Sec.  3.1 ; (3)  Neutral hydrogen mass (the gas mass used in the BTFR is M bar = M ⋆ + 1.4 ​ M HI subscript 𝑀 bar subscript 𝑀 ⋆ 1.4 subscript 𝑀 HI M_{\rm bar}=M_{\star}+1.4M_{\rm HI} ); (4)  Stellar mass; (5)  H  i  size defined as the radius where Σ HI subscript Σ HI \Sigma_{\rm HI} drops to 1 ​ M ⊙ ​ pc − 2 1 subscript M direct-product superscript pc 2 1\,{\rm M}_{\odot}\,{\rm pc}^{-2} , measured from the 0 th superscript 0 th 0^{\rm th}  moment maps; (6)  Maximum circular velocity ( V circ = G M ( < r ) / r V_{\rm circ}=\sqrt{GM(<r)/r} ); (7)  Circular velocity at 2 ​ kpc 2 kpc 2\,{\rm kpc} ; (8)  H  i  gas azimuthal velocity at 2 ​ kpc 2 kpc 2\,{\rm kpc} ; (9)  Rotation velocity at 2 ​ kpc 2 kpc 2\,{\rm kpc} as fit using 3D barolo .

Symbol		$M_{\mathrm{HI}}$	$M_{\star }$	$R_{\mathrm{HI}}$	$\max \left(V_{\mathrm{circ}}\right)$	$V_{\mathrm{circ}}\left(2\mathrm{kpc}\right)$	$V_{\mathrm{gas}}\left(2\mathrm{kpc}\right)$	$V_{\mathrm{rot}}\left(2\mathrm{kpc}\right)$
Number	Object	$\left(M_{\odot }\right)$	$\left(M_{\odot }\right)$	$\left(\mathrm{kpc}\right)$	$\left(\mathrm{km}{s}^{-1}\right)$	$\left(\mathrm{km}{s}^{-1}\right)$	$\left(\mathrm{km}{s}^{-1}\right)$	$\left(\mathrm{km}{s}^{-1}\right)$
$\mathrm{𝟏}$	AP-L1-V1-4-0	$2.6\times {10}^9$	$2.6\times {10}^9$	$23.8$	$91.0$	$68.2$	$60.1$	$51.4$
$\mathrm{𝟐}$	AP-L1-V1-6-0	$1.2\times {10}^8$	$8.8\times {10}^8$	$3.3$	$60.2$	$45.5$	$45.4$	$32.8$
$\mathrm{𝟑}$	AP-L1-V1-7-0	$3.0\times {10}^8$	$6.7\times {10}^8$	$4.2$	$72.0$	$61.6$	$59.4$	$41.0$
$\mathrm{𝟒}$	AP-L1-V1-8-0	$2.9\times {10}^8$	$7.7\times {10}^8$	$3.4$	$68.2$	$52.0$	$52.2$	$50.8$
$5$	AP-L1-V4-6-0	$4.4\times {10}^8$	$1.4\times {10}^9$	$7.5$	$86.4$	$74.8$	$62.6$	$25.2$
$\mathrm{𝟔}$	AP-L1-V4-8-0	$6.1\times {10}^8$	$8.1\times {10}^8$	$8.7$	$69.1$	$48.0$	$41.4$	$29.5$
$7$	AP-L1-V4-10-0	$2.2\times {10}^9$	$4.9\times {10}^8$	$17.5$	$66.2$	$45.5$	$25.3$	$25.0$
$8$	AP-L1-V4-13-0	$3.9\times {10}^8$	$2.7\times {10}^8$	$8.5$	$64.8$	$53.3$	$35.3$	$31.2$
$\mathrm{𝟗}$	AP-L1-V4-14-0	$1.0\times {10}^9$	$4.2\times {10}^8$	$11.8$	$60.2$	$44.7$	$40.5$	$31.4$
$10$	AP-L1-V6-5-0	$1.3\times {10}^9$	$2.3\times {10}^9$	$12.4$	$89.5$	$72.4$	$55.0$	$71.1$
$\mathrm{𝟏𝟏}$	AP-L1-V6-6-0	$2.1\times {10}^9$	$7.6\times {10}^8$	$10.3$	$68.0$	$46.4$	$39.7$	$31.1$
$\mathrm{𝟏𝟐}$	AP-L1-V6-7-0	$1.5\times {10}^9$	$5.1\times {10}^8$	$12.7$	$70.2$	$43.9$	$48.7$	$29.9$
$\mathrm{𝟏𝟑}$	AP-L1-V6-8-0	$2.8\times {10}^9$	$8.4\times {10}^8$	$21.2$	$75.9$	$55.6$	$48.0$	$43.3$
$14$	AP-L1-V6-11-0	$3.1\times {10}^8$	$7.5\times {10}^8$	$5.6$	$60.3$	$49.8$	$33.0$	$16.8$
$15$	AP-L1-V6-12-0	$1.7\times {10}^9$	$1.5\times {10}^9$	$16.3$	$76.6$	$61.8$	$51.3$	$38.8$
$16$	AP-L1-V6-15-0	$6.3\times {10}^8$	$3.4\times {10}^8$	$8.4$	$61.6$	$46.0$	$32.2$	$52.5$
$17$	AP-L1-V6-16-0	$4.1\times {10}^8$	$5.3\times {10}^8$	$7.9$	$65.5$	$55.0$	$44.1$	$58.8$
$\mathrm{𝟏𝟖}$	AP-L1-V6-18-0	$1.9\times {10}^8$	$3.0\times {10}^8$	$4.4$	$61.9$	$56.6$	$53.5$	$57.2$
$\mathrm{𝟏𝟗}$	AP-L1-V6-19-0	$6.5\times {10}^8$	$5.2\times {10}^8$	$9.2$	$61.0$	$46.7$	$50.7$	$26.0$
$\mathrm{𝟐𝟎}$	AP-L1-V6-20-0	$2.3\times {10}^7$	$3.8\times {10}^8$	$1.7$	$67.7$	$59.9$	$62.3$	$73.5$
$\mathrm{𝟐𝟏}$	AP-L1-V10-5-0	$2.0\times {10}^9$	$6.1\times {10}^9$	$12.9$	$108.8$	$83.3$	$80.5$	$68.6$
$22$	AP-L1-V10-6-0	$2.1\times {10}^9$	$4.0\times {10}^9$	$13.3$	$103.9$	$71.2$	$59.9$	$30.8$
$23$	AP-L1-V10-13-0	$6.0\times {10}^8$	$2.0\times {10}^9$	$5.6$	$83.6$	$74.9$	$59.6$	$50.1$
$24$	AP-L1-V10-14-0	$1.4\times {10}^9$	$9.7\times {10}^8$	$13.6$	$66.3$	$51.0$	$43.1$	$28.5$
$25$	AP-L1-V10-16-0	$1.1\times {10}^9$	$9.6\times {10}^8$	$12.5$	$75.5$	$51.6$	$35.5$	$43.2$
$26$	AP-L1-V10-17-0	$4.4\times {10}^8$	$7.6\times {10}^8$	$6.4$	$67.4$	$48.5$	$38.5$	$20.5$
$27$	AP-L1-V10-19-0	$4.9\times {10}^8$	$4.9\times {10}^8$	$9.3$	$67.2$	$48.2$	$26.9$	$20.9$
$28$	AP-L1-V10-20-0	$3.0\times {10}^8$	$7.0\times {10}^8$	$1.9$	$73.5$	$62.3$	$43.0$	$47.4$
$29$	AP-L1-V10-22-0	$3.9\times {10}^8$	$7.6\times {10}^8$	$7.2$	$65.5$	$48.0$	$31.2$	$22.4$
$\mathrm{𝟑𝟎}$	AP-L1-V10-30-0	$6.4\times {10}^8$	$3.6\times {10}^8$	$9.4$	$61.4$	$49.5$	$45.2$	$22.4$
$31$	AP-L1-V11-3-0	$4.6\times {10}^9$	$9.5\times {10}^9$	$29.1$	$118.4$	$94.6$	$79.5$	$81.3$
$\mathrm{𝟑𝟐}$	AP-L1-V11-5-0	$4.7\times {10}^9$	$3.2\times {10}^9$	$23.8$	$91.1$	$65.5$	$60.4$	$37.2$
$33$	AP-L1-V11-6-0	$4.1\times {10}^9$	$1.4\times {10}^9$	$28.6$	$88.5$	$67.6$	$47.3$	$43.2$

Table 3. Table 4: Parameters used to configure the 3D barolo  software, omitting parameters that have no impact on the result of the calculation (e.g., file paths, flags to enable optional outputs, etc.). Our parameter choices where they differ from the default values are inspired by the choices made by Iorio et al. ( 2017 ) . These authors model the LITTLE THINGS galaxies; given that we mimic the observing setup of this survey and our simulated galaxies of interest are broadly similar to those in this sample, many of their parameter choices are applicable here. The most significant change we make is to restrict the inclination and position angles somewhat closer to their ‘true’ values.

Parameter	Value	Units	Description & comments
CHECKCHANNELS	FALSE	–	Check for bad channels in the data cube?
FLAGROBUSTSTATS	TRUE	–	Use robust statistics?
FLAGSEARCH	FALSE	–	Search for sources in the data cube?
FLAGRING	FALSE	–	Fit velocity field with a ring model?
SMOOTH	FALSE	–	Smooth the data cube?
GALFIT	TRUE	–	Fit a 3D model to the data cube?
BOX	NONE	$\mathrm{px}$	Select a sub-region of the cube?
NRADII	varies	–	Number of rings to use; we use enough rings to reach the radius enclosing 90 per cent of the H i mass of the galaxy.
RADSEP	$14.101$	$\mathrm{arcsec}$	Separation of rings.
XPOS	varies	$\mathrm{px}$	Centre of rings; set individually for each galaxy at the projected centre of the stellar light distribution.
YPOS	varies	$\mathrm{px}$	As XPOS.
VSYS	$257.4528$	$\mathrm{km}{s}^{-1}$	Systemic velocity.
VROT	$30$	$\mathrm{km}{s}^{-1}$	Initial guess for rotation velocity.
VDISP	$8$	$\mathrm{km}{s}^{-1}$	Initial guess for velocity dispersion.
INC	$60$	$\mathrm{degrees}$	Initial guess for inclination.
DELTAINC	$15$	$\mathrm{degrees}$	Allowed deviation of inclination from initial guess.
PA	$270$	$\mathrm{degrees}$	Initial guess for position angle.
DELTAPA	$20$	$\mathrm{degrees}$	Allowed deviation of position angle from initial guess.
Z0	$2.136$	$\mathrm{arcsec}$	Disc scale height.
DENS	$-1$	$\mathrm{atoms}{\mathrm{cm}}^{-2}$	Global column density of gas (unused when NORM LOCAL is set).
FREE	VROT VDISP INC PA	–	Parameters to fit for each ring.
MASK	SMOOTH	–	The data cube is smoothed by a factor of $2$ and a signal-to-noise cut is used to define a mask.
BLANKCUT	$2.5$	–	Signal-to-noise threshold for mask construction.
SIDE	B	–	The entire galaxy is modelled; the approaching and receding sides can also be fit separately.
NORM	LOCAL	–	The model is normalized pixel by pixel, i.e. the surface brightness is not explicitly fit.
LTYPE	$1$	–	Layer type along z is Gaussian.
FTYPE	$2$	–	Minimization function is $\left|\mathrm{model}-\mathrm{observed}\right|$.
WFUNC	$1$	–	Azimuthal weighting function is $\left|\cos \theta \right|$.
TOL	$0.001$	–	Minimization tolerance.
TWOSTAGE	TRUE	–	Two stage fitting, i.e. geometric parameters are regularized and rotation velocity is fit again.
POLYN	bezier	–	Degree of polynomial fitting INC and PA.
FLAGERRORS	varies	–	Errors are estimated only for those galaxies we use for illustrative purposes.
BWEIGHT	$1$	–	Exponent of weight for blank pixels.
LINEAR	$0.424$	$\mathrm{channels}$	Instrumental spectral broadening (standard deviation).
CDENS	$10$	–	Number of clouds to use in building the ring model.
NV	$200$	–	Number of sub-clouds used within each cloud (see CDENS) to populate the spectral axis of the ring model.

Equations

$$V_{\mathrm{LoS}}(\varphi )=V_{\mathrm{sys}}+V_0\cos (\varphi -{\varphi }_0)$$

$$V_{\mathrm{circ}}^2=V_{\mathrm{rot}}^2-{\sigma }^2\frac{\mathrm{d}\log ({\Sigma }_{\mathrm{HI}}{\sigma }^2)}{\mathrm{d}\log \mathrm{R}}$$

$$\frac{{\Sigma }_{\mathrm{HI}}{\sigma }^2}{({\Sigma }_{\mathrm{HI}}{\sigma }^2{)}_0}=\frac{(R_0+1)}{R_0+e^{\alpha R}}$$

$$\sigma =\sqrt{\frac{k_{\mathrm{B}}T}{\mu m_{\mathrm{p}}}+\frac{1}{3}\left({\sigma }_{\varphi }^2+{\sigma }_r^2+{\sigma }_z^2\right)}$$

$$V(\varphi )=\sum _m{V}_m\cos [m(\varphi -{\varphi }_m))]$$

$$V_{\mathrm{LoS}}(\varphi )=\sin i\cos \varphi \sum _m{V}_m\cos [m(\varphi -{\varphi }_m)]$$

$$\begin{gathered} V_{\mathrm{LoS}}(\varphi )=\sin i\sum _m\frac{V_m}{2}[\cos ((m-1)\varphi -m{\varphi }_m) \\ +\cos ((m+1)\varphi -m{\varphi }_m)] \end{gathered}$$

$$V(\varphi )=V_0\sin (i)\cos (\varphi )\left[1+\frac{V_2}{V_0}\cos (2(\varphi -{\varphi }_2))\right]$$

$$\Delta V_{\mathrm{LoS}}(\varphi )=V_2\sin (i)\cos (\varphi )\cos (2(\varphi -{\varphi }_2)),$$

$$\Delta V_{\mathrm{LoS}}(\varphi )=\frac{V_2}{2}\sin (i)\left[\cos (3\varphi -2{\varphi }_2)+\cos (\varphi -2{\varphi }_2)\right].$$

Full text

1 Introduction

The ‘cusp-core problem’ is a long-standing controversy that arises when contrasting the steep central density profiles of cold dark matter haloes (‘cusps’) predicted by N-body simulations (Navarro et al., 1996b, 1997) with the mass profiles inferred from disc galaxy rotation curves after subtracting the contribution of the stellar and gaseous (baryonic) components (Flores & Primack, 1994; Moore, 1994, and see de Blok, 2010 for a review). The comparison usually involves fitting a power-law density profile to the dark matter contribution in the innermost resolved region of the rotation curve. The power-law slope is then compared with that of simulated CDM haloes at similar distances from the centre.

Although legitimate in principle, this procedure is in practice fraught with difficulties. One difficulty is that rotation speeds rapidly approach zero near the centre, which implies that inferences about the dark matter cusp are made by fitting small rotation velocities at small radii, a regime where even small errors can have a large influence on the results. A further difficulty is that, in order to recover the dark mass profile, one must remove the baryonic contribution to the circular velocity, which involves making assumptions about relatively poorly known parameters, such as, for example, the mass-to-light ratio of the stars, and the $X_{\rm CO}$ parameter used to infer the molecular hydrogen distribution from CO observations. However, even with extreme assumptions for the values of these parameters, the recovered innermost slopes are in general shallower than predicted for CDM haloes (e.g. de Blok & McGaugh, 1997; Swaters, 1999; Oh et al., 2011).

The reason for this emphasis on the innermost slope may be traced to early N-body simulation work (Dubinski & Carlberg, 1991; Navarro et al., 1996b, 1997), which reported that the the central cusp of CDM halo density profiles was at least as steep as $\rho\propto r^{-1}$. Any slope measured to be shallower than that could then argued to be in conflict with CDM, an idea that has guided many studies of inner rotation curves since.

Our understanding of this issue, however, has evolved, largely as a result of improved simulations and of better constraints on the cosmological parameters. Indeed, observations of the cosmic microwave background and of large-scale galaxy clustering have now constrained the cosmological parameters to great precision (Planck Collaboration, 2016). At the same time, cosmological simulations have improved to the point that the halo scaling parameters, their dependence on mass, as well as their scatter, are now well understood (e.g. Ludlow et al., 2014).

Because of these advances, the dark matter contribution to a rotation curve can now be predicted at essentially all relevant radii once a single parameter is specified for a halo, such as the maximum circular velocity, $V_{\rm max}$. Since baryons can dissipate and add mass to the inner regions, this dark matter contribution may be regarded as a minimum mass (or, equivalently, a minimum circular velocity) at each radius. The advantage of focussing on this minimum velocity is that it allows the theoretical predictions to be directly confronted with observations at radii not too close to the centre, where rotation curves are less prone to uncertainty and where the theoretical predictions are less vulnerable to artefact.

This is the approach we adopted in an earlier study (Oman et al., 2015), where we proposed to reformulate the cusp-core problem as an ‘inner mass deficit’ problem that afflicts galaxies where the inner circular velocities fall below the minimum expected from the dark matter alone. At a fiducial radius111We note that there is nothing exceptional about this choice; it is just a compromise radius that both simulations and observations resolve well for a wide range of galaxy masses. of $2\,{\rm kpc}$, a number of galaxies have rotation speeds much lower than predicted, given their $V_{\rm max}$. As discussed in that paper, this mass deficit does not affect all galaxies but it does affect galaxies of all masses, from dwarfs to massive discs, and varies widely from galaxy to galaxy at fixed $V_{\rm max}$.

In some cases, the deficit at $2\,{\rm kpc}$ is so pronounced that it far exceeds the total baryonic mass of the system. This argues (at least in those extreme examples) against the idea that the deficit is caused by dark matter ‘cores’ produced by the baryonic assembly of the galaxy (see, e.g., Navarro et al., 1996a; Read & Gilmore, 2005; Mashchenko et al., 2008; Governato et al., 2012; Pontzen & Governato, 2012; Brooks & Zolotov, 2014; Oñorbe et al., 2015; Tollet et al., 2016, see also the review of Pontzen & Governato, 2014). Indeed, these baryon-induced modifications are in general modest, and the effect restricted to a small range of galaxy masses (Di Cintio et al., 2014; Chan et al., 2015). In addition, we note that baryon-induced cores are not a general prediction, but rather a result of some implementations of star formation and feedback in galaxy formation simulations. Indeed, simulations like those from the EAGLE (Schaye et al., 2015), APOSTLE (Sawala et al., 2016) and Illustris (Vogelsberger et al., 2014) projects are able to reproduce most properties of the galaxy population without producing any such cores.

An alternative is that gas rotation curves do not faithfully trace the circular velocity in the inner regions of some galaxies and have therefore been erroneously interpreted as evidence for ‘cores’. This possibility has been repeatedly raised in the past: early observations were subject to concerns around beam smearing in the case of H i data (Swaters et al., 2009, and references therein), and centering, alignment, and seeing in the case of H $\alpha$ slit spectroscopy (Swaters et al., 2003; Spekkens et al., 2005). Such worries have largely been laid to rest with the advent of new, high-resolution datasets that often combine H $\alpha$, H i, and CO maps to yield 2D gas velocity fields of excellent angular resolution (e.g., Kuzio de Naray et al., 2006; Walter et al., 2008; Hunter et al., 2012; Adams et al., 2014).

Other concerns, however, remain. Observed velocity data must be processed through a model before they can be contrasted with theoretical predictions, and a number of modelling issues have yet to be properly understood. The main purpose of modelling the data is to infer the speed of a hypothetical circular orbit, $V_{\rm circ}(r)$, which can then be directly compared with the mass distribution interior to radius $r$ predicted by the simulations. Observations, however, can at best only constrain the average azimuthal speed of the gas at each radius, usually referred to as the ‘rotation speed’, $V_{\rm rot}(r)$. In general, $V_{\rm rot}\neq V_{\rm circ}$, and corrections must be applied.

One such correction concerns the support provided by pressure gradients and velocity dispersion of the gas. This depends on the gas surface density profile, as well as on the gas velocity dispersion and its radial gradient, in a manner akin to the familiar ‘asymmetric drift’ that affects the average rotation speed of stars in a disc (e.g. Valenzuela et al., 2007). Although the corrections are approximate, for galaxies with V_{\rm max}\mathrel{\raise 1.50696pt\hbox{\scriptstyle>}\kern-6.00006pt\lower 1.72218pt\hbox{{\scriptstyle\sim}}}30\,{\rm km}\,{\rm s}^{-1} the changes they imply are usually too small to compromise the results (e.g. Oh et al., 2011).

Non-circular motions in the gas, on the other hand, are a greater concern. Although these are likely ubiquitous at some level, they are seldom considered in the modelling. One reason for this is that there is no simple and general way of assessing the effect of non-circular motions, which may affect both the estimates of $V_{\rm rot}$ and the translation of $V_{\rm rot}$ into $V_{\rm circ}$. As a result, $V_{\rm rot}(r)$ is often used as a direct measure of $V_{\rm circ}(r)$ without further correction. Although this practice may sometimes be acceptable, such as in the case of the 19 galaxy discs from the THINGS survey studied in detail by Trachternach et al. (2008), it can also lead to erroneous conclusions (as we shall see below) and must be carefully scrutinized for each individual galaxy.

The case of NGC 2976 provides a sobering example. Obvious asymmetries in the velocity field led Simon et al. (2003) to use a harmonic decomposition of the velocity field, where circular gas ‘rings’ were allowed to have non-zero radial velocity (i.e., they may be expanding or contracting), in addition to the usual rotation speed. With this assumption, a ‘tilted-ring’ model (Rogstad et al., 1974) can reproduce the observed gas velocity field quite accurately, yielding a well-defined mean azimuthal velocity as a function of radius, $V_{\rm rot}(r)$. (And, of course, non-zero radial velocities too.)

A visual inspection of the velocity field of NGC 2976 shows that it differs significantly in contiguous quadrants, but is roughly antisymmetric in diagonally opposite ones. This is a clear signature of eccentric, rather than circular, gas orbits (Hayashi & Navarro, 2006), which led Spekkens & Sellwood (2007) to model this galaxy assuming the presence of a radially-coherent bisymmetric ($m=2$) velocity pattern, as expected in a barred galaxy, or when the halo potential is triaxial. This model also fits the data quite well, but yields a rather different radial dependence for $V_{\rm rot}(r)$, especially near the centre.

In addition to this model degeneracy, in neither case does the derived ‘rotation curve’ $V_{\rm rot}(r)$ – the mean azimuthal rotation of the gas – trace the circular velocity, $V_{\rm circ}(r)$. Translating $V_{\rm rot}$ into $V_{\rm circ}$ in such cases requires a dynamical gas flow model of the whole galaxy, something that can only be accomplished when specific assumptions are made about the ellipticity of the gravitational potential, its dependence on radius, and/or the radial dependence of its phase angle. Given these complexities, it is not surprising that non-circular motions are usually not modelled in detail and instead treated as a source of error.

We explore these issues here by analysing synthetic H i observations of $33$ simulated central galaxies with $60<V_{\rm max}/{\rm km}\,{\rm s}^{-1}<120$ selected from the APOSTLE suite of cosmological hydrodynamical simulations (Fattahi et al., 2016; Sawala et al., 2016). We use 3Dbarolo, a tilted-ring processing tool that has recently been used to derive the rotation curves of galaxies in the LITTLE THINGS survey (Iorio et al., 2017). We begin in Sec. 2 with a brief review of earlier work on kinematic models of galactic discs. In Sec. 3 we briefly describe the simulations and how we construct our synthetic observations. We compare these with measurements of real galaxies in Sec. 3.4 to show that the two are similar in terms of several kinematic and symmetry metrics. In Sec. 4 we describe the process of fitting kinematic models to our synthetic observations. In Sec. 5 we present our main results: in Sec. 5.2 we demonstrate that the recovered rotation curve depends sensitively on the orientation of non-circular motions in the disc with respect to the line of sight, and in Sec. 5.6 we show that the mixing of different annuli and layers in the disc along the line of sight is a further source of error. In Sec. 5.5 we show that these non-circular motions can lead to ‘inner mass deficits’ comparable to those reported by Oman et al. (2015), and compare with observed galaxies in Sec. 6. Finally, we discuss the implications of our findings and summarize our conclusions in Sec. 7.

2 Prior work

We are hardly the first to suggest that non-circular motions can substantially impact the rotation curve modelling process. The pioneering work of Teuben (1991) was followed up by Franx et al. (1994) and Schoenmakers et al. (1997) to extract the signature, in projection, of perturbations to circular orbits. This work showed that, when expanded into harmonic modes, each mode of order $m$ gives rise to patterns of order $m\pm 1$ in projection.

Projection effects thus introduce degeneracies in the interpretation of non-circular motions, which add to others introduced by errors in geometric parameters such as the systemic velocity, centroid, inclination, and position angle of the kinematic principal axes. Some of the degeneracies may be lifted if the amplitude of the perturbations is small, which allows the epicyclic approximation to be used to constrain the various amplitudes and phases, but this approach is of little help when perturbations are a substantial fraction of the mean azimuthal velocity, as is often the case near the galaxy centre.

Rhee et al. (2004) discuss the kinematic modelling of a simulated barred galaxy (their ‘Model I’). In this case, the bar is strong enough to drive down the mean rotation velocity substantially. The authors find that the rotation curve they measure for the system depends sensitively on the orientation of the bar, with the apparent rotation falling far below the circular velocity of the system when the bar is aligned along the major axis of the galaxy in projection, just the case when asymmetries in the projected velocity field are hardest to detect. They also note that even very small systematic errors in the velocity, on the order of $10$ per cent, are enough to cause large changes in the inferred dark matter density profile slope.

A similar cautionary tale is told by Valenzuela et al. (2007), who analyse synthetic observations of simulations of isolated galaxies set up initially in equilibrium. They argue that non-circular motions, coupled with the extra support provided by gas pressure, are enough to explain the slowly rising rotation curves of NGC 3109 and NGC 6822, two galaxies where the inferred dark matter profiles depart substantially from that expected from $\Lambda$CDM. Their conclusion is that cuspy profiles are actually consistent with the data once these effects are taken into account.

Spekkens & Sellwood (2007) extended earlier work to account for perturbations that are large compared with the mean azimuthal velocity, a regime where the epicyclic approximation fails. Their model applies to bisymmetric ($m=2$) perturbations of radially-constant phase, as in the case of a bar or a triaxial dark halo. Applied to a real galaxy, the method enables estimates of the mean orbital speed from a velocity map, even when strongly non-axisymmetric. These authors also caution that this is just a measure of the average azimuthal speed around a circle, and not a precise indicator of centrifugal balance. In other words, $V_{\rm rot}\neq V_{\rm circ}$ when non-circular motions are not negligible. In this case, the only way to estimate $V_{\rm circ}(r)$ is to find a non-axisymmetric model that produces a fluid dynamical flow pattern matching the observed one.

Given these complexities, it is an illuminating exercise to ‘observe’ simulated galaxies and analyse them on an even footing with analogous observed data. This is the approach we adopt in this paper. Previous comparisons between synthetic rotation curves from gas dynamic simulations and real data include the recent work by Read et al. (2016), who analysed simulations of $4$ idealized galaxies which have dark matter cores, and by Pineda et al. (2017), who analysed a set of $6$ simulated isolated galaxies which retain their initial dark matter cusps. We draw a sample of galaxies from the APOSTLE suite of cosmological hydrodynamical simulations (Sawala et al., 2016; Fattahi et al., 2016), construct synthetic ‘observations’ of their H i gas kinematics, and apply the same ‘tilted-ring’ modelling commonly adopted to analyse the kinematics of nearby discs. None of the APOSTLE galaxies have ‘cores’, which simplifies the interpretation: any possible deviations between the recovered $V_{\rm rot}(r)$ and the known $V_{\rm circ}(r)$ are either due to the fact that the gas is not truly in centrifugal equilibrium, or to inadequacies in the modelling of the data.

3 Simulations

3.1 The APOSTLE simulations

The APOSTLE222A Project Of Simulations of The Local Environment. simulation suite comprises $12$ volumes selected from a large cosmological volume and resimulated using the zoom-in technique (Frenk et al., 1996; Power et al., 2003; Jenkins, 2013) with the full hydrodynamics and galaxy formation treatment of the ‘Ref’ model of the EAGLE project (Schaye et al., 2015; Crain et al., 2015). The regions are selected to resemble the Local Group of galaxies in terms of the mass, separation and kinematics of two haloes analogous to the Milky Way and M 31, and relative isolation from more massive systems. Full details of the simulation setup and target selection are available in Fattahi et al. (2016) and Sawala et al. (2016); we summarize a few key points here.

EAGLE, and by extension APOSTLE, use the pressure-entropy formulation of smoothed particle hydrodynamics (Hopkins, 2013) and the numerical methods from the anarchy module [Dalla Vecchia et al. (in preparation); see Schaye et al., 2015 for a short summary]. The galaxy formation model includes subgrid recipes for radiative cooling (Wiersma et al., 2009a), star formation (Schaye, 2004; Schaye & Dalla Vecchia, 2008), stellar and chemical enrichment (Wiersma et al., 2009b), energetic stellar feedback (Dalla Vecchia & Schaye, 2012), and cosmic reionization (Haardt & Madau, 2001; Wiersma et al., 2009b). The stellar feedback is calibrated to reproduce approximately the galaxy stellar mass function and the sizes of $M_{\star}>10^{8}\,{\rm M}_{\odot}$ galaxies at $z=0$. EAGLE also includes black holes and AGN feedback. In APOSTLE these processes are neglected, which is reasonable for the mass scale of interest here (Crain et al., 2015; Bower et al., 2017).

The APOSTLE volumes are simulated at three resolution levels, labelled AP-L3 (similar to the fiducial resolution of the EAGLE suite), AP-L2 (similar to the ‘high resolution’ EAGLE runs) and an even higher resolution level, AP-L1. Each resolution level represents an increase by a factor of $\sim 12$ in mass and $\sim 2$ in force softening over the next lowest level. Typical values of particle mass, gravitational softening, and other numerical parameters vary slightly from volume to volume; representative values are shown in Table 1. All $12$ volumes have been simulated at AP-L2 and AP-L3 resolution, but only $5$ volumes: V1, V4, V6, V10 & V11 have thus far been simulated at AP-L1. APOSTLE assumes the WMAP-7 cosmological parameters (Komatsu et al., 2011): $\Omega_{\rm m}=0.2727$, $\Omega_{\Lambda}=0.728$, $\Omega_{\rm b}=0.04557$, $h=0.702$, $\sigma_{8}=0.807$.

The subfind algorithm (Springel et al., 2001; Dolag et al., 2009) is used to identify structures and galaxies in the APOSTLE volumes. Particles are first grouped into friend-of-friends (FoF) haloes by iteratively linking particles separated by at most $0.2\times$ the mean interparticle separation (Davis et al., 1985); gas and star particles are attached to the same FoF halo as their nearest dark matter particle. Saddle points in the density distribution are then used to separate substructures, and particles which are not gravitationally bound to substructures are removed. The end result is a collection of groups each containing at least one ‘subhalo’; the most massive subhalo in each group is referred to as ‘central’; others are ‘satellites’. In this analysis we focus exclusively on central objects as satellites are subject to additional dynamical processes which complicate their interpretation.

We label our simulated galaxies according to the resolution level, volume number, FoF group and subgroup, so for instance AP-L1-V1-8-0 corresponds to resolution AP-L1, volume V1, FoF group 8 and subgroup 0 (the central object). We focus primarily on the AP-L1 resolution. At this resolution level the circular velocity curves of our galaxies of interest (defined below) are numerically converged at all radii \mathrel{\raise 1.50696pt\hbox{\scriptstyle>}\kern-6.00006pt\lower 1.72218pt\hbox{{\scriptstyle\sim}}}700\,{\rm pc}, as defined by the criterion of Power et al. (2003, for further details pertaining to the numerical convergence of the APOSTLE simulations see , ).

3.2 Galaxy sample selection

We select $33$ galaxies from the APOSTLE simulations for further consideration based on two criteria. First, as noted above, we restrict ourselves to the highest AP-L1 resolution level so that the central regions of the galaxies, which are of particular interest in the context of the cusp-core problem, are sufficiently well resolved. Second, we choose galaxies in the interval $60<V_{\rm max}/{\rm km}\,{\rm s}^{-1}<120$, where $V_{\rm max}={\rm max}(V_{\rm circ}(R))$. The lower bound ensures that the gas distribution of the galaxies is well-sampled (\mathrel{\raise 1.50696pt\hbox{\scriptstyle>}\kern-6.00006pt\lower 1.72218pt\hbox{{\scriptstyle\sim}}}10^{4} gas particles contribute to the H i distribution of each galaxy). The upper bound is chosen to exclude massive galaxies were the baryonic component dominates the kinematics. Indeed, the kinematics of all galaxies in our sample are largely dictated by their dark matter haloes.

APOSTLE galaxies have realistic masses, sizes, and velocities. This is shown in Fig. 1, where we compare the simulated sample (large black points) with observational data from the Spitzer Photometry and Accurate Rotation Curves (SPARC) database (Lelli et al., 2016) and the THINGS (Walter et al., 2008) and LITTLE THINGS (Hunter et al., 2012) surveys. For the SPARC galaxies, we plot only those with the highest quality flag ($Q=1$). For the THINGS and LITTLE THINGS surveys, we plot only those galaxies which were selected for kinematic modelling by the survey teams (THINGS: de Blok et al., 2008; Oh et al., 2011; LITTLE THINGS: Oh et al., 2015).

APOSTLE galaxies comfortably match three key scaling relations. The left panel of Fig. 1 shows the baryonic Tully-Fisher relation (BTFR). The quantity plotted on the horizontal axis varies by dataset: for APOSTLE galaxies we show the maximum of the circular velocity curve, $V_{\rm circ}(r)=\sqrt{GM(<r)/r}$, for SPARC galaxies we show the asymptotically flat rotation velocity, whereas for THINGS & LITTLE THINGS galaxies we show the maximum of the rotation curve. The baryonic masses are, in all cases, calculated as $M_{\rm bar}=M_{\star}+1.4M_{\rm HI}$ (e.g. McGaugh, 2012, and see Sec. 3.3 for the method used to calculate the H i masses). Our selection in $V_{\rm max}$ is highlighted by the shaded vertical band. It is clear from this panel that the BTFR of APOSTLE galaxies is in good agreement with the observed scaling, provided that the observed velocities trace the maximum circular velocity of the halo (see; e.g., Sales et al., 2017; Oman et al., 2016, for a more in-depth discussion of this point).

The middle panel shows the H i mass – stellar mass relation. The simulated galaxies once again lie comfortably within the scatter of the observed relation. The right panel shows the H i mass – size relation, where the size is defined as the radius at which the H i surface density, $\Sigma_{\rm HI}$, drops below $1\,{\rm M}_{\odot}\,{\rm pc}^{-2}\,(\approx 10^{20}\,{\rm atoms}\,{\rm cm}^{-2})$. APOSTLE galaxies seem to have, at fixed H i mass, slightly larger sizes ($\sim 0.2\,{\rm dex}$) than observed. The offset in size shown in the right-hand panel of Fig. 1 should be of little consequence to our analysis.

3.3 Synthetic H i data cubes

For each simulated galaxy in our sample we carry out a synthetic H i observation, as follows. First, we compute an H i mass fraction for each gas particle in the central galaxy, following the prescription of Rahmati et al. (2013) for self-shielding from the metagalactic ionizing background radiation, and including an empirical pressure-dependent correction for the molecular gas fraction, as detailed in Blitz & Rosolowsky (2006).

Second, we adopt a coordinate system centred on the potential minimum of the galaxy, and choose a $z$-axis aligned with the direction of $\vec{L}_{\rm HI}$, the specific angular momentum vector of the H i gas disc. The velocity coordinate frame is chosen such that the average (linear) momentum of the H i gas in the central $500\,{\rm pc}$ is zero. A viewing angle inclined by $60^{\circ}$ relative to the $z$-axis is adopted, with random azimuthal orientation. Each galaxy is placed in the Hubble flow at a nominal distance of $3.7\,{\rm Mpc}$, the median distance of galaxies in the LITTLE THINGS sample (Hunter et al., 2012). We choose an arbitrary position on the sky at $(0^{\rm h}\,0^{\rm m}\,0.0^{\rm s},+10^{\circ}\,0^{\prime}\,0.0^{\prime\prime})$ and adopt an ‘observing setup’ similar to that used in the LITTLE THINGS survey, with a $6\,{\rm arcsec}$ circular Gaussian beam and $1024^{2}$ pixels spaced $3\,{\rm arcsec}$ apart. This yields an effective physical resolution (FWHM) of $\sim 110\,{\rm pc}$. We use a velocity channel spacing of $4\,{\rm km}\,{\rm s}^{-1}$ and enough channels to accommodate comfortably all of the galactic H i emission.

The gas particles are spatially smoothed with the $C^{2}$ Wendland (1995) smoothing kernel used in the EAGLE model. The integral of the kernel over each pixel is approximated by the value at the pixel centre. Provided the pixel size is $\leq\frac{1}{2}$ the smoothing length, this approximation is accurate to better than 1 per cent; we explicitly verify that this condition is satisfied. We also verified that omitting this smoothing step does not significantly change our main results.

In the velocity direction, the $21\,{\rm cm}$ emission is modelled with a Gaussian line profile centred at the particle velocity and a fixed width of $7\,{\rm km}\,{\rm s}^{-1}$, which models the (unresolved) thermal broadening of the H i line (e.g. Pineda et al., 2017). Our main results are insensitive to the precise width we choose for the line, provided it is \mathrel{\raise 1.50696pt\hbox{\scriptstyle<}\kern-6.00006pt\lower 1.72218pt\hbox{{\scriptstyle\sim}}}12\,{\rm km}\,{\rm s}^{-1}, because then the integrated H i profile is dominated by the dispersion in the particle velocities. Each particle contributes flux proportionally to its H i mass, i.e. the gas is assumed to be optically thin. Finally, the synthetic data cube is convolved along the spatial axes with the ‘beam’, implemented as a $6\,{\rm arcsec}$ circular Gaussian kernel. The completed cube is saved in the fits format (Pence et al., 2010) with appropriate header information.

In Fig. 2 we illustrate the synthetic observations of three of our simulated galaxies. The left column shows the surface density ($0^{\rm th}$ moment) maps. The red contour marks the $\log_{10}(\Sigma_{\rm HI}/{\rm atoms}\,{\rm cm}^{-2})=19.5$ isodensity contour. This is about $0.5\,{\rm dex}$ deeper than the typical limiting depth of observations in the THINGS and LITTLE THINGS surveys of $\sim 10^{20}\,{\rm atoms}\,{\rm cm}^{-2}$. We choose this because galaxies in our sample are slightly larger than observed ones, by roughly $\sim 0.2\,{\rm dex}$ in $M_{\rm HI}$–$R_{\rm HI}$. In light of this, a slightly deeper nominal limiting column density allows for more reasonable comparisons than a strict cut at $10^{20}\,{\rm atoms}\,{\rm cm}^{-2}$. The central column shows the line-of-sight velocity ($1^{\rm st}$ moment) maps333We show intensity weighted mean (IWM) velocity fields. The choice of velocity field type can have a significant impact on the fit rotation curve for techniques that model the velocity field directly (de Blok et al., 2008). For our purposes, however, the choice of velocity field impacts only the visualization of the data because our model of choice, 3Dbarolo, models the full data cube., and the right column the velocity dispersion ($2^{\rm nd}$ moment) maps.

3.4 Kinematics properties of simulated and observed galaxies

Are the kinematic properties of simulated galaxies broadly consistent with observed ones? We have already seen in Fig. 1 that APOSTLE galaxies have structural parameters that follow scaling laws similar to observed discs, but it is important to check that they also resemble observations in their internal kinematics. We explore this using three simple metrics that we can apply both to the publicly available moment maps444We use the ‘robust weighted’, not the ‘natural weighted’, maps (de Blok et al., 2008), though both give very similar results. of observed galaxies as well as to our simulated data cubes with minimal extra processing. The observational maps are provided cleaned of noise, with low signal-to-noise pixels masked out. We approximate this by masking in the simulated maps all pixels where the H i column density drops below $10^{19.5}\,{\rm atoms}\,{\rm cm}^{-2}$ (see Fig. 2).

The first metric is the median velocity dispersion along the line of sight (i.e. the median of all unmasked pixels in the $2^{\rm nd}$ moment map) as a function of H i mass, which we show in the left panel of Fig. 3. APOSTLE galaxies are shown by grey/black symbols, whereas galaxies from the THINGS and LITTLE THINGS surveys are shown with blue squares and diamonds (note that we plot only those galaxies regular enough to have been selected for mass modelling by the survey teams). Each symbol has a number that identifies the galaxy as listed in Tables A and A. At given H i mass the simulated galaxies have slightly larger velocity dispersions than observed galaxies, but the difference is less than a factor of two on average.

The second metric estimates the symmetry of the $1^{\rm st}$ moment maps (i.e., the line-of-sight velocity field). This is computed by rotating a map by $180^{\circ}$ about the galaxy centre and subtracting it from the unrotated fields (with a change of sign so that in the perfectly symmetric case the residual would be zero everywhere). The mean of the residual map indicates whether there is an offset in the average velocity of the approaching and receding sides of the galaxy; its rms is a crude estimate of the lack of circular symmetry of the velocity field. As shown by the middle panel of Fig. 3, APOSTLE and observed galaxies seem to deviate from perfect axisymmetry by similar amounts.

Finally, the third metric uses a measure of the residuals produced by subtracting a very simple kinematic model from the $1^{\rm st}$ moment map. Assuming a single inclination, position angle and systemic velocity for each galaxy (as listed in Table A), we fit the function:

[TABLE]

to a series of concentric, inclined ‘rings’ (ellipses in projection). We use the same ring spacings as de Blok et al. (2008); Oh et al. (2011, 2015), typically about $130\,{\rm pc}$. $V_{0}$ and $\phi_{0}$ are free555The freedom in $\phi_{0}$ means that, strictly speaking, we are not removing a pure rotation field. We recall that the purpose of this measurement is to compare synthetic and real data cubes, and the measurement is made identically in both cases. parameters fit to each ring independently. The residual map is then analysed as for the preceding metric: its rms is shown as a function of H i mass in the right-hand panel of Fig. 3. As in the other cases, the simulated and observed galaxies are nearly indistinguishable according to these metrics.

These results, together with those shown earlier in Fig. 1, give us confidence that the kinematics of the simulated galaxies are, to zeroth order, similar enough to those of their observed counterparts to warrant applying similar analysis tools.

4 Kinematic modelling

4.1 Tilted-ring model

The standard tool for kinematic modelling of disc galaxy velocity fields is known as a ‘tilted-ring’ model (Rogstad et al., 1974). In such a model, a disc is represented as a series of rings of increasing size. The properties of each ring are described by a set of parameters which can be categorized as geometric (radius, width, thickness, centroid, inclination, position angle, systemic velocity) and physical (surface density, rotation velocity, velocity dispersion). A number of publicly-available tilted-ring models exist; we use here the 3Dbarolo666http://editeodoro.github.io/Bbarolo/, we used the latest version available at the time of writing: 1.3 (github commit d54e901). software package (for a detailed description see Di Teodoro & Fraternali, 2015).

Whereas most older versions of tilted-ring models only use the first few moments of the kinematics – the surface density and velocity fields, and in some cases the velocity dispersion field – 3Dbarolo belongs to a class of more recent tools that model the full data cube directly, and therefore nominally utilize all available kinematic information. The software has many configurable parameters; we discuss our choices for several of the most important ones below, and in Table 4 we summarize the full configuration used.

4.2 Parameter choices

The most important parameters of the model are those that define the handling of the geometric parameters of each ring. When applied to projections of APOSTLE galaxies, and in order to facilitate convergence, we provide 3Dbarolo with a ‘correct’ guess of $i=60^{\circ}$ for the inclination angle (and allow it to deviate by no more than $15^{\circ}$ from this value). We also initialize the software with the ‘correct’ guess for the position angle of the rings ($270^{\circ}$ counter-clockwise from North), and allow deviations of no more than $20^{\circ}$. Providing reasonably accurate initial guesses (within $\sim 15^{\circ}$) for these two parameters is, unfortunately, necessary for the fitting procedure to converge to a correct solution (Di Teodoro & Fraternali, 2015). For real galaxies, these must be estimated from the geometry of either the gas or stellar distribution. The inclination and position angles that would be derived from the gas isodensity contours for our sample of APOSTLE galaxies typically differ from the ‘true’ values by less than the maximum variations we allow in the fitting routine. The ring widths are fixed at $14.1\,{\rm arcsec}$, corresponding to a physical separation of $250\,{\rm pc}$ at the distance of $3.7\,{\rm Mpc}$ chosen for our synthetic observations.

We fix the centre of each ring to the density peak of the projected stellar distribution of the galaxy. This coincides, within a few pixels ($<3\,{\rm px}\sim 10\,{\rm arcsec}$), with the minimum potential centre returned by the subfind algorithm (Springel et al., 2001; Dolag et al., 2009). For simplicity, the systemic velocity is fixed at $257\,{\rm km}\,{\rm s}^{-1}$, determined from the distance as $V_{\rm sys}=H_{0}D$. The initial guesses for the rotation speed and velocity dispersion of each ring are set to $30$ and $8\,{\rm km}\,{\rm s}^{-1}$, respectively. These initial guesses have little impact on the final fits to the rotation curve and velocity dispersion profile.

We fix the thickness of the rings at $2\,{\rm arcsec}=40\,{\rm pc}$. This is much thinner than the actual thicknesses of the simulated gas discs, where the half-mass height can reach $\sim 1\,{\rm kpc}$. Modelling thick discs is a well-known limitation of tilted ring models. Future codes may be able to capture better the vertical structure of discs (e.g. Iorio et al., 2017), but for the present we are bound by the limitations of current implementations.

We model each galaxy out to the radius enclosing 90 per cent of its H i mass. This roughly coincides with the $\log_{10}(\Sigma_{\rm HI}/{\rm atoms}\,{\rm cm}^{-2})=19.5$ isodensity contour, and is, in all cases, extended enough to reach the asymptotically flat (maximum) portion of the circular velocity curve.

4.3 Fitting procedure

Using the parameter choices outlined above (see also Table 4), the tilted ring model is fit to each galaxy in two stages (e.g. Iorio et al., 2017). In the first stage the free parameters are the rotation speed, velocity dispersion, inclination and position angle of each ring (in 3Dbarolo’s ‘locally normalized’ mode the surface brightness is not explicitly fit). The inclination and position angle profiles are then smoothed with a low-order polynomial fit and, in a second stage, the rotation speeds and velocity dispersions of the rings are fit again with the geometric parameters held fixed at their smoothed values.

4.4 Correction for pressure support

The procedure above yields the mean azimuthal velocity of the galaxy as a function of radius, $V_{\rm rot}(r)$. This is usually smaller than the true circular velocity because the gas may be partially supported by ‘pressure’ forces. We therefore correct the rotation speeds as in, e.g., Valenzuela et al. (2007):

[TABLE]

where $\Sigma_{\rm HI}$ is the surface density of the H i gas and $\sigma$ is the component of the velocity dispersion along the line of sight. This formulation of the pressure support correction is the one most commonly employed in the rotation curve literature. It is often called the ‘asymmetric drift’ correction because its formulation is analogous to the familiar correction that applies to (collisionless) stellar discs, although the two corrections have different physical origins (see, e.g., Pineda et al., 2017, for a discussion). This correction is not, strictly speaking, correct, as it assumes a single gas phase and that no bulk flows are present in the disc. Neither assumption holds exactly, of course, but this formula is enough to assess whether pressure forces make an important contribution to the disc kinematics.

We measure the surface density along the (projection of) each of the best-fitting rings directly from the synthetic data cubes. In practice, we measure the gradient of the ‘pressure’ profile $\Sigma_{\rm HI}\sigma^{2}$ using the following fitting function ($\alpha$, $(\Sigma_{\rm HI}\sigma^{2})_{0}$ and $R_{0}$ are free parameters):

[TABLE]

This is the same functional form used in recent analyses of the THINGS and LITTLE THINGS galaxies777de Blok et al. (2008) make no mention of pressure support corrections in their analysis of THINGS galaxies, though for the majority of the galaxies in their sample the correction would be expected to be very small. (Oh et al., 2011, 2015; Iorio et al., 2017).

5 Results

5.1 Gas rotation velocities

Before discussing the application of the tilted-ring model described in the previous section to APOSTLE galaxies, we begin by comparing the mean azimuthal speed of the gas in the disc plane, $V_{\rm rot}(r)$, with the true circular velocity of the system, $V_{\rm circ}(r)$. The purpose of this exercise is to weed out cases where the gas is patently out of equilibrium, since our main goal is to examine the possible shortcomings of the tilted-ring model for galaxies where the disc is close to equilibrium. This is, very roughly, analogous to the common practice of omitting galaxies with obvious kinematic irregularities (mergers, strong bars or tidal features) from rotation curve observing campaigns or analyses. We stress that the ‘equilibrium’ criterion used here is indicative only, and cannot be replicated in observed galaxies, where the true circular velocity profile is unknown. The distinction between equilibrium and non-equilibrium galaxies is only adopted in order to simplify the interpretation of our analysis, and not to compare with observations. In particular, we note that several of the galaxies which we discard as out-of-equilibrium would very likely be included in observational samples of relaxed galaxies.

The rotation profile was measured using the H i mass-weighted mean azimuthal velocity of gas particles in a series of $2\,{\rm kpc}$ thick, $500\,{\rm pc}$ wide cylindrical shells aligned along the disc plane. The velocity dispersion profile was measured using the same series of rings. The 1D line-of-sight gas velocity dispersion, $\sigma$, results from the contribution from the (isotropic) thermal pressure plus that of the ‘bulk’ motion of the gas; i.e.

[TABLE]

where $k_{\rm B}$ is Boltzmann’s constant, $T$ is the particle temperature, $\mu$ is its mean molecular weight, $m_{\rm p}$ is the proton mass, and $\sigma_{\phi}$, $\sigma_{R}$ and $\sigma_{z}$ are the azimuthal, radial and vertical components of the gas particle velocity dispersion. Both components are reflected in the synthetic data cubes (Sec. 3.3), though in practice the ‘bulk’ component always dominates by a factor $>2$.

We show three examples from APOSTLE in Fig. 4, where each column refers to a different galaxy, and, from top to bottom, each panel shows, respectively, the rotation speed, the 1D velocity dispersion, and the H i surface density profiles. The thick black curve in the top panels denotes $V_{\rm circ}(r)$; the thin grey curve $V_{\rm rot}(r)$; and the thick grey curve the ‘pressure-corrected’ rotation speed, as in Eq. 2. Note that, as anticipated in Sec. 1, the pressure corrections are usually small.

If we focus on the inner (rising) part of the rotation curves shown in Fig. 4, we see that the mean gas rotation speed closely traces the circular velocity in two of the three galaxies. The gas rotation curve of the galaxy in the rightmost column, on the other hand, deviates quite strongly from $V_{\rm circ}(r)$ in the inner regions, indicating that this galaxy has likely undergone a recent perturbation that has pushed the gas component temporarily out of equilibrium. Galaxies like the latter are highlighted in Fig. 3 by a lighter shade of grey and are excluded from the analysis that follows, leaving $15$ galaxies. (The actual criterion adopted is that the pressure-corrected $V_{\rm rot}$ differs from $V_{\rm circ}$ by more than $15$ per cent at a fiducial radius of $2\,{\rm kpc}$.)

5.2 Orientation and tilted-ring rotation curves

Having excluded galaxies where the inner gas disc is clearly out of equilibrium, we proceed to model the remaining galaxies using 3Dbarolo. Although the inclination is fixed at $60^{\circ}$ in all synthetic observations, there is still freedom to choose a second angle to define the line of sight. Fig. 5 shows two 3Dbarolo fits to the ‘equilibrium’ galaxies in Fig. 4, as obtained for two different line-of-sight orientations. These were not chosen at random, but have instead been selected to demonstrate the importance of orientation effects on the rotation curves of seemingly ‘equilibrium’ galaxies in APOSTLE.

The tilted-ring modelling returns rotation curves that at times underestimate significantly the mean azimuthal speed of the gas (see blue curves), and, consequently, its circular velocity. The situation changes when the galaxy is rotated by $90^{\circ}$, keeping the same inclination: in this case (shown in red) the inferred rotation speeds are substantially higher, and at times even exceed the true circular velocity of the system. Note that the difference between the red and blue rotation curves is much greater than the ‘errors’ that the model assigns to the recovered $V_{\rm rot}(r)$; these are shown by the shaded area888Errors shown are as estimated by 3Dbarolo: the model parameters are resampled around the best-fitting values to determine the variations required to change the model residual by 5 per cent. This yields an error similar to what might be derived from differences between the approaching and receding sides of the galaxy (Di Teodoro & Fraternali, 2015). around the rotation curves. (Shaded areas are only shown for a couple of curves for clarity.)

The differences in the recovered rotation curves cannot be ascribed to variations in the inclination (a difference of $10^{\circ}$ at $i=60^{\circ}$ only changes $V_{\rm rot}$ by $\sim 10$ per cent), or in the velocity dispersion (differences of \mathrel{\raise 1.50696pt\hbox{\scriptstyle<}\kern-6.00006pt\lower 1.72218pt\hbox{{\scriptstyle\sim}}}4\,{\rm km}\,{\rm s}^{-1}), or in the pressure correction inferred by the model, as may be seen from the other panels in Fig. 5. The model actually recovers these parameters quite well, which implies that the orientation dependence must be due to the presence of large-scale, coherent non-circular motions in the plane of the disc.

5.3 Non-circular motions and orientation effects

The presence of large-scale non-circular motions is illustrated in Fig. 6, where the top panels show, for each galaxy, the residual azimuthal motions in the discs after subtracting the mean rotation at each radius, where the mean is measured directly from the simulation particle information. The line of nodes (i.e., the projected major axis) of the two projection axes shown in Fig. 5 are illustrated by the lines, with corresponding colours. Note the presence of a clear radially coherent bisymmetric pattern in the residual velocities, which explains the results obtained by 3Dbarolo in projection. The bisymmetric perturbation (which resembles that of a slowly-rotating bar-like pattern) is caused by the triaxial nature of the dark matter halo that hosts the galaxy (Hayashi et al., 2004), as we discuss in detail in a companion paper (Marasco et al., 2018).

When the projected major axis slices through the minima of the pattern (blue lines) the recovered rotation velocities underestimate the true rotation speed; the opposite happens when the major axis slices through the two maxima of the residual map (red lines). This is because, in projection, most of the information about the rotation velocity is contained in sight lines near the major axis – gas rotating faster or slower than the average on the projected major axis drives the rotation curve up or down, respectively.

The bottom panels of Fig. 6 further illustrate the non-circular motion pattern. The points indicate the rotation speed as a function of azimuthal angle at a radius $R=5\,{\rm kpc}$ (innermost grey ring in the upper panels). We fit the $m=0$ and $m=2$ terms of a Fourier series:

[TABLE]

with amplitudes $V_{m}$ and phases $\phi_{m}$, to these points, and plot the two terms separately with dashed line styles. In both cases there is a strong $m=2$ component. The maxima of this mode align with the projection axis drawn in red in the upper panels (red vertical lines in lower panels); the minima align with the direction drawn in blue.

In principle there may also be an $m=1$ component in the non-circular motions; however, at any given radius its amplitude is degenerate with the assumed systemic velocity. Here we have adopted a systemic velocity that minimizes the $m=1$ term in the harmonic expansion at $5\,{\rm kpc}$ in order to focus on the bisymmetric component.

Fig. 7 confirms unambiguously the effect of this $m=2$ pattern on the rotation curve recovered by 3Dbarolo. Here we show the rotation speed recovered by the tilted-ring model at two different radii ($R=2$ and $10\,{\rm kpc}$) as a function of the orientation of the line of sight (keeping the inclination always fixed at $i=60^{\circ}$). Clearly, as the orientation varies the inferred rotation speed varies as expected from a dominant $m=2$ mode (i.e., two maxima and two minima as the galaxy is spun by $360^{\circ}$). Note that the phase of the modulation varies between the two radii, slightly in the case of AP-L1-V1-8-0 and more strongly for AP-L1-V4-8-0. This indicates that the phase of the $m=2$ mode shifts gradually with radius, as may be corroborated by visual inspection of the residual maps in Fig. 6.

5.4 Non-circular motions in projection

Harmonic modulations of the velocity field (of order $m$) are not always easily discernible in projection, where they are mapped into a combination of $m\pm 1$ modes. The analogue in projection along the line of sight of Eq. 5 is:

[TABLE]

where, as usual, $i$ is the inclination angle and we have assumed that the position angle of the major axis of the projected circle (ellipse) is at $\phi=0$. This may also999$V_{\rm LoS}(\phi)=\sin i\cos\phi\sum_{m}a_{m}\sin(m\phi)+b_{m}\cos(m\phi)$, as used by, e.g., Schoenmakers et al. (1997), is also equivalent. be written as:

[TABLE]

For example, when projecting an inclined circle of radius $R$ with average azimuthal velocity $V_{0}$, perturbed by an $m=2$ pattern with amplitude $V_{2}$ and phase $\phi_{2}$, the line-of-sight velocity along the resulting ellipse will be

[TABLE]

where we have assumed that the position angle of the major axis of the ellipse is at $\phi=0$.

When the maxima of the mode are aligned with the major axis of the projection (i.e., $\phi_{2}=0^{\circ}$ or $180^{\circ}$), then the modulation increases the inferred rotation velocity. When $\phi_{2}=90^{\circ}$ or $270^{\circ}$, on the other hand, the projected kinematic major axis lines up with the minima and the inferred rotation velocities decrease.

More generally, the velocity fluctuation about what is expected from uniform circular motion may be expressed by subtracting $V_{0}\sin(i)\cos(\phi)$ from Eq. 8,

[TABLE]

which may also be written as

[TABLE]

In other words, in projection, a bisymmetric perturbation would be seen as simultaneous $m=1$ and $m=3$ perturbation to the line-of-sight velocities. The more visually obvious of the two is the three-peaked $m=3$ component. Such a three-peaked modulation may be quite difficult to detect in residual maps from tilted-ring models, as we show in Fig. 8. Here we show the same projections of the two galaxies from Fig. 5. The left column shows the line-of-sight velocity map, the middle column the $1^{\rm st}$ moment of the 3Dbarolo map, and the residuals from the difference between these two are shown in the rightmost column.

The expected three-peaked pattern in the residuals is not readily apparent in any of the four cases illustrated. There are a number of reasons for this. First, the amplitude of the pattern is not very large ($\sim 15\,{\rm km}\,{\rm s}^{-1}$ on average according to Fig. 6) and therefore only comparable to the rotation speed near the centre. Second, other residuals, not necessarily caused by the bisymmetric mode, dominate in the outer regions, obscuring the effect. Third, there are other harmonic modes in the velocity field which hinder a straightforward interpretation of the line-of-sight velocity field. For instance, there is a $m=2$ symmetric modulation of the radial velocities, as expected for gas orbits in a bar-like potential (e.g. Spekkens & Sellwood, 2007), which tends to partially cancel the projected signature of the $m=2$ term in the azimuthal harmonic expansion. Finally, the tilted ring model attempts to provide a ‘best fit’ by varying all available parameters so as to minimize the residuals. Given the number of parameters available (each ring has, in principle, independent velocity, inclination, dispersion, and position angle), the resulting residuals are quite small, masking the expected three-peaked pattern, except perhaps in the most obvious cases.

We have attempted to measure the amplitudes and phases of the $m=2$ azimuthal harmonic perturbations seen in Fig. 6 from the line-of-sight velocity fields. However, all harmonic terms contribute to the line of sight velocities and must therefore in principle be modeled. Even a harmonic expansion fit up to only order $m=2$ is a nine-parameter problem prone to degeneracies: radial and azimuthal amplitudes and phases for $m=1$ and $m=2$, and the $m=0$ amplitude (circular velocity) must be determined, even when assuming (probably incorrectly) that vertical motions in the disc are negligible. In particular, the $m=1$ amplitude is degenerate with a combination of the systemic velocity and centroid; the radial and azimuthal terms of the same order are also degenerate given freedom in the phase; and the $m=2$ amplitude is partially degenerate with the inclination.

Breaking these degeneracies requires strong assumptions, e.g. regarding the relative phases of the various terms in the expansion, or that the gas orbits form closed loops, which are difficult to justify based on ’observable’ information. Ultimately, we find that we are unable to accurately recover the harmonic modes present in the discs of APOSTLE galaxies from their line of sight velocity maps without recourse to information which would be unavailable for real galaxies.

In spite of this complication, the azimuthal $m=2$ term clearly dominates in the case of APOSTLE galaxies, and is largely responsible for setting the recovered rotation velocities. Although the $m=1$ term in the azimuthal expansion sometimes has an amplitude comparable to the $m=2$ term, its degeneracy with the systemic velocity implies that in practice the effect on the recovered rotation curve is dominated by the $m=2$ harmonic. This is seen in Fig. 7, where the inferred rotation velocity fluctuates with orientation angle as expected from an $m=2$ modulation, with the same phase as that shown in Fig. 6. We have verified this empirically by using 3Dbarolo to fit simple analytic models of differentially rotating discs with harmonic perturbations to their velocity fields. For independent $m=1$ and $m=2$ perturbations of the same amplitude, the $m=2$ symmetric perturbation always has a much stronger effect on the recovered rotation velocities, as measured by the average amplitude of their variation with phase angle.

5.5 Non-circular motions and the inner mass deficit problem

The discussion above shows that non-circular motions can substantially affect the rotation curves of APOSTLE galaxies inferred from tilted-ring models. Because of their ubiquity, these motions can affect the inferred inner matter content of APOSTLE galaxies in a manner relevant to the ‘inner mass deficits’ discussed in Sec. 1.

Following Oman et al. (2015), we estimate the inner mass deficit (a proxy for the importance of a putative ‘core’) of a galaxy by comparing the observed and simulated relations between $V_{\rm circ}(2\,{\rm kpc})$ and the maximum circular velocity of the system, $V_{\rm max}$. The latter can usually be accurately determined because it is reached in the outer regions, which are less affected101010The overall inclination of the system may in some cases be the ultimate impediment for accurate estimates of the rotation speed, especially when derived from kinematics alone (see, e.g., Oman et al., 2016, for a discussion of a few examples). by observational uncertainties.

A lower limit for $V_{\rm circ}(2\,{\rm kpc})$ can be derived solely from the dark matter enclosed within $2\,{\rm kpc}$ for a $\Lambda$CDM halo of given $V_{\rm max}$, and is indicated by the grey band in Fig. 9. In the same figure, the red band indicates the $V_{\rm circ}(2\,{\rm kpc})$–$V_{\rm max}$ relation for all galaxies in the APOSTLE and EAGLE suites of cosmological hydrodynamical simulations, including scatter. At low circular velocities the red and grey bands overlap, indicating that most low-mass EAGLE and APOSTLE galaxies are dark matter dominated111111Note that at the very low velocity end, V_{\rm max}\mathrel{\raise 1.50696pt\hbox{\scriptstyle<}\kern-6.00006pt\lower 1.72218pt\hbox{{\scriptstyle\sim}}}30\,{\rm km}\,{\rm s}^{-1}, the maximum circular velocity is reached at radii close to $2$ kpc, so that $V_{\rm max}\approx V_{\rm circ}$(2 kpc).. Note the small scatter in the $V_{\rm max}$–$V_{\rm circ}(2\,{\rm kpc})$ relation expected from these simulations.

The red circles indicate the ‘equilibrium’ APOSTLE galaxies selected for the present study, where $V_{\rm circ}(2\,{\rm kpc})$ and $V_{\rm max}$ are estimated directly from the mass profile as derived from the particle data. Green circles in Fig. 9 indicate the average azimuthal velocities at $2\,{\rm kpc}$ of the gas for the same APOSTLE sample. The good agreement between red and green circles is not surprising and simply reflects our definition of ‘equilibrium’ as cases where the difference between average rotation speed and circular velocity at $2\,{\rm kpc}$ is smaller than $15$ per cent.

Blue symbols indicate the results for THINGS and LITTLE THINGS galaxies collated directly from the literature (de Blok et al., 2008; Oh et al., 2011, 2015). Note that a number of galaxies lie below the red band; these are galaxies with an apparent ‘inner deficit’ of mass at that radius (i.e., ‘cores’) compared with the cold dark matter prediction. The broken lines are labelled by an estimate of this mass deficit, expressed in solar masses. Galaxies like DDO 47 (blue square ‘4’) and DDO 87 (blue square ‘9’) lie well below the expected relation; they are clear examples of galaxies with the slowly rising rotation curves traditionally associated with ‘cores’ in the dark matter (Fig. 18).

The black circles in Fig. 9, finally, indicate rotation speeds at $2\,{\rm kpc}$ for all ‘equilibrium’ APOSTLE galaxies, derived using 3Dbarolo for a fixed inclination ($i=60^{\circ}$) and a single random orientation per galaxy. The obvious difference between black and green circles highlights two important conclusions. One is that non-circular motions in the inner regions are important enough to produce, at times, deviations from the expected relation as large as measured in observed galaxies. The second is that $V_{\rm rot}$ generally underestimates $V_{\rm circ}$ at $2\,{\rm kpc}$. Overestimates also occur in some cases, but these are rare and usually milder.

This indicates that the discrepancy between inner rotation and circular speeds is not solely a result of a bisymmetric modulation of the velocity field, where overestimates should occur as frequently as underestimates. The systematic underestimate of the circular velocity must arise from other effects, such as (i) the non-negligible thickness of the gas disc (which causes gas at different radii and heights to fall along the line of sight; see below); (ii) morphological irregularities that may push the gas temporarily out of equilibrium (e.g. H i bubbles, see also Read et al., 2016; Verbeke et al., 2017); and (iii) underestimated ‘pressure’ support from random motions in the gas (Pineda et al., 2017). Our main conclusion is that tilted-ring modelling of APOSTLE galaxies results in a diversity of inner rotation curve shapes and apparent inner mass deficits that are comparable to those of observed galaxies, mainly due to non-circular motions in the gas.

5.6 Disc thickness and projection effects

Of the three possible explanations for the systematic underestimate of the circular velocity at $2\,{\rm kpc}$, the thickness of the gaseous discs seems to be the leading factor in APOSTLE galaxies. We show this in Fig. 10, using as an example AP-L1-V1-8-0, whose 3Dbarolo rotation curves for two orthogonal orientations are shown in the top-left panel of Fig. 5. As discussed above (Sec. 5.3), the blue curve ($\Phi=255^{\circ}$)121212We show here the orientation offset $180^{\circ}$ from that in Fig. 5 – the effect of the bisymmetric perturbation is much the same, but by chance this orientation illustrates the effect of the thick disc more clearly than the $\Phi=75^{\circ}$ orientation. subtantially underestimates $V_{\rm circ}(r)$ because the kinematic major axis of the projection coincides with the minima of the $m=2$ perturbation pattern shown in Fig. 6. Why doesn’t then the red curve ($\Phi=165^{\circ}$), where the kinematic major axis traces the maxima of the $m=2$ pattern, overestimate $V_{circ}$ by a similar amount?

The answer may be gleaned from Fig. 10, where we show ‘position-velocity’ (PV) diagrams (i.e., the distribution of line-of-sight velocities along the kinematic major axis) for this galaxy. Red and blue curves are the rotation velocites estimated by 3Dbarolo; black curves indicate the true circular velocity on the plane of the disc. The left panels show the PV diagrams including all the H i gas in the galaxy; the panels on the right, on the other hand, only include gas very close to the disc plane (i.e., $|z|<200\,{\rm pc}$)131313This height was chosen as the minimum height for which a synthetic observation could be constructed without becoming unduly dominated by shot noise in the simulation; the precise value is not otherwise significant.. The right-hand panels exclude extra-planar gas, which tends to rotate more slowly and to lower the average speed at given $R$. When only gas near the plane is considered, the inferred rotation velocities do indeed under and overestimate the true circular velocity by similar amounts, as expected for a bisymmetric perturbation. The effect of the extraplanar gas is to bring down average speeds in both cases, effectively cancelling the overestimate in the case of $\Phi=165^{\circ}$, and leading to systematic underestimation of the inner circular velocities when averaged over all orientations.

This is shown explicitly in Fig. 11 (see also inset panels in Fig. 10), where we plot the distribution of 3Dbarolo inferred rotation velocities at $2\,{\rm kpc}$ for $72$ different orientations of the galaxy, each separated by $5^{\circ}$ in azimuth. The filled histogram correspond to all H i gas, the open histogram only to gas with $|z|<200\,{\rm pc}$. The combined effects of non-circular motions and extraplanar gas lead to a relatively large dispersion in the estimated values (of order $8.3\,{\rm km}\,{\rm s}^{-1}$ and $7.7\,{\rm km}\,{\rm s}^{-1}$ for the open and filled histograms, respectively), as well as a substantial shift in the median value; from $52.8\,{\rm km}\,{\rm s}^{-1}$ to $46.7\,{\rm km}\,{\rm s}^{-1}$ (for reference, the circular velocity at this radius is $52.0\,{\rm km}\,{\rm s}^{-1}$).

We end this discussion by presenting, in Fig. 12, the thickness of APOSTLE H i discs. Here we show, as a function of cylindrical radius, the height ${\rm HWHM}_{z}$ containing half the H i mass. Black and grey thin lines indicate our sample of ‘equilibrium’ and ‘non-equilibrium’ galaxies, respectively. The thick red line and shaded areas highlight the median and interquartile half-mass heights for the ‘equilibrium’ sample. Note that ‘equilibrium’ galaxies are slightly thinner, on average; their typical half-HI height is $\sim 500$ pc at $R=2\,{\rm kpc}$. APOSTLE discs also flare outwards (Benítez-Llambay et al., 2018); the typical H i half-height climbs to $\sim 1\,{\rm kpc}$ at $R=5\,{\rm kpc}$.

6 Inner mass deficits of observed galaxies

We explore next whether non-circular motions and thick, differentially rotating discs seen in projection might be at least partially responsible for the observed ‘inner mass deficits’ seen in Fig. 9. This is a lengthy topic that requires a detailed individual study of each galaxy, which we plan to pursue in future work. Our main purpose here is to provide a ‘proof of concept’ that at least some of the galaxies with large inferred deficits/cores might indeed have rotation curves that are unduly affected by model shortcomings similar to those seen in the analysis of APOSTLE galaxies.

6.1 Non-circular motions

We motivate the potential applicability of our findings from Sec. 5.3 to real galaxies by considering two example systems. DDO 47 and DDO 87 have reported inner rotation speeds which are well below those expected for CDM haloes (see blue squares labelled ‘4’ and ‘9’ in Fig. 9). These galaxies have been modelled with 3Dbarolo by Iorio et al. (2017), and we adopt here the same parameters as in that study. The 3Dbarolo rotation curves are consistent (within the reported errors) with those of Oh et al. (2015), despite significant differences in the algorithms implemented in the codes used in each case.

The $1^{\rm st}$ moment maps of DDO 47 and DDO 87 are shown in the left column of Fig. 13. The 3Dbarolo velocity fields are shown in the middle, and the residuals in the panels on the right. Although, as discussed in Sec. 5.4, we are unable to accurately constrain the parameters of a harmonic expansion from the line of sight velocity maps, the clear three-peaked (‘clover leaf’) pattern in the residuals (with amplitudes of order $\sim 10\,{\rm km}\,{\rm s}^{-1}$; about half the inferred rotation velocity at $2\,{\rm kpc}$) strongly suggests the presence of bisymmetric kinematic modulations near the centres of these galaxies. The phase of the residual patterns – with maxima (red) lying approximately along the major axis on the approaching side of the disc (blue in the left column; see Fig. 19 for a schematic) – further suggests that the recovered rotation velocities might very well underestimate the true inner circular velocity. Although a definitive verdict on whether these galaxies are actually consistent with cuspy CDM haloes must await a more detailed analysis, we regard the evidence for non-circular motions shown in the right-hand panels of Fig. 13 as strong enough to call into question the conclusion that large ‘cores’ are present in these galaxies.

6.2 Disc thickness and projection effects

We argued in Sec. 5.6 that the finite thickness of APOSTLE discs is mainly responsible for the systematic underestimation of inner circular velocities shown in Fig. 9. Could this explanation also hold for observed galaxies?

We begin by noting that assessing the impact of disc thickness on the inferred rotation curve is challenging, since it depends not only on the gas scale height and its radial dependence, but also on the gas velocity gradient in both the radial and vertical direction, none of which are entirely straightforward to constrain in real galaxies. We defer a more detailed analysis to a future study, and will only attempt here a preliminary comparison between the thickness of observed and APOSTLE discs to motivate the fact that our conclusion in Sec. 5.6 may indeed be applicable to real galaxies.

A simple (but statistical) measure of disc thickness is provided by the distribution of apparent axis ratios (i.e., projected shapes) of H i discs. This has been reported for FIGGS141414Faint Irregulars Galaxies GMRT Survey. (Begum et al., 2008) galaxies by Roychowdhury et al. (2010), whose results we reproduce in Fig. 14 (shaded histogram). The symbols with error bars in the same figure correspond to a sample of APOSTLE discs chosen to match the median and width of the H i mass distribution of FIGGS galaxies. (Note that this is not the same sample we have used for the kinematic analysis in previous sections.) The two distributions are in very good agreement, suggesting that FIGGS and APOSTLE discs have similar thicknesses, at least according to this statistic.

Direct estimates of the half-mass thickness of H i discs are available only for small samples of galaxies, and are subject to selection biases and observational uncertainties that hinder a proper comparison with the APOSTLE results presented in Fig. 12. For example, edge-on galaxies are good targets in principle, but they are typically identified from samples of thin galaxies, thereby possibly incurring selection biases difficult to quantify and correct for.

O’Brien et al. (2010) analyze a sample of 8 H i-rich, approximately edge-on late-type galaxies and report ${\rm HWHM}_{z}$ in the range $0.3$–$1.0\,{\rm kpc}$ at a radius of $5\,{\rm kpc}$ (see their fig. 24; note that they show ${\rm FWHM}_{z}=2\times{\rm HWHM}_{z}$). Peters et al. (2017) re-analyzed a subset of the same galaxies and also report estimates of disc thickness in the central regions that are consistent with those from O’Brien et al. (2010). These scale heights are on average slightly thinner than the APOSTLE equilibrium sample, where we find a half-mass H i height of $1.0\pm 0.2$ at $5\,{\rm kpc}$, but this is perhaps not unexpected given the caveat expressed in the preceding paragraph.

Another set of estimates of disc thickness for dwarf galaxies comparable to our APOSTLE sample (DDO 154, Ho II, IC 2574 and NGC 2366), is provided by Banerjee et al. (2011). These authors use a dynamical model that requires assumptions about geometry and about the contribution of stars, gas, and dark matter to the gravitational potential. They report ${\rm HWHM}_{z}(5\,{\rm kpc})$ in the range $0.2$–$0.8\,{\rm kpc}$, or a factor of $\sim 2$ thinner than the ‘equilibrium’ APOSTLE galaxies at the same radius.

Finally, an estimate of the importance of extra-planar gas can also be made by integrating the amount of ‘kinematically anomalous’ gas in a galaxy, which refers to the ‘non-Gaussian’ tail of the H i velocity distribution at each pixel (Fraternali et al., 2002; Sancisi et al., 2008). When applied to our projected APOSTLE velocity fields, we find that between $10$ and $20$ per cent of the H i flux of a galaxy may be classified as ‘kinemaically anomalous’. This is actually in good agreement with the results of the same procedure applied to LITTLE THINGS galaxies, although a quantitative comparison is not straightforward because it depends on observational issues such as the treatment of noise, smoothing, masking, etc., which are difficult to replicate in simulations. Interestingly, we find that, in APOSTLE, a fair fraction of the ‘kinematically anomalous’ gas is actually not extra-planar, but, rather, gas close to the mid-plane that has been disturbed by recent episodes of star formation, or other kinematic perturbations. Again, a detailed comparison between observation and simulation requires a more meticulous study than the preliminary exploration we attempt here, but the overall agreement between observation and simulation for this measure is reassuring.

In conclusion, and taking these various measurements at face value, it appears that our sample of APOSTLE H i discs are typically somewhat thicker at R\mathrel{\raise 1.50696pt\hbox{\scriptstyle<}\kern-6.00006pt\lower 1.72218pt\hbox{{\scriptstyle\sim}}}5\,{\rm kpc} than their observed counterparts. Whether this slight offset is enough to invalidate our conclusion that the finite thickness of gas may be responsible for the systematic ‘inner mass deficits’ shown in Fig. 9 is unclear. Reaching a more definitive conclusion demands a more detailed analysis than warranted by the scope of the present paper.

6.3 Non-circular motions, surface brightness, and rotation curve slopes

Finally, we consider whether non-circular motions are consistent with the existence of relatively tight correlations betwen the shape of the rotation cuves and structural properties of the galaxy, such as the inner surface density of gas and stars. Indeed, it has long been appreciated that high-surface brightness galaxies typically have faster rising inner rotation curves than lower surface brightness systems (de Blok et al., 1996). Could these correlations exist if the inner parts of rotation curves are substantially affected by non-circular motions?

We address this question by discussing the correlation between the central rotation curve slope, which we parametrize using the ratio between the radius at which the circular velocity curve reaches $90$ per cent of its maximum value and the velocity at that radius151515This measure has the advantage of being independent of galaxy size; APOSTLE dwarfs have systematically larger H i (Fig. 1) and stellar (Campbell et al., 2017) sizes than observed., $V(R_{90})/R_{90}$, and the central surface density of the stellar component, $\Sigma_{\star}(0)$.

We show the observed correlation based on measurements reported in the recent compilation of late type galaxies by Lelli et al. (2013, note that these authors use a different definition of central rotation curve slope) and compare it with APOSTLE results in Fig. 15. We make two measurements of $V(R_{90})/R_{90}$ for APOSTLE galaxies; one using the true circular velocities of the galaxies (red points), and a second one using the rotation curves recovered using 3Dbarolo (black points), for the same random orientations as used in Fig. 9. $\Sigma_{\star}(0)$ is measured, as in the observations, by extrapolating a fit to the exponential portion of the radial stellar mass profile.

Fig. 15 shows good agreement between observations and APOSTLE. For the surface brightness range sampled by our sample of APOSTLE galaxies rotation curves rise in general fairly slowly, with typical slopes of order $\sim 10\,{\rm km}\,{\rm s}^{-1}\,{\rm kpc}^{-1}$. This is quite well reproduced in the APOSTLE sample, regardless of whether the true circular velocity curves or the 3Dbarolo rotation curves are used. Fig. 15 suggests that there is no conflict between the slope-surface brightness correlations and our results, at least in the range explored by the APOSTLE sample.

7 Summary and Conclusions

We analyse synthetic H i observations of $33$ galaxies from the APOSTLE suite of $\Lambda$CDM cosmological hydrodynamical simulations of galaxy formation. Our sample includes all systems with maximum circular velocities in the range $60$–$120\,{\rm km}\,{\rm s}^{-1}$, simulated at the highest (L1) resolution available in APOSTLE. The H i properties of these galaxies follow closely the scaling laws relating the size, mass and kinematics of observed galaxy discs.

We derive H i rotation curves using the same kind of tilted-ring modelling usually applied to interferometric H i observations. We find that in many cases the rotation curves inferred by this modelling depend sensitively on the orientation of the chosen projection. Particularly affected are the rotation speeds in the inner regions of individual galaxies, which fluctuate systematically as the orientation angle of the line of sight is varied, at fixed inclination. The fluctuations are due to non-circular motions of the gas velocity field in the disc plane. Bisymmetric ($m=2$) variations in the azimuthal velocity of the gas have the most obvious effect, but other asymmetric patterns are also present.

Maximal deviations from the actual circular velocity are obtained when the line of nodes of the projection coincides with the principal axes of the $m=2$ pattern. In particular, rotation speeds can severely underestimate the circular velocity when the line of nodes is aligned with the minima of the modulation pattern. Conversely, alignment of the line of nodes with the maxima can result in a substantial overestimate of the rotation speed. The case of an underestimate is further exacerbated by lower velocity extra-planar gas along the line of sight, whereas in the case of an overestimate the two effects partially compensate each other. Those cases which result in a slowly-rising rotation curves would be erroneously interpreted as signaling a severe ‘inner mass deficit’ compared with the $\Lambda$CDM predictions.

Slowly-rising rotation curves are often interpreted as evidence for a ‘core’ (rather than a cusp) in the dark matter central density profile. This raises the possibility that some galaxies with ‘cores’, and especially those where the inner mass deficit is extremely large, are galaxies whose inner rotation curves have been unduly affected by modelling uncertainties. In fact, all APOSTLE galaxies have central dark matter cusps that have been minimally affected, if at all, by the baryonic assembly of the galaxy.

Bisymmetric patterns map, in projection, into one- and three-peaked patterns in the residual velocity fields, but these are often difficult to detect in tilted ring model residuals, as they may be effectively masked by the numerous free parameters of model fits and by degeneracies which arise in projected coordinates. They are easily appreciated, however, in some galaxies that show extreme deviations from the steeply-rising rotation curves expected from the cuspy CDM halo density profiles, such as DDO 47 and DDO 87. The residual pattern in these galaxies suggests the presence of bisymmetric gas flows strong enough to substantially affect the kinematic modelling. Not only is the true average azimuthal speed of the gas, $V_{\rm rot}(R)$, very difficult to extract in this case, but even if that were possible, it may still differ substantially from the true circular velocity of the system. Galaxies such as these may well have a cuspy dark matter halo despite the apparent slow rise of their inner rotation curves.

Our results thus suggest a possible resolution of the rotation curve diversity problem highlighted by Oman et al. (2015). The problem could, at least in part, result from the failure of tilted-ring models applied to discs of finite thickness with non-negligible non-circular motions. In this context, the systematic underestimation of inner circular velocities may just reflect the finite thickness of gas discs.

It is therefore important to design simple and model-independent diagnostics of non-circular motions to identify the regions where tilted-ring model results are suspect and where uncertainties from such techniques are almost certainly underestimated. Non-circular motions should not be ‘assumed’ to be unimportant, as is often done, but demonstrated to be so through careful analysis, keeping in mind that they may appear weaker than they are. Tighter constraints on the scale heights of the H i discs of dwarfs would also be particularly useful.

Much remains to be done to explore fully these ideas. One issue is the origin of the bisymmetric perturbations, which seem ubiquitous in APOSTLE, as well as their amplitude and radial dependence: are they radially-coherent modes of constant phase angle, as in a bar or in a triaxial halo potential, or of varying pitch angle, as in two-armed spirals? (See Marasco et al., 2018, for a first attempt.) Are they associated with disc instabilities, halo triaxiality, accretion events, or perhaps even triggered and sustained by the stochastic nature of star formation and feedback? Are non-circular motions in observed dwarfs as strong as they seem to be in our simulated galaxies? What are the implications for $\Lambda$CDM models?

Another issue that remains pending is a case-by-case demonstration that non-circular motions provide a viable explanation for galaxies with alleged ‘cores’. Although difficult to detect, we are hopeful that our models will provide enough guidance to devise effective methods for finding and characterizing them. Finally, what are the statistics of these modulations? What fraction of galaxies do they affect, and by how much? Is this consistent with available data? For simplicity, we have fixed the inclination of all galaxies in our sample to $60^{\circ}$; what is the result of adopting a more realistic distribution of inclination angles?

These are all important considerations that we plan to address in future contributions, where we hope to clarify whether the inner mass deficit problem, and the related cusp-core problem, results in whole or in part from an underestimate of the impact of non-circular motions on the kinematic modelling of disc galaxies, or else is a true predicament for the otherwise extremely successful $\Lambda$CDM paradigm.

Acknowledgements

KO and AM thank E. di Teodoro for technical assistance with 3Dbarolo and useful discussions. We thank S.-H. Oh and E. de Blok for data contributions, and the THINGS, LITTLE THINGS and SPARC survey teams for making their data publicly available. We thank the anonymous referee whose comments helped to improve substantially the manuscript. This work was supported by the Science and Technology Facilities Council (grant number ST/F001166/1). CSF and ABL acknowledge support by the Science and Technology Facilities Council (grant number ST/L00075X/1). JS acknowledges ERC grant agreement 278594-GasAroundGalaxies and the Netherlands Organisation for Scientific Research (NWO) through VICI grant 639.043.409. This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grant ST/H008519/1, and STFC DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure. This research has made use of NASA’s Astrophysics Data System.

Appendix A Properties of simulated and observed galaxies

In Tables A & A we collect some key properties of our sample of simulated galaxies and the THINGS (Walter et al., 2008) and LITTLE THINGS (Hunter et al., 2012) galaxies, respectively, including the key to the numeric labels used in Figs. 1, 3 & 9.

Not all galaxies in these surveys are amenable to kinematic analysis. We select the THINGS galaxies analysed in de Blok et al. (2008) and Oh et al. (2011). In cases where the same galaxy is analysed in both studies, we use the more recent analysis. We exclude some of the galaxies:

• •

NGC 4826 has two counter-rotating discs, making the interpretation of the rotation curve as a circular velocity curve highly uncertain.

• •

Ho II is the same galaxy as DDO 50, which appears in the LITTLE THINGS sample.

• •

We omit DDO 53, DDO 154 and NGC 2366, which are included as part of the LITTLE THINGS sample.

• •

Ho I has a particularly low inclination of \mathrel{\raise 1.50696pt\hbox{\scriptstyle<}\kern-6.00006pt\lower 1.72218pt\hbox{{\scriptstyle\sim}}}14^{\circ} (Oh et al., 2011), which makes the necessary inclination correction unacceptably large and uncertain.

We select all of the LITTLE THINGS galaxies analysed in Oh et al. (2015).