# The Odd Dark Matter Halos of Isolated Gas-rich Ultradiffuse Galaxies

Demao Kong<EMAIL_ADDRESS>Department of Physics and Astronomy, Tufts
University, Medford, MA 02155, USA  Department of Physics and Astronomy,
University of California, Irvine, CA, 92697, USA  Manoj Kaplinghat
<EMAIL_ADDRESS>Department of Physics and Astronomy, University of
California, Irvine, CA, 92697, USA  Hai-Bo Yu<EMAIL_ADDRESS>Department of
Physics and Astronomy, University of California, Riverside, CA 92521, USA
Filippo Fraternali Kapteyn Astronomical Institute, University of Groningen,
Landleven 12, 9747 AD, Groningen, The Netherlands  Pavel E. Mancera Piña
Kapteyn Astronomical Institute, University of Groningen, Landleven 12, 9747
AD, Groningen, The Netherlands  ASTRON, Netherlands Institute for Radio
Astronomy, Postbus, 7900 AA, Dwingeloo, The Netherlands

###### Abstract

We analyze circular velocity profiles of seven ultradiffuse galaxies (UDGs)
that are isolated and gas-rich. Assuming that the dark matter halos of these
UDGs have a Navarro-Frenk-White (NFW) density profile or a Read density
profile (which allows for constant-density cores), the inferred halo
concentrations are systematically lower than the cosmological median, even as
low as $-0.6$ dex (about $5\sigma$ away) in some cases. Alternatively, similar
fits can be obtained with a density profile that scales roughly as $1/r^{2}$
for radii larger than a few kiloparsecs. Both solutions require the radius
where the halo circular velocity peaks ($R_{\rm max}$) to be much larger than
the median expectation. Surprisingly, we find an overabundance of such
large-$R_{\rm max}$ halos in the IllustrisTNG dark-matter-only simulations
compared to the Gaussian expectation. These halos form late and have higher
spins compared to median halos of similar masses. The inner densities of the
most extreme among these late-forming halos are higher than their NFW
counterparts, leading to a $\sim 1/r^{2}$ density profile. However, the two
well-resolved UDGs in our sample strongly prefer lower dark matter densities
in the center than the simulated ones. Comparing to IllustrisTNG
hydrodynamical simulations, we also find a tension in getting both low enough
circular velocities and high enough halo mass to accommodate the measurements.
Our results indicate that the gas-rich UDGs present a significant challenge
for galaxy formation models.

Low surface brightness galaxies (940), Dark matter (353)

##  1 Introduction

Ultradiffuse galaxies (UDGs) are a class of galaxies that have an extremely
low luminosity. Compared to usual low surface brightness galaxies, UDGs are
more extended in the light distribution and hence fainter. Recent observations
have revealed a large number of UDGs in galaxy clusters (van Dokkum et al.,
2015; Koda et al., 2015; Yagi et al., 2016; Janssens et al., 2017; Wittmann et
al., 2017; Mancera Piña et al., 2018; Zaritsky et al., 2019; Iodice et al.,
2020; Lee et al., 2020); see also Binggeli et al. (1984); Impey et al. (1988)
for earlier studies and Conselice (2018) for related discussions. UDGs have
also been found in the in groups or in the field (Martínez-Delgado et al.,
2016; Bennet et al., 2017; Leisman et al., 2017; Román & Trujillo, 2017; van
der Burg et al., 2017; He et al., 2019; Janowiecki et al., 2019; Prole et al.,
2019; Román et al., 2019; Guo et al., 2020; Mancera Piña et al., 2020; Prole
et al., 2021; Shi et al., 2021).

The existence of these extreme galaxies in both high- and low-density
environments has motivated many studies on their formation mechanisms. For
UDGs in clusters and groups, environmental effects, such as ram pressure
stripping and tidal stripping, likely play an important role (van Dokkum et
al., 2015; Yozin & Bekki, 2015; Conselice, 2018; Ogiya, 2018; Carleton et al.,
2019; Jiang et al., 2019; Collins et al., 2020; Yang et al., 2020; Tremmel et
al., 2020; Sales et al., 2020; Benavides et al., 2021; Moreno et al., 2022).
For UDGs in the field, they could form in dwarf halos on the distribution tail
of high angular momentum (Amorisco & Loeb, 2016; Román & Trujillo, 2017; Liao
et al., 2019; Wright et al., 2021). In addition, gas outflows driven by
baryonic feedback could change the gravitational potential and expel stars to
more external orbits, producing extended and faint galaxies (Di Cintio et al.,
2016; Chan et al., 2018). It is possible that UDGs are produced owing to a
combination of multiple formation channels. For example, both feedback and
environmental effects could be important for UDGs in high-density environments
(Martin et al., 2019; Jackson et al., 2021).

Kinematic measurements are essential for understanding mass distributions in
UDGs and further testing their formation mechanisms, although they are often
difficult owning to the low surface brightness nature (Greco et al., 2018;
Martín-Navarro et al., 2019; Gannon et al., 2020). Most isolated UDGs in the
field are gas-rich and recent measurements from their neutral gas Hi emissions
indicate that their circular velocities are low (Leisman et al., 2017; Mancera
Piña et al., 2020). In fact, Shi et al. (2021) and Mancera Piña et al. (2022)
showed, using high-resolution data, that the host halos of two gas-rich UDGs
have low concentrations.

In this work, we focus on a sample of seven gas-rich UDGs in the field: five
of them have low-resolution kinematic data from Mancera Piña et al. (2020),
and two have high-resolution ones from Mancera Piña et al. (2022); Shi et al.
(2021). We propose a unified model to describe the surface mass density of Hi
gas for both high- and low-resolution UDGs, fit their circular velocity data,
and infer properties of dark matter halos. We will show that the host halos
are “odd” in the sense that they are extremely low concentration or that their
profiles depart dramatically from the Navarro–Frenk–White (NFW) profile in the
inner and outer regions or both. We search the IllustrisTNG simulations
(Nelson et al., 2019) and find an overabundance of such halos, much higher
than expected from a Gaussian tail. We will further discuss characteristic
properties of the simulated halos and challenges in fully accommodating the
measurements of the gas-rich UDGs.

The rest of the paper is organized as follows. In Sec. 2, we discuss mass
modeling. We provide fits to the circular velocity profiles in Sec. 3, and we
infer halo parameters and compare them to the IllustrisTNG simulations in Sec
4. We study inner dark matter densities of the halos in Sec 5. We discuss
properties of the IllustrisTNG halos in Sec. 6. We highlight tensions between
the field UDGs and their IllustrisTNG analogs and discuss possible solutions
in Sec. 7. We conclude in Sec. 8. In the Appendix, we provide additional
information on the mass modeling, fits, and density profiles of the simulated
halos.

##  2 Mass Modeling

Figure 1: Gas surface mass density profiles (solid orange curves) for resolved
AGC 114905 (top panel) and AGC 242019 (bottom panel), reconstructed using the
gas model in equation 1. The measured profiles are shown for comparison (gray
circles), taken from Mancera Piña et al. (2022) and Shi et al. (2021). The
helium correction is included.

We consider a sample of seven gas-rich isolated UDGs with measured circular
velocity profiles based on their Hi kinematics. Two of them have high-
resolution data and their gas surface mass density profiles that are well
resolved: AGC 114905 (Mancera Piña et al., 2022) and AGC 242019 (Shi et al.,
2021). Five are low resolution (two resolution elements): AGC 122966, AGC
219533, AGC 248945, AGC 334315 and AGC 749290 (Mancera Piña et al., 2019,
2020). For these low-resolution UDGs, we propose the following model to
describe their Hi gas surface mass density:

$\displaystyle\Sigma_{\rm\footnotesize\textsc{Hi}}(R)={\Sigma_{\rm\footnotesize\textsc{Hi}}(0)}{\left[1+(bR/R_{\rm\footnotesize\textsc{Hi}})^{m}\right]^{-n/m}}$
(1)

where $\Sigma_{\rm\footnotesize\textsc{Hi}}(0)$ is the central surface
density, $R_{\rm\footnotesize\textsc{Hi}}$ is the radius where
$\Sigma_{\rm\footnotesize\textsc{Hi}}(R_{\rm\footnotesize\textsc{Hi}})=1~{}{\rm
M_{\odot}/pc^{2}}$, the parameter $b$ is given by the condition
$b=\\{[\Sigma_{\rm\footnotesize\textsc{Hi}}(0)/{{\rm(M_{\odot}/pc^{2})}}]^{m/n}-1\\}^{1/m}$,
and we determine numerical factors $n$ and $m$ by comparing the model with the
high-resolution UDGs.

For $m=n=4$, we find that the model can reproduce the observed Hi surface mass
density of the two UDGs with high-resolution data. In this case, the total
enclosed Hi gas mass within the radius $R_{\rm out}$ can be calculated as

$\displaystyle M_{\rm\footnotesize\textsc{Hi}}(R_{\rm
out})=\pi\Sigma_{\rm\footnotesize\textsc{Hi}}(0)\left[\frac{R_{\rm\footnotesize\textsc{Hi}}}{b}\right]^{2}{\rm
Arctan}\left[\left(\frac{bR_{\rm
out}}{R_{\rm\footnotesize\textsc{Hi}}}\right)^{2}\right].$ (2)

Fig. 1 shows our gas surface mass density profiles for AGC 114905 and AGC
242019 (solid orange curves), compared to measured ones (gray circles), after
including helium correction $\Sigma_{\rm
gas}=1.33\Sigma_{\rm\footnotesize\textsc{Hi}}$. Our model well reproduces the
observations of the two galaxies. For AGC 114905, we fix
$R_{\rm\footnotesize\textsc{Hi}}=11.4~{}{\rm kpc}$, which is $1\sigma$ higher
than the median reported in Gault et al. (2021). For AGC 242019,
$R_{\rm\textsc{Hi}}=10.4~{}{\rm kpc}$. We have also checked that the enclosed
Hi mass calculated using Eq. 2 agrees well with the measured one for the two
UDGs. Our model slightly overestimates the gas surface mass density in the
inner region $R\lesssim 2~{}{\rm kpc}$ for AGC 114905. Since this region is
dominated by the stellar mass, the discrepancy in $\Sigma_{\rm gas}$ has a
negligible effect on the fit.

We apply the model ($m=n=4$) to the five low-resolution UDGs in our sample. We
take two approaches to determine $R_{\rm\footnotesize\textsc{Hi}}$ in our
fits: fix $R_{\rm\footnotesize\textsc{Hi}}$ to be $1\sigma$ higher than the
median for individual UDGs as listed in Table 3 of Gault et al. (2021), see
our Table 1, motivated by the comparison shown in Fig. 1, and treat
$R_{\rm\footnotesize\textsc{Hi}}$ as a varying parameter. We find that the
inferred halo properties are almost the same in both approaches. In this
paper, we will report the results from the former. We obtain $\Sigma_{\rm
gas}(0)$ using the condition in Eq. 2, where we fix $R_{\rm out}$ for each of
the low-resolution UDGs as listed in Table 1, while varying
$M_{\rm\footnotesize\textsc{Hi}}$ (see Sec. 3 for details). Fig. 11 in the
Appendix shows detailed $\Sigma_{\rm\footnotesize\textsc{Hi}}(R)$ profiles for
the low-resolution UDGs after taking representative mode parameters, where we
also include profiles with $m=4$ and $n=2.5$ for comparison.

For the stellar component, we use a thin-disk model as (Freeman, 1970)

$\Sigma_{\star}(R)=\Sigma_{\star}(0)\mathrm{e}^{-R/R_{\star}},$ (3)

where $\Sigma_{\star}(0)$ is the central surface density and $R_{\star}$ is
the scale length of the stellar disk. The total stellar mass can be calculated
as $M_{\star}=2\pi\Sigma_{\star}(0)R^{2}_{\star}$. For the five low-resolution
UDGs, we fix $M_{\star}$ and $R_{\star}$ to their corresponding median values;
see Table 1. For the high-resolution ones, Mancera Piña et al. (2022) and Shi
et al. (2021) have obtained profiles of baryonic circular velocities based on
the measured stellar and gas distributions, and we directly use their profiles
in our fits.

For the halo, we first use a restricted form of the Read profile (Read et al.,
2016). The enclosed mass within a radius $r$ is given by

$M(r)=4\pi\rho_{s}r^{3}_{s}\left[\ln\left(1+\frac{r}{r_{s}}\right)-\frac{r}{r+r_{s}}\right]f(r)^{\delta},$
(4)

where $r_{s}$ is the scale radius, $\rho_{s}$ is the scale density, $\delta$
is a numerical factor, and $f(r)=\tanh(r/r_{c})$ with $r_{c}$ being the core
radius. For $r\gg r_{c}$, the function $f(r)$ tends to $1$ and we recover the
NFW profile with the same $r_{s}$ and $\rho_{s}$. Mancera Piña et al. (2022)
analyzed AGC 114905 with $\delta$ and $r_{c}$ being varied according to the
relations found in simulations Read et al. (2016), and found that $\delta$ is
effectively $0$ if the concentration is low. Since those relations are not
calibrated for low-concentration halos and current measurements of circular
velocities are unlikely to be sensitive to both $r_{c}$ and $\delta$, we relax
the constraints and set $\delta=1$ to ensure a complete transition to a cored
profile, while leaving the core radius $r_{c}$ as a free parameter. This
choice allows the fits to explore a broad range of profiles from NFW-like to
large core sizes of order $r_{s}$.

For completeness, we perform additional fits with the Read profile in two
extreme limits. By setting $\delta=0$, we consider an exact NFW profile.
Moreover, we will also treat both $\delta$ and $r_{c}$ as free parameters,
i.e., a general Read profile. The goal is to demonstrate that the low halo
concentration inferred from the gas-rich UDGs is driven by the high Hi gas
mass, and it is robust to the variation of the inner density profile. The
results of these additional fits will be shown in the Appendix.

One of the benefits of using the Read profile is that we can directly relate
$r_{s}$ and $\rho_{s}$ to constraints from the cosmological concentration-mass
$c_{200}\textup{--}M_{200}$ relation (see Diemer & Joyce 2019), which is often
based on fitting simulated halos with an NFW profile. In this work, we will
also use the maximal circular velocity $V_{\rm max}$ and its corresponding
radius $R_{\rm max}$ to specify a halo. For an NFW-like profile like the Read
one, $c_{200}=r_{200}/r_{s}$ with $r_{200}$ being the virial radius of the
halo; $V_{\rm max}=1.64r_{s}\sqrt{G\rho_{s}}$, where $G$ is Newton’s constant;
and and $R_{\rm max}=2.16r_{s}$.

The assumption of an NFW halo profile in the outer regions may not be a good
approximation for halos far below the median concentration. Hence, we also
consider a more general double-power-law (DPL) density profile, defined as

$\rho(r)=\rho_{s}\left(\frac{r}{r_{s}}\right)^{-\gamma}\left(1+\frac{r}{r_{s}}\right)^{\gamma-\beta},$
(5)

where $\gamma$ and $\beta$ are numerical factors. For $\gamma=1$ and
$\beta=3$, we recover the NFW profile. For $\beta\sim 2$, $\rho(r)\propto
1/r^{2}$ as $r>r_{s}$ and the resulting profile of halo circular velocities is
nearly flat. We further generalize the calculation of the concentration
parameter to $c_{200}=r_{200}/r_{-2}$, where
$r_{-2}=(2-\gamma)r_{s}/(-2+\beta)$ is the radius at which the logarithmic
slope of the density profile is $-2$. For the NFW profile with $\gamma=1$ and
$\beta=3$, we have the well-known relations $r_{-2}=r_{s}$ and
$c_{200}=r_{200}/r_{s}$. However, for $\beta\sim 2$, $r_{-2}$ can be much
larger than $r_{s}$.

As we will discuss later, for many individual “low-concentration” halos in the
IllustrisTNG simulations that could potentially host the UDGs, $R_{\rm max}$
is close to $r_{200}$ and their profiles of circular velocities are
surprisingly flat. Thus, the DPL profile with $\beta\sim 2$ is particularly
interesting, and we will explicitly show that it provides a good fit to the
IllustrisTNG halos in Sec. 7. Note that for those simulated halos the density
drops sharply around the virial radius $r_{200}$ and hence the halo mass
$M_{200}=M(r_{200})$ is well defined. Overall, the Read and DPL ($\beta\sim
2$) profiles together allow us to test the robustness of our inferences about
UDG halos to variations in both the inner and outer density profiles.

AGC ID | $\left[\frac{M_{\star}}{\rm M_{\odot}}\right]$ | $\frac{R_{\star}}{\rm kpc}$ | $\left[\frac{M_{\rm\footnotesize\textsc{Hi}}}{\rm M_{\odot}}\right]$ | $\frac{R_{\rm\footnotesize\textsc{Hi}}}{\rm kpc}$ | $\frac{R_{\rm out}}{\rm kpc}$ | $\left[\frac{\rho_{s}}{\rm M_{\odot}kpc^{3}}\right]$ | $\frac{r_{s}}{\rm kpc}$ | $\frac{r_{c}}{\rm kpc}$ | $\left[\frac{M_{\rm\footnotesize\textsc{Hi}}}{\rm M_{\odot}}\right]$ | $\left[\frac{M_{200}}{\rm M_{\odot}}\right]$ | $\frac{r_{200}}{r_{s}}$
---|---|---|---|---|---|---|---|---|---|---|---
114905 | $8.12$ | $1.79$ | $8.99_{-0.07}^{+0.07}$ | $11.4$ | $12$ | $5.24_{-0.12}^{+0.25}$ | $23.82_{-6.13}^{+7.40}$ | $6.74_{-3.22}^{+2.29}$ | - | $10.14_{-0.11}^{+0.24}$ | $2.11_{-0.28}^{+0.70}$
122966 | $7.45$ | $4.15$ | $9.03_{-0.05}^{+0.05}$ | $11.2$ | $14.3$ | $6.11_{-0.72}^{+0.66}$ | $9.29_{-5.01}^{+15.75}$ | $0.59_{-0.41}^{+1.42}$ | $8.97_{-0.07}^{+0.11}$ | $10.08_{-0.19}^{+0.38}$ | $5.76_{-3.16}^{+6.06}$
219533 | $7.82$ | $2.35$ | $9.24_{-0.06}^{+0.06}$ | $14.2$ | $12.9$ | $5.97_{-0.71}^{+0.62}$ | $12.78_{-6.89}^{+25.76}$ | $0.57_{-0.40}^{+1.35}$ | $9.14_{-0.06}^{+0.10}$ | $10.32_{-0.23}^{+0.51}$ | $4.83_{-2.63}^{+4.44}$
248945 | $8.35$ | $2.08$ | $8.78_{-0.06}^{+0.06}$ | $7.6$ | $11.4$ | $5.48_{-0.28}^{+0.50}$ | $15.73_{-7.23}^{+9.60}$ | $0.85_{-0.65}^{+3.15}$ | $8.7_{-0.08}^{+0.12}$ | $9.92_{-0.18}^{+0.32}$ | $2.84_{-0.82}^{+2.31}$
334315 | $7.65$ | $3.76$ | $9.22_{-0.05}^{+0.05}$ | $14.4$ | $11.2$ | $5.32_{-0.19}^{+0.33}$ | $22.35_{-7.62}^{+9.38}$ | $0.62_{-0.44}^{+1.90}$ | $9.1_{-0.02}^{+0.05}$ | $10.16_{-0.14}^{+0.26}$ | $2.31_{-0.45}^{+1.04}$
749290 | $8.11$ | $2.38$ | $8.95_{-0.05}^{+0.05}$ | $9.7$ | $11.2$ | $5.40_{-0.22}^{+0.40}$ | $18.15_{-7.09}^{+9.33}$ | $0.91_{-0.71}^{+3.60}$ | $8.85_{-0.04}^{+0.08}$ | $10.00_{-0.17}^{+0.30}$ | $2.54_{-0.59}^{+1.50}$
242019 | $8.14$ | - | $8.93_{-0.02}^{+0.02}$ | 10.4 | 15.7 | $5.08_{-0.13}^{+0.21}$ | $47.12_{-15.96}^{+14.77}$ | $0.26_{-0.12}^{+0.27}$ | - | $10.76_{-0.22}^{+0.17}$ | $1.75_{-0.25}^{+0.48}$
242019† | $8.14$ | - | $8.93_{-0.02}^{+0.02}$ | 10.4 | 15.7 | $5.47_{-0.43}^{+0.61}$ | $24.14_{-14.42}^{+8.98}$ | $3.15_{-1.23}^{+1.67}$ | - | $10.39_{-0.35}^{0.44}$ | $3.18_{-1.46}^{+0.44}$
Table 1: Parameters used in mass modeling and those inferred from the fits with the Read density profile in Eq. 4. Columns from left to right: galaxy AGC ID, stellar mass ($\log_{10}M_{\star}$), scale radius of stellar disks ($R_{\star}$), Hi gas mass ($\log_{10}M_{\rm\footnotesize\textsc{Hi}}$), scale radius of Hi surface mass densities ($R_{\rm\footnotesize\textsc{Hi}}$), the farthest radius at which the total Hi mass is measured ($R_{\rm out}$); inferred halo scale density ($\log_{10}\rho_{s}$), scale radius ($r_{s}$), Hi gas mass ($\log_{10}M_{\rm\footnotesize\textsc{Hi}}$), halo mass ($\log_{10}M_{200}$) and concentration ($c_{200}=r_{200}/r_{s}$). For AGC 114905, the measurement data shown in the table are from Mancera Piña et al. (2022) and we take their baryonic contribution to the total circular velocity, see the left panel of their Fig. 5 (distance $D=67~{}{\rm Mpc}$, inclination $i=26^{\rm o}$). For AGC 242019, the data are compiled from Shi et al. (2021) and we directly fit the circular velocity profile of the dark matter halo listed in their Table 2. For the low-resolution UDGs, the Hi measurement and the stellar data are compiled from Mancera Piña et al. (2019, 2020); Gault et al. (2021), and we assume the distance and inclination for individual UDGs as in Mancera Piña et al. (2020). We re-calibrated the stellar mass using the mass-to-light-color relation in Du et al. (2020), yielding slightly smaller $M_{\star}$ values than those reported in Mancera Piña et al. (2020).  AGC ID | $\left[\frac{\rho_{s}}{\rm M_{\odot}kpc^{3}}\right]$ | $\frac{r_{s}}{\rm kpc}$ | ${\beta}$ | ${\gamma}$ | $\left[\frac{M_{\rm\footnotesize\textsc{Hi}}}{\rm M_{\odot}}\right]$ | $\left[\frac{M_{200}}{\rm M_{\odot}}\right]$ | $\frac{r_{200}}{r_{s}}$ | $\frac{r_{200}}{r_{-2}}$
---|---|---|---|---|---|---|---|---
114905 | $6.25_{-0.07}^{+0.1}$ | $4.15_{-0.51}^{+0.46}$ | $2.05_{-0.03}^{+0.07}$ | $0.09_{-0.07}^{+0.13}$ | - | $9.98_{-0.06}^{+0.12}$ | $10.84_{-0.63}^{+1.44}$ | $0.28_{-0.21}^{+0.41}$
122966 | $7.34_{-0.71}^{+0.59}$ | $1.34_{-0.69}^{+1.9}$ | $2.17_{-0.08}^{+0.09}$ | $0.85_{-0.53}^{+0.45}$ | $8.98_{-0.07}^{+0.1}$ | $10.12_{-0.16}^{+0.22}$ | $37.59_{-20.67}^{+35.58}$ | $5.11_{-2.85}^{+6.19}$
219533 | $7.29_{-0.69}^{+0.6}$ | $1.61_{-0.84}^{+2.25}$ | $2.17_{-0.08}^{+0.09}$ | $0.85_{-0.53}^{+0.44}$ | $9.15_{-0.06}^{+0.09}$ | $10.28_{-0.16}^{+0.23}$ | $35.66_{-19.29}^{+34.84}$ | $4.83_{-2.61}^{+5.39}$
248945 | $6.85_{-0.47}^{+0.74}$ | $1.94_{-1.16}^{+1.67}$ | $2.15_{-0.07}^{+0.09}$ | $0.59_{-0.42}^{+0.53}$ | $8.69_{-0.06}^{+0.11}$ | $9.86_{-0.14}^{+0.23}$ | $21.54_{-9.24}^{+29.2}$ | $2.59_{-1.43}^{+3.64}$
334315 | $6.63_{-0.31}^{+0.54}$ | $3.01_{-1.49}^{+1.48}$ | $2.14_{-0.06}^{+0.1}$ | $0.5_{-0.34}^{+0.46}$ | $9.1_{-0.02}^{+0.04}$ | $10.09_{-0.09}^{+0.18}$ | $16.6_{-5.13}^{+15.04}$ | $1.79_{-0.87}^{+1.6}$
749290 | $6.68_{-0.34}^{+0.63}$ | $2.52_{-1.39}^{+1.49}$ | $2.14_{-0.07}^{+0.09}$ | $0.49_{-0.35}^{+0.5}$ | $8.85_{-0.04}^{+0.07}$ | $9.93_{-0.12}^{+0.23}$ | $17.56_{-5.95}^{+19.37}$ | $1.95_{-1.0}^{+2.19}$
242019 | $6.44_{-0.19}^{+0.29}$ | $5.24_{-1.83}^{+1.76}$ | $2.11_{-0.08}^{+0.11}$ | $0.5_{-0.22}^{+0.16}$ | - | $10.57_{-0.16}^{+0.18}$ | $13.51_{-2.64}^{+5.7}$ | $1.05_{-0.74}^{+0.91}$
242019† | $6.51_{-0.2}^{+0.2}$ | $5.73_{-1.75}^{+2.21}$ | $2.13_{-0.09}^{+0.11}$ | $0.24_{-0.16}^{+0.19}$ | - | $10.75_{-0.19}^{+0.22}$ | $14.19_{-2.91}^{+4.26}$ | $1.07_{-0.73}^{+0.74}$
Table 2: Model parameters inferred from the fits with the DPL density profile
in Eq. 5, including the halo scale density $\log_{10}\rho_{s}$, scale radius
($r_{s}$), outer logarithmic slope of the profile ($\beta$), inner logarithmic
slope ($\gamma$), Hi gas mass ($\log_{10}M_{\rm\footnotesize\textsc{Hi}}$),
halo mass ($\log_{10}M_{200}$), ratio of $r_{200}$ to $r_{s}$, and
concentration ($c_{200}=r_{200}/r_{s}$). The parameters for modeling stellar
and gas distributions are the same as those in Table 1. The inferences for
$\beta$ reflect the right prior on it $2.05<\beta<2.3$.  Figure 2: Circular
velocity profiles of the best fit (solid black), including halo (solid red)
and baryonic (dashed magenta) contributions, compared to observational data
(black circles), based on the Read profile in Eq. 4. For AGC 242019†, the fit
is performed without including the innermost datapoint. The horizontal lines
denote the halo core size $r_{c}$ within the $68$th percentile inferred from
the fits (red). The observational data are taken from Mancera Piña et al.
(2022, 2020) and Shi et al. (2021). For the fits based on the DPL profile in
Eq. 5, see Fig. 17 in the Appendix.

##  3 Fits to the measured circular velocities

We calculate the total predicted circular velocity as $V^{\rm mod}_{\rm
circ}=\sqrt{V^{2}_{\rm star}+V^{2}_{\rm gas}+V^{2}_{\rm halo}}$, where $V_{\rm
star}$, $V_{\rm gas}$ and $V_{\rm halo}$ are contributions from the stars,
gas, and dark matter halo, respectively. For AGC 114905, we take the profile
of baryonic circular velocities from Mancera Piña et al. (2022); see the left
panel of their Fig. 5 ($D=67~{}{\rm Mpc}$, $i=26^{\rm o}$). For AGC 242019, we
directly fit to the halo circular velocities after subtraction of baryonic
contributions, as listed in Table 2 of Shi et al. (2021). For the low-
resolution UDGs, we calculate $V_{\rm gas}$ and $V_{\rm star}$ using the thin-
disk approximation with the surface density profiles in Eqs. 1 and 3,
respectively (Binney & Tremaine 2008).

For the contribution from gas,

$\displaystyle V^{2}_{\rm
gas}=-4G\int^{R}_{0}\frac{ada}{\sqrt{R^{2}-a^{2}}}\frac{d}{da}\int^{\infty}_{a}\frac{R^{\prime}\Sigma_{\rm
gas}(R^{\prime})dR^{\prime}}{\sqrt{R^{\prime 2}-a^{2}}}.$ (6)

In practice, we convert Eq. 6 into a dimensionless form and generate numerical
templates for the fits. For the stellar component, $V_{\rm star}$ has an
analytical form

$\displaystyle V^{2}_{\rm star}=4\pi
G\Sigma_{\star}(0)R_{\star}y^{2}[I_{0}(y)K_{0}(y)-I_{1}(y)K_{1}(y)],$ (7)

where $y\equiv R/2R_{\star}$, and $I_{0,1}$ and $K_{0,1}$ are modified Bessel
functions.

We use an Markov Chain Monte Carlo (MCMC) routine emcee (Foreman-Mackey et
al., 2013) to fully explore the halo parameters and find the best fit that
minimizes the residuals of $V^{\rm mod}_{\rm circ}$. The likelihood function
is $\exp(-0.5\chi^{2})$, where $\chi^{2}={(V^{\rm mod}_{\rm circ}-V^{\rm
obs}_{\rm circ})^{2}/\delta V^{2}_{\rm circ}}$. We impose the following flat
priors on the relevant parameters. For the Read profile, we impose
$7\leq\log_{10}(M_{200}/{\rm M_{\odot}})\leq 12$, $0\leq\log_{10}c_{200}\leq
1.7$, $-1\leq\log_{10}(r_{c}/{\rm kpc})\leq 1$, and $r_{c}/r_{s}<1$. For the
DPL profile, we have the same prior on $M_{200}$, while for the other
parameters $1\leq\log_{10}(r_{200}/r_{s})\leq 2$, $0<\gamma<1.5$, and
$2.05<\beta<2.3$. Both the $r_{s}$ and $\beta$ priors are chosen to recover a
profile of $\rho\propto 1/r^{2}$ behavior as evident within $10~{}{\rm kpc}$,
where we have the measurement data. The prior on the inner slope $\gamma$ is
generous and allows the profile to be cored or cuspy.

For the low-resolution UDGs, we impose a flat prior on
$\log(M_{\rm\footnotesize\textsc{Hi}}/{\rm M_{\odot}})$ within the $\pm
3\sigma$ range. In addition, we require that $M_{200}$ is large enough such
that the cosmological limit $(M_{\rm
gas}+M_{\star})/M_{200}\leq\Omega_{b}/(\Omega_{m}-\Omega_{b})=0.187$ is
satisfied (Planck Collaboration et al., 2020). This puts a conservative lower
limit on the halo mass for a given total baryon mass. We summarize the
parameters used for mass modeling and those inferred from our fits in Tables 1
and 2.

Fig. 2 shows circular velocity profiles of the total best fit (black),
including halo (red) and baryonic (magenta) contributions, based on the Read
profile in Eq. 4. Compared to the measured kinematic data (black circles), our
constructed model works reasonably well. We notice that except for AGC 242019,
the fits overestimate the circular velocity at large radii because of
significant contributions from the gas disk. In particular, for AGC 114905,
AGC 248945, and AGC 749290, the gas and halo contributions are comparable at
the radius where the last datapoint was measured (close to
$R_{\rm\footnotesize\textsc{Hi}}$). This reflects the fact that these UDGs
have a high baryon content (Mancera Piña et al., 2019). We note that $M_{200}$
is close to the lower limit set by the cosmological baryon fraction prior
because the circular velocity data push the dark matter contribution to be
small. The best fits with the DPL profile are similar to the fits with the
Read profile (Fig. 17 in Appendix). The inferred halo masses are similar as
well, because of the strong influence of the lower limit on $M_{200}/(M_{\rm
gas}+M_{\star})$.

In Fig. 2, we also show the core size $r_{c}$ inferred from the fits (red
horizontal line). In all low-resolution cases, core sizes of a few kiloparsecs
are allowed. But it cannot be larger because circular velocities between $2$
and $8~{}{\rm kpc}$ do not change significantly. The core sizes can also be
zero and hence consistent with an NFW profile. AGC 114905 allows for a core
size of $r_{c}\sim 7~{}{\rm kpc}$, but the other high-resolution one, AGC
242019, prefers effectively zero core size (Shi et al., 2021). The preference
is largely driven by the innermost datapoint of AGC 242019, which has a small
error of $1.2~{}\rm kms^{-1}$. When we remove this datapoint (denoted by
$242019^{\dagger}$ in the legend), there is mild evidence for a nonzero core
size of $r_{c}\sim 4~{}{\rm kpc}$. Removing the innermost datapoint also
brings AGC 242019 in agreement with the other UDGs in terms of the preferred
$V_{\rm max}$ and $M_{200}$ values as we discuss below. This highlights the
importance of accurate measurements at radii below $5~{}{\rm kpc}$.

##  4 Inferred halo parameters compared to simulations

Figure 3: Left: $R_{\rm max}\textup{--}V_{\rm max}$ distributions inferred
from the Read fits (colored), together with those from the TNG50-1-Dark
simulations (gray circles). The solid black lines show the median relation
(Diemer & Joyce, 2019), as well as $0.3$ dex and $0.6$ dex larger $R_{\rm
max}$ (“lower concentration”) compared to the median. The dashed black lines
denote the halo masses $\log_{10}(M_{200}/{\rm M_{\odot}})=10,~{}10.5$ and
$11$. The IllustrisTNG halos within the red box are selected for further study
in the right panel. Right: the halo circular velocity at radii $8~{}{\rm kpc}$
(top panel) and $2~{}{\rm kpc}$ (bottom panel) vs. maximum circular velocity
for the observed UDGs (symbols colored according to the legend in the panel on
the left), compared to the selected TNG50-1-Dark halos (gray circles).

Figure 4: Left: $R_{\rm max}\textup{--}V_{\rm max}$ distributions inferred
from the DPL fits (colored), together with those from TNG50-1-Dark simulation
(gray circles). The solid and dashed black lines and the red box are the same
as in Fig. 3. Right: the halo circular velocity at radii $8~{}{\rm kpc}$ (top
panel) and $2~{}{\rm kpc}$ (bottom panel) vs. maximum circular velocity for
the observed UDGs (symbols colored according to the legend in the panel on the
left), compared to the TNG50-1-Dark halos (gray circles) within the red box.

We will compare the halo parameters inferred from our fits with those from the
IllustrisTNG project in Nelson et al. (2019), which implements a comprehensive
model for galaxy formation physics. We take the snapshot $99$ ($z=0$) in the
TNG50-1-Dark dataset based on dark-matter-only simulations, which is publicly
available, and choose simulated halos with the mass $\log_{10}(M_{200}/{\rm
M_{\odot}})=9.3\textup{--}11$ (Primary Flag=1). As we will show later, this
mass range is appropriate for the UDGs we consider. In total, there are
$37411$ halos, and we retrieve their properties, including concentrations,
density profiles, formation histories, and spins.

We choose the TNG50-1-Dark simulated halos for comparison becuase they provide
a clean benchmark sample where baryonic effects on the halo properties are
absent. For the IllustrisTNG hydrodynamical simulations, supernova feedback is
“weak” in the sense that it does not lead to a lowering of the inner dark
matter densities (Lovell et al., 2018). Thus, the results from the dark-
matter-only simulations should provide a good approximation for the purposes
of the current work. We will leave a detailed comparison with the full
hydrodynamical TNG50 simulations for future work.

Fig. 3 (left) shows $R_{\rm max}$ and $V_{\rm max}$ of the halos of the UDGs
inferred from the Read fits (colored). If $R_{\rm max}$ is larger than
$r_{200}$, we use $r_{200}$ instead. For the UDGs we consider, the values of
their halo parameters are $V_{\rm max}\sim 30\textup{--}50~{}{\rm km/s}$,
$R_{\rm max}\sim 20\textup{--}100~{}{\rm kpc}$, and $M_{200}\sim
10^{10}\textup{--}10^{11}~{}{\rm M_{\odot}}$. For comparison, we plot the
$R_{\rm max}\textup{--}V_{\rm max}$ median relation predicted in cosmological
simulations (Diemer & Joyce, 2019), as well as the lines with $R_{\rm max}$
values $0.3$ dex and $0.6$ dex larger than the median for given $V_{\rm max}$
(solid black). Since the $1\sigma$ scatter for the lognormal
$c_{200}\textup{--}M_{200}$ distribution is $0.11$ dex, the concentration of
the UDG halos would be $\sim 3\sigma\textup{--}5\sigma$ below the cosmological
median when if we were to extrapolate the lognormal distribution to the low
end. Thus, we conclude that the UDG halos must have extremely low
concentrations (large $R_{\rm max}$) if they are hosted in NFW-like halos; see
also Mancera Piña et al. (2022); Shi et al. (2021) for the two high-resolution
UDGs. We reach the same conclusion if we adopt the NFW profile or the general
Read profile by varying both $r_{c}$ and $\delta$ simultaneously; see the
Appendix.

We have further checked that the conclusion holds for systematically denser
gas surface densities, for example, with $n=2.5$ in Eq. 1 (see Fig. 12 in the
Appendix), and it holds even when only including the outmost measured
datapoint in the fits. The preference for a large $R_{\rm max}$ solution does
not depend on whether the profile is cored or cuspy, as we have marginalized
over the core radius in the Read profile. We have checked that the core radius
is not strongly correlated with other model parameters. From Fig. 3 (left), we
also see that removing the innermost datapoint brings AGC 242019 more in line
with the other UDGs in the $R_{\rm max}\textup{--}V_{\rm max}$ plane.

If a Gaussian tail of the $c_{200}\textup{--}M_{200}$ distribution is extended
to the low-concentration end, we would not expect to find halos with
concentrations as low as $0.5$ dex below the median. However, the Gaussian
extrapolation is unwarranted. Fig. 3 (left) shows results from the
TNG50-1-Dark simulations (gray circles), where we have imposed the selection
condition $\log_{10}(M_{200}/{\rm M_{\odot}})=9.3\textup{--}11$ corresponding
to $V_{\rm max}\sim 20\textup{--}80~{}{\rm km/s}$. There is clearly a non-
Gaussian tail to low concentrations that dominates about $0.2$ dex below the
median; see also Fig. 6 (left). In the $50h^{-1}\rm Mpc$ box, we find $956$
halos within this mass range with concentrations of $0.3$ dex below the median
and $195$ halos with concentrations $0.5$ dex below the median.

The gas-rich UDGs strongly favor halos with anomalously large $R_{\rm max}$,
and in many cases inferred $R_{\rm max}$ is close to the virial radius
$r_{200}$. This brings up the possibility that the halo density profile could
be close to $1/r^{2}$, which would lead to a flat profile of circular
velocities. To investigate this possibility, we have intended to choose the
prior for the outer slope in a narrow range of $2.05<\beta<2.3$ for the DPL
fits, a significant difference from the NFW one, $\beta=3$. We further take
$1\leq\log(r_{200}/r_{s})\leq 2$, again with the intention of picking
solutions leading to flat halo $V_{\rm circ}$. Fig. 4 (left) shows $R_{\rm
max}$ and $V_{\rm max}$ of the halos, based on the DPL fits (colored), where
$r_{200}$ is reported if $R_{\rm max}>r_{200}$. We see that the inferred
$V_{\rm max}\textup{--}R_{\rm max}$ values are consistent with those from the
fits with the Read profile.

It is also useful to compare the concentration values inferred from the Read
and DPL fits, as summarized in Tables 1 and 2. Both fits show that the halo
concentration has to be order unity, and they follow a similar variation
pattern for individual UDGs. The $r_{s}$ values from the DPL fits are much
smaller than those from the Read fits. Both fits give rise to similar halo
masses owing to the constraint from the cosmological baryon abundance, and the
DPL fits require smaller $r_{s}$ for given $M_{200}$ as the density drops more
slower for $r>r_{s}$, compared to the Read fits.

Based on the fits above, a consistent picture emerges, i.e., the gas-rich UDGs
are hosted by halos with $R_{\rm max}$ much larger than the median. This would
imply anomalously low concentrations if the halo profile is NFW-like, allowing
for a flat density core. On the other hand, equally good fits can be obtained
with halo density profiles that are approximately $1/r^{2}$ and transition to
another power law in the inner regions. In both cases, it is clear that the
UDG hosts are very different from the median halos in this mass range.

##  5 Inferred internal densities compared to simulations

To study the large-$R_{\rm max}$, low-concentration IllustrisTNG halos that
are overabundant, we selected the simulated halos within the red box of Fig. 3
(left), which are similar to those hosting the UDGs, and examined their dark
matter distributions. The bottom edge of the red box corresponds to a
concentration of $0.2$ dex below the median (assuming an NFW profile). Fig. 3
(right) shows the circular velocities at radii $8$ and $2~{}{\rm kpc}$ versus
$V_{\rm max}$ for the halos inferred from the fits (colored), compared to the
simulated ones that are calculated using the mass profiles of the IllustrisTNG
halos (gray circles). The simulated halos shown in Fig. 3 are all well
resolved at $2~{}\rm kpc$ using the Power et al. (2003) criteria. About $4\%$
of the halos have resolution radii close to $2~{}{\rm kpc}$, and in all cases
the resolution radius is smaller than $2.3~{}{\rm kpc}$.

We see that while the simulated halos get close to the required densities at
$8~{}\rm kpc$, they are too dense at $2~{}{\rm kpc}$ compared to the Read
fits. The preference for a low density is mainly because the host halos have a
low concentration, while a cored profile could also contribute to lowering the
central density. All of the UDGs are consistent with core sizes of a few
kiloparsecs, except for AGC 242019. For this UDG, the circular velocity at
$2~{}{\rm kpc}$ is about a factor of $2$ smaller than that expected for the
simulated halos despite being consistent with zero core size. As discussed in
Sec. 2, AGC 242019 allows to have a core with $r_{c}\sim 6~{}{\rm kpc}$,
without including the innermost datapoint at $0.5~{}{\rm kpc}$. Indeed, AGC
242019† has lower halo $V_{\rm circ}$ at $2~{}{\rm kpc}$, and $V_{\rm max}$
becomes smaller as well. Thus, the tension is reduced, but it still remains
compared to the simulated halos.

Fig. 4 (right) shows the inferences from the DPL fits. We remind the reader
that the fits impose a strong prior on the outer logarithmic density slope
$2.05<\beta<2.3$, leading to a flat profile of halo circular velocities. The
DPL profile has the freedom to be steeper than the NFW profile in the inner
region through the parameter $\gamma$; however, we do not see a systematic
shift to $\gamma>1$ in those fits (see Table 2). The inferred halo circular
velocities at $2$ and $8~{}{\rm kpc}$ are higher than those obtained with the
Read fits, resulting in better agreement with the simulated halos. However,
overall the simulated halos are still too dense. For the five low-resolution
UDGs, their halo circular velocities at $2~{}{\rm kpc}$ are systematically
lower than the simulated ones, although the tension is reduced compared to the
Read fits. The high-resolution AGC 242019 and AGC 114905 are still highly
discrepant with the TNG50-1-Dark halos. We expect that the five low-resolution
UDGs would follow a similar trend if future measurements could resolve the
rising part of their inner circular velocity profile.

For both Read and DPL fits, we have imposed a conservative prior $M_{\rm
baryon}/M_{200}\leq 0.187$, the cosmological baryon fraction. This prior has a
large influence on the inferences for the halo parameters. If the halo mass
were allowed to be lower, then $V_{\rm max}$ would be lower and more in line
with the TNG50-1-Dark simulations shown in the right panels of Figs. 3 and 4.
Thus, the tension between the inferred UDG halos and the simulated ones in the
halo $V_{\rm circ}\textup{--}V_{\rm max}$ plane could be related to the prior
on the halo mass from the cosmological baryon fraction. In Sec. 7, we will
show examples after relaxing this constraint.

To amplify this point regarding the baryon mass fraction, we check the full
hydrodynamic TNG50-1 (highest-resolution TNG50 box) simulations to compare
$M_{200}\textup{--}M_{\rm baryon}$ distributions of the simulated galaxies
with those from our fits; see Fig. 5 (gray and red squares). To focus on
extremely gas-rich systems similar to many of the observed field UDGs, where
the ratio of total gas-to-stellar masses ranges from a factor of $3$ to $50$,
we have selected the simulated galaxies in a similar mass range to the UDGs
with a further cut of $M_{\rm gas}/M_{\rm baryon}>0.9$. About $76\%$ of the
TNG50-1 galaxies in the queried mass range satisfy the cut, and their baryon
masses are similar to those plotted in Fig. 5. In those TNG50-1 galaxies, Hi
is the dominant gas component (Diemer et al., 2018).

As shown in Fig. 5, for a given total baryon mass, the inferred $M_{200}$
values of six UDGs are close to the lower end in the $M_{200}\textup{--}M_{\rm
baryon}$ plane, and the simulated galaxies have a systematically higher halo
mass. Given the uncertainty of the halo mass from the fits, we see that some
of the simulated galaxies can be compatible with the observed ones. AGC 242019
favors a high halo mass owing to its sharp rising circular velocities, and
again it becomes more aligned with the other UDGs, as well as the simulated
ones, if the innermost datapoint is not included. The red squares in Fig. 5
denote TNG50-1 halos with the concentration $0.3$ dex below the median. It
seems that there is no correlation between halo concentration and gas content.
Overall, we find that the IllustrisTNG simulations could produce field UDG
analogs in low-concentration halos with high gas content, along with similar
halo and gas masses to those of our sample. But we need an additional
mechanism to reduce the inner densities to be consistent with the
measurements.

Another important factor contributing to the tension between the TNG50-1-Dark
simulations and the inferred circular velocities at $2~{}{\rm kpc}$ (for both
DPL and Read fits) is that some of the simulated halos have a steep density
profile in the inner regions. We find that this is the case for halos with the
largest $R_{\rm max}$, relative to the cosmological median; see Fig. 18 in the
Appendix for illustrative examples. For these halos, the circular velocity at
$2~{}{\rm kpc}$ is higher than it would be for an NFW profile with the same
$V_{\rm max}$ and $R_{\rm max}$. It is likely that the strong deviation from
the NFW profile for these outlier halos is due to the fact that they have not
yet relaxed. As we look at halos less than $-0.3~{}\rm dex$ from the median,
we find that the halo profiles start to look more NFW-like. We will further
discuss these points in the next section.

In short, the internal densities of the UDGs measured at $2~{}{\rm kpc}$ and
$8~{}{\rm kpc}$ are inconsistent with the TNG50-1-Dark halo density profiles.
In particular, the halos of two high-resolution UDGs are significantly
underdense compared to the simulated halos, and the other low-resolution UDGs
also show the same behavior to varying degrees. Thus, in addition to these
inferred halo profiles being clear outliers in terms of their $R_{\rm max}$
values, our results indicate that the TNG50-1-Dark halos are also
systematically overdense within the inner $2~{}{\rm kpc}$ compared to this
sample of field UDGs. This leaves open the question of whether baryonic
feedback or dark matter physics like self-interactions can alleviate the
discrepancy sufficiently.

Figure 5: Inferred $M_{200}$ vs. $M_{\rm baryon}$ for the UDGs (colored with
error bars), compared to selected TNG50-1 galaxies satisfying $M_{\rm
gas}/M_{\rm baryon}>0.9$ in a similar mass range (full hydro; gray and red
squares). The dashed line denotes the lower limit on the halo mass from the
cosmological ratio of baryonic to dark matter masses. The red squares denote
low-concentration halos satisfying $\log_{10}(c_{200}/c^{\rm
median}_{200})<-0.3$, see Fig. 6.

Figure 6: Left: The histogram of “halo concentration” of the TNG50-1-Dark
halos with $\log_{10}(M_{200}/{\rm M_{\odot}})=9.3\textup{--}11$ inferred
assuming an NFW profile, compared to the extrapolation with a lognormal
distribution (Diemer & Joyce, 2019) (dashed black). The red and blue colors
highlight two subsets of the halos satisfying $\log_{10}(c_{200}/c^{\rm
median}_{200})<-0.3$ and $|\log_{10}(c_{200}/c^{\rm median}_{200})|<0.1$,
respectively. Note that many of the systems in the tail have density profiles
that fall much slower than $1/r^{3}$ in the outer regions leading to larger
$R_{\rm max}$ and hence lower $c_{200}$. Middle: the histogram of the redshift
at which half of the halo mass is assembled, for the low (red) and median
(blue) concentration halos. Right: the histogram of the total halo spin
parameter for the low (red) and median (blue) concentration halos.

##  6 Late-forming halos

We have seen that halos with large $R_{\rm max}$ required to host the UDGs are
predicted in hierarchical structure formation theory, and their abundance is
much higher than expected from a Gaussian tail. In this section, we look
further into these IllustrisTNG halos. We use the TNG50-1-Dark dataset since
our focus will be on the concentration, spin, and assembly history of these
halos. In Fig. 6 (left), we show the number fraction of TNG50-1-Dark halos
versus the normalized concentration (red), compared to a Gaussian
extrapolation (dashed black). Indeed, the number of halos is not exponentially
suppressed toward low $c_{200}$. Instead, it follows a power-law behavior in
$|\log(c_{200}/c^{\rm median}_{200})|$ to low concentrations. For
$\log(c_{200}/c^{\rm median}_{200})=-0.4$, the Gaussian extrapolation,
corresponding to $\sim 4\sigma$ below the median, underestimates the halo
population by two orders of magnitude. The tail in $c_{200}$ values is a
direct reflection for the tail in the $R_{\rm max}$ distribution shown in Fig.
3 (left).

Before we delve further, a few comments on interpreting $c_{200}$ may be
useful. When plotting the TNG50-1-Dark halos in Fig. 6 (left), we have
computed their $c_{200}$ and $M_{200}$ from $V_{\rm max}$ and $R_{\rm max}$
assuming an NFW profile. Thus, $c_{200}$ and $M_{200}$ should be seen as
convenient proxies for $V_{\rm max}$ and $R_{\rm max}$ from the simulations.
There are other ways of computing the concentration that are profile
independent, for example, $\Delta_{V/2}$, the average density within the
radius where the circular velocity is $V_{\rm max}/2$ in units of the critical
density (Alam et al., 2002). Given the fact that the halos with the lowest
$c_{200}$ seem to have a density profile steeper than an NFW one in the inner
regions, a measure like $\Delta_{V/2}$ may not identify them as being
outliers. On the other hand, they are outliers because the NFW density profile
is a poor fit, unlike the median halos. In this work, we have chosen to retain
the simplicity of $c_{200}$ for quantitatively identifying outliers.

Fig. 6 (middle) shows the histograms of redshifts at which half of the halo
mass is assembled for the TNG50-1-Dark halos with $\log_{10}(c_{200}/c^{\rm
median}_{200})<-0.3$ (red) and $|\log_{10}(c_{200}/c^{\rm median}_{200})|<0.1$
(blue). We see that low-concentration (large $R_{\rm max}$) halos assemble
systematically later. The formation redshift peaks around $z\sim 0.5$ and
$1.5$ for low and median concentration halos, respectively. Thus, the large
deviations from the NFW profile apparent in density profiles of many of these
halos may be tied to the late formation history of these halos and the fact
that these halos may be less relaxed than median halos that formed earlier.
Our fits to the circular velocity data support the hypothesis that the field
UDGs may form in these halos with larger-than-expected $R_{\rm max}$ values.
If so, they would also be forming later than the galaxies hosted in the median
halos of the same mass. This opens up the possibility that their star
formation histories (compared to more abundant galaxies) could shed light on
this issue, although the small numbers of UDGs and the inherent spread in the
formation times seen in Fig. 6 (middle) may make it difficult to use this as a
discriminant.

The late-forming halos also have a high spin, as demonstrated in Fig. 6
(right). We calculate the spin parameter as $\lambda_{\rm
halo}=J/(\sqrt{2}M_{200}V_{200}R_{200})$, where $J$ is the total angular
momentum and $V_{200}=\sqrt{GM_{200}/R_{200}}$ (Bullock et al., 2001). It is
clear that the low-concentration halos have higher spin parameters on average
than the median halos (see also Neto et al. (2007); Macciò et al. (2007). This
may be relevant for creating the extended gas distribution of UDGs (Amorisco &
Loeb, 2016), although AGC 114905 and AGC 242019 seem to have a normal baryonic
specific angular momentum (Mancera Piña et al., 2021). It would be interesting
to take the late-forming IllustrisTNG halos we have identified and perform a
detailed study on halo and baryonic angular momenta for the gas-rich UDGs. We
note that the histograms for the spin parameters are wide with a significant
overlap. To ensure that they are indeed two distinct halo populations, we
conducted a Kolmogorov–Smirnov test and found that the p-value is $\sim
10^{-10}$.

Figure 7: The fraction of the dark matter mass within the virial radius that
is contained in the halo (Primary Flag=1) vs normalized halo concentration for
the TNG50-1-Dark halos. The red and blue colors denote the halos satisfying
$\log_{10}(c_{200}/c^{\rm median}_{200})<-0.3$ and
$\left|\log_{10}(c_{200}/c^{\rm median}_{200})\right|<0.1$, respectively.

For the gas-rich UDGs we consider, there is little evidence that they
experienced a recent major merger. To check the merger history of the
TNG50-1-Dark halos (Primary Flag=1), we obtain the fraction of the halo mass
over the total group mass within the virial radius (Anbajagane et al., 2021),
as shown in Fig. 7. Overall, the median (blue) and low-concentration (red)
halos are similar in the halo mass fraction. We note that all halos with
$\log_{10}(c_{200}/c^{\rm median}_{200})\lesssim-0.5$ have a mass fraction
less than $1$. Thus, there is a tendency that the extremely low concentration
halos are more likely to have a relatively massive neighbor, an indication of
late major mergers. However, the majority of the low-concentration halos we
have identified have a mass fraction larger than $\sim 0.8\textup{--}0.9$, and
those halos could be “quiet” enough to host the gas-rich UDGs discussed in
this work.

To summarize this section, the IllustrisTNG halos with concentrations below
about $0.3~{}\rm dex$ from the median are promising candidates for hosting the
UDGs. Their low formation redshifts and high spins also seem to be relevant
for setting the properties of the extended gas and stellar distributions of
the UDGs. However, the internal densities in the simulated halos are
systematically higher than what is observed in the seven field UDGs that we
have analyzed.

##  7 Possible solutions

Figure 8: Circular velocity profiles (solid black) constructed from example
halos in TNG50-1-Dark simulations (solid red) and the baryonic contribution
(dashed magenta) compared to observational data (black circles). Each panel
also shows the halo-to-baryon mass ratio in the upper right corner. All other
details are the same as in Fig. 2. Figure 9: Dark matter density profiles of
the TNG50-1-Dark halos adapted in Fig. 8 (solid red) vs the fits using the DPL
density profile in Eq. 5 (dashed black). The parameter values of the DPL
profile are shown in each panel.

Our analysis indicates that finding a consistent solution for the UDGs will
require some modifications of the inner regions of the simulated halos. Such
modifications are constrained by the fact that these galaxies are gas-rich and
the ratio of baryon to halo masses is constrained in our fits to be close to
the cosmological baryon fraction. Here we discuss possible ways in which the
TNG50-1-Dark halos we have identified could be modified to provide a better
fit to the circular velocity profiles of the field UDGs.

To further highlight the tension and modifications that are needed, we have
tried to find TNG50-1-Dark halos with circular velocity profiles that are
close to those of the UDGs, and the results are shown in Fig. 8. Note that we
do not perform a detailed fit here. Instead, we identify example TNG50-1-Dark
halos that seem consistent with the observed sample in terms of the required
$V_{\rm max}$. This exercise reveals critical issues exemplified by the
comparisons to the two well-resolved circular velocity profiles. For AGC
242019, the ratio of simulated halo to observed baryon masses is
$M_{200}/M_{\rm baryon}\approx 8$, above the lower limit expected from the
cosmological baryon fraction. This is not surprising, as this UDG is dark
matter dominated. However, the inner density of the simulated halo is too high
for AGC 242019.

AGC 114905 is even more difficult to explain, as the gas content within the
inner $10~{}{\rm kpc}$ dominates over dark matter and stellar disk. In order
to get low enough $V_{\rm max}$, we had to pick a halo with
$M_{200}=1.74M_{\rm baryon}$, a factor of $3$ below the expected halo mass
from the cosmological baryon fraction, which is already very conservative.
Even with such a small halo mass, the inner density needs to be further
reduced. This is also clear from the inferred core size $r_{c}\approx
6.7~{}{\rm kpc}$ and the inner logarithmic slope $\gamma\approx 0.09$ from the
Read and DPL fits, respectively. It should be noted that a low dark matter
density could be in tension with disk stability requirements unless the
inclination is lower or the gas dispersion is higher (Mancera Piña et al.,
2022; Sellwood & Sanders, 2022). Deeper Hi observations and more refined
simulations are needed to further address the stability issue.

For the five low-resolution UDGs, the circular velocity profiles of the
simulated halos are in good agreement with the observed ones. However, since
they have only two datapoints, it is difficult to make an accurate assessment.
Furthermore, for AGC 334315, AGC 749290, and AGC 248945, the simulated halo
masses are too low to be consistent with the cosmological limit, as in the
case of AGC 114905. It would be fascinating to see what will be revealed when
higher-resolution measurement data become available for these UDGs.

Fig. 8 also shows that the circular velocity profiles of the simulated halos
are surprisingly flat for $r\gtrsim 2~{}{\rm kpc}$, indicating that their
density profiles significantly deviate from the NFW profile. In fact, the DPL
profile in Eq. 5 can provide a good fit as shown in Fig. 9. For the sample of
five simulated halos we consider, we find the transition radius $r_{s}\approx
1\textup{--}2~{}{\rm kpc}$, $\gamma\sim 1$, and $\beta\approx 2.3$,
drastically different from an NFW halo with a low concentration. The result
further indicates that “low-concentration” TNG50-1-Dark halos we have
identified are not those expected from a simple extrapolation of the
concentration-mass relation with an NFW profile. Thus, our application of the
DPL profile to fit the UDGs is well motivated and justified. For the five
halos, the fraction of the halo mass over the group mass within $r_{200}$
varies in the range of $0.8\textup{--}0.95$. Thus, we expect that they did not
experience late major mergers.

Figure 10: Measured circular velocity profiles of seven isolated UDGs studied
in this work (color-coded with error bars), compared to those predicted in
NIHAO simulations (gray curves), based on the total mass profiles of simulated
UDGs in Di Cintio et al. (2016), where $\log_{10}(M_{\rm halo}/{\rm
M_{\odot}})=10.22\textup{--}10.85$ and
$\log_{10}(M_{\rm\footnotesize\textsc{Hi}}/{\rm
M_{\odot}})=7.22\textup{--}9.24$.

Di Cintio et al. (2016) suggested that field UDGs could form in more normal
halos using the simulations from project NIHAO (Numerical Investigation of a
Hundred Astrophysical Objects); see also Chan et al. (2018). These simulations
produce extended stellar distributions, and some of the simulated galaxies are
gas-rich as well. However, we find that the circular velocity profiles of the
simulated UDGs in NIHAO do not match the observed UDGs in our sample as shown
in Fig. 10. 111Preliminary measurements in Leisman et al. (2017) suggested
that AGC 122966, AGC 219533, and AGC 334315 had similar rotation velocities as
NIHAO UDGs, but more detailed and accurate kinematic modeling in Mancera Piña
et al. (2019, 2020) reveals that their circular velocities are much lower. In
addition, the simulated NIHAO UDGs have overall lower gas fractions, compared
to the observed gas-rich UDGs we consider. Brook et al. (2021) used a similar
DPL profile to Eq. 5 and highlighted that AGC 242019 could be hosted by a
“median” halo that has been modified by strong feedback. However, our fits
show that just modifying the central density is not sufficient; large $R_{\rm
max}$ is also required to get the correct $M_{\rm baryon}/M_{200}$ ratio. In
fact, their best-fit model has $\beta\approx 2.15$, resulting in $R_{\rm
max}\approx 261~{}{\rm kpc}$ if the density profile is extrapolated,
consistent with our results.

If the gas-rich UDGs do indeed form in halos that are far from the median,
then it is interesting to ask what aspects of strong feedback might change
compared to the simulations based on median halos (Di Cintio et al., 2016;
Chan et al., 2018). For example, the dynamical time $t_{\rm dyn}=2\pi
r_{s}/V_{\rm circ}(r_{s})$ is a factor of $2\textup{--}5$ larger for
$10^{10}\textup{--}10^{10.5}~{}{\rm M_{\odot}}$ halos with a concentration of
$0.3\textup{--}0.6$ dex below, compared to the median halos. Thus, the low-
concentration halos have a weaker potential, and it may take longer time to
replenish gas necessary for bursty start formation, and the efficiency of
feedback could be decreased. In addition, these low-concentration halos formed
at $z\lesssim 1$ (see Fig. 6), which reduces the time for star formation.
Based on these arguments, we expect that feedback will have a lesser impact on
the field UDGs compared to galaxies hosted by median concentration halos.
Whether the reduced feedback efficiency is sufficient to explain the circular
velocity profiles of the observed field UDGs, with gas fractions close to the
cosmological value, remains to be seen.

Another possibility to consider is self-interacting dark matter (SIDM; see
Tulin & Yu (2018) for a review). The self-interactions thermalize the inner
halo and heat up dark matter particles, thereby reducing the inner halo
density. Studies have shown that SIDM can explain the diverse rotation curves
over a wide mass range of spiral galaxies (Ren et al., 2019; Kaplinghat et
al., 2020) and the properties of UDGs in clusters and groups (Carleton et al.,
2019; Yang et al., 2020). For field UDGs, self-interactions would need to
create constant-density core sizes of few kiloparsecs to be able to lower the
circular velocity at $\sim 2~{}{\rm kpc}$; see Fig. 3 (right).

The low “concentrations” inferred in the UDGs imply that the radius out to
which their host halos are thermalized (Kaplinghat et al., 2016) would be a
small fraction of $r_{s}$ for self-interacting cross sections over mass of
order $1~{}\rm cm^{2}~{}g^{-1}$, which have been used to explain the rotation
curves of spiral galaxies (Ren et al., 2019). Conversely, we expect core size
constraints in the UDGs to provide significant constraints on SIDM models. In
this regard, the two high-resolution UDGs AGC 114905 and AGC 242019 present
different challenges. AGC 114905 is completely dominated by gas, which will
need to be folded in to determine the SIDM density profile (Kaplinghat et al.,
2014), while the core size inference in AGC 242019 is highly dependent on the
innermost data point. More well-resolved circular velocity curves for the UDGs
are required to make further progress on this front.

An interesting connection is that some of the field UDGs are close to the
higher end of the velocity range relevant for the Milky Way dwarf spheroidal
satellites, approximately $5\textup{--}30~{}\rm km~{}s^{-1}$. Recent studies
have argued for cross sections much larger than $1~{}\rm cm^{2}~{}g^{-1}$ at
those velocities based on the observed diversity in the central dark matter
densities of the satellites (Nishikawa et al., 2020; Zavala et al., 2019;
Kahlhoefer et al., 2019; Sameie et al., 2020; Correa, 2021; Jiang et al.,
2021). Thus, in conjunction with the Milky Way satellites, we expect the field
gas-rich UDGs to be a strong test of SIDM models.

##  8 Conclusions

We have analyzed the circular velocity profiles of seven gas-rich UDGs in the
field and inferred properties of their host dark matter halos. The preferred
halo mass is $M_{200}\sim 10^{10}\textup{--}10^{11}~{}{\rm M_{\odot}}$
($V_{\rm max}\sim 30\textup{--}50~{}{\rm km~{}s^{-1}}$), and the resulting
ratio of baryon-to-halo masses is $\sim 0.07$. The baryons in these UDGs are
mostly in the form of Hi gas, with the ratio of gas to stellar masses in the
range of $3\textup{--}50$.

We have argued in this paper that simulated halos with $c_{200}$ values about
$0.3$ dex or lower than the median seem to be the best possible hosts for the
field UDGs we have analyzed. Remarkably, we find that such low-concentration
(large $R_{\rm max}$) halos are formed in cold dark matter models (and by
extension other models with a similar initial power spectrum). We have
demonstrated this using the IllustrisTNG simulations, in particular the dark-
matter-only run, and showed that the distribution of halo concentrations has a
large non-Gaussian tail to low concentrations. We have identified that this
population of simulated hosts, on average, forms later than the median halos
and has higher spins. It is possible that these features may be tied to the
unusual properties of the gas-rich UDGs.

The gas-rich field UDGs we have analyzed are outliers in terms of their
inferred halo properties. We emphasize that the host halos are low-
concentration in the sense of having large $R_{\rm max}$, but their density
profiles could depart drastically from the NFW profile in both inner and outer
regions. The low dark matter densities and the odd halo properties of these
gas-rich field UDGs provide an exciting opportunity to test both feedback
models and dark matter models. To take advantage of this opportunity, we need
to discover more field UDGs and obtain well-resolved circular velocity
profiles for them.

## Acknowledgements

We thank Volker Springel and collaborators for the IllustrisTNG simulation
data and the user-friendly data access interface, and Salvador Cardona Barrero
for clarifications on the NIHAO data. M.K. is supported by the NSF under award
No. 1915005. H.-B.Y. is supported by the U.S. Department of Energy under grant
No. de-sc0008541, the John Templeton Foundation under grant ID #61884, and
NASA under grant No. 80NSSC20K0566. F.F. and P.E.M.P. are supported by the
Netherlands Research School for Astronomy (Nederlandse Onderzoekschool voor
Astronomie, NOVA), Network 1, Project 10.1.5.6. The opinions expressed in this
publication are those of the authors and do not necessarily reflect the views
of the John Templeton Foundation.

AGC ID |  $V_{\rm sys}$ (km/s) |  $D$ (Mpc) |  $i$ (deg) | References
---|---|---|---|---
114905 | $5435$ | $67\pm 1$ | $26.1^{+0.2}_{-0.1}$ | Leisman et al. (2017); Mancera Piña et al. (2020, 2022)
122966 | $6509$ | $90\pm 5$ | $34\pm 5$ | Leisman et al. (2017); Mancera Piña et al. (2020)
219533 | $6384$ | $96\pm 5$ | $42\pm 5$ | Leisman et al. (2017); Mancera Piña et al. (2020)
248945 | $5703$ | $84\pm 5$ | $66\pm 5$ | Leisman et al. (2017); Mancera Piña et al. (2020)
334315 | $5107$ | $73\pm 5$ | $45\pm 5$ | Leisman et al. (2017); Mancera Piña et al. (2020)
749290 | $6516$ | $97\pm 5$ | $39\pm 5$ | Leisman et al. (2017); Mancera Piña et al. (2020)
242019 | $1840.4$ | $30.8\pm 1.5$ | $73.0\textup{--}70.2$ | Shi et al. (2021)
Table 3:  Additional information about properties of the gas-rich UDGs We
Consider. From left to right columns: galaxy AGC ID, systemic velocity,
distance, inclination, source of the data. For AGC 114905, Mancera Piña et al.
(2022) reports another set of distance and inclination values $D=76\pm
5~{}{\rm Mpc}$ and $i=32\pm 3^{\rm o}$. In our fits, we take the values listed
in the table for AGC 114905, see Mancera Piña et al. (2022) for details about
measuring the inclination. For AGC 242019, Shi et al. (2021) reports that the
inclination varies from $i=73.0^{\rm o}$ to $70.2^{\rm o}$ for the distance in
the range $0.67\textup{--}8.75~{}{\rm kpc}$. Figure 11: Hi gas surface mass
density profiles of the low-resolution UDGs for two sets of model parameters
$n=4$ (solid) and $2.5$ (dashed), where $m=4$ in both cases; see Eq. 1. The
central surface density $\Sigma_{\rm\footnotesize\textsc{Hi}}(0)$ is set by
demanding the predicted Hi mass within $R_{\rm out}$ to be the median value of
Hi mass in Gault et al. (2021), as denoted in each panel. Our MCMC routine
assumes a flat prior for $M_{\rm\footnotesize\textsc{Hi}}$, and thus
$\Sigma_{\rm\footnotesize\textsc{Hi}}(0)$ varies accordingly. Figure 12:
Similar to the left panel of Fig. 3 for the low-resolution UDGs, but with a
gas surface density profile of $n=2.5$ and $m=4$ (dashed error bars). Compared
to the case with $n=4$ and $m=4$ in Fig. 3, the inferred $R_{\rm max}$ values
are somewhat smaller.  Figure 13: Circular velocity profiles of the best fit
model (solid black), including halo (solid red) and baryonic (dashed magenta)
contributions compared to observational data (black dots), where we assume the
NFW density profile. For AGC 242019†, the fit is performed without including
the innermost datapoint. The observational data are taken from Mancera Piña et
al. (2022, 2020); Shi et al. (2021).

Figure 14: Left: similar to the left panel of Fig. 3, but for fits with the
NFW density profile. The inferred $V_{\rm max}$ and $R_{\rm max}$ values are
similar to those from the Read fits. For AGC 114905 and AGC 242019† (overlaps
with AGC 242019 in the $R_{\rm max}\textup{--}V_{\rm max}$ plane), the
difference is relatively large, because the halo is not cored in this case.
Right: $c_{200}\textup{--}M_{200}$ distributions from the NFW fits,
complementary to the $R_{\rm max}\textup{--}V_{\rm max}$ distributions shown
in the left panel. Figure 15: Circular velocity profiles of the best fit model
(solid black), including halo (solid red) and baryonic (dashed magenta)
contributions compared to observational data (black circles), where we assume
the general Read profile. For AGC 242019†, the fit is performed without
including the innermost datapoint. The observational data are taken from
Mancera Piña et al. (2022, 2020) and Shi et al. (2021).

Figure 16: Left: inferred $r_{c}/r_{s}$ and $\delta$ distributions from the
fits with the general Read profile. Right: similar to the left panel of Fig.
3, but with the general Read profile. Figure 17: Circular velocity profiles of
the best-fit model (solid black), including halo (solid red) and baryonic
(dashed magenta) contributions compared to observational data (black circles).
The halo profile is the DPL profile in Eq. 5. For AGC 242019†, the fit is
performed without including the innermost datapoint. The observational data
are taken from Mancera Piña et al. (2022, 2020) and Shi et al. (2021). For the
inferred parameter errors, see Table 2.

Figure 18: Left: dark matter density profiles of example $M_{200}\approx
3.2\times 10^{10}~{}{\rm M_{\odot}}$ TNG50-1-Dark halos with a median
concentration (blue) and $0.3$ dex (gray) and $0.6$ dex (red) below. The solid
curves are computed from the particle data, while the dashed curves are NFW
profiles using the catalog halo parameters $(V_{\rm max},~{}R_{\rm max})$:
$(56.9~{}{\rm km~{}s^{-1}},~{}11.9~{}{\rm kpc})$ (blue), $(49.1~{}{\rm
km~{}s^{-1}},~{}24.3~{}{\rm kpc})$ (gray), and $(46.2~{}{\rm
km~{}s^{-1}},~{}48.6~{}{\rm kpc})$ (red). As expected, the density profile of
the median halo is NFW-like, but it becomes steeper for the low concentration
halos. The halo with the lowest concentration (red) has a logarithmic slope of
$-1.4$ within $\sim 8~{}{\rm kpc}$, which is approximately
$R_{\rm\footnotesize\textsc{Hi}}$ of the UDGs. The shaded region indicates the
resolution limit. Right: the corresponding circular velocity profiles for the
halos shown in the left panel computed from their simulated mass distribution
(solid), compared to the calculated ones assuming an NFW profile (dashed).

In this appendix, we provide additional information about properties of the
gas-rich UDGs, the Hi gas surface density for the low-resolution UDGs, the
fits using the NFW, general Read and DPL profiles, inferred $R_{\rm
max}\textup{--}V_{\rm max}$ distributions using a different Hi gas surface
density profile, and dark matter distributions of three representative
IllustrisTNG halos.

Table 3 provides addition information about the gas-rich UDGs considered in
this work, including their systemic velocity, distance, and inclination.

Fig. 11 shows Hi gas surface mass density profiles of the low-resolution UDGs
for two sets of model parameters $n=4$ and $2.5$. We fix $m=4$ in both cases.

Fig. 12 shows the inferred $R_{\rm max}\textup{--}V_{\rm max}$ distributions
for the low-resolution UDGs, similar to the left panel of Fig. 3, but using a
gas surface density profile with $n=2.5$ and $m=4$. See Fig. 11 for the
difference in the gas surface density for $n=4$ and $2.5$.

Fig. 13 shows circular velocity profiles of the best fit model for the NFW
profile. The inferred $R_{\rm max}\textup{--}V_{\rm max}$ and
$c_{200}\textup{--}M_{200}$ distributions are shown in Fig. 14.

Fig. 15 shows circular velocity profiles of the best fit model for the general
Read profile. The inferred $r_{c}/r_{s}$ and $\delta$ distributions, as well
as the $R_{\rm max}\textup{--}V_{\rm max}$ distributions, are shown in the
left and right panels of Fig. 16, respectively.

Fig. 17 shows circular velocity profiles of the best-fit model for the DPL
profile in Eq. 5. The model parameters are summarized in Table 2, and the
inferred $R_{\rm max}\textup{--}V_{\rm max}$ distributions are shown in Fig.
4.

Fig. 18 shows dark matter density profiles of three TNG50-1-Dark halos with
three different concentrations (left panel) and their corresponding profiles
of circular velocities (right panel).

## References

  * Alam et al. (2002) Alam, S. M. K., Bullock, J. S., & Weinberg, D. H. 2002, Astrophys. J., 572, 34, doi: 10.1086/340190
  * Amorisco & Loeb (2016) Amorisco, N. C., & Loeb, A. 2016, MNRAS, 459, L51, doi: 10.1093/mnrasl/slw055
  * Anbajagane et al. (2021) Anbajagane, D., Evrard, A. E., & Farahi, A. 2021, Mon. Not. Roy. Astron. Soc., 509, 3441, doi: 10.1093/mnras/stab3177
  * Benavides et al. (2021) Benavides, J. A., Sales, L. V., Abadi, M. G., et al. 2021, Nature Astronomy, 5, 1255, doi: 10.1038/s41550-021-01458-1
  * Bennet et al. (2017) Bennet, P., Sand, D. J., Crnojević, D., et al. 2017, ApJ, 850, 109, doi: 10.3847/1538-4357/aa9180
  * Binggeli et al. (1984) Binggeli, B., Sandage, A., & Tarenghi, M. 1984, AJ, 89, 64, doi: 10.1086/113484
  * Binney & Tremaine (2008) Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition
  * Brook et al. (2021) Brook, C. B., Di Cintio, A., Maccio, A. V., & Blank, M. 2021, Astrophys. J. Lett., 919, L1, doi: 10.3847/2041-8213/ac236a
  * Bullock et al. (2001) Bullock, J. S., Dekel, A., Kolatt, T. S., et al. 2001, ApJ, 555, 240, doi: 10.1086/321477
  * Carleton et al. (2019) Carleton, T., Errani, R., Cooper, M., et al. 2019, MNRAS, 485, 382, doi: 10.1093/mnras/stz383
  * Chan et al. (2018) Chan, T. K., Kereš, D., Wetzel, A., et al. 2018, MNRAS, 478, 906, doi: 10.1093/mnras/sty1153
  * Collins et al. (2020) Collins, M. L. M., Tollerud, E. J., Rich, R. M., et al. 2020, MNRAS, 491, 3496, doi: 10.1093/mnras/stz3252
  * Conselice (2018) Conselice, C. J. 2018, Research Notes of the American Astronomical Society, 2, 43, doi: 10.3847/2515-5172/aab7f6
  * Correa (2021) Correa, C. A. 2021, MNRAS, 503, 920, doi: 10.1093/mnras/stab506
  * Di Cintio et al. (2016) Di Cintio, A., Brook, C. B., Dutton, A. A., et al. 2016, Monthly Notices of the Royal Astronomical Society: Letters, 466, L1, doi: 10.1093/mnrasl/slw210
  * Diemer & Joyce (2019) Diemer, B., & Joyce, M. 2019, ApJ, 871, 168, doi: 10.3847/1538-4357/aafad6
  * Diemer et al. (2018) Diemer, B., Stevens, A. R. H., Forbes, J. C., et al. 2018, ApJS, 238, 33, doi: 10.3847/1538-4365/aae387
  * Du et al. (2020) Du, W., Cheng, C., Zheng, Z., & Wu, H. 2020, The Astronomical Journal, 159, 138, doi: 10.3847/1538-3881/ab6efb
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Freeman (1970) Freeman, K. C. 1970, ApJ, 160, 811, doi: 10.1086/150474
  * Gannon et al. (2020) Gannon, J. S., Forbes, D. A., Romanowsky, A. J., et al. 2020, MNRAS, 495, 2582, doi: 10.1093/mnras/staa1282
  * Gault et al. (2021) Gault, L., Leisman, L., Adams, E. A. K., et al. 2021, ApJ, 909, 19, doi: 10.3847/1538-4357/abd79d
  * Greco et al. (2018) Greco, J. P., Goulding, A. D., Greene, J. E., et al. 2018, ApJ, 866, 112, doi: 10.3847/1538-4357/aae0f4
  * Guo et al. (2020) Guo, Q., Hu, H., Zheng, Z., et al. 2020, Nature Astronomy, 4, 246, doi: 10.1038/s41550-019-0930-9
  * He et al. (2019) He, M., Wu, H., Du, W., et al. 2019, ApJ, 880, 30, doi: 10.3847/1538-4357/ab2710
  * Impey et al. (1988) Impey, C., Bothun, G., & Malin, D. 1988, ApJ, 330, 634, doi: 10.1086/166500
  * Iodice et al. (2020) Iodice, E., Cantiello, M., Hilker, M., et al. 2020, A&A, 642, A48, doi: 10.1051/0004-6361/202038523
  * Jackson et al. (2021) Jackson, R. A., Martin, G., Kaviraj, S., et al. 2021, MNRAS, 502, 4262, doi: 10.1093/mnras/stab077
  * Janowiecki et al. (2019) Janowiecki, S., Jones, M. G., Leisman, L., & Webb, A. 2019, MNRAS, 490, 566, doi: 10.1093/mnras/stz1868
  * Janssens et al. (2017) Janssens, S., Abraham, R., Brodie, J., et al. 2017, ApJ, 839, L17, doi: 10.3847/2041-8213/aa667d
  * Jiang et al. (2019) Jiang, F., Dekel, A., Freundlich, J., et al. 2019, MNRAS, 487, 5272, doi: 10.1093/mnras/stz1499
  * Jiang et al. (2021) Jiang, F., Kaplinghat, M., Lisanti, M., & Slone, O. 2021.  https://arxiv.org/abs/2108.03243
  * Kahlhoefer et al. (2019) Kahlhoefer, F., Kaplinghat, M., Slatyer, T. R., & Wu, C.-L. 2019, JCAP, 12, 010, doi: 10.1088/1475-7516/2019/12/010
  * Kaplinghat et al. (2014) Kaplinghat, M., Keeley, R. E., Linden, T., & Yu, H.-B. 2014, Phys. Rev. Lett., 113, 021302, doi: 10.1103/PhysRevLett.113.021302
  * Kaplinghat et al. (2020) Kaplinghat, M., Ren, T., & Yu, H.-B. 2020, J. Cosmology Astropart. Phys, 2020, 027, doi: 10.1088/1475-7516/2020/06/027
  * Kaplinghat et al. (2016) Kaplinghat, M., Tulin, S., & Yu, H.-B. 2016, Phys. Rev. Lett., 116, 041302, doi: 10.1103/PhysRevLett.116.041302
  * Koda et al. (2015) Koda, J., Yagi, M., Yamanoi, H., & Komiyama, Y. 2015, ApJ, 807, L2, doi: 10.1088/2041-8205/807/1/L2
  * Lee et al. (2020) Lee, J. H., Kang, J., Lee, M. G., & Jang, I. S. 2020, ApJ, 894, 75, doi: 10.3847/1538-4357/ab8632
  * Leisman et al. (2017) Leisman, L., Haynes, M. P., Janowiecki, S., et al. 2017, ApJ, 842, 133, doi: 10.3847/1538-4357/aa7575
  * Liao et al. (2019) Liao, S., Gao, L., Frenk, C. S., et al. 2019, MNRAS, 490, 5182, doi: 10.1093/mnras/stz2969
  * Lovell et al. (2018) Lovell, M. R., Pillepich, A., Genel, S., et al. 2018, MNRAS, 481, 1950, doi: 10.1093/mnras/sty2339
  * Macciò et al. (2007) Macciò, A. V., Dutton, A. A., van den Bosch, F. C., et al. 2007, MNRAS, 378, 55, doi: 10.1111/j.1365-2966.2007.11720.x
  * Mancera Piña et al. (2022) Mancera Piña, P. E., Fraternali, F., Oosterloo, T., et al. 2022, MNRAS, 512, 3230, doi: 10.1093/mnras/stab3491
  * Mancera Piña et al. (2018) Mancera Piña, P. E., Peletier, R. F., Aguerri, J. A. L., et al. 2018, MNRAS, 481, 4381, doi: 10.1093/mnras/sty2574
  * Mancera Piña et al. (2021) Mancera Piña, P. E., Posti, L., Pezzulli, G., et al. 2021, A&A, 651, L15, doi: 10.1051/0004-6361/202141574
  * Mancera Piña et al. (2019) Mancera Piña, P. E., Fraternali, F., Adams, E. A. K., et al. 2019, ApJ, 883, L33, doi: 10.3847/2041-8213/ab40c7
  * Mancera Piña et al. (2020) Mancera Piña, P. E., Fraternali, F., Oman, K. A., et al. 2020, MNRAS, 495, 3636, doi: 10.1093/mnras/staa1256
  * Martin et al. (2019) Martin, G., Kaviraj, S., Laigle, C., et al. 2019, MNRAS, 485, 796, doi: 10.1093/mnras/stz356
  * Martín-Navarro et al. (2019) Martín-Navarro, I., Romanowsky, A. J., Brodie, J. P., et al. 2019, MNRAS, 484, 3425, doi: 10.1093/mnras/stz252
  * Martínez-Delgado et al. (2016) Martínez-Delgado, D., Läsker, R., Sharina, M., et al. 2016, AJ, 151, 96, doi: 10.3847/0004-6256/151/4/96
  * Moreno et al. (2022) Moreno, J., Danieli, S., Bullock, J. S., et al. 2022, Nature Astronomy, 6, 496, doi: 10.1038/s41550-021-01598-4
  * Nelson et al. (2019) Nelson, D., Springel, V., Pillepich, A., et al. 2019, Computational Astrophysics and Cosmology, 6, 2, doi: 10.1186/s40668-019-0028-x
  * Neto et al. (2007) Neto, A. F., Gao, L., Bett, P., et al. 2007, Mon. Not. Roy. Astron. Soc., 381, 1450, doi: 10.1111/j.1365-2966.2007.12381.x
  * Nishikawa et al. (2020) Nishikawa, H., Boddy, K. K., & Kaplinghat, M. 2020, Phys. Rev. D, 101, 063009, doi: 10.1103/PhysRevD.101.063009
  * Ogiya (2018) Ogiya, G. 2018, MNRAS, 480, L106, doi: 10.1093/mnrasl/sly138
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6, doi: 10.1051/0004-6361/201833910
  * Power et al. (2003) Power, C., Navarro, J. F., Jenkins, A., et al. 2003, MNRAS, 338, 14, doi: 10.1046/j.1365-8711.2003.05925.x
  * Prole et al. (2019) Prole, D. J., van der Burg, R. F. J., Hilker, M., & Davies, J. I. 2019, MNRAS, 488, 2143, doi: 10.1093/mnras/stz1843
  * Prole et al. (2021) Prole, D. J., van der Burg, R. F. J., Hilker, M., & Spitler, L. R. 2021, MNRAS, 500, 2049, doi: 10.1093/mnras/staa3296
  * Read et al. (2016) Read, J. I., Agertz, O., & Collins, M. L. M. 2016, MNRAS, 459, 2573, doi: 10.1093/mnras/stw713
  * Ren et al. (2019) Ren, T., Kwa, A., Kaplinghat, M., & Yu, H.-B. 2019, Physical Review X, 9, 031020, doi: 10.1103/PhysRevX.9.031020
  * Román et al. (2019) Román, J., Beasley, M. A., Ruiz-Lara, T., & Valls-Gabaud, D. 2019, MNRAS, 486, 823, doi: 10.1093/mnras/stz835
  * Román & Trujillo (2017) Román, J., & Trujillo, I. 2017, MNRAS, 468, 4039, doi: 10.1093/mnras/stx694
  * Sales et al. (2020) Sales, L. V., Navarro, J. F., Peñafiel, L., et al. 2020, MNRAS, 494, 1848, doi: 10.1093/mnras/staa854
  * Sameie et al. (2020) Sameie, O., Yu, H.-B., Sales, L. V., Vogelsberger, M., & Zavala, J. 2020, Phys. Rev. Lett., 124, 141102, doi: 10.1103/PhysRevLett.124.141102
  * Sellwood & Sanders (2022) Sellwood, J. A., & Sanders, R. H. 2022, arXiv e-prints, arXiv:2202.08678.  https://arxiv.org/abs/2202.08678
  * Shi et al. (2021) Shi, Y., Zhang, Z.-Y., Wang, J., et al. 2021, ApJ, 909, 20, doi: 10.3847/1538-4357/abd777
  * Tremmel et al. (2020) Tremmel, M., Wright, A. C., Brooks, A. M., et al. 2020, MNRAS, 497, 2786, doi: 10.1093/mnras/staa2015
  * Tulin & Yu (2018) Tulin, S., & Yu, H.-B. 2018, Phys. Rep., 730, 1, doi: 10.1016/j.physrep.2017.11.004
  * van der Burg et al. (2017) van der Burg, R. F. J., Hoekstra, H., Muzzin, A., et al. 2017, A&A, 607, A79, doi: 10.1051/0004-6361/201731335
  * van Dokkum et al. (2015) van Dokkum, P. G., Abraham, R., Merritt, A., et al. 2015, ApJ, 798, L45, doi: 10.1088/2041-8205/798/2/L45
  * Wittmann et al. (2017) Wittmann, C., Lisker, T., Ambachew Tilahun, L., et al. 2017, MNRAS, 470, 1512, doi: 10.1093/mnras/stx1229
  * Wright et al. (2021) Wright, A. C., Tremmel, M., Brooks, A. M., et al. 2021, MNRAS, 502, 5370, doi: 10.1093/mnras/stab081
  * Yagi et al. (2016) Yagi, M., Koda, J., Komiyama, Y., & Yamanoi, H. 2016, ApJS, 225, 11, doi: 10.3847/0067-0049/225/1/11
  * Yang et al. (2020) Yang, D., Yu, H.-B., & An, H. 2020, Phys. Rev. Lett., 125, 111105, doi: 10.1103/PhysRevLett.125.111105
  * Yozin & Bekki (2015) Yozin, C., & Bekki, K. 2015, MNRAS, 452, 937, doi: 10.1093/mnras/stv1073
  * Zaritsky et al. (2019) Zaritsky, D., Donnerstein, R., Dey, A., et al. 2019, ApJS, 240, 1, doi: 10.3847/1538-4365/aaefe9
  * Zavala et al. (2019) Zavala, J., Lovell, M. R., Vogelsberger, M., & Burger, J. D. 2019, Phys. Rev. D, 100, 063007, doi: 10.1103/PhysRevD.100.063007