NIHAO XVII: The diversity of dwarf galaxy kinematics and implications for the HI velocity function
Aaron A. Dutton Aura Obreja Andrea V. Macciò


Phonon coherences reveal the polaronic character of excitons in two-dimensional lead-halide perovskites
    Ajay Ram Srimath Kandada
    December 30, 2023
=======================================================================================================






We use 85 pairs of high resolution LCDM cosmological simulations from
the NIHAO project to investigate why in dwarf galaxies neutral
hydrogen (HI) linewidths measured at 50% of the peak flux W_50/2
(from the hydrodynamical simulations) tend to underpredict the maximum
circular velocity  (from the corresponding dark matter only
simulations).  There are two main contributing processes.  1) Lower
mass galaxies are less rotationally supported. This is confirmed
observationally from the skewness of linewidths in bins of HI mass in
both ALFALFA and HIPASS observations.  2) The HI distributions are
less extended (relative to the dark matter halo) in dwarf galaxies.
Coupled to the lower baryon-to-halo ratio results in rotation curves
that are still rising at the last measured data point, in agreement
with observations from SPARC.  Combining these two effects, in both
simulations and observations lower mass galaxies have, on average,
lower W_50/W_20.  Additionally,  mass loss driven by supernovae
and projection effects (dwarf galaxies are in general not thin disks)
further reduce the linewidths.   The implied HI linewidth velocity
function from NIHAO is in good agreement with observations in the
nearby Universe of dwarf galaxies: 10 < W_50/2 < 80. The
dark matter only slope of ≈ -2.9 is reduced to ≈ -1.0
in the hydro simulations.  Future radio observations of unbiased
samples with higher spatial resolution will enable stricter tests of
the simulations, and thus of the LCDM model.



cosmology: theory – dark matter – galaxies: formation – galaxies:
kinematics and dynamics – galaxies: structure – methods: numerical






§ INTRODUCTION


The Lambda Cold Dark Matter (LCDM) model is very successful in
reproducing the large scale structure of the Universe
<cit.> and the anisotropies in the Cosmic Microwave
Background <cit.>.  In a LCDM universe there is a larger
number of low mass structures compared to more massive ones,  or in
other words a steeply declining halo mass function N(M) ∝
M^-1 or halo velocity function N(V)∝ V^-3
<cit.>.  This prediction is naively at odds with
observational data around galaxies, i.e. the satellite abundance
<cit.>.  This “missing satellite” problem has
natural baryonic solution: galaxy formation becomes increasingly
inefficient in low-mass dark matter haloes, due to the ionizing
background, supernovae (SN) explosions and gas removal due to ram
pressure <cit.>.  The latest
cosmological hydrodynamical simulations of the Local Group are now
consistent with the observed satellite stellar mass and velocity
functions <cit.>.

The velocity function of galaxies using HI linewidths has been
measured in the nearby universe (within ∼ 200 Mpc) from the
HIPASS <cit.> and ALFALFA surveys <cit.>, and
within the Local Volume (within 10 Mpc) by <cit.>. These
and other authors <cit.>  have shown that
LCDM provides very good estimates of the number of galaxies with
circular velocities around and above 80.  With a careful
treatment of various systematics including survey selection effects
<cit.> showed agreement between LCDM down to circular
velocities of ∼ 60.  However, LCDM appears to fail  quite
dramatically at lower circular velocities, overestimating by a factor
of ≃ 3 the number of dwarf galaxies at velocity scale of
50, and by a factor of ≃ 5 at velocity scale of
30 <cit.>.   Galaxies at these velocity scales are
essentially insensitive to the ionization background and, by not being
satellites, they are not affected by gas depletion via ram pressure or
by stellar stripping. This makes the mismatch between the observed
linewidth function and the halo velocity function a more serious
challenge to the LCDM paradigm. Furthermore, simple modifications to
the LCDM paradigm, such as warm dark matter, are unable to resolve
this issue as candidates with a warm enough particle to significantly
reduce the number density of dark matter haloes result in a
characteristic scale in the velocity function that is not observed
<cit.>. The normalization of the velocity function is
dependent on key cosmological parameters such as the amplitude of
primordial perturbations (σ_8) and the total matter density
() <cit.>. For example, the difference between a
Planck <cit.> and WMAP7 <cit.> cosmology is a factor
of 1.43  in number density at fixed velocity <cit.>. 


The solution to this discrepancy is that in dwarf galaxies HI
linewidths measured at 50% of the peak flux, W_50, are a poor
tracer of the maximum circular velocity of the host dark matter halo
<cit.>.  This should not be a
surprise, since it is generally known that lower mass and younger
galaxies are less rotationally supported than higher mass, and older
galaxies. This is seen in both observations <cit.>, and hydrodynamical galaxy formation
simulations <cit.>.  Furthermore, as
discussed in <cit.>, unlike resolved rotation curves,
linewidths do not necessarily sample outer disks effectively since
the surface density drops rapidly beyond the optical radius.


In this paper we investigate the cause of the difference between
W_50/2 and  further using galaxy formation simulations
from the NIHAO project <cit.>.  Possible explanations we
consider include: gas mass loss from the halo – resulting in lower
circular velocity in the hydro than dark matter only simulation; HI in
the simulations is not extended enough – so that the maximum rotation
velocity does not reach the maximum circular velocity; Non-circular
motions – resulting in the rotation curve underestimating the
circular velocity; Projection effects – due to dwarf galaxies not
being thin rotating disks.

The outline of the paper is as follows. Section <ref>
describes the cosmological hydrodynamical simulations and the various
velocity definitions we use. In Section <ref> we present the
transformation from maximum circular velocity to HI linewidth. In
section <ref> we present the implied HI velocity function from
our simulations. In section <ref> we discuss what causes the
linewidths to underestimate the maximum halo velocity, and present
observational tests of these effects. A summary in given in section
<ref>.











§ SIMULATIONS



Here we give a brief overview of the NIHAO simulations. We refer the
reader to <cit.> for a more complete discussion.  NIHAO  is a
sample of ∼ 90 hydrodynamical cosmological zoom-in simulations
using the SPH code gasoline2 <cit.>.

Haloes are selected at redshift z=0 from parent dissipationless
simulations of size 60, 20, & 15 Mpc/h, presented in
<cit.> which adopt a flat ΛCDM cosmology with
parameters from the <cit.>: Hubble parameter H_0= 67.1
Mpc^-1, matter density =0.3175, dark energy density
Ω_Λ=1-=0.6825, baryon density =0.0490,
power spectrum normalization σ_8 = 0.8344, power spectrum slope
n=0.9624.  Haloes are selected uniformly in log halo mass from ∼
10 to ∼ 12 without reference to the halo merger history,
concentration or spin parameter.  Star formation and feedback is
implemented as described in <cit.>.  Mass and
force softening are chosen to resolve the mass profile at 1%
the virial radius, which results in ∼ 10^6 dark matter particles
inside the virial radius of all haloes at z=0. The motivation of
this choice is to ensure that the simulations resolve the galaxy
dynamics on the scale of the half-light radii, which are typically
∼ 1.5% of the virial radius <cit.>. Hydro particles
have force softenings a factor of 2.34 smaller than the dark matter
particles, and range from ≃ 75  pc in the lowest mass
haloes to ≃ 400  pc in the most massive haloes.

Each hydro simulation has a corresponding dark matter only (DMO)
simulation of the same resolution. These simulations have been started
using the identical initial conditions, replacing baryonic particles
with dark matter particles.  We remove the four  most massive haloes
(g1.12e12, g1.77e12, g1.92e12, g2.79e12), as these have formed too
many stars, in particular near the galaxy centers, resulting in
strongly peaked central circular velocity profiles <cit.>.
The final sample used in this work consists of 85 simulations.

Haloes in NIHAO zoom-in simulations were identified using the
MPI+OpenMP hybrid halo finder
[http://popia.ft.uam.es/AMIGA] <cit.>.  locates local over-densities in an
adaptively smoothed density field as prospective halo centers. The
virial masses of the haloes are defined as the masses within a sphere
whose average density is 200 times the cosmic critical matter density,
=3H_0^2/8π G.  The virial mass, size and circular velocity
of the hydro simulations are denoted: M_200, R_200, V_200.
The corresponding properties for the dark matter only simulations are
denoted with a superscript, DMO.  For the baryons we calculate
masses enclosed within spheres of radius r_ gal=0.2R_200,
which corresponds to ∼ 10 to ∼ 50 kpc.  The stellar mass
inside r_ gal is M_ star, the neutral hydrogen, inside
r_ gal is computed following <cit.> as described in
<cit.>. In principle the neutral hydrogen should be
separated further into atomic (Hi), and molecular, (H_2),
gas. In practice this has only marginal impact on the derived linewidths. This is expected because in dwarf galaxies the ratio
between molecular and atomic gas is very low, while in high mass
galaxies both atomic and molecular gas trace the same flat rotation
velocity, albeit from gas at different radii.  As a fiducial choice we
consider atomic gas to be neutral gas with a density less than
10  cm^-3 (which is also the star formation threshold used
in our simulations).  We measure galaxy velocities using a number of
definitions as discussed below.


The NIHAO simulations are the largest set of cosmological zoom-ins
covering the halo mass range 10^10 to 10^12. Their
uniqueness is in the combination of high spatial resolution coupled to
a statistical sample of haloes. In the context of LCDM they form the
“right” amount of stars both today and at earlier times
<cit.>. Their cold gas masses and sizes are consistent with
observations <cit.>, they follow the gas,
stellar, and baryonic Tully-Fisher relations <cit.>, and
they reproduce the diversity of dwarf galaxy rotation curve shapes
<cit.>.  As such they provide a good template with
which to connect galaxy observables (such as HI linewidths) to
intrinsic properties of the host dark matter halo (such as maximum
circular velocity).   









 §.§ Velocity definitions

We use a number of velocity definitions in this paper. We define them
with the help of some examples as shown in Fig. <ref>.  These
galaxies were chosen to represent both rotationally supported galaxies
(left panels) and those with more pressure support (right panels).


  
  *  – maximum circular velocity of the DMO simulation.
We start with the spherical circular velocity profile of the dark
matter only simulation (solid black line), which is equivalent to the
cumulative mass profile  M(<r) = r  V(r)^2/G.  The maximum spherical
circular velocity is shown as a black circle,
and occurs at 15% to 25% of the virial radius, R_200 (upper axis scale).


  * V_ max – maximum circular velocity of hydro simulation.
We next consider the circular velocity profile of the hydrodynamical
simulation (blue lines). The solid line shows the spherical circular
velocity, while the dashed line shows that derived from the potential
in the disk plane, V^2_ pot=-R ∂Φ/∂ R.
We define ≡ V_ circ(R_
  max^ DMO) as the spherical circular velocity of the hydro
simulation at the radius where the maximum circular velocity of the
DMO simulation occurs. We use this definition as it explicitly shows
when there is mass loss from the system (<), and where
there is baryon dissipation (>). In these examples,
there is mass loss from the galaxies in the right panels in Fig. <ref>,
and no change for the galaxies in the left panels.

Notice that the potential based circular velocity is lower than the
spherical circular velocity at small radii, and higher at large
radii. This is a well known feature of the differences between a
flattened and spherical mass distribution <cit.>






  * V_ HI - circular velocity at HI radius.  We define
  the HI radius, R_ HI, as enclosing 90% of the HI mass (in
  the face-on projection). This approximates the edge of the
  observable HI rotation curve, and occurs close to the radius where
  the projected HI surface density profile reaches 1 
    pc^-2, but it has the advantage of being free from projection
  effects and is measurable in every galaxy (that contains HI).  The
  potential based circular velocity at the HI radius is then, V_
    HI≡ V_ circ(R_ HI).


  * V_ rot(R), σ(R) – rotation and dispersion
  profile.  We calculate the rotation curve (red line), and
  dispersion profile (pink line) as follows. We rotate the galaxy to a
  face-on view, using the angular momentum of the cold (T<15000 K)
  gas particles inside 10% of the virial radius. We then divide the
  galaxy into a series of annuli of equal width. In each ring we
  calculate the mean velocity (V_ rot) and velocity dispersion
  (σ) in the azimuthal (ϕ) direction. Note that this
  dispersion only takes into account the bulk motions of the gas.  In
  the galaxies in the left panels, the rotation velocity closely traces the
  potential based circular velocity, and the dispersion is low, ∼
  5-10. In the galaxies in the right panels, the rotation velocity
  systematically underpredicts the circular velocity, and the velocity
  dispersion is higher, ∼ 10-20.


  * W_50, W_20 - HI linewidths.  We next consider the HI
  linewidths, measured at 50% and 20% of the peak flux (see
  Fig. <ref> for examples from our reference galaxies).  In an
  update to the calculation in <cit.> here we include the
  thermal motions of the HI gas. At large linewidths the changes are
  minor, but for small linewidths the temperature adds a floor to the
  linewidth of ∼ 20, corresponding to σ∼ 8.5.  For each gas particle we represent the line-of-sight
  velocity distribution with a Gaussian with mean corresponding to the
  line-of-sight velocity of the particle, a flux corresponding to the
  HI mass of the particle, and a dispersion computed from the
  temperature of the particle assuming σ=√(k_ BT/m_
      H)=9.09 [] √(T/10000  K).  We sum up the
  Gaussians to give the flux vs velocity histogram.
  
  As is commonly done in observations, we define the peak flux as the
  average of the peak flux on each velocity wing (i.e., positive and
  negative). See Fig. <ref> for examples. The peak fluxes on
  each velocity wing are marked with dotted horizontal lines. The
  average of these is shown with solid horizontal lines.  By
  definition W_20≥ W_50.  Because W_20 requires higher
  signal to noise, and thus can be measured from fewer galaxies,
  W_50 is the more commonly used definition in large HI surveys,
  so will be our default definition.

  We calculate the linewidths from 100 random projections (uniformly
  selected from the surface of a unit sphere), and measure the
  maximum, median, and minimum values. The maximum values typically
  correspond to the values derived from an edge-on view of the
  galaxy. They are shown as cyan (W_20^ max) and green
  (W_50^ max) horizontal lines in Fig. <ref>.  The
  median W_50/2 is shown with an orange horizontal line. The
  minimum linewidth (purple horizontal lines) can be used to
  approximate the gas velocity dispersion of the system, since for a
  Gaussian, σ = W_50/2.35. 




For each projection where we measure a linewidth we also measure the minor
to major axis ratio (b/a) of the HI gas using moments as follows.  For
each gas particle we have its projected coordinates, (x, y),
distance from the galaxy center R^2=x^2+y^2, and HI mass, m. Using
these quantities we calculate the moments weighted by the HI mass:


    Sxx=(∑ m x^2/R^2 ) / ∑ m,
    
      Syy=(∑ m y^2/R^2 ) / ∑ m,
    
      Sxy=(∑ m xy/R^2 ) / ∑ m.

The axis ratio is then given by

    (b/a)=1-√(Q^2 + U^2)/1+√(Q^2 + U^2),

where Q=Sxx-Syy and U=2Sxy.













§ FROM MAXIMUM CIRCULAR VELOCITY TO HI LINEWIDTH



Fig. <ref> shows the relations between half-linewidths
W_50/2 and W_20/2 with .  The left panel is very
similar to Fig.2 of <cit.>, but here we use an updated
calculation of the linewidth which introduces a floor to the linewidth
of W_50∼ 20.  Galaxies lie systematically below the
one-to-one relation (dotted lines) indicating the linewidth is a
biased tracer of the maximum halo velocity.  Linear fits to the
relations (using all the data points) are given by

    log(W_50/2)= 1.614 +1.557(log -1.904),

with a scatter of σ=0.16, and

    log(W_20/2)= 1.751 +1.419(log -1.904),

with a scatter of σ=0.12.  The slope of these fits is greater
than unity, showing that linewidths are a poorer tracer of halo
velocity in lower velocity haloes. The scatters show that W_50 is
a poorer tracer of the halo velocity than W_20.













Fig. <ref> shows the relations between various velocity
definitions and . Linear fits (in log-log space) are
given in Table. <ref>.  Notice that the normalization,
slope, and scatter systematically change at each step. At anyone of
these steps the changes are not large, but they accumulate to a
relation between linewidth and halo velocity that is very different
from a one-to-one correspondence.  For galaxies in Milky Way mass
haloes ≃ 180 the average projected half-linewidth
is W_50/2=141, so that up to the expected project effects, the
linewidth is a good tracer of the maximum halo velocity. However, for
dwarf galaxies at ≃ 50 the average
half-linewdith of  W_50^ med/2 ≃ 20, is a factor of 2.5
times lower than the halo velocity.

There are a number of physical processes that play a role in the
conversion of maximum circular velocity into projected HI linewidths:
inflows and outflows; HI extent and rotation curve shape; non-circular
motions, and projection. These are related, but not uniquely, to  four
ratios, which we show in Fig. <ref>. 



 §.§ Inflows and Outflows

The upper left panel shows /, which depends on the net
inflow and outflow of baryons.  At low velocities the ratio ≃
0.85 indicating net outflows, while at high velocities the ratio is
greater than 1, indicating net inflow. As a reference, if a system
lost or did not accrete any of its baryons we would naively expect the
ratio of √(()1-)≃ 0.92, where the cosmic baryon
fraction, =Ω_ b/Ω_ m≃ 0.15.  However,
the reduction in halo velocity is larger than this because the loss of
baryons at early times (as a result of the stellar feedback) reduces
the accretion rate of dark matter, and hence the halo mass and maximum
velocity of the halo by redshift z=0.  Similar results are seem in
other cosmological hydrodynamical simulations,
<cit.>, so the effect appears unrelated to the
details of the sub-grid model for star formation and feedback.



 §.§ HI extent and rotation curve shape

The upper right panel shows the / ratio vs .  In
all galaxies  is measured at a smaller radius than , so
the ratio depends on the shape of the circular velocity profile, and
the extent of the HI gas. In high velocity haloes the ratio is greater
than 1 because the inflow of baryons. In low mass haloes the ratio is
less than 1 because the HI does not extend to the flat part of the
velocity profile. At a halo velocity of 50 km/s the average ratio is
0.8, but there are some dwarf galaxies with a ratio of unity and
others with a ratio of 0.6.  This is already pointing to diversity of
HI extent and/or rotation curve shapes.



 §.§ Non-circular motions

The lower left panel shows the W_50^ max/2 ratio. The
linewidth is a convolution of the rotation curve with the HI
distribution. In galaxies that are rotationally supported with a flat
rotation curve this ratio is close to unity. When the ratio is less
than unity this signals a rising rotation curve and/or significant
non-circular motions.  These can be dispersion in the gas, but also
non-axisymmetric features such as bars, spiral arms, and warps, and
out of equilibrium gas flows due to e.g., supernova driven winds or
mergers.  At a halo velocity of 50 km/s the average ratio is 0.7, but
again there is significant scatter, with some dwarf galaxies with a
ratio close to unity.



 §.§ Projection

The lower right panel shows the ratio between the median and maximum
W_50 linewidth.  This ratio tells us about projection effects. The
inclination angle, i, where i=0 is face-on, and i=90 is edge-on,
is distributed uniformly in cos(i) from 0 to 1. For a thin rotating
disk with negligible velocity dispersion, the linewidth at
inclination, i, is related to the linewidth at inclination i=90 by
W_i=sin(i)W_90. The ratio between the median linewidth (i=60)
and the maximum (i=90) is thus sin(60)=0.866.  In the figure we
see a few galaxies at all halo velocities near this value (dotted
line). However, the majority of galaxies lie below this line
indicating the galaxies are in general not thin rotating disks.  This
is another manifestation of the role of non-circular motions.  At a
halo velocity of 50  the average ratio is 0.8.  At the lowest
velocities the ratio is close to unity, indicating these systems (or
at least the HI lines) have very little rotational support.









§ IMPLIED VELOCITY FUNCTIONS


We now discuss how the various velocity definitions result in velocity
functions with different slopes.  First we show results with
power-laws, as the transformations are simple and analytic. Then we
consider the impact of Gaussian scatter, and finally the actual
distribution of linewidths from our simulations corresponding to
non-power-laws and non-Gaussian scatter.



 §.§ Analytic considerations for the velocity function

For a cumulative velocity function with a power-law form

    N(>V) = AV^α,

the differential velocity function has the same slope,
α, but with a normalization that is different by the slope:

    dN/ dln V= α  N(>V).

If two different velocities are related by
V_2=a V_1^b, then the cumulative velocity function is simply
translated:

    N(>V_2) = A (V_2/a)^α/b.

The differential velocity function is changed due to the dln V, and since

    dln V_1 = (1/b)  dln V_2,

we have the new differential velocity function:

    dN/dln V_2 = (α/b)  N(>V_2) = (α/b)  A (V_2/a)^α/b


The differential velocity function for  for Cold Dark Matter
for a Planck cosmology <cit.> is given by <cit.>:

    dN/dlog_10V = 0.186 (/100)^-2.9.

We can thus use the equations in Table <ref> to translate
this into the velocity functions of the various velocity
definitions. We first do this ignoring the scatter to show the
magnitude of the effect, later we will include the scatter.  The
result is shown in Fig. <ref>, and Table <ref> gives
the slopes and normalizations at V=50. The majority of the
differences come from four conversions: /≃ 0.63,
/≃ 0.53, W_50^ max/2≃ 0.54, and
W_50^ med/W_50^ max≃ 0.67.  The overall shift in
normalization is a factor of 8.7. We stress that this is entirely due
to a change in the velocities, rather than any change in the number
densities of the objects. We note that the strong impact of the
velocity definition on the slope of the velocity function has been
shown previously by  <cit.> using halo abundance matching to
link baryonic masses to halo masses, and then various Baryonic
Tully-Fisher relations to relate baryonic masses to velocities. 




















 §.§ Impact of Gaussian scatter

Fig. <ref> shows the impact of Gaussian (in log(V)) scatter
on the velocity function.  The black line shows the velocity function
using , sampled between a velocity of 30 to 200 . The
red dotted line shows the velocity function using W_50/2. The other
lines show the impact of log-normal scatter (0.1, 0.2, and 0.3 dex) in
the relation between W_50/2 and .  The effect of scatter
is to increase the normalization of the velocity function.  Galaxies
get preferentially scattered to higher velocity because there are many
more low velocity haloes than high velocity haloes. The overall effect
is not large, just 0.1 dex change for a scatter of 0.2 dex.








 §.§ Deriving the velocity function from linewidth vs maximum circular velocity


We now go into more detail and show how the distribution of galaxies
in the W_50/2 vs  plane can be converted into the
velocity function.  We set up a grid in x=log() and y=log(W_50/2) of width dx and
dy and count the number of galaxies in each cell, n(x,y).  From
the cumulative DMO velocity function we straightforwardly know the
number of dark matter haloes in a given x-bin in a given volume of
space, dN_ CDM(x). Comparing the actual number of data points
N(x) (number of haloes in each x-bin multiplied by number of
projections per halo) to dN_ CDM(x) we get the weight w(x)
for each halo in bin x:

    w(x)=dN_ CDM(x)/N(x)

We then go to each y-bin and count up the
number of galaxies using the weights:

    N(y)=Σ_x n(x,y) w(x)/Σ_x w(x).


As fiducial bin widths we adopt dx=0.05 and dy=0.10.  The
resulting velocity functions for W_50 (left) and W_20 (right)
are shown with red points in Fig. <ref>, and tabulated in
Table <ref>. The error bars (e_y) are calculated as
e_y=y/√(N), where N is the number of simulations that
contribute to each y-bin.  At low velocities (W_50/2 < 80)
the NIHAO W_50/2 function has a normalization and shallow slope
∼-1 in agreement with observations from the Local Volume
<cit.> and ALFALFA <cit.>.  At high
velocities (100) the NIHAO simulations underpredict the
observed number densities, or alternatively the linewidths are too
low.  We return to a possible cause of this in section
 <ref> below.








In our simulations W_20/2 is a better predictor of the halo
velocity than W_50/2 (See Fig. <ref>) so this is our
preferred definition for future observations.  Currently observations
of the W_20 function are not available, so we shift the
observations using the approximation W_20=W_50+25
<cit.>, and a minimum of W_50/W_20=0.63.
In this case the simulations are closer to 
the observed velocity function at high velocities, while
maintaining the agreement at low velocities.








§ TESTING THE SIMULATIONS


We have shown that the discrepancy at low velocities between the
observed HI linewidth function and the  function of LCDM is
resolved by the NIHAO simulations. We now discuss in more detail why
this is the case, and observational ways in which the simulations can
be tested. We focus our discussion on six (mostly dimensionless)
parameters that are correlated with the variation in HI linewidth at
fixed maximum dark halo velocity.



  * Rotation-to-dispersion ratio, V_ rot/σ_ HI. The
  HI dispersion is calculated from the minimum HI linewidth,
  σ_ HI=W_50^ min/2.35.  For the rotation velocity
  we use the maximum of the rotation curve, which tends to occur near
  the HI radius (see Fig. <ref> for examples). This definition is not
  directly observable, but similar definitions are.

  * γ_1, Skewness of the W_50 distribution, see below for
  definition.  Can be measured on samples of galaxies.

  * HI line profile shape, W_50/W_20. Directly observable.

  * Outer circular velocity curve slope, Δlog V /Δlog R. For our simulations is measured between 0.5 and 1.0 R_
    HI. Can be measured for individual galaxies. 

  * HI-to-virial radius, R_ HI/R_200. Dimensionless, but
  not directly observable. HI sizes are observable.

  * HI disk thickness. Computed as the minimum HI axis ratio from
  all projections. Can be measured from samples of galaxies.









 §.§ Testing rotational support with projection effects



One of the key results from our simulations is that galaxies, and
dwarf galaxies in particular, do not always have HI kinematics as
expected from idealized thin rotating disks. An observationally
testable consequence of this is the effect of projection.

For a thin rotationally supported axisymmetric disk the minor-to-major
axis ratio (b/a) is uniformly distributed between 0 and 1. In terms
of the disk inclination, i, where i=0 is face-on and i=90 is
edge-on, b/a=cos(i), the observed linewidth then varies as
sin(i). Fig. <ref> shows the linewidth vs axis ratio
for our four test galaxies. The axis ratio is computed for each
projection using moments as described in section 2.  The upper
galaxies have maximum halo velocity ≃ 55, while
the lower two galaxies have ∼ 100. The corresponding
stellar masses are ∼ 2× 10^8 and ∼ 10^9,
respectively.

The two galaxies on the left panels of Fig. <ref> have
close to uniform axis ratio distributions (upper histograms) and
strongly skewed linewidth distributions (right histograms) close to
that predicted for thin disks rotating at the maximum velocity of the
dark matter halo (solid red lines).  On the other hand, the two
galaxies in the right panels of Fig. <ref> strongly
deviate from this prediction with a lower mean linewidths and a more
symmetric distribution of linewidths. The axis ratio distribution is
also more centrally concentrated.


A parameter that is related to how rotationally supported a
galaxy is, is the skewness of the linewidth distribution. Here we
adopt Pearson's moment coefficient of skewness, which 
is defined as

    γ_1 = μ_3 / σ^3.

Here μ_3=Σ (W-W)^3 / N is the third central
moment, σ is the standard deviation, and W is the
mean. A value of γ_1=0 corresponds to a symmetric distribution,
γ_1<0 corresponds to a tail to low values, while γ_1>0
corresponds to a tail to high values. A thin rotating disk has a
strongly negatively skewed distribution of linewidths with
γ_1=-1.13.

The galaxies in the left panels of Fig. <ref>  have
γ_1=-0.97 and -1.06, i.e., they behave like thin rotating
disks. In contrast the galaxies on the right have γ_1=-0.10,
and γ_1=-0.73, and are thus less rotationally supported
systems.  These inferences agree with the actual rotation to
dispersion ratios (V_ rot/σ_ HI). 

Another  parameter that is related to the rotational support of a
galaxy is minimum minor-to-major axis ratio, q_ min (i.e., the
disk thickness). The galaxies on the left have smaller minimum axis
ratios (corresponding to thinner edge-on disks): q_ min=0.11
and 0.06 compared to 0.28 and 0.25.

The ratio between linewidths W_50/W_20, also depends on the
rotational support, as well as the circular velocity curve shape, and
the extent of the HI.  The galaxies on the left have steeper line
profiles, and more extended and shallower outer circular velocity
curve slopes.



 §.§ What causes linewidths to underestimate maximum circular velocity?


Fig. <ref> shows the impact of several (mostly
dimensionless) parameters on the relation between median HI linewidth
and maximum dark halo circular velocity. The points are color coded by
the parameter in question.  The four numbers in each panel indicate
the mean value of each quartile, the solid line is a fit to each
quartile of points. We discuss each of them in turn below.



  * Outer circular velocity curve slope (top left). A positive value
  indicates that the curve is still rising in the outer parts, while a
  value close to zero indicates a flat velocity curve.  Galaxies with
  more rising circular velocity profiles (red points) tend to have
  lower linewidth-to-halo velocity ratios. 


  * HI-to-virial radius (top middle). Smaller galaxies (blue
  triangles) have lower linewidths.  At least partially because if the
  HI does not reach the flat part of the circular velocity curve, then
  the kinematics will progressively underestimate the maximum circular
  velocity. 


  * HI axis ratio (top right). The thinnest galaxies
  (red points, b/a∼ 0.1) have linewidths that trace the halo
  velocity, while thicker galaxies have progressively lower
  linewidths.
  

  * γ_1, Skewness of the W_50 distribution (top
  left). Galaxies with strong negative skew (blue points) have
  linewidths that trace the maximum halo velocity. Larger values of
  γ_1 result in lower linewidths.
  

  * Rotation-to-dispersion ratio (bottom middle).  Rotationally
  supported systems (red circles) tend to have linewidths that trace
  the halo velocity, while pressure supported systems (blue triangles)
  tend to have low linewidths.


  * HI line profile shape (bottom left). Steeper line profiles
  W_50/W_20∼ 0.9 (red points) have linewidths close to halo
  velocity, while shallow line profiles W_50/W_20∼ 0.6 (blue
  triangles) have linewidths significantly below the halo velocity. The
  profile shape depends on a combination of the rotation curve shape and
  extent, and the rotation to dispersion ratio.




The number in the top left indicates the correlation coefficient
between the specific parameter and W_50^ med/2.
Line profile shape shows the strongest correlation (ρ=0.89)
closely followed by the skewness (ρ=-0.83), and HI rotation to
dispersion ratio  (ρ=0.79).  At a given halo velocity all six
parameter show a correlation with linewidth. To see which of these
parameters is most likely to explain the trend with halo mass we  show
the dependence on halo velocity in Fig. <ref>. The straight
and dashed lines show a fit with 1σ scatter. The number gives
the correlation coefficient. We see clear trends but also significant
scatter.  Structurally, we see that lower velocity haloes have steeper
outer circular velocity curve slopes and relatively smaller HI sizes,
and marginally thicker HI disks.  Kinematically, we see that lower
velocity haloes have less negatively skewed HI linewidth
distributions, lower rotation-to-dispersion ratios, and shallower line
profiles.

In summary we see that the dependence of W_50 with  is
driven primarily by two effects: The degree of rotational support, and
the shape of the rotation curve. Galaxies in lower mass halos have
less rotational support, and less extended HI. Below we show current
observations that show evidence of these two effects.










 §.§ Line profile shape


The shape of the HI line profile depends on the amount of rotational
support and the shape of the rotation curve.  A simple way to
parameterize this is the ratio between  linewidths measured at the
50% and 20% level of peak flux, W_50 and W_20. The ratio is
close to 1 for a rotating disk with a flat rotation curve whose
velocity dispersion is small compared to the line-of-sight component
of the rotation velocity.  As the outer rotation curve slope becomes
more positive the ratio decreases, since a larger fraction of the HI
flux is coming from gas rotating lower than the maximum. Pressure
support also causes ratio to decrease. In the limit V/σ=0 the
line profile is simply that due to random motions of the gas, and is
independent of the rotation curve shape.  For a Gaussian
W_50/W_20=√(ln(0.5)/ln(0.2))≃ 0.66.  The
observational advantage of this ratio is it does not require spatially
resolved HI observations, which are currently not feasible for large
samples of galaxies.

A comparison between NIHAO simulations and observations was previously
shown in <cit.>. In both observations and simulations the
W_50/W_20 decreases for lower linewidth galaxies. This result
was also shown by <cit.> using a smaller sample of galaxies
from the MaGICC project (the precursor to NIHAO), and by
<cit.> using galaxies from the FIRE project.

Fig. <ref> shows the relation between W_50 and W_20
from NIHAO simulations compared to observations from HIPASS
<cit.> and THINGS <cit.>.  In both simulations
and observations the ratio varies from  ∼ 0.95 at high velocities
to ∼ 0.65 at low velocities, and the amount of scatter is also
similar.  However, when looking at the median relations (lower right
panel) the NIHAO simulations underpredict the HIPASS observations for
W_50∼ 100. This is the same scale where the simulations
underpredict the observed velocity function (Fig. <ref>).  The
upper right panel shows the distributions of W_50/W_20 in five
bins of W_50/2. Again this shows an excess of galaxies in NIHAO
with low linewidth ratios at intermediate velocities. In principle
this could be due to there being too much non-circular motions and/or
the HI disks being too small at this velocity scale. Since the NIHAO
simulations are in good agreement with observations of HI sizes (see
below), we think the low linewidth ratio is more likely due to too an
excess of non-circular motions, likely driven by feedback. Due to the
large intrinsic variation in W_50/W_20 and the relatively small
sample size, another contributing cause could be sampling effects,
i.e., we by chance happen to simulate more galaxies with lower
ratios. A similar discrepancy in W_50/W_20 at this velocity
scale is seen in the FIRE simulations <cit.> and likely
has a common origin, which points to feedback driven turbulence.

In summary, the line profile shape is a powerful, yet simple test
of the simulations, and thus motivates future blind surveys
obtaining deep enough data to accurately measure both W_50 and
W_20.













 §.§ Skewness of the linewidth distribution

For observed galaxies we only get a single projection per
galaxy. However, for a sample of galaxies, we can expect a random
distribution of projections, provided we select galaxies by a
parameter that is independent of the projection angle, such as the HI
mass.

Fig. <ref> shows the distributions of HI linewidth vs HI
mass.  Single dish observations are from the ALFALFA survey
<cit.> and HIPASS <cit.>.  The lower right
panel shows the dependence of γ_1 measured in bins of HI mass.
The two observational datasets and simulations find consistent
results, for the skewness of both W_50 (solid lines) and W_20
(dashed lines).  Namely high mass galaxies are strongly negatively
skewed, close to the value expected for thin rotating disks.  Lower
mass galaxies have higher γ_1.  This is an observational
signature that lower mass galaxies are more kinematically disordered.
For the lowest mass galaxies the distribution of linewidths is
positively skewed, which is a signature of a floor in the linewidth
from thermal broadening of the HI line.

The distribution of linewidths in bins of HI mass can also be used to
estimate the velocity dispersion of the HI gas, σ_
  HI. This is because for a galaxy disk that is face-on, there will
be zero rotational component projected into the line-of-sight, and
thus the observed linewidth will be a reflection of the disordered
motions of the gas. Assuming the line profiles are Gaussian
σ_ HI=W_50/2.355.  The minimum linewidth is subject to
observational measurement errors, and may not be representative of the
population, so we take the average of the lowest 5% of linewidths.
This will give an upper limit to σ_ HI because for a thin
rotating disk the lowest 5% of inclinations result in a line-of-sight
rotation that is ∼23% of the edge-on value. 

The results of this exercise are shown with the solid points in
Fig. <ref>. Overall there is good agreement between the
observations and simulations.  In detail there are some differences.
Comparing the observations, in the low mass bins HIPASS results in
∼ 3 higher σ_ HI than ALFALFA, likely due to the
lower spectral resolution of HIPASS (indicated with horizontal
lines). Except for the highest mass bin, the NIHAO simulations have
slightly higher σ_ HI than ALFALFA. For NIHAO  the red
open squares show the thermal velocity dispersion.  The difference
between the solid red and open red points is thus the contribution of
turbulence (and projected rotation).  Since projection affects both
simulations and observations the slightly higher dispersions in NIHAO
are thus caused by turbulence, likely as a result of the feedback
being too strong.  

The black and grey points show measurements from individual galaxies
using spatially resolved observtions from THINGS <cit.>
and LITTLE THINGS <cit.>. The solid lines show linear fits
with 3σ clipping. The slopes are similar, but there is a small
offset of ≃ 2, likely reflecting a systematic in the
measurement techniques. In the mass range 10^7 < M_ HI  < 10^9 the resolved measurements are in good agreement with the
spatially unresolved linewidth based measurements. Thus in the future
our technique can be used to infer the velocity dispersions of large
samples of spatially unresolved observations.  At higher masses
(M_ HI> 10^9) the resolved measurements give lower
dispersions than the unresolved measurements. As mentioned above,
projected rotation likely contributes to this difference. In addition,
the linewidth based dispersion includes deviations from a perfectly
axisymmetric disk. For example, a warp will prevent the projected
rotation from being zero at any viewing angle, increasing the implied
dispersion.  Another possible contributing effect is the role of
sample selection. The THINGS survey is not a complete sample of
galaxies, it is possibly biased towards that more rotation dominated
galaxies at these masses than the galaxy population.

In summary, the distribution of linewidths in bins of mass shows
that in both simulations and observations lower mass galaxies are less
rotationally supported than more massive galaxies.  This is driven by
a reduction in the rotation velocity in low mass galaxies rather than
an increase in the gas dispersion, which also decreases in lower mass
galaxies.
  







 §.§ HI sizes

We showed previously in <cit.> that the NIHAO galaxies match
two HI size vs HI mass relations (the exponential scale length, and
the radius where the surface density reaches 1 
  pc^-2). We are interested in the relation between HI size and
halo mass, but halo masses are not directly observable.  In the NIHAO
simulations, stellar mass is the galaxy observable that is most
strongly correlated with halo mass (correlation coefficient of
0.98). The scatter in the stellar mass vs halo mass relation is small
(∼ 0.2 dex) both in theory <cit.>, and
observations <cit.>. Here we use galaxies
from the SPARC survey <cit.> as this is the largest
compilation of resolved HI rotation curves with 3.6 μ m Spitzer
imaging (which gives the most reliable stellar mass tracer).  

Fig. <ref> shows HI size vs stellar mass relations and a
comparison between various HI size definitions.  The upper panels show
that R_1 (the radius where the HI surface density drops to 1) is closely related to R_90 in the NIHAO simulations and
R_ last in the SPARC observations.  The lower left panel shows
a very good agreement between the R_1 vs stellar mass relations in
NIHAO (red) and SPARC (blue).  Fits to the size-mass relations are
given in Table <ref>.

While the observations are not a complete sample of nearby galaxies
the comparison is a good place to start.  It contradicts the claim by
<cit.> that the HI distribution in NIHAO galaxies is
too compact compared to observations. This was suggested by
<cit.> as a reason for the low
W_50/2 ratio in NIHAO dwarf galaxies. As shown in
(Fig. <ref>), the HI in NIHAO dwarf galaxies is less
extended relative to the dark matter halo than in more massive
galaxies. So this does contribute to the low W_50/2 ratios
in dwarfs. However, rather than being an artifact of 
NIHAO, this appears to be a feature of real galaxies.






 §.§ Outer HI rotation curve slopes


Fig. <ref> shows the logarithmic slope (Δlog V /
Δlog R) of the outer rotation curve vs the stellar mass (left)
and the HI radius (right). Points show individual galaxies, while
symbols with error bars show the mean and error on the mean in bins of
mass or radius. For the observations the slope is measured between 0.5
and 1.0 R_ last  on the rotation curve, for the simulations the
fiducial slope is measured between 0.5 and 1.0 R_ HI on the
circular velocity curve (red filled circles), and on the rotation
curve (pink open circles).

The simulations are in good agreement with the observations, except
that the observed relations have slightly larger scatter, plausibly
due to larger measurement errors.  We see that high mass or large
galaxies have, on average, flat outer rotation curves at the HI
radius. As we go to lower stellar masses and smaller sizes the outer
rotation curves are progressively rising.  Notice that even though
both observations and simulations include cored dark matter density
profiles in dwarf galaxies, at large radii the typical outer
rotation curve slope is ∼ 0.5 corresponding to a ρ∝
r^-1 density profile.  When the rotation curve is rising, the
half-linewidth W_50/2 will neccesarily underestimate the maximum
rotation velocity.





 §.§ Axis ratios

Fig. <ref> shows the relations between projected axis ratio
and HI mass. Each galaxy has 100 projections shown. For clarity we
have randomly shifted the HI masses by a small amount for each
projection. There is a weak trend, similar to that between minimum
axis ratio and halo velocity shown in Fig. <ref>, with a
scatter of 0.20.  The average projected axis ratio is about 0.6,
compared to 0.5 for an idealized thin disk.  The clearest testable
prediction is the minimum axis ratio is larger in lower mass
galaxies. The most massive galaxies can have disks as thin as
0.05. Below a mass of M_ HI∼ 10^9 the minimum axis
ratio steadily increases reaching 0.4 at M_ HI∼ 10^6
The lower envelope (dotted line in Fig. <ref>) is given by

    (b/a)_ min = 0.06  + 0.106( log_10(M_ HI) -9.4)

for log_10(M_ HI) ≤ 9.4, and   (b/a)_ min = 0.06 otherwise.

The cyan circles show observed galaxies from the nearby Universe from
<cit.>.  This study combines data from several different
projects so is representative of a wide variety of galaxy types and
environments: Atlas3D <cit.>, Bluedisk
<cit.>, Faint Irregular Galaxies GMRT Survey (FIGGS)
<cit.>.  LITTLE THINGS <cit.>, Local Volume HI
Survey (LVHIS) <cit.>,  starbursting dwarf galaxies
<cit.>, The HI Nearby Galaxy Survey (THINGS)
<cit.>, Ursa Major cluster <cit.>, Void Galaxy
Survey (VGS) <cit.>, VLA Imaging of Virgo Spirals in Atomic
Gas (VIVA) <cit.>, Westerbork HI survey of spiral and
irregular galaxies (WHISP) <cit.>.




For the vast majority of galaxies (88%) the axis ratio  is determined
by <cit.> from the maps, based on the second order moments
of the pixel distributions where Σ_ HI > 1. For
the galaxies in VIVA (36) and Atlas3D (8) axis ratios are from  tilted
ring fits to the velocity fields.  Galaxies that are poorly resolved
(HI radius less than the beam major axis) have been excluded, though
this does not guarantee the axis ratios are well resolved in all
galaxies.

The observations and simulations have similar distributions of axis
ratios, in particular there are no observed galaxies in the lower left
corner where our simulations predict there to be none.  This appears
to confirm the lack of dwarf galaxies with thin HI disks.  Future
observations of a complete sample of galaxies with higher spatial
resolution are needed to conclusively confirm the lack of thin dwarf
galaxies.




§ SUMMARY



We use a sample of 85 galaxies simulated in a LCDM cosmology from the
NIHAO project to investigate why HI linewidths systematically
underpredict the maximum dark halo circular velocities in dwarf
galaxies (Fig. <ref>).  We trace this to two primary effects.



  * Lower mass galaxies are less rotationally supported. This is
  confirmed observationally from the skewness of linewidths in bins of HI mass in both
  ALFALFA and HIPASS observations (Fig. <ref>).


  * The HI distributions are less extended (relative to the dark
  matter halo) in dwarf galaxies, so that the rotation curves are
  still rising at the last measured data point, in
  agreement with observations (Fig. <ref>). 


The HI profile shape parameterized by W_50/W_20 decreases in
lower mass galaxies (Fig. <ref>) consistent with both these two
effects. An additional observational test is in the distribution of HI
axis ratios. In particular, in the NIHAO simulations the minimum axis
ratio is larger in lower mass galaxies (Fig. <ref>).

In our simulations the HI kinematics are an inhomogeneous
population. There is a significant range of rotational support and HI
extent at any given halo or galaxy mass. This variation drives the
variation in W_50/2 (Fig. <ref>).  Thin and
extended rotating HI disks exist in our simulations at all halo
masses, but they are not a fair sample.  This implies that one cannot
use a sample of well ordered rotating disks to interpret the HI
linewidths of large unbiased samples of galaxies, as is sometimes done
<cit.>.

The implied linewidth velocity function from the NIHAO simulations has
a shallow slope (≃ -1) at low velocities (10 < W_50/2 <
80), consistent with observations (Fig. <ref>). Thus
we conclude that the apparent discrepancy between the
predictions of the LCDM cosmological model and observations as
highlighted by previous authors <cit.> is due
to their incorrect assumption that W_50/2 ≃. 

We look forward to the next generation of blind HI line surveys, APERTIF
<cit.> and ASKAP <cit.> that will provide
higher spatial resolution data with which to further
test our simulations and thus the LCDM model. 



§ ACKNOWLEDGMENTS

We thank the anonymous referee for a detailed and constructive report.
This research was carried out on the High Performance Computing
resources at New York University Abu Dhabi; on the  theo cluster
of the Max-Planck-Institut für Astronomie and on the hydra
clusters at the Rechenzentrum in Garching.

The authors gratefully acknowledge the Gauss Centre for Supercomputing
e.V. (www.gauss-centre.eu) for funding this project by providing
computing time on the GCS Supercomputer SuperMUC at Leibniz
Supercomputing Centre (www.lrz.de).
AO is funded by the Deutsche Forschungsgemeinschaft (DFG, German  
Research Foundation) – MO 2979/1-1.






99


[Begum et al.(2008)]Begum08 Begum, A., Chengalur, J. N., Karachentsev, I. D., Sharina, M. E., & Kaisin, S. S. 2008, , 386, 1667 
  

[Benson et al.(2002)]Benson02 Benson, A. J., Frenk, C. S., Lacey, C. G., Baugh, C. M., & Cole, S. 2002, , 333, 177 
  

[Binney & Tremaine(1987)]Binney87 Binney, J., & Tremaine, S. 1987, Princeton, NJ, Princeton University Press, 1987, 747 p.,  


[Brook & Di Cintio(2015)]Brook15 Brook, C. B., & Di Cintio, A. 2015, , 453, 2133 

 
[Brook & Shankar(2016)]Brook16a Brook, C. B., & Shankar, F. 2016, , 455, 3841 


[Brook et al.(2016)]Brook16b Brook, C. B., Santos-Santos, I., & Stinson, G. 2016, , 459, 638  
  

[Brooks et al.(2017)]Brooks17 Brooks, A. M., Papastergis, E., Christensen, C. R., et al. 2017, , 850, 97 


[Buck et al.(2018)]Buck18 Buck, T., Macciò, A. V., Dutton, A. A., Obreja, A., & Frings, J. 2018, MNRAS in press, arXiv:1804.04667 

 
[Bullock et al.(2000)]Bullock00 Bullock, J. S., Kravtsov, A. V., & Weinberg, D. H. 2000, , 539, 517 


[Courteau(1997)]Courteau97 Courteau, S. 1997, , 114, 2402 


[Chung et al.(2009)]Chung09 Chung, A., van Gorkom, J. H., Kenney, J. D. P., Crowl, H., & Vollmer, B. 2009, , 138, 1741 


[Ceverino et al.(2017)]Ceverino17 Ceverino, D., Primack, J., Dekel, A., & Kassin, S. A. 2017, , 467, 2664 


[Dutton & Macciò(2014)]Dutton14 Dutton,
  A. A., & Macciò, A. V. 2014, , 441, 3359


[Dutton et al.(2017)]Dutton17 Dutton, A. A., Obreja, A., Wang, L., et al. 2017, , 467, 4937 


[El-Badry et al.(2018)]El-Badry18a El-Badry, K., Quataert, E., Wetzel, A., et al. 2018a, , 473, 1930 
    
  
[El-Badry et al.(2018)]El-Badry18b El-Badry, K., Bradford, J., Quataert, E., et al. 2018b, , 477, 1536 


[Gill et al.(2004)]Gill04 Gill, S. P. D.,
  Knebe, A., & Gibson, B. K. 2004, , 351, 399


[Gutcke et al.(2017)]Gutcke17 Gutcke, T. A., Stinson, G. S., Macciò, A. V., Wang, L., & Dutton, A. A. 2017, , 464, 2796 


[Haynes et al.(2018)]Haynes18 Haynes, M. P., Giovanelli, R., Kent, B. R., et al. 2018, , 861, 49 


[Hunter et al.(2012)]Hunter12 Hunter, D. A., Ficut-Vicas, D., Ashley, T., et al. 2012, , 144, 134 


[Ianjamasimanana et al.(2012)]Ianja12 Ianjamasimanana, R., de Blok, W. J. G., Walter, F., & Heald, G. H. 2012, , 144, 96 


[Iorio et al.(2017)]Iorio17 Iorio, G., Fraternali, F., Nipoti, C., et al. 2017, , 466, 4159 


[Jarosik et al.(2011)]WMAP7 Jarosik, N., Bennett, C. L., Dunkley, J., et al. 2011, , 192, 14 


[Johnston et al.(2008)]Johnston08 Johnston, S., Taylor, R., Bailes, M., et al. 2008, Experimental Astronomy, 22, 151 


[Kassin et al.(2012)]Kassin12 Kassin, S. A., Weiner, B. J., Faber, S. M., et al. 2012, , 758, 106 
  

[Knollmann & Knebe(2009)]Knollmann09 Knollmann,
  S. R., & Knebe, A. 2009, , 182, 608


[Koribalski et al.(2004)]Koribalski04 Koribalski, B. S., Staveley-Smith, L., Kilborn, V. A., et al. 2004, , 128, 16 


[Koribalski et al.(2018)]Koribalski18 Koribalski, B. S., Wang, J., Kamphuis, P., et al. 2018, ,  
  

[Klypin et al.(1999)]Klypin99 Klypin, A., Kravtsov, A. V., Valenzuela, O., & Prada, F. 1999, , 522, 82 
  

[Klypin et al.(2015)]Klypin15 Klypin, A., Karachentsev, I., Makarov, D., & Nasonova, O. 2015, , 454, 1798 


[Kravtsov et al.(2004)]Kravtsov04 Kravtsov, A. V., Gnedin, O. Y., & Klypin, A. A. 2004, , 609, 482 
  

[Kravtsov(2013)]Kravtsov13 Kravtsov, A. V. 2013, ,  764, L31


[Kreckel et al.(2012)]Kreckel12 Kreckel, K., Platen, E., Aragón-Calvo, M. A., et al. 2012, , 144, 16 


[Lelli et al.(2014)]Lelli14 Lelli, F., Verheijen, M., & Fraternali, F. 2014, , 566, A71 
 

[Lelli et al.(2016)]Lelli16 Lelli, F., McGaugh, S. S., & Schombert, J. M. 2016, , 152, 157 


[Macciò et al.(2010)]Maccio10 Macciò, A. V., Kang, X., Fontanot, F., et al. 2010, , 402, 1995 
  

[Macciò et al.(2016)]Maccio16 Macciò, A. V., Udrescu, S. M., Dutton, A. A., et al. 2016, , 463, L69 


[Matthee et al.(2017)]Matthee17 Matthee, J., Schaye, J., Crain, R. A., et al. 2017, , 465, 2381 
  

[Moore et al.(1999)]Moore99 Moore, B., Ghigna, S., Governato, F., et al. 1999, , 524, L19 



[More et al.(2011)]More11 More, S., van den Bosch, F. C., Cacciato, M., et al. 2011, , 410, 210

  
[Moster et al.(2018)]Moster18 Moster, B. P., Naab, T., & White, S. D. M. 2018, , 477, 1822 


[Newman et al.(2013)]Newman13 Newman, S. F., Genzel, R., Förster Schreiber, N. M., et al. 2013, , 767, 104 


[Noordermeer et al.(2005)]Noordermeer05 Noordermeer, E., van der Hulst, J. M., Sancisi, R., Swaters, R. A., & van Albada, T. S. 2005, , 442, 137 
  

[Obreschkow et al.(2013)]Obreschkow13 Obreschkow, D., Ma, X., Meyer, M., et al. 2013, , 766, 137 
  

[Papastergis et al.(2011)]Papastergis11 Papastergis, E., Martin, A. M., Giovanelli, R., & Haynes, M. P. 2011, , 739, 38


[Papastergis & Shankar(2016)]Papastergis16 Papastergis, E., & Shankar, F. 2016, , 591, A58 
 

[Papastergis & Ponomareva(2017)]Papastergis17 Papastergis, E., & Ponomareva, A. A. 2017, , 601, A1 
  

[Planck Collaboration et al.(2014)]Planck14 Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2014, , 571, A16


[Rahmati et al.(2013)]Rahmati13 Rahmati, A., Schaye, J., Pawlik, A. H., & Raičevic̀, M. 2013, , 431, 2261 


[Reddick et al.(2013)]Reddick13 Reddick, R. M., 
  Wechsler, R. H., Tinker, J. L., & Behroozi, P. S. 2013, , 771, 30
  

[Santos-Santos et al.(2018)]Santos-Santos18 Santos-Santos, I. M., Di Cintio, A., Brook, C. B., et al. 2018, , 473, 4392 


[Sawala et al.(2016)]Sawala16 Sawala, T., Frenk, C. S., Fattahi, A., et al. 2016, , 457, 1931 
  

[Schneider et al.(2014)]Schneider14 Schneider, A., Anderhalden, D., Macciò, A. V., & Diemand, J. 2014, , 441, L6 


[Schneider & Trujillo-Gomez(2018)]Schneider18 Schneider, A., & Trujillo-Gomez, S. 2018, , 475, 4809 
  

[Serra et al.(2012)]Serra12 Serra, P., Oosterloo, T., Morganti, R., et al. 2012, , 422, 1835 


[Serra et al.(2014)]Serra14 Serra, P., Oser, L., Krajnović, D., et al. 2014, , 444, 3388 


[Springel et al.(2005)]Springel05 Springel, V., White, S. D. M., Jenkins, A., et al. 2005, , 435, 629 



[Stinson et al.(2006)]Stinson06 Stinson, G., Seth, A., 
Katz, N., et al. 2006, , 373, 1074 


[Stinson et al.(2013)]Stinson13 Stinson, G. S., Brook, 
C., Macciò, A. V., et al. 2013, , 428, 129 


[Stinson et al.(2015)]Stinson15 Stinson, G. S., Dutton, A. A., Wang, L., et al. 2015, , 454, 1105 


[Swaters et al.(2002)]Swaters02 Swaters, R. A., van Albada, T. S., van der Hulst, J. M., & Sancisi, R. 2002, , 390, 829 
  

[Trujillo-Gomez et al.(2011)]Trujillo-Gomez11 Trujillo-Gomez, S., Klypin, A., Primack, J., & Romanowsky, A. J. 2011, , 742, 16 


[Trujillo-Gomez et al.(2018)]Trujillo-Gomez18 Trujillo-Gomez, S., Schneider, A., Papastergis, E., Reed, D. S., & Lake, G. 2018, , 475, 4825 


[Verbeke et al.(2017)]Verbeke17 Verbeke, R., Papastergis, E., Ponomareva, A. A., Rathi, S., & De Rijcke, S. 2017, , 607, A13 


[Verheijen & Sancisi(2001)]Verheijen01 Verheijen, M. A. W., & Sancisi, R. 2001, , 370, 765 


[Verheijen et al.(2008)]Verheijen08 Verheijen, M. A. W., Oosterloo, T. A., van Cappellen, W. A., et al. 2008, The Evolution of Galaxies Through the Neutral Hydrogen Window, 1035, 265 


[Walter et al.(2008)]Walter08 Walter, F., Brinks, E., de Blok, W. J. G., et al. 2008, , 136, 2563 


[Wang et al.(2013)]Wang13 Wang, J., Kauffmann, G., Józsa, G. I. G., et al. 2013, , 433, 270 
  

[Wang et al.(2015)]Wang15 Wang, L., Dutton, A. A., Stinson, G. S., et al. 2015, , 454, 83 


[Wang et al.(2016)]Wang16 Wang, J., Koribalski, B. S., Serra, P., et al. 2016, , 460, 2143 

  
[Wadsley et al.(2017)]Wadsley17 Wadsley, J. W., Keller, B. W., & Quinn, T. R. 2017, , 471, 2357 


[Wisnioski et al.(2015)]Wisnioski15 Wisnioski, E., Förster Schreiber, N. M., Wuyts, S., et al. 2015, , 799, 209 
  

[Zavala et al.(2009)]Zavala09 Zavala, J., Jing, Y. P., Faltenbacher, A., et al. 2009, , 700, 1779 


[Zwaan et al.(2010)]Zwaan10 Zwaan, M. A., Meyer, M. J., & Staveley-Smith, L. 2010, , 403, 1969 
  


