How to Reconcile the Observed Velocity Function of Galaxies with Theory

TL;DR

This study uses high-resolution cosmological simulations and mock HI data to reconcile the discrepancy between observed dwarf galaxy velocities and the predictions of the LCDM model, highlighting observational biases and dark matter profile effects.

Contribution

It demonstrates that accounting for observational methods and dark matter profiles resolves the velocity function discrepancy in LCDM.

Findings

Observed velocity functions align with LCDM when using mock observations.

Underestimation of halo velocities by HI linewidth measurements explains part of the discrepancy.

Faint, low-mass halos are often undetected in current surveys.

Abstract

Within a Lambda Cold Dark Matter (LCDM) scenario, we use high resolution cosmological simulations spanning over four orders of magnitude in galaxy mass to understand the deficit of dwarf galaxies in observed velocity functions. We measure velocities in as similar a way as possible to observations, including generating mock HI data cubes for our simulated galaxies. We demonstrate that this apples-to-apples comparison yields an "observed" velocity function in agreement with observations, reconciling the large number of low-mass halos expected in a LCDM cosmological model with the low number of observed dwarfs at a given velocity. We then explore the source of the discrepancy between observations and theory, and conclude that the dearth of observed dwarf galaxies is primarily explained by two effects. The first effect is that galactic rotational velocities derived from the HI linewidth severely underestimate the maximum halo velocity. The second effect is that a large fraction of halos at the lowest masses are too faint to be detected by current galaxy surveys. We find that cored dark matter density profiles can contribute to the lower observed velocity of galaxies, but only for galaxies in which the velocity is measured interior to the size of the core (~3 kpc).

Figures

Tables

Table 1. Table 1 Properties of the Simulated Galaxies

Simulation	$v_{max,dmo}$ range	Mstar range	mDM,part	mstar,part	Softening	Overdensity	N${}_{\text{DM}}$
	km s-1	$M_{\odot }$	$M_{\odot }$	$M_{\odot }$	pc	$\Delta \rho /\rho$	within Rvir
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
Fields 1 – 6	30-150	2$\times$105-1010	1.6$\times$105	8 $\times$103	174	-0.15 to 1.35	0.03-3.4$\times$106
Field 7	43-56	2$\times$107-3$\times$108	2$\times$104	103	85	-0.02	0.05-2$\times$106
Field 8	38	108	6$\times$103	4.2$\times$102	64	0.01	2$\times$106

Table 2. Table 2 Characteristic Velocity Definitions in the Text

symbol	definition
$v_{circ}$	circular velocity; $v_{circ}=\sqrt{GM/r}$ where $M$ is the mass enclosed within radius $r$
$v_{max,dmo}$	the maximum value of $v_{circ}$ for a dark matter-only simulated halo; 2$v_{max,dmo}$ is the theoretical counterpart to $w_{50}^e$
$2v_{max,dmo}\sin i$	twice the maximum value of $v_{circ}$ measured for a dark matter-only simulated halos multiplied by the $sin$ of the
	observational inclination angle $i$; the theoretical counterpart to $w_{50}$
$w_{50}^e$	for galaxies with measurable HI: the full width of a galaxy’s HI line profile, measured at 50% of the profile peak
	height when the galaxy is viewed edge-on (inclination $i$ = 90∘); for galaxies with no measurable HI: twice the
	stellar velocity dispersion; the observational counterpart to 2$v_{max,dmo}$
$w_{50}$	$w_{50}^e$$\times$ sin($i$); $w_{50}$ for a galaxy viewed at a random inclination angle $i$; observational counterpart to $2v_{max,dmo}\sin i$
$w_{20}^e$	similar to $w_{50}^e$ but measured at 20% of the HI profile peak height
$V_f$	velocity of a galaxy measured on the flat part of the rotation curve
$v_{max,sph}$	the maximum value of $v_{circ}$ for a galaxy halo in a baryonic simulation
$v_{out}$	$v_{circ}$ measured at Rout, the radius at which a galaxy’s HI surface density falls below 1 M⊙/pc2
$v_{out,dmo+b}$	$v_{circ}$ for a dark matter-only halo (reduced by a velocity consistent with removing the cosmic baryon fraction) + $v_{circ}$
	for only the baryons in the counterpart baryonic simulation. Measured at Rout, where Rout is determined from the
	simulated baryonic counterpart

Equations

$${\varphi }_h(v_{max,dmo})=\frac{d{N}_h}{dVd{\log }_{10}(v_{max,dmo})}.$$

$${\varphi }_{h,det}=f_{det}(v_{max,dmo})\times {\varphi }_h(v_{max,dmo}).$$

$${\varphi }_{h,det}(v_{max,dmo})\to {\varphi }_{h,det}(w_{50}).$$

Full text

††slugcomment: Accepted for publication in ApJ

1. Introduction

The velocity function of galaxies is indicative of the number of galactic halos that exist as a function of mass, and is therefore a powerful test of our cosmological galaxy formation model. For galaxies with velocities above $\sim$100 km s*-1* that are primarily dispersion-dominated, the observed velocity function (VF) is generally in agreement with theoretical expectations within a $\Lambda$ Cold Dark Matter (CDM) cosmology, as long as the effects of baryons are included (Gonzalez et al., 2000; Sheth et al., 2003; Chae, 2010; Obreschkow et al., 2013). To probe to lower galaxy masses which are more likely to be rotation-dominated, HI rotation data is ideal. Zwaan et al. (2010) and Trujillo-Gomez et al. (2011) were some of the first to combine the early type galaxy VF from SDSS with the HIVF from HIPASS to probe to lower galaxy masses, and found that theory predicted more dwarfs below $\sim$80 km s*-1* than observed (see also Abramson et al., 2014; Bekeraitė et al., 2016).

The HIVF has since been updated, thanks in large part to data from the ALFALFA HI survey (Giovanelli et al., 2005), and from systematic optical searches for neighboring galaxies (Karachentsev et al., 2013). Papastergis et al. (2011) used early data from the ALFALFA survey (at 40% of its eventual sample size) to confirm that there is a deficit of low mass observed galaxies compared to that expected in a $\Lambda$CDM cosmology. In galaxies with rotational velocities, $v_{rot}$, $\sim$25 km s*-1*, the ALFALFA HIVF shows nearly an order of magnitude fewer galaxies than expected based on straightforward $\Lambda$CDM estimates (e.g., that each dark matter halo contains one luminous galaxy). Klypin et al. (2015) made a separate measurement of the VF using the catalog of Local Volume galaxies out to 10 Mpc (Karachentsev et al., 2013). They derived velocities for galaxies as faint as $M_{B}=-10$. Despite probing to these low masses and correcting for completeness, they still found a dearth of low velocity galaxies compared to the number expected in CDM. Likely due to the fact that they could include gas-poor faint galaxies, the discrepancy is not as large as seen in the ALFALFA HIVF sample, but they confirm the nearly factor of 10 discrepancy between theory and observation at $v_{rot}\sim 25$ km s*-1*.

This missing dwarf problem is reminiscent of the missing satellites problem (Moore et al., 1999; Klypin et al., 1999), but now extends into the field, beyond the virial radius of more massive galaxies. This means that solutions that rely on the tidal field of the host galaxy to reduce the numbers and masses of satellite dwarfs (Zolotov et al., 2012; Brooks et al., 2013; Arraki et al., 2014; Brooks & Zolotov, 2014; Wetzel et al., 2016) should not apply, and a new mechanism to reduce the number of field galaxies needs to be invoked. One long-standing solution to the missing dwarf problem is warm dark matter (WDM, e.g., Bode et al., 2001; Polisensky & Ricotti, 2011; Menci et al., 2012; Lovell et al., 2012; Nierenberg et al., 2013), in which the thermal relic mass of the dark matter particle is $\gtrsim 2$ keV. However, Klypin et al. (2015) and Brook & Di Cintio (2015b) showed that WDM more massive than 1.5 keV doesn’t suppress enough structure at low masses to be compatible with the observed VF (though see Schneider et al., 2017). Lighter masses have already been ruled out based on the small scale structure observed in the Lyman-$\alpha$ forest at high redshift (Viel et al., 2006; Seljak et al., 2006; Viel et al., 2008, 2013). Hence, WDM is difficult to make compatible with all available observational constraints.

Another interpretation of the observed VF is not that there are dwarfs missing, but that low mass galaxies display lower velocities than anticipated. This may be due to complications related to the way rotational velocities are measured observationally (Brook & Shankar, 2016; Macciò et al., 2016; Yaryura et al., 2016), baryonic physics (e.g., Brook & Di Cintio, 2015a, b), or to dark matter physics if the dark matter is self-interacting (Spergel & Steinhardt, 2000; Loeb & Weiner, 2011; Vogelsberger et al., 2012; Zavala et al., 2013; Elbert et al., 2015; Fry et al., 2015). However, assigning galaxies with low HI rotational velocities to relatively large halos in order to reproduce the observed VF has also proven challenging in the $\Lambda$CDM context. This is because the internal kinematics of dwarfs seem to indicate low-mass hosts. Ferrero et al. (2012) demonstrated that galaxies with stellar masses in the 106-108 Mstar range appear to be hosted by much smaller halos than predicted by abundance matching, based on their observed rotation velocities. Papastergis et al. (2015) extended this to a much larger sample, but confirmed that galaxies with HI rotation velocities below $\sim$25 km s*-1* were incompatible with residing in the more massive halos that abundance matching predicts. Garrison-Kimmel et al. (2014) used the observed densities of local group dwarf irregulars to derive the halo masses that they reside in. All of the galaxies seemed to be in halos of similar mass, but it was a much lower mass than predicted by abundance matching. They concluded that it does not seem possible to simultaneously reproduce the measured velocity function (i.e., satisfy abundance matching) and the observed densities of galaxies, an issue referred to as the too big to fail problem.

However, recent work by Brook & Di Cintio (2015a, b) used results from simulations in which stellar feedback processes alter the dark matter content of dwarf galaxies to show that they can simultaneously match the densities and velocities of observed dwarfs. In this scenario, feedback from stars and supernovae create bursty star formation histories in dwarf galaxies that fluctuate the gravitational potential well at the center of the dwarf (Pontzen & Governato, 2012; Teyssier et al., 2013; Pontzen & Governato, 2014; Oñorbe et al., 2015; Chan et al., 2015; Dutton et al., 2016). Dark matter core creation leads to a better match between theory and observed rotation curves (Katz et al., 2017; Santos-Santos et al., 2017). Feedback is particularly effective in dwarf galaxies with halo masses of a few 1010 M*⊙* (Governato et al., 2012; Di Cintio et al., 2014b), where it can transform an initially steep inner dark matter density profile into a flatter “cored” profile. At lower halo masses there is less star formation, leading to less energy injection and lower core formation efficiencies (Peñarrubia et al., 2012; Maxwell et al., 2015). At higher masses, the deeper potential wells of galaxies make core formation increasingly difficult (Di Cintio et al., 2014b; Pontzen & Governato, 2014), at least if an additional source of feedback is neglected, such as AGN (Martizzi et al., 2013). In this model with baryonic feedback, it is possible to assign dwarf galaxies to relatively massive halos, despite the low rotational velocities measured from their spatially resolved stellar kinematics. This is because baryonic feedback can push dark matter out of the central regions, lowering the enclosed mass at the radii that stellar kinematics probe (but without affecting the total halo mass111modulo a slight reduction in halo mass caused by the loss of baryons or preventive feedback (Munshi et al., 2013; Sawala et al., 2013)). Hence, the densities are lowered, and the apparent velocities of the galaxies, reconciling the observations with theory.

Based on such simulations, Di Cintio et al. (2014a) derived an analytic model for the dark matter density profile that varies with stellar-to-halo mass ratio. Brook & Di Cintio (2015a) and Brook & Di Cintio (2015b) used this analytic model to derive galaxy trends that they claim reconcile the halo densities and the observed VF. In this work, we use simulations directly. These simulations also create dark matter cores (Governato et al., 2012; Zolotov et al., 2012; Christensen et al., 2014; Brooks & Zolotov, 2014), following very similar trends to those in Di Cintio et al. (2014a). However, because we use the simulations directly, we do not have to resort to analytic models for the baryon distribution in the galaxies. Macciò et al. (2016) also recently used simulations directly to show that baryonic simulations can be reconciled with observations. However, they did not investigate the role of dark matter cores in their results. We show that accounting for the gas distribution is important and not straightforward. Unlike Brook & Di Cintio (2015b), we do not find that dark matter core creation consistently has a large impact on observed velocities of galaxies, yet we do find that we can reproduce the observed VF.

This paper is organized as follows: In Section 2 we present information about the simulations. In Section 3 we demonstrate that deriving velocities from baryons yields a substantially lower velocity in dwarf galaxies than expected from theoretical results that rely on dark matter-only simulations. In Section 3.1 we explore how completeness (i.e., the number of detectable halos at low velocity) affects the observed VF. In Section 3.2 we describe our method to mimic observations and derive velocities in as close a way as possible to the observations. In Section 3.3 we re-derive the expected VF given our completeness results and mock observed velocities. In Section 3.4 we demonstrate that our simulations match other essential scaling relations. We systematically explore the importance of various effects in reducing the observed velocities relative to to theoretical velocities in Section 4. In Section 5 we explore the role of dark matter cores on the reduced observed velocities. We find that cores are only important in galaxies where the velocity is measured interior to the size of the core. We compare are results to previous work in Section 6, and conclude in Section 7.

2. The Simulations

The high-resolution simulations used in this work were run with pkdgrav (Stadel, 2001) and its baryonic (SPH) version gasoline (Wadsley et al., 2004), using a $\Lambda$ Cold Dark Matter ($\Lambda$CDM) cosmology with $\Omega_{m}$ = 0.24, $\Omega_{\Lambda}$ = 0.76, $H_{0}$ = 73 km s*-1*, $\sigma_{8}$=0.77, and n=0.96. The galaxies were originally selected from two uniform dark matter-only simulations of 25 and 50 comoving Mpc per side. From these volumes, eight field–like regions where selected, each centered on a galaxy with halo mass222The virial radius is defined relative to critical density, $\rho_{c}$, where the mean density enclosed is $\rho/\rho_{c}\approx 100$ at $z=0$. ranging from 1010 to 1012 M*⊙. Each field was then resimulated using the “zoom-in” volume renormalization technique (Katz & White, 1993), which simulates a region out to roughly 1 Mpc of the primary halo at the highest resolution, while fully preserving the surrounding large scale structure that builds angular momentum in tidal torque theory (Peebles, 1969; Barnes & Efstathiou, 1987). These simulations were run from approximately $z=150$ to $z=0$. A uniform UV background turns on at $z=9$, mimicking cosmic reionization following a modified version of Haardt & Madau (2001). The rms mass fluctuation relative to the cosmic average, $\delta\rho$/$\rho$, for each chosen field ranges from -0.15 to1.35 when measured on a scale of 8h-1* Mpc (see Table 1). Five of the fields fall within 0.05 standard deviations of the cosmic mean density.

The spline force softening, $\epsilon$, ranges from 64 pc to 174 pc in the high resolution regions (see Table 1), and is kept fixed in physical pc at $z<10$. The dark matter (DM) and stellar mass resolutions are listed in Table 1. The gas smoothing length is allowed to shrink as small as 0.1$\epsilon$ in very dense regions (0.5$\epsilon$ is typical) to ensure that hydro forces dominate at very small scales. The main galaxy in every zoomed region contains several millions of DM particles within its virial radius.

The high resolution of these cosmological simulations allows us to identify the high density peaks where H2 can form. We track the non-equilibrium formation and destruction of H2, following both a gas-phase and a dust (and hence metallicity) dependent scheme that traces the Lyman-Werner radiation field and allows for gas and dust self-shielding (Gnedin et al., 2009; Christensen et al., 2012). We include cooling from both metal lines and H2 (Shen et al., 2010; Christensen et al., 2012). Metal cooling, H2 fractions, and self-shielding of high density gas from local radiation play an important role in determining the structure of the interstellar medium and where star formation can occur (Kennicutt, 1998; Krumholz & McKee, 2008; Bigiel et al., 2008; Blanc et al., 2009; Bigiel et al., 2010; Schruba et al., 2011; Gnedin & Kravtsov, 2011; Narayanan et al., 2012). With this approach, we link the local star formation efficiency directly to the local H2 abundance. As described in Christensen et al. (2012), the efficiency of star formation, $c^{*}$, is tied to the H2 fraction, $X_{\rm H_{2}}$. The resulting star formation rate (SFR) depends on the local gas density such that SFR $\propto$ $c^{*}$X${}_{\rm H_{2}}(\rho_{\rm gas})^{1.5}$, with $c^{*}=0.1$. This value of $c^{*}$ gives the correct normalization of the Kennicutt-Schmidt relation.333Note that the efficiency of star formation in any given region is actually much lower than the implied 10%, due to the fact that feedback from newly formed stars quickly disrupts gas, shuts off cooling, and lowers the overall efficiency (e.g., Agertz et al., 2013). Because star formation is restricted to occurring in the presence of H2, stars naturally form in high density regions ($>$ 100 amu cm*-3*), with no star formation density threshold imposed.

Star particles represent a simple stellar population born with a Kroupa et al. (1993) initial mass function. The star particles lose mass through stellar winds and supernovae (SN Ia and SN II). Supernovae deposit 1051 ergs of thermal energy into the surrounding gas following the “blastwave” scheme described in Stinson et al. (2006). No velocity “kicks” are given to the surrounding gas particles, but thermal energy is deposited and cooling is turned off within a “blastwave” radius and adiabatic expansion phase calculated following Ostriker & McKee (1988). The thermal energy deposition from supernovae can lead to bubbles of hot gas that expand, driving winds from the galaxies. Unlike other “sub-grid” schemes, the gas stays hydrodynamically coupled while in galactic outflows. Despite its reliance on supernovae, this model should be interpreted as a scheme to model the effect of energy deposited in the local interstellar medium by all processes related to young stars, including UV radiation from massive stars (Hopkins et al., 2011; Wise et al., 2012; Agertz & Kravtsov, 2015). The rate of ejected mass in winds in these simulations is dependent on galaxy mass, ranging from less than the current SFR in Milky Way-mass galaxies, to typically a few times the current SFR in galaxies with $v_{circ}\sim 50$ km s*-1*, to more than 10 times the current SFR in galaxies with $v_{circ}\sim 20$ km s*-1*(Christensen et al., 2016). These ejection rates are similar to what is observed in real galaxies over a range of redshifts (Martin, 1998; Kirby et al., 2011; Kornei et al., 2012). Additionally, Munshi et al. (2013) demonstrated that these simulations match the observed stellar mass to halo mass relation (Moster et al., 2013), by creating a more realistic star formation efficiency as a function of galaxy mass.

The star formation and feedback in these simulations leads to important trends in the resulting galaxies that are important for the present study. First, feedback strongly suppresses star formation, but the amount of suppression scales with galaxy mass (Brooks et al., 2007). In the deeper potential wells of massive galaxies, high densities make it easier for the gas to cool quickly after being heated by supernovae. The lower densities in dwarf galaxies are more susceptible to heating, driving the star formation efficiencies even lower in dwarfs. Hence, even though the simulated dwarf galaxies may lose much of their gas in winds (Christensen et al., 2016), the gas that stays behind is very inefficient at forming stars, so that the dwarfs are very gas rich (Munshi et al., 2013). Second, when star formation is tied directly to high density regions with H2, subsequent feedback causes these cold, dense regions to become massively over-pressurized. This leads to very bursty star formation histories in dwarf galaxies (Domínguez et al., 2015; van der Wel et al., 2011; Kauffmann, 2014). Bursty star formation creates fluctuations in the galaxy potential well, particularly in halos with masses a few 1010 M*⊙* (Governato et al., 2012; Di Cintio et al., 2013, 2014b), which causes initially cuspy dark matter density profiles to transform into flatter “cores”. In Section 5 we examine whether this core formation lowers the measured $v_{rot}$ of halos from that predicted in DM-only simulations.

3. The Impact of Baryons on the VF

The goal of this study is to identify whether baryonic processes can reconcile the VF expected theoretically in a $\Lambda$CDM universe with observations. The first attempts to compare the theoretical and observational VFs were based on the results of DM-only cosmological simulations (e.g., Zavala et al., 2009; Zwaan et al., 2010; Papastergis et al., 2011; Trujillo-Gomez et al., 2011). However, this simple approach neglects several important baryonic effects that are important for making a fair comparison between theory and observations. Here, we analyze the impact of baryonic effects on the theoretical VF by using simulations run both with baryons and as DM-only. Our aim is to perform “mock observations” of our baryonic simulations, and derive a theoretical VF in as similar a way as possible to current observational determinations.

In what follows, various definitions of velocity arise. In order to compare results from observations to results from theory, one must define a characteristic galaxy velocity that can be compared. In practice, each defined characteristic velocity is slightly different, being derived in a slightly different way. Below, we explore in detail the results of various definitions of characteristic velocity. To minimize confusion for the reader, in Table 2 we define each velocity that we use in the remainder of this paper.

We focus our comparison on the VF measured in the Local Volume ($D\lesssim 10$ Mpc) by Klypin et al. (2015), based on the catalog of nearby galaxies of Karachentsev et al. (2013). The catalog is optically selected, and probes with reasonable completeness galaxies as faint as $M_{B}=-10$. The majority ($\sim$80%) of galaxies in this Local Volume catalog have measurements of their rotational velocity based on the width of their HI profile, $w_{50}$. Some fraction of galaxies lack HI data, either because they are intrinsically gas-poor (e.g., satellites of nearby massive galaxies), or because they have not been targeted by HI observations. These galaxies are assigned rotational velocities based on stellar kinematic measurements when available, or otherwise according to an empirical luminosity-velocity relation. Note that the Klypin et al. (2015) VF is consistent with other independent observational measurements of the VF, such as the one performed by the ALFALFA blind HI survey (Papastergis et al., 2011; Papastergis & Shankar, 2016).

In order to make an appropriate comparison with the observational VF measured by Klypin et al. (2015), we need to model two key observational effects using our baryonic simulations. First, we need to replicate the completeness limitations of the Karachentsev et al. (2013) catalog at low luminosities. Faint galaxies tend to have low rotational velocities. Hence, if a halo is not detectable in current surveys, the density of observed galaxies at the low velocity end of the VF will be surpressed relative to theoretical expectations that populate each halo with a detectable galaxy. Thus completeness can significantly impact the measurement of the low-velocity end of the VF and must be accounted for. Second, we need to compute the theoretical VF in terms of the rotational velocity measured observationally, $w_{50}$. This entails deriving realistic estimates of the HI linewidths for our baryonic halos.

3.1. Detectability of halos

Each “zoomed” simulation contains a high resolution region centered on a halo, ranging in virial mass from 1010 M*⊙* to 1012 M*⊙*. In addition to the central halo, every zoomed region contains smaller galaxies that we also include in our analysis. Because the theoretical VF is traditionally derived using results from DM-only simulations, for every simulated baryonic halo we identify its counterpart in the DM-only run in order to assign a $v_{max,dmo}$ value, the maximum circular velocity in the DM-only runs. Because the DM particles are identical in both the baryonic and DM-only initial conditions, identifying a counterpart is relatively straightforward. For all halos in the DM-only run with more than 64 particles, we identify the DM particles that make up each halo in the DM-only run, then find those same particles in the baryonic run and note the halo444Halos are identified with AHF, AMIGA’s Halo Finder (Gill et al., 2004; Knollmann & Knebe, 2009). AHF is available for download at http://popia.ft.uam.es/AMIGA/. that most of those particles belong to. We find a matching counterpart for 6271 halos and subhalos.

We use this sample to compute the fraction of halos hosting simulated galaxies with $M_{\ast}>10^{6}\;M_{\odot}$ in the baryonic runs, as a function of the maximum circular velocity in the DM-only runs, $f_{det}(v_{max,dmo})$. The $M_{\ast}>10^{6}\;M_{\odot}$ cutoff is chosen because it corresponds to the typical stellar mass of galaxies with $M_{B}=-10$, which define the faint limit of the Klypin et al. (2015) measurement. The result is shown in Figure 1. As the figure shows, virtually all halos with $v_{max,dmo}$ $\gtrsim 35$ km s*-1* host detectable galaxies, and thus are expected to be included in the VF measurement of Klypin et al. (2015). On the other hand, the detectable fraction drops precipitously at lower values of $v_{max,dmo}$, falling below the 5% level at $v_{max,dmo}$ $\lesssim 25$ km s*-1*. As shown in §3.3, this sharp drop in the fraction of detectable galaxies at low values of $v_{max,dmo}$ has important consequences for the measurement of the low-velocity end of the VF.

Keep in mind that the value of $v_{max,dmo}$ where the dramatic drop in detectability takes place is dependent on the depth of the galaxy catalog used to measure the VF. If a deeper census of Local Volume galaxies were available, the minimum detectable stellar mass would be lower than $\sim 10^{6}\;M_{\odot}$, and the drop in detectability would consequently appear at lower values of $v_{max,dmo}$ than shown in Figure 1. Eventually, a physical effect will limit galaxy formation in halos with very low values of $v_{max,dmo}$, namely reionization feedback (e.g., Okamoto et al., 2008; Sawala et al., 2015).

3.2. Mock “observed” rotational velocities

Most of the galaxies ($\sim$80%) in Klypin et al. (2015) have rotational velocities derived from HI. For this reason, we analyze the HI content of our baryonic halos and derive observationally motivated rotational velocities for our simulated galaxies that contain enough HI mass to fall into the Karachentsev et al. (2013) catalog. Klypin et al. (2015) also include dispersion-supported galaxies with no measurable HI down to $M_{B}=-10$. In this section, we describe our selection criteria to mimic this sample and derive mock observational velocities.

To restrict our sample to halos with enough baryonic material to fall into the Klypin et al. (2015) sample, we identify all halos in the DM-only zoomed runs that have $v_{max,dmo}\geq 15$ km s*-1* at $z=0$ and their counterparts in the baryonic zoomed runs. This yields an initial sample of 57 halos. From this initial sample, we identify those with an HI mass, MHI, greater than 106 M*⊙*, corresponding to the HI mass of the faintest galaxies in the Karachentsev et al. (2013) catalog that have HI linewidth data. This yields a sample of 42 galaxies with enough HI mass to generate mock HI data cubes (described below). Of the remaining gas-poor galaxies, we keep only those with $r$-band magnitudes brighter than -10 in the baryonic runs, to approximately mimic the $M_{B}=-10$ limit of the Karachentsev et al. (2013) catalog. Five out of the initial sample of 57 halos are fainter than this $r$-band cutoff, and are therefore not included in the subsequent analysis.

Ten gas-poor halos with HI masses below our adopted cutoff, $M_{HI}<10^{6}\;M_{\odot}$, remain in the sample. Four of these dispersion-supported halos with no HI are satellites of a Milky Way-mass galaxy, and we adopt for them the stellar velocity dispersion as the mock “observed” velocity. For the other six faint galaxies without HI data cubes, we adopted the procedure of Klypin et al. (2015), who assigned a velocity dispersion of 10 km s*-1* to all halos with $M_{K}$ fainter than -15.5. Hence, we assign these halos a velocity dispersion of 10 km s*-1*. However, whether we use a fixed 10 km s*-1* or the stellar velocity dispersion measured directly from the simulation makes no change to our results, as the simulated velocity dispersions are on the order of 10 km s*-1*, similar to the observational data.

The HI mass fraction of every gas particle in the baryonic runs is calculated based on the particle’s temperature, density, and the cosmic UV background radiation flux, while including a prescription for self-shielding of H2 and dust shielding in both HI and H2 (Christensen et al., 2012). This allows for the straightforward calculation of the total HI mass of each simulated galaxy. We create mock HI data cubes only for the 42 halos that contain M${}_{\rm HI}>$ 106 M*⊙. Specifically, we create mock data cubes that mimic ALFALFA observations (Haynes et al., 2011). After specifying a viewing angle (see below), our code considers the line-of-sight velocity of each gas particle. The velocity of each particle is tracked in the simulation by solving Newton’s equations of motion, but any turbulent velocity of the gas is not taken into account. Velocity dispersions in dwarf galaxies can be on the order of the rotational velocity, $\sim$10-15 km s-1* (e.g., Stanimirović et al., 2004; Tamburro et al., 2009; Oh et al., 2015). Dispersions are thought to be driven at least partially by thermal velocities or supernovae (Tamburro et al., 2009; Stilp et al., 2013b, a). In our simulations, supernovae inject thermal energy, and the thermal state of the HI gas needs to be considered in the mock HI linewidth for a realistic comparison to observations. To account for the thermal velocity, the HI mass of each gas particle is assumed to be distributed along the line-of-sight in a Gaussian distribution with a standard deviation given by the thermal velocity dispersion, $\sigma=\sqrt{kT/m_{HI}}$, where $T$ is the temperature of the gas particle. After this thermal broadening is calculated, a mock HI data cube can be generated by specifying the spatial and velocity resolution. For all of our galaxies, we adopted a spatial resolution of 54 pixels across 2$R_{vir}$. In practice, this corresponds to $\sim$1kpc resolution in our lowest mass galaxies to up to $\sim$9kpc resolution in our most massive galaxies. However, the spatial resolution plays no role in our study, since measurements of the VF are based on spatially unresolved HI data. For the velocity resolution, we match the ALFALFA specification of 11.2 km s*-1* (two-channel boxcar smoothed).

For each of the 42 galaxies with $M_{HI}>10^{6}\;M_{\odot}$, we create two HI data cubes. In the first case, we orient each galaxy to be viewed edge-on, i.e., such that the HI angular momentum vector is lying in the image plane. This generates HI data cubes without inclination effects. In the second case, we pick a random orientation of each simulated galaxy (the x-axis of the simulation volume in all cases) and generate HI data cubes that capture inclination effects. In both cases, we measure the width of the HI profile at 50% of the peak height. Hereafter, we denote the edge-on velocity width by $w_{50}^{e}$, while we denote the velocity width projected at a random inclination angle by $w_{50}$. The latter projected velocity width, $w_{50}$, is the one that can be directly measured observationally. For the 10 gas-poor, dispersion-dominated galaxies, we define both $w_{50}^{e}$ and $w_{50}$ to be twice the stellar velocity dispersion.

Example HI line profiles for three simulated galaxies spanning a large range of mass are shown in Figure 2. The HI profiles are derived at random inclination angles, which are indicated in each panel. The figure demonstrates how the HI rotational velocity can differ from the simple theoretical expectation based on the DM-only runs. In particular, we compare the measured $w_{50}$ of the simulated galaxies to its simplest theoretical equivalent555The form of the “theoretical velocity width”, $2\,v_{max,dmo}\,\sin i$, follows from the fact that the HI profiles plotted in Fig. 2 are projected on a viewing angle of inclination $i$, and include emission from both the approaching and receding sides of the HI disk. In the text we will generally use $2\,v_{max,dmo}$ to compare DM-only velocities with edge-on HI velocity widths, $w_{50}^{e}$, and $2\,v_{max,dmo}\,\sin i$ to compare with projected HI velocity widths, $w_{50}$., $2\,v_{max,dmo}\,\sin i$. The example simulated galaxies shown in Figure 2 demonstrate a trend that has a profound impact on the computation of the theoretical VF. In particular, the HI velocity for massive galaxies is larger than the DM-only velocity, w_{50}$$>$$2\,v_{max,dmo}\,\sin i. This shift to higher velocities in the baryonic run is attributed to the cooling of baryons onto the central halo in massive galaxies (e.g., Trujillo-Gomez et al., 2011; Dutton et al., 2011). On the contrary, low-mass simulated galaxies display the opposite effect. The single-peaked shape of their HI profile leads to a measured value of $w_{50}$ that is significantly smaller than $2\,v_{max,dmo}\,\sin i$.

Figure 3 shows the relation between the mock observational and DM-only rotational velocities for all our baryonic halos. More specifically, we compare the edge-on velocity widths, $w_{50}^{e}$, with the equivalent edge-on DM-only widths, $2\,v_{max,dmo}$. This is done in order to facilitate a direct comparison that neglects inclination effects. The red points show the average relation in bins containing four (in the highest velocity bins) to nine (in the lowest velocity bins) data points, depending on the density of the data. Error bars reflect the 1$\sigma$ standard deviation about the average. The dashed line in both panels shows a one-to-one relation between the baryon and DM-only results. It is obvious from this plot that galaxies with 2v${}_{max,dmo}\gtrsim 150$ km s*-1* show higher velocities in the baryonic runs than the DM-only runs, while the trend is reversed at lower masses. The dotted line in the right panel shows the decrease expected in velocity from the DM-only runs if all of the baryons had been lost from the halo. The lowest mass galaxies show a much larger change than can be explained due to baryon loss alone.666Note that even if a simulated dwarf galaxy loses a large percentage of the cosmic baryon fraction it remains gas-rich at $z=0$ due to the fact that the gas that remains behind is inefficient at forming stars, unless it is a satellite and has had its gas stripped. We dissect the reasons for this lower-than-expected velocity in Section 4.

Twenty of the 52 halos plotted in Figure 3 are subhalos (denoted by squares) of larger halos. As seen in this figure and those that follow, the simulated galaxies hosted by subhalos follow similar kinematic trends to those hosted by central halos.

3.3. Re-deriving the Expected VF

Based on the results of §3.1 and §3.2, we can now compute a realistic expectation for the VF of galaxies in a $\Lambda$CDM universe. The process is illustrated in Figure 4. In particular, we start from the VF of halos in a $\Lambda$CDM universe with Planck cosmological parameters (Planck Collaboration et al., 2014). This DM-only VF is plotted as a black dashed line in Fig. 4, and is obtained from the BolshoiP dissipationless cosmological simulation (Rodríguez-Puebla et al., 2016). The halo VF represents the number density of halos as a function of their maximum circular velocity $v_{max,dmo}$. We denote the theoretical DM-only VF by

[TABLE]

In the equation above, $dN_{H}$ is the number of halos contained in a representative volume element $dV$ of the universe, that have rotational velocities within the logarithmic velocity bin $d\log_{10}(v_{max,dmo})$.

Second, we correct the plotted DM-only halo VF to take into account the detectability of halos as a function of $v_{max,dmo}$. We perform this correction based on the result of Figure 1. In particular,

[TABLE]

The corrected DM-only VF is plotted in Fig. 4 as thin grey lines. The bundles of lines represent the uncertainty due to the number of simulated halos used to make Figure 1.

Lastly, we compute the change in the theoretical VF that is due to the difference between the theoretical and observational measures of rotational velocity,

[TABLE]

This is done by first generating a large number of $v_{max,dmo}$ values according to the DM-only halo VF corrected for halo detectability (grey lines in Fig. 4). We then assign to each generated halo an edge-on velocity width value, $w_{50}^{e}$, based on the mean and scatter of the $w_{50}^{e}$-$2\,v_{max,dmo}$ relation shown in Figure 3. Lastly, we calculate the projected HI velocity width as $w_{50}=w_{50}^{e}\times\sin i$. Inclination values, $i$, are drawn assuming random orientations, i.e., such that $\cos i$ is uniformly distributed in the $[0,1]$ interval. The final results for the baryonic VF expected in a $\Lambda$CDM cosmology according to our simulations are shown by the blue lines in Figure 4. The bundles again represent the uncertainty due to the combined uncertainties introduced by the number of simulated halos used to calculate detectability and the number of galaxies in each of the bins in Figure 3. This distribution is directly comparable to the observational VF measured in the Local Volume by Klypin et al. (2015).

Figure 4 clearly demonstrates that taking into account both the “observed” velocities and the luminous fraction of halos has a dramatic effect on the theoretical VF. At the high velocity end, the baryonic VF displays a higher normalization than the DM-only distribution, which is caused by the fact that the HI velocity width, $w_{50}$, is larger than $2\,v_{max,dmo}\,\sin i$ for massive halos (refer to Figs. 2 & 3; though note the effect appears less strong in Figure 4 because it shows 2$v_{max,dmo}$ instead of $2\,v_{max,dmo}\,\sin i$). However, baryonic effects have their largest impact on the low-velocity end of the theoretical VF. In particular, the fact that low-mass halos have $w_{50}^{e}$ values significantly smaller than than 2$v_{max,dmo}$ means that the theoretical VF systematically “shifts” towards lower velocities in the dwarf regime. This translates into a substantial reduction of the VF normalization at $w_{50}\lesssim 100$ km s*-1*.

At even lower velocities, $w_{50}\lesssim 40$ km s*-1*, the very low detectability of small halos further suppresses the normalization of the baryonic VF. Together, the effects of the baryonic velocity shift and of halo detectability lead to a dramatic decrease on the number of low-velocity galaxies expected in $\Lambda$CDM, compared to the simplistic DM-only estimate. As Fig. 4 shows, the difference is more than an order of magnitude already at $w_{50}=50$ km s*-1*. This huge suppression in the number density at low velocities brings our theoretical VF in agreement with the observational measurements, and shows no signs of the overproduction of dwarf galaxies typically encountered in $\Lambda$CDM.

3.4. Validation Against Other Scaling Relations

A key point regarding the results of Fig. 4 is that reproducing the observational VF in a simulation is not physically meaningful unless the typical HI disk sizes in dwarf galaxies are also reproduced correctly. This is because the ratio between $w_{50}^{e}$ and $2\,v_{max,dmo}$ can be made arbitrarily small in dwarf galaxies by producing simulated galaxies with very small HI disks. Because the innermost portion of the rotation curve is rapidly rising, it could be possible to reproduce the observed VF but not accurately reproduce observed disk sizes. Figure 5 compares the sizes of HI disks in our simulated galaxies with the observed sizes in the sample of galaxies with interferometric HI observations compiled by Papastergis et al. (2015). The observational datapoints show the outermost radius where the HI rotational velocity can be measured by the interferometric observations for each galaxy. One complication here is that the outermost HI radius for the galaxies in the Papastergis et al. (2015) sample is not defined in a consistent way, but depends on the depth of each interferometric observation and the quality of each galaxy’s kinematics. For the simulations, we derive “outermost” HI radii where the HI surface density profiles of our simulated galaxies fall below 1 M*⊙/pc2. The adopted HI surface density cutoff corresponds to the value probed by typical interferometric HI observations. We examined the results using different definitions of “outermost” HI radius for our simulated galaxies, and found that the results were generally consistent but that this definition produces the least scatter. This is not surprising, because we have also verified that our simulated galaxies follow the observed HI mass – radius relation from Wang et al. (2016), where the HI radius is again defined at the 1 M⊙*/pc2 isophote. The observed relation has remarkably low scatter, so it is reassuring that using a similar definition for the simulations also produces the smallest scatter. Overall, Fig. 5 shows that our simulated galaxies have HI disk sizes that are in agreement with observations, indicating that the mock observational velocities computed in §3.2 are realistic.

Similarly, the fraction of detectable halos computed in §3.1 is not physically meaningful unless our simulations reproduce the baryonic content of real galaxies. In Figure 6 we show the baryonic (cold gas plus stellar mass) Tully-Fisher relation for the simulated galaxies used in this work (black points, top panel). We restrict ourselves to central galaxies only (excluding subhalos) for comparison to the observational data, which is taken from McGaugh & Schombert (2015). The line in both panels is the baryonic Tully Fisher relation fit to observed galaxies in McGaugh & Schombert (2015, their figure 6 and table 5), $log(M_{b})=1.61+4.04V_{f}$. Since the McGaugh & Schombert (2015) measurement refers to the flat outer velocity of galactic rotation curves, we adopt for the simulations the circular velocity of the baryonic runs measured at 4 disk scale lengths as $V_{f}$. The bottom panel is for the cold gas mass (1.33MHI* in the simulations) only. The simulations have been divided into a gas-rich (Mstar/1.33M${}_{HI}>2.0$, blue points) and gas-poor (Mstar/1.33M${}_{HI}<2.0$, red points) sample. Like the observational data, gas-rich galaxies follow the observed baryonic Tully Fisher relation, while gas-poor galaxies lie below the relation (the dotted line shows a reduction of the relation by a factor of 5). This plot demonstrates that our simulated galaxies match the stellar and HI masses of galaxies as a function of velocity.

Overall, Figures 5 & 6 give us confidence that the theoretical VF computed in §3.3 is physically well motivated. Consequently, moving from predictions based on DM-only runs to baryonic simulations may be the key to reconciling the theoretical expectation of the VF with the observational measurements.

4. Velocity Changes in the Presence of Baryons

In this section, we examine the baryonic effects that lead to dwarf halos being observed at lower velocities than predicted based on DM-only simulations, and that help to reconcile the theory with the observations.

It is well known that the rotation curves of many dwarf galaxies are still rising at their outermost measured point (Catinella et al., 2006; de Blok et al., 2008; Swaters et al., 2009; Oh et al., 2011), suggesting that the true $v_{max}$ of the halo is higher than HI measures. In this section we use the velocity at the outermost HI data point in our baryonic simulations in order to determine how much of a role this plays in the lowered velocities we see in the dwarf simulations compared to their DM-only vmax,dmo values. Recall that in Figure 5 we defined the outermost HI data point, Rout, in our simulations to be the point at which the HI surface density falls below 1 M*⊙/pc2. In what follows, we refer to the circular velocity at Rout* as $v_{out}$.

In our dwarf galaxies, the radius of the outermost HI data is generally still on the rising part of the rotation curve. We quantify this in the top panel of Figure 7, where we compare $v_{out}$ to the maximum value of the circular velocity in the baryonic run, $v_{max,sph}$.777Note that up until now the $v_{max}$ we have been dealing with comes from the DM-only runs, $v_{max,dmo}$. $v_{max,sph}$ will differ from from $v_{max,dmo}$ due to processes like baryonic contraction at high masses, or loss of most of the baryons from the smallest mass halos. We wish to quantify how well HI traces the rotation velocity after these other factors have had their influence, and ultimately determine how well $w_{50}$ is tracing the outermost HI rotation velocity. Hence, we switch to $v_{max,sph}$ in Figure 7. In the more massive galaxies, $v_{out}$ is indeed capturing the maximum value of the rotation curve. However, in galaxies below $\sim$50 km s*-1*, the outermost HI rotation velocity systematically underestimates $v_{max,sph}$.

Next we wish to know if $w_{50}^{e}$ is tracing the outermost velocity, $v_{out}$. The second panel in Figure 7 shows the ratio of the two. In the four most massive galaxies, $w_{50}^{e}$ traces a slightly higher velocity than the outermost HI rotation velocity, due to the fact that these galaxies have large bulges and higher velocities near their center. More importantly for interpreting dwarf galaxy data, $w_{50}^{e}$ is systematically smaller than $v_{out}$. In other words, $v_{out}$ is already under-measuring the maximum rotational velocity of the galaxy because it is on the rising rotation curve, but $w_{50}^{e}$ is measuring an even lower velocity. This suggests that $w_{50}^{e}$ may be measuring a velocity even closer to the center than $v_{out}$.

Evidence that this is the case is found in the third panel of Figure 7, where we compare $w^{e}_{20}$ to $v_{out}$ instead. $w^{e}_{20}$ measures the width of the HI profile at 20% of the peak height rather than 50%. While it slightly overestimates $v_{out}$ in galaxies above $\sim$50 km s*-1*, it does a much better job of capturing $v_{out}$ in the lower mass galaxies. In summary, it seems that $w^{e}_{20}$ is a more reliable indicator of the outermost measurable rotation velocity in the dwarf galaxies.

Finally, the bottom panel of Figure 7 shows the ratio between our $w^{e}_{50}$ and $w^{e}_{20}$ measurements, and demonstrates that $w^{e}_{20}$ can measure a much larger velocity in the dwarfs than $w^{e}_{50}$, up to a factor of two larger in the lowest mass galaxies. This difference has been noted previously. Using ALFALFA data, Bradford et al. (2015) showed that the difference between the two velocities is well described by the relation $w^{e}_{20}=w^{e}_{50}+25$ km s*-1* (see also Koribalski et al., 2004). This relation is shown as the black line in the bottom panel of Figure 7. Brook et al. (2016) showed that the discrepancy between $w^{e}_{20}$ and $w^{e}_{50}$ could lead to substantial differences in the slope of the baryonic Tully Fisher relation, while Brook & Shankar (2016) demonstrated that the use of $w^{e}_{50}$ instead of $v_{max,dmo}$ could fully explain the difference in the theoretical VF compared to observations. We note that almost all observational measurements of the VF are based on $w^{e}_{50}$ rather than $w^{e}_{20}$ (Zwaan et al., 2010; Papastergis et al., 2011; Klypin et al., 2015) due to the fact that it can be hard to measure the line width at 20% of the peak height due to spectrum noise at typical signal-to-noise ratios.

The change between $w^{e}_{20}$ and $w^{e}_{50}$ is likely due to the shape of the HI profile as a function of mass. As was seen in Figure 2, more massive galaxies exhibit a double-horned profile. The horns are built up due to the piling up of velocity along the flat part of the rotation curve in large spirals. However, lower mass galaxies are usually still rising at the outermost HI data point, as discussed above. This leads to an HI profile that is more Gaussian. The drop-off at the edges of the double-horned profile is rapid, so that the difference between $w^{e}_{20}$ and $w^{e}_{50}$ is small. However, the Gaussian shape in the dwarfs ensures that this is no longer true. Measuring lower in the HI profile can lead to a much larger velocity width. These higher velocities must come from further out on the rotation curve.

In summary, the maximum rotational velocity traced by HI does not generally trace the full $v_{max,sph}$ for dwarf galaxies below $\sim$50 km s*-1*. This is due to the fact that the outermost HI is still on the rising part of the rotation curve. Additionally, $w^{e}_{50}$ does not measure the the outermost HI rotation velocity in dwarf galaxies, compounding the problem further. The combination of these two effects leads to the shift in velocities measured between the baryonic and DM-only simulations seen in Figure 3.

5. Does dark matter core creation matter?

Recent high resolution cosmological simulations of galaxies, including those used in this study, have shown that feedback from young stars and supernovae can create dark matter cores in galaxies (e.g., Governato et al., 2010; Teyssier et al., 2013; Chan et al., 2015; Dutton et al., 2016). Governato et al. (2012) and Di Cintio et al. (2014b) showed that this result varies with stellar mass (and thus also with halo mass, given that there is a stellar-to-halo mass relation). The shallow potential wells of dwarf galaxies at $M_{vir}\sim 10^{10}M_{\odot}$ are particularly susceptible to core creation, but the deeper potential wells of MW-mass galaxies are less so, and galaxies have a harder time creating large cores in lower mass halos that form less stars and therefore inject less energy (Maxwell et al., 2015; Read et al., 2016a).

In this section we explore whether the change in the dark matter profile in dwarf galaxies has any impact on the observed VF. Work by Brook & Di Cintio (2015b) concluded that measuring a theoretical velocity at the radius which reproduces $w_{50}^{e}$ is not enough to match observed velocities in models that retain a cuspy, NFW dark matter density profile. Instead, they showed that additionally considering dark matter core creation could lower the theoretical velocities enough to bring them in line with observations. We demonstrate here that this is true only for galaxies which have R${}_{out}\lesssim 3$kpc.

Assessing the impact of core creation is not simple because the densities in the baryonic simulations may also be subjected to some level of contraction due to the presence of the baryons, and disentangling the two effects is not straightforward. Note that this contraction does not have to be adiabatic contraction of the dark matter, and in fact adiabatic contraction of the dark matter is unlikely to occur in the dwarf regime that we are exploring here. However, as we demonstrate below, the fact that gas can cool to the center of the galaxy can increase the rotation velocity in the inner regions, even in dwarf galaxies, in the baryonic simulations. This effect must be accounted for before a direct comparison can be made between the velocities in the baryonic runs and the DM-only runs. If it is not accounted for, a comparison between the baryonic and DM-only velocities would minimize the impact of dark matter core creation.

To overcome this, we develop a proxy for a contracted model without core creation by adding together the velocity profile in a DM-only run888A DM-only simulation contains the cosmic density of all matter, $\Omega_{baryon}$ and $\Omega_{DM}$. Here, we scale down the velocity profile of the DM-only run by an amount consistent with removing the cosmic baryonic fraction, so that we can add the baryonic contribution from the SPH runs instead. with the velocity profile of only the baryonic component in its counterpart SPH run. This effectively “contracts” the profile due to the presence of baryons, but does not include dark matter cores since the DM-only runs do not experience core creation. We measure the velocity from this combined model at the outermost HI radius, Rout, determined from the SPH runs, and label it $v_{out,dmo+b}$.

In Figure 8 we compare $v_{out}$ measured in the SPH runs to $v_{out,dmo+b}$ as a function of Rout. The ratio $v_{out}/v_{out,dmo+b}$ gives us an estimate of how much core creation alone has supressed the rotation curve in the baryonic runs. The data points are color coded based on the slope of their dark matter density profile, measured between 300-700pc, labeled $\alpha_{500pc}$. We include in this plot all simulated galaxies with HI. The lowest mass galaxies have star formation efficiencies too low to create substantial dark matter cores. Core creation is not the only mechanism that can suppress the rotation curve, as loss of baryons alone can lower the baryonic rotation curve relative to the DM-only case. The ratio expected for pure baryonic mass loss is shown by the dotted line in Figure 8. If core creation is important we would expect to see that the strongly cored galaxies lie systematically lower than other galaxies. We find this is only true for galaxies with R${}_{out}<$ 2-3 kpc.

For galaxies with R${}_{out}<3$ kpc (corresponding to $v_{out}<50$ km s*-1*), the galaxies with dark matter cores generally occupy the lowest velocity ratios. This suggests that core creation contributes to velocity suppression in this regime. The velocities can be lower by up to 40%, comparable to the reduction from measuring on the rising part of the rotation curve alone (see top panel of Figure 7). Thus, cores do seem to substantially contribute to lowered velocities for galaxies with R${}_{out}<3$ kpc.

For galaxies with R${}_{out}>3$ kpc, there are strongly cored galaxies that do not show any signs of having their velocities reduced. In Figure 9, we provide an example of why a galaxy with a strong dark matter core may not have a lower velocity. Figure 9 shows the rotation curve for one of our dwarf galaxies that undergoes significant dark matter core creation. At $z=0$ this halo has a dark matter density slope of -0.3. This profile causes the rotation curve to rise much more slowly in the baryonic run (red solid line) compared to the combined DM-only/baryonic model (solid black line) or the DM-only run (black dashed line). The red dashed line shows the DM contribution to the baryonic run’s total $v_{circ}$, to emphasize the presence of the dark matter core.

It can be seen that the baryonic run has a lower rotational velocity than the combined DM-only/baryonic model interior to $\sim$2 kpc. It is clear from Figure 9 that if the HI is tracing velocity interior to $\sim$2kpc, then core creation would reduce the measured velocity in this galaxy. However, this galaxy has HI gas that extends out to roughly 5 kpc, where it is tracing the flat part of the rotation curve, and is an excellent measure of $v_{max,sph}$.

This galaxy also highlights another subtle point. The DM-only run reaches a $v_{max,dmo}=55.8$ km s*-1* at 27 kpc. The baryonic run reaches $v_{max,sph}=58.3$ km s*-1* at 7.5 kpc. The velocity from the HI profile, $w^{e}_{50}$, is 55 km s*-1*, comparable to the $v_{max,dmo}$ measured in the DM-only run. Thus, there is almost no change in $v_{max}$ between the two runs, i.e., this halo does not undergo adiabatic contraction in the usual sense. It is simply that the radius at which $v_{max}$ occurs is quite different. In the baryonic run, the fact that gas can cool leads to the mass being more centralized than in the DM-only run, without increasing $v_{max}$ overall. Likewise, the “contracted” model combining the DM-only profile with the baryonic profile is not adiabatically contracted, but simply reaches $v_{max}$ at a smaller radius. The cold gas increases the central velocity relative to the DM-only run despite the fact that this dwarf is dark matter dominated overall, with a baryon ratio (cold gas and stellar mass to total DM mass) of only 2% at $z=0$ (but remains gas-rich due to the fact that star formation is inefficient).

There are a total of four galaxies in our sample where the DM-only counterpart has $v_{max,dmo}\sim$ 55 km s*-1*, like the galaxy shown in Figure 9. All of these galaxies have stellar masses between 1.5-3$\times$108 M*⊙, and all have a cored dark matter density profile, but their HI masses vary by an order of magnitude. Two of them have R${}_{out}\sim$ 1.5 kpc, while two have R${}_{out}\sim$ 5 kpc. As expected, the two with small Rout* have substantially lower $w^{e}_{50}$ values compared to $v_{max,dmo}$. Hence, scatter in the HI content at a given halo mass leads to scatter in the role of dark matter cores.

From these examples, we learn that if core creation is to impact the measured velocity in a galaxy, the HI must not extend significantly further than the size of the dark matter core. A similar conclusion was found by Papastergis & Shankar (2016) by analyzing observational dwarf data. In simulated galaxies with efficient core creation, the dark matter cores are often 1-2 kpc. From Figure 8, we see that the strongly cored galaxies with R${}_{out}<2$ kpc do indeed tend to show a lower rotation velocity in the baryonic run than their DM-only counterpart.

6. Comparison with Previous Works

In this section we discuss how our results compare to previous works on this topic. First we focus on the ability to reproduce the VF, then specifically on the impact of dark matter cores.

6.1. Velocities

Brook & Shankar (2016) were the first to show explicitly that using $w^{e}_{50}$ instead of $v_{max,dmo}$ could reconcile the theoretical VF with the observed VF. Their approach was semi-empirical, using abundance matching (a relationship between baryonic mass and halo mass) convolved with a relation between between baryonic mass and velocity. They showed the impact of using various definitions of velocity, with only $w^{e}_{50}$ recovering the observed VF.

Like the work presented in this paper, Macciò et al. (2016) also used cosmological zoomed simulations, the NIHAO suite, to make mock HI profiles, and showed that their measured $w^{e}_{50}$ could reproduce the observed VF. Macciò et al. (2016) followed a similar analysis as in this paper, and both works use galaxies simulated with the code Gasoline, but the simulations vary in terms of details. A slightly lower resolution in most of the NIHAO galaxies prevents the use of H2-based star formation as used here, but NIHAO includes a prescription for early stellar feedback (feedback from young massive stars that is deposited prior to the first SNII from any given star particle). A detailed comparison of mock observed velocities at a given $v_{max,dmo}$ shows that the mock velocities in NIHAO are lower than in this work. Perhaps because of this, Macciò et al. (2016) need not consider completeness in order to reproduce the observed VF; the lower velocities of $w^{e}_{50}$ alone are enough to allow the NIHAO galaxies to match the data (and may even slightly over-reduce the velocities in the lowest halos; see their figure 3).

Thus, this work and both Brook & Shankar (2016) and Macciò et al. (2016) have concluded that the difference between $v_{max,dmo}$ and $w^{e}_{50}$ is the primary reason for the disagreement between theory and observations. An apples-to-apples comparison between models and real galaxies alleviates the tension.

On the other hand, Trujillo-Gomez et al. (2016) attempt to correct observed velocities to their underling $v_{max}$. Using a sample of galaxies with resolved HI rotation curves from Papastergis & Shankar (2016), they fit $v_{out}$ to both NFW and cored rotation curve models in order to infer the true $v_{max}$ of each galaxy. This correction can then be applied to galaxies with unresolved HI velocities of similar baryonic mass. However, they conclude that there is not enough of a shift to resolve the discrepancy between the theoretical and observed VF, even when the effects of dark matter core creation are taken into account.

To reconcile the work of Trujillo-Gomez et al. (2016) with the conclusions in this paper, Brook & Shankar (2016), and Macciò et al. (2016), the correction from observed $v_{out}$ to $v_{max}$ must fail. Mock resolved HI rotation curves of the simulated galaxies should, in principal, be able to address this question. However, results so far are inconclusive. Read et al. (2016b) made mock HI rotation curves of two simulated dwarf galaxies and tested the conditions under which they could reliably recover the model halo masses. They found that starburst and post-starburst dwarf galaxies have large HI bubbles that push the rotation curve out of equilibrium, and that galaxies viewed near face-on also presented problems, but could otherwise recover their model inputs (as long as they used a model with a dark matter core). They concluded that a carefully selected sample should allow for a reliable recovery of true halo masses. In Read et al. (2017) they applied their method to 19 observed galaxies, and derived a stellar mass-to-halo mass relation in agreement with abundance matching results for field galaxies, concluding that there are no dwarf galaxy problems in CDM. On the other hand, Verbeke et al. (2017) failed to recover the true $v_{circ}$ of any of their 10 dwarf galaxies (from the Moria simulation suite) when producing mock HI rotation curves. They conclude that the disks of dwarfs are simply too thick, combined with feedback causing significant structure and disequilibrium so that the HI rotation curve fails to be a good measure of the underlying gravitational potential. Given the mixed results, more work in this area is required.

6.2. Dark Matter Cores

A recent analysis by Brook & Di Cintio (2015b) also examined the effects of dark matter core creation on the observed galaxy VF, comparing to the Local Volume VF derived in Klypin et al. (2015). The top panel of Figure 10 shows the measured velocity dispersion (for HI poor galaxies) or edge-on $w^{e}_{50}$/2 (for HI rich galaxies) versus the stellar mass in the simulated galaxies. The flattening of $w_{50}^{e}$/2 below $\sim$107 M*⊙* is attributed to core creation in Brook & Di Cintio (2015b). This flattening is not reproduced in their models with an NFW profile (they examine galaxies down to 106 M*⊙* in stellar mass). Only their model that includes dark matter core creation reproduces this flattening. Although we also find this flattening to occur at $\sim$107 M*⊙*, the bottom panel of Figure 10 demonstrates that this flattening in velocity cannot be due to core creation, as the trend is found in DM-only runs as well.

While $w_{50}^{e}$/2 is a quantity derived from the baryonic simulations, the bottom panel of Figure 10 shows $v_{max,dmo}$ plotted against the stellar mass of the galaxies in the baryonic version of the runs. Recall that $v_{max,dmo}$ is a quantity derived from the DM-only versions of the galaxies. DM core creation requires the presence of baryons, and hence cores cannot form in the DM-only runs. The galaxies in the DM-only runs retain a steep, cuspy DM density profile. Despite the steep inner profile, the flattening of the trend at low stellar masses persists in each panel. Hence, core creation cannot be responsible for the flattening. This is contrary to the conclusions in Brook & Di Cintio (2015b). Reinforcing this conclusion, the data points in Figure 10 are again color coded by the slope of the DM density profile in the baryonic version of the runs. While cored galaxies tend to cluster in a given stellar mass range, they do not appear to play a role in the flattening of the trend below stellar masses of 107 M*⊙*.

We note that the flattening below stellar masses of 107 M*⊙* is consistent with observational data (see figure 1 of Klypin et al., 2015), which show a roughly constant velocity of $\sim$10 km s*-1* and tend to be in dispersion supported galaxies. As previously discussed in Section 3.2, our faintest simulated galaxies have velocity dispersions $\sim$10 km s*-1*, consistent with the observations. This is a more direct comparison to the observations than presented in Brook & Di Cintio (2015b), where they measured $v_{circ}$ of a model galaxy at the radius that best reproduced $w_{50}$ values.

We offer a different interpretation for the flattening of velocities at low galaxy masses: the steep relation between Mstar and $v_{max,dmo}$ (or, equivalently, between Mstar and $M_{halo}$) at low halo masses. As has been noted by previous authors (e.g., Ferrero et al., 2012), the steep relation at low halo masses suggests that galaxies over a wide range of stellar masses (106-108 M*⊙) reside in nearly the same host halo mass ($\sim$1010* M*⊙*, though see Read et al., 2017). The bottom panel of Figure 10 confirms that this trend also occurs in our simulations. All of the low stellar mass galaxies reside in a narrow range of $v_{max,dmo}$ (or equivalently, $M_{halo}$). This will lead them to have similar observed velocities as well, as seen in the top panel.

We note that unlike Brook & Di Cintio (2015b), our results are not in conflict with the conclusions in Brook & Di Cintio (2015a). In that paper, the role of core creation on the Too Big to Fail Problem (Boylan-Kolchin et al., 2011; Garrison-Kimmel et al., 2014) was explored using an analytic model (not simulations) for galaxies that had halo masses determined using stellar velocity dispersions at their half light radii. In all cases, the half light radius is $\lesssim$1 kpc. As we have seen, core creation can reduce the velocities of dwarfs interior to 1 kpc, and will thus alter the derived masses (though the magnitude of the reduction may not be as significant as Brook & Di Cintio (2015a) predicted for dwarf Irregulars, since they neglected gas in the inner kpc). HI, on the other hand, can extend to much larger radii than typical half light radii, and eliminates any impact of dark matter cores on the measured velocity (see also Papastergis & Ponomareva, 2017).

7. Conclusions

In this work, we have used high resolution cosmological simulations of individual galaxies in order to resolve the discrepancy between the observed galaxy VF and the predicted VF within $\Lambda$CDM. In particular we study the apparent dearth of observed low velocity galaxies.

To ensure that the simulated galaxies have realistic sizes and gas contents, and thus can be used to interpret observations, we verified that the simulated galaxies with baryons match observed scaling relations. In particular, the simulations match the HI sizes of galaxies as a function of velocity and the baryonic Tully-Fisher relation.

We use these realistic galaxies to generate mock “observed” velocities. For galaxies with M${}_{HI}>10^{6}$ M*⊙*, we produce mock HI datacubes and derive a characteristic velocity using the width of the HI profile at 50% of the peak height, $w_{50}$. This is the velocity commonly used to generate the observed galaxy VF (Zwaan et al., 2010; Papastergis et al., 2011; Klypin et al., 2015). For gas poor galaxies, we follow the procedure of Klypin et al. (2015) and use stellar velocity dispersion. When the “observed” velocities from baryonic simulations are compared to theoretical velocities (derived from the maximum circular velocity of matched counterpart halos in dark matter-only simulations), we find that there is a systematic shift in dwarf galaxies to lower velocities (see Figure 3). The magnitude of this velocity shift, combined with a proper accounting of luminous halos, reconciles the observed VF with the theoretical VF (Figure 4).

Thus, there are two primary considerations necessary to bring the theoretical VF into agreement with the observed VF. First, to match the observed VF at velocities below $w_{50}$ $\sim$40 km s*-1*, the fraction of luminous halos must be accounted for. If a halo does not host a luminous galaxy, it will remain undetected in current surveys, lowering the observed number of galaxies at low velocities compared to theoretical expectations that allow all halos to host a detectable galaxy. Here, we calculated the luminous fraction for halos with M${}_{*}>10^{6}$ M*⊙, which corresponds to the lower luminosity limit used to calculate the observed VF in Klypin et al. (2015). The fraction of luminous halos drops precipitously below 40 km s-1*. Without considering this effect, the velocity difference alone between our mock observations and theoretical velocities is not sufficient to reproduce the observed VF at the low velocity end. We note that previous work on this subject did not explicitly consider the fraction of luminous halos (e.g., Brook & Di Cintio, 2015b; Macciò et al., 2016).

Second, to match the observed VF it is necessary to derive a relationship between observed characteristic velocities of galaxies and theoretical velocities for halos. We have demonstrated here that this relationship shifts the predicted VF into agreement with the current observations. The source of the velocity shift in dwarf galaxies is a combination of factors:

(1) The primary shift that makes observed velocities lower than theoretical velocities in dwarf galaxies is due to the fact that the velocity tracer (typically HI) does not trace the full potential wells of dwarfs. That is, the outermost HI is still on the rising part of the rotation curve (Catinella et al., 2006; de Blok et al., 2008; Swaters et al., 2009; Oh et al., 2011). We demonstrate this in the top panel of Figure 7, where we explicitly found the circular velocity at the radius that the HI surface density dropped below 1 M*⊙*/pc2, i.e., the outermost observable rotation velocity, $v_{out}$, in simulated galaxies with baryons. The top panel of Figure 7 shows that $v_{out}$ underpredicts the maximum value of the circular velocity in dwarf galaxies.

(2) For galaxies that derive a characteristic velocity using $w_{50}$, there is an additional reduction in observed velocity. We demonstrate this in the second panel of Figure 7, where we compare the velocity results derived from $w_{50}$ to $v_{out}$. Although $v_{out}$ was already lower than the true maximum velocity of a galaxy’s dark matter halo, $w_{50}$ can be an additional 50% lower than $v_{out}$. This is because the HI line profile shape in dwarfs tends to be gaussian. Measuring at a lower peak height, 20%, instead agrees with $v_{out}$ (third panel of Figure 7). To date, essentially all observed VF measurements have been made with $w_{50}$ rather than $w_{20}$ because typical signal-to-noise ratios generally prevent a reliable measurement of $w_{20}$.

(3) For galaxies with HI sizes under $\sim$3 kpc, an additional reduction in velocity can occur if the galaxy has a dark matter core. Typical core sizes found in simulations are on the order of 1-2 kpc, and we demonstrate in Figure 9 that core creation reduces the overall circular velocity in the very center of cored galaxies. However, if the characteristic rotational velocity is derived at a larger radius, then the measured circular velocity is usually comparable to the expected theoretical velocity. We attempted to quantify the contribution of core creation to the reduction in velocity in Figure 8. There, we compared the circular velocity of halos in both the baryonic and contracted (dark matter-only + baryons) models at $R_{out}$. This removes the contribution to the reduction in velocity due to being on the rising part of the rotation curve, and avoids the reduction due to $w_{50}$. Figure 8 shows that galaxies with HI sizes $<$ 3 kpc typically have lower circular velocities than the contracted dark matter models, by up to 40%. This reduction is comparable to the reduction in velocity from measuring on the rising part of the rotation curve alone (see top panel of Figure 7). Hence, core creation leads to a further reduction in observed velocities for galaxies with R${}_{out}<3$ kpc.

Overall, we have demonstrated in this paper that we can start with the abundance of dwarf galaxies predicted in $\Lambda$CDM and reconcile the theoretical predictions with the observed VF. We do this by properly accounting for the relation between characteristic velocities derived from observations and the characteristic velocities typically derived from theory, and by accounting for the fraction of observable halos detectable in current VF studies. We conclude that there is no missing dwarf problem in $\Lambda$CDM.

AB acknowledges support from National Science Foundation (NSF) grant AST-1411399. EP is supported by a NOVA postdoctoral fellowship to the Kapteyn Astronomical Institute. FG was partially supported by NSF grant AST-1410012, HST theory grant AR-14281, and NASA grant NNX15AB17G. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. This work was performed in part at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1066293.