The Dark Matter Distributions in Low-mass Disk Galaxies. II. The Inner Density Profiles

The rotation curves derived in Paper I are a result of the gravitational pull of both the baryonic matter and dark matter in each galaxy. To isolate the effect of the baryons from that of the dark matter, we need to determine how much the stars and gas contribute to the rotation curve. We will discuss the effects due to the stars only in this section, while the discussion of gas follows in Section 3.

To estimate the potential generated by the stellar distribution we use the package presented by Cappellari (2002), which gives an accurate and robust algorithm for determining multi-Gaussian expansion (MGE) fits to galaxy images. We use images taken in the r-band for all galaxies in our sample, in addition to infrared (4.5 μm) images for the 18 galaxies that have archival data from Channel 2 of Spitzer's IRAC (details about the imaging data can be found in Paper I). When both are available, we preferentially use the infrared images, as they more closely trace the stellar mass distribution. The center, PA, and ellipticity of each galaxy are kept fixed to the values determined in Section 3.2 of Paper I (see below for exceptions). Note that our sample was restricted to galaxies with inclinations between 30° and 70°, as galaxies close to edge-on would not allow us to see the internal motions, and those close to face-on would prevent us from measuring the rotation. The MGE package fits a series of Gaussians to the image in order to produce a surface brightness profile and a luminosity. The stellar contributions to the rotation curves are then found using Jeans Anisotropic Modeling, as described in Cappellari (2008), to calculate the circular velocity in the potential generated by the stars. We model the stellar disk as an oblate ellipsoid with thickness c/a = 0.14 (Kregel et al. 2002).

As mentioned, during the MGE fitting process, the galaxy geometry is kept fixed to values determined in Paper I. For NGC 959 and NGC 7320, we observed a monotonic decrease in ellipticity toward the center, potentially indicating the presence of a bulge. Since this could modestly influence the rotation curve, in these two cases we allowed the ellipticity to vary with radius for the MGE fits.

The above procedure traces the luminous stellar distributions and therefore the shape of the stellar contribution to the rotation curve. However, stellar population synthesis (SPS) models are needed to estimate the stellar mass, and therefore the amplitude of the stellar contribution to the rotation curve. In this section, we use SPS models to give us a plausible estimate of the stellar mass-to-light ratio ϒ_* for each galaxy. Due to the uncertainty in these models and the possible dependence of the outcome on model choice, we will explore the effect of different stellar mass-to-light ratios on the inferred dark matter density profile in Section 5.1.

The NASA-Sloan Atlas (NSA; Blanton et al. 2011) contains 17 of the galaxies in our sample. The information in this catalog allows us to investigate the mass-to-light ratio with homogeneous, matched-aperture photometry in the ugriz bands.

Fitting and Assessment of Synthetic Templates (FAST), introduced by Kriek et al. (2009), is a program that can fit SPS templates to given photometric information. We use FAST to fit models to the 17 galaxies with data in the NSA catalog using the Bruzual & Charlot (2003) stellar population synthesis model, the Calzetti et al. (2000) model for dust attenuation, a Chabrier (2003) initial mass function (IMF), and an exponential star formation history. We then integrate the best-fit model spectra to derive mass-to-light ratios in the r-band and in the IRAC 4.5 μm band, which we denote ϒ_*,r and ${{\rm{\Upsilon }}}_{* ,4.5}$, respectively.

We compare the results from FAST against a second method. The program kcorrect, as presented in Blanton & Roweis (2007), is also based on Bruzual & Charlot (2003) models and a Chabrier (2003) IMF. However, instead it fits the photometry as a linear combination of spectral energy distribution model templates. We ran this procedure on the same 17 galaxies.

The 4.5 μm data were not used to constrain the fits and so they provide a test of the models' ability to predict the near-infrared luminosities. The mean difference between the IRAC 4.5 μm flux and the FAST-predicted flux at the same wavelength is 0.19 dex, compared to a mean difference of −0.02 dex for kcorrect. The models from kcorrect are both an accurate fit to the input optical data as well as a better predictor of the 4.5 μm luminosities, and thus we chose to proceed using the values from this method. Comparison plots of the FAST and kcorrect models alongside the data can be found in Figure 1.

Zoom In Zoom Out Reset image size

Histograms of our results from kcorrect for each wavelength band can be found in Figure 2. The mass-to-light values for 4.5 μm luminosities were closely grouped across our sample (scatter is 20% of mean), leading us to adopt the mean value of ${{\rm{\Upsilon }}}_{* ,4.5}=0.21{M}_{\odot }/{L}_{\odot }$ for all galaxies with data in this band. The values of ϒ_*,r, however, show more variation within the sample (scatter is 37% of mean). This is consistent with our expectation that the mass-to-light ratio is less sensitive to age, dust, and metallicity in the near-infrared than at optical wavelengths (McGaugh & Schombert 2014).

Zoom In Zoom Out Reset image size

There is a wide range in SPS model predictions in the near-infrared, as reviewed by McGaugh & Schombert (2014). The most important model dependences are the treatment of TP-AGB stars and the assumed IMF. Dynamical constraints are needed to calibrate the absolute mass scale. The DiskMass Survey (Martinsson et al. 2013) consider a sample of 30 spiral galaxies with masses M_* ∼ 10^10.3 $\,{\text{}}{M}_{\odot }$. They use measurements of the vertical stellar velocity dispersion to find an estimate of ϒ_*,K that is independent of stellar population synthesis models, finding ϒ_*,K = 0.31 ± 0.07M_☉/L_☉. Converting with ϒ_*,4.5 = 0.91 × ϒ_*,K − 0.08 (Oh et al. 2008) gives a value of ${{\rm{\Upsilon }}}_{* ,4.5}^{\mathrm{DiskMass}}\,=0.20{M}_{\odot }/{L}_{\odot }$, in excellent agreement with our SPS-based estimate of ϒ_*,4.5 = 0.21M_☉/L_☉.

Although this agreement is reassuring, we note that other methods have delivered different results. McGaugh & Schombert (2015) calibrated the tight relation between rotation speed and baryonic mass (the baryonic Tully–Fisher relation) with gas-rich galaxies; they then applied this relation to normal galaxies to estimate ϒ_* and inferred a typical value about twice that of the DiskMass Survey. They conclude that, at present, there is a factor of two systematic uncertainty in the mass-to-light ratio. We will explore a wide range of ϒ_* in Section 5.1 to evaluate the effect of this uncertainty on our dark matter profile measurements.

While a single value of ϒ_*,4.5 suffices for the infrared data, we find a much larger scatter in the r-band (see Figure 2). This presents a problem for estimating ϒ_*,r for the nine galaxies in our sample that are not present in the NSA catalog. The g−r color has been shown to correlate well with the optical mass-to-light ratio (e.g., Bell et al. 2003), which motivates us to explore this correlation using the 17 NSA galaxies. We plotted their g−r magnitudes against the r-band mass-to-light ratios found with kcorrect, as shown in Figure 3, and we fit a line to the data, enabling us to predict the value of ϒ_*,r using only the g−r color of a galaxy.

Zoom In Zoom Out Reset image size

The fit has a linear scatter of 0.17 $\,{\text{}}{M}_{\odot }$/L_☉. Bell et al. (2003) provides a relation between mass-to-light ratio and color, finding ${\mathrm{log}}_{10}{{\rm{\Upsilon }}}_{* ,r}=1.14(g-r)-0.54$. As a comparison, we plot this line alongside our linear fit in Figure 3 and find excellent agreement between the two. For consistency, we will use the mass-to-light ratio predicted by our linear fit for all galaxies.

In order to predict ϒ_*,r for the galaxies not in the NSA, we require r-band magnitudes, which we already have for all galaxies in our sample (see Paper I), and g-band magnitudes. Of the galaxies not in the NSA, five have g-band imaging taken with SPICAM on the 3.5 m telescope at Apache Point Observatory on 2013 October 6 and 2014 November 11. We use the MGE process described in Section 2.1 to determine the g-band magnitudes. The g-band magnitude for NGC 2976 is provided by the Spitzer Local Volume Legacy Survey (Cook et al. 2014). Three galaxies did not have g-band data available, so we instead calculate the B−V color and convert it to g−r color following Jester et al. (2005). NGC 6503 has B- and V-band data from Cook et al. (2014), while NGC 7320 and UGC 3371 have the same from de Vaucouleurs et al. (1991).

Table 1 gives the values of ϒ_* as well as the wavelength used for mass models (r-band or 4.5 μm).

Table 1.  Stellar Mass-to-light Ratio Estimates

Galaxy	${{\rm{\Upsilon }}}_{* }\pm \sigma$	$\mathrm{log}{M}_{* }$	ϒ*	λ
	(SPS est.)	($\,{\text{}}{M}_{\odot }$)	(Max. Disk.)
NGC 746	0.87 ± 0.15	8.9	0.67	r-band
NGC 853	0.67 ± 0.15	9.3	1.59	r-band
NGC 949	1.09 ± 0.15	9.3	1.35	r-band
NGC 959	0.21 ± 0.04	8.8	0.53	4.5 μm
NGC 1012	0.99 ± 0.15	9.4	0.18	r-band
NGC 1035	0.21 ± 0.04	9.6	0.22	4.5 μm
NGC 2644	0.21 ± 0.04	9.7	0.47	4.5 μm
NGC 2976	0.21 ± 0.04	8.9	0.28	4.5 μm
NGC 3622	0.21 ± 0.04	9.2	0.06	4.5 μm
NGC 4376	0.21 ± 0.04	9.1	0.53	4.5 μm
NGC 4396	0.21 ± 0.04	9.2	0.42	4.5 μm
NGC 4451	0.21 ± 0.04	9.6	0.27	4.5 μm
NGC 4632	0.21 ± 0.04	9.5	0.41	4.5 μm
NGC 5303	0.21 ± 0.04	9.7	0.24	4.5 μm
NGC 5692	0.83 ± 0.15	9.5	1.06	r-band
NGC 5949	0.21 ± 0.04	9.3	0.50	4.5 μm
NGC 6106	0.21 ± 0.04	9.7	0.22	4.5 μm
NGC 6207	0.21 ± 0.04	9.6	0.31	4.5 μm
NGC 6503	0.21 ± 0.04	9.4	0.22	4.5 μm
NGC 7320	0.21 ± 0.04	9.0	0.53	4.5 μm
UGC 1104	0.21 ± 0.04	8.1	0.35	4.5 μm
UGC 3371	0.41 ± 0.15	8.5	0.89	r-band
UGC 4169	0.21 ± 0.04	9.4	0.26	4.5 μm
UGC 8516	0.21 ± 0.04	9.1	0.14	4.5 μm
UGC 11891	0.62 ± 0.15	8.9	1.55	r-band
UGC 12009	0.67 ± 0.15	9.1	0.45	r-band

Note. Stellar masses are calculated using the SPS-predicted value of ϒ_*. The maximum disk ϒ_* is found through scaling the stellar contribution to the rotation curve; see Section 2.3 for more detail.

Download table as:  ASCIITypeset image

In order to test the sensitivity of our derived dark matter profiles to our estimates of the stellar distribution, we wish to explore the full range of plausible ϒ_*. In addition to the estimates based on SPS models described in the previous section, we will also consider limits on ϒ_* that are based only on kinematic constraints. The "minimum disk" and "maximum disk" hypotheses bracket the kinematically allowed values of ϒ_*.

The minimum value the mass-to-light ratio can take is when there is hypothetically no contribution to the gravitational potential from the stars. In this case, ϒ_* = 0 and the rotation curve is a dark-matter-only model (we will also assume gas is negligible, see Section 3 for further discussion). A common procedure for estimating the maximum ϒ_* value is to increase the amplitude of the model stellar rotation curve until it exceeds the data at some radius. This simple approach has some shortcomings, however, which necessitates minor modifications.

First, the radius at which the stellar rotation curve exceeds the data is generally found to be the innermost observed radius, which often leads to an unphysically small estimate of the maximum ϒ_*. This is likely due to random uncertainties in the rotation curve that this procedure does not account for, as well as any non-circular motions that might affect the innermost point. As a more conservative estimate of the maximum ϒ_*, we instead consider the value that produces a stellar rotation curve that exceeds the total rotation curve at two radial bins.

While an improvement over surpassing one radial bin, there remain several special cases we must address. For four galaxies (NGC 959, NGC 4376, NGC 7320, UGC 11891), surpassing two radial bins leads instead to an unphysically large value of ϒ_*, motivating us to set an upper limit of 2.5 times the value of ϒ_* predicted from SPS models. There were also seven cases where the maximum disk estimate of ϒ_* was smaller than the SPS-predicted value. For four of these cases, studying the maximum disk case is inconsequential because we ultimately find that these galaxies are not suitable for measuring the dark matter profile (see Section 4.2); however, three require addressing. The SPS-predicted ϒ_* value for UGC 12009 already surpasses two (nearly three) radial bins out of a total of five (as seen in Figure 4). As the maximum disk value could not reasonably be larger, we keep the original estimate, despite being lower than the SPS-predicted value. For NGC 2976 and NGC 6503, there is a small peak in the inner region of the stellar model rotation curve that causes the low estimate of ϒ_* when surpassing two radial bins. We therefore define the maximum disk ϒ_* as the value that exceeds the first radial bin that falls outside the range affected by the peak. Such peaks arise from a bright nucleus, which is likely to have a different mass-to-light ratio than the rest of the galaxy, and so their omission from our maximum disk estimate is justified.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Although these criteria are ad hoc, they are tailored to our goal of studying the inner mass distribution. We found it important to have a maximum disk definition that ensured physical results within ∼1 kpc. For comparison, we also considered the definition of the maximum disk presented by Sackett (1997), in which the disk produces 85% of the velocity at 2.2 times the disk scale radius. This produced a higher maximum disk ϒ_* than our previous method in 12 of 18 galaxies, and in almost all of those cases, the stellar contribution would significantly exceed the total rotation throughout all or most of the inner kiloparsec. This may be related to the fact that the Sackett (1997) maximum disk estimate is based on observations of much more massive galaxies, where the variation of the stellar mass fraction with radius is likely different.

The final values of the maximum disk ϒ_* can be found in Table 1. We explore the effects of these values on our mass models in Section 5.1.

While most of our galaxies do not have resolved H i data, total H i fluxes exist in the literature for all of the galaxies in our sample. We are able to use these total flux values to estimate the H i mass surface density, which can then be used to estimate the contribution to the rotation curve.

Some galaxies have only one source giving a value for the H i flux, but for the cases with multiple sources, we preferentially select those that had been corrected for beam attenuation, pointing offsets, and H i self-absorption. Data for 13 of our galaxies are found in Springob et al. (2005), which includes all of the aforementioned corrections, 6 are found in Huchtmeier & Richter (1989), Doyle et al. (2005) and Popping & Braun (2011) each contain 2, while Schneider et al. (1992), Greisen et al. (2009), and Walter et al. (2008) contain 1 each. The H i fluxes were converted to a total mass using ${M}_{\mathrm{HI}}\,=(2.36\times {10}^{5}){D}^{2}f$, where D is the distance in Mpc, f is the flux in Jy km s⁻¹, and M_HI is in solar masses (Roberts 1962).

We follow the method in Martinsson et al. (2013) to use these values to estimate the H i mass distribution. They find that the radial H i mass surface density profile is well fit with a Gaussian as follows:

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{HI}}(R)={{\rm{\Sigma }}}_{\mathrm{HI}}^{\max }\exp \left[-\displaystyle \frac{{\left(R-{R}_{{\rm{\Sigma }},\max }\right)}^{2}}{2{\sigma }_{{\rm{\Sigma }}}^{2}}\right]\end{eqnarray} \tag{ 1 }$$

where ${{\rm{\Sigma }}}_{\mathrm{HI}}^{\max }$ is the peak density, ${R}_{{\rm{\Sigma }},\max }$ is the radius at which the peak occurs, and σ_Σ is the width of the profile. Since ${{\rm{\Sigma }}}_{\mathrm{HI}}^{\max }$ is a normalization constant, it is determined explicitly by integrating over the disk. The remaining values are found using R_HI, defined as the radius for which ${{\rm{\Sigma }}}_{\mathrm{HI}}=1{M}_{\odot }\,{\mathrm{pc}}^{-2}$. There is a relationship between R_HI and the total H i mass, though it varies somewhat depending on the mass range of the sample. The galaxies used in Martinsson et al. (2013) are more massive than ours, so we use the relation from Swaters et al. (2002) as their mass range is closer to ours. They find:

$$\begin{eqnarray}&&\mathrm{log}{M}_{\mathrm{HI}}=1.86\mathrm{log}(2{R}_{\mathrm{HI}})+6.6,\end{eqnarray} \tag{ 2 }$$

where M_HI is measured in units of $\,{\text{}}{M}_{\odot }$ and R_HI is in kiloparscecs. Taking the radial H i surface density profile (averaged over all of the dwarf galaxies in the Swaters et al. 2002 sample), we fit a Gaussian to determine ${R}_{{\rm{\Sigma }},\max }=0.20{R}_{\mathrm{HI}}$ and σ_Σ = 0.44R_HI. Finally, we multiply Σ_HI by 1.4 to account for helium and metals, Σ_atomic = 1.4Σ_HI (Martinsson et al. 2013). The values used in this calculation for each galaxy can be found in Table 2.

Table 2.  H i Parameters

Galaxy	H i Flux	MHI	RHI
	(Jy km s−1)	(108 $\,{\text{}}{M}_{\odot }$)	(kpc)
NGC 746	18.6	6.3	7.6
NGC 853	6.0	5.7	7.2
NGC 949	16.4	3.9	5.9
NGC 959	16.7	4.0	5.9
NGC 1012	55.1	25.5	16.1
NGC 1035	14.2	8.6	9.0
NGC 2644	3.7	7.9	8.6
NGC 2976	45.4	1.4	3.4
NGC 3622	13.0	16.2	12.7
NGC 4376	3.7	5.0	6.7
NGC 4396	21.4	12.9	11.2
NGC 4451	3.2	5.0	6.8
NGC 4632	55.9	25.8	16.2
NGC 5303	11.3	16.7	12.8
NGC 5692	3.8	6.5	7.7
NGC 5949	5.8	2.3	4.4
NGC 6106	33.0	44.9	21.9
NGC 6207	34.4	20.8	14.5
NGC 6503	205.0	19.2	13.9
NGC 7320	8.3	3.8	5.8
UGC 1104	9.8	2.3	4.4
UGC 3371	31.6	16.8	12.9
UGC 4169	29.0	61.7	25.9
UGC 8516	3.9	4.4	6.3
UGC 11891	88.5	12.4	10.9
UGC 12009	9.8	9.2	9.3

Note. See Section 3.1 for derivation details.

Download table as:  ASCIITypeset image

The H i surface densities are then used to estimate the atomic gas contribution to the rotation curve. This is done analytically, following the method in Section 2.6.3 of Binney & Tremaine (1987), which derives an expression for v_circ for a thin disk given the surface density Σ(R).

We have resolved CO data for the 11 galaxies in our sample that were studied in Truong et al. (2017), which allows us to determine the distribution of molecular gas in these galaxies. Using the CO moment-0 maps provided by Truong et al. (2017), we convert to H₂ mass surface density maps using:

$$\begin{eqnarray}\begin{array}{rcl}\left[\displaystyle \frac{{{\rm{\Sigma }}}_{{{\rm{H}}}_{2}}}{\,{\text{}}{M}_{\odot }\,{\mathrm{pc}}^{-2}}\right] & = & 1.6\left[\displaystyle \frac{{I}_{\mathrm{CO}}{\rm{\Delta }}V}{{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right]\\ & & \times \,\left[\displaystyle \frac{{X}_{\mathrm{CO}}}{{10}^{20}{\mathrm{cm}}^{-2}{\left({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1}\right)}^{-1}}\right]\cos i\end{array}\end{eqnarray} \tag{ 3 }$$

where i is the inclination of the galaxy and X_CO is the CO to H₂ conversion factor. We use a value of ${X}_{\mathrm{CO}}=2.0\,\times {10}^{20}\,{\mathrm{cm}}^{-2}{({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1})}^{-1}$, following the recommendation in Bolatto et al. (2013). The resulting H₂ mass surface density maps are then multiplied by 1.4 to account for helium and metals. We average the maps over a series of elliptical bins to produce a radial molecular mass surface density profile for each galaxy.

We again analytically derive the rotation curves by assuming a thin disk and following Section 2.6.3 of Binney & Tremaine (1987). Before we can include them in our mass model, we first need to extend the profiles, as several of the CO maps do not extend as far out as our Hα data. We do this by fitting an exponential model, which is then scaled to agree with the outermost point of the CO data to ensure continuity. When necessary, we use this function to extrapolate the curve to match the extent of the Hα data. This did not produce changes in the rotation curves within the inner 1 kpc for any of our galaxies, so this extrapolation does not have a significant effect on our results.

The rotation curves derived in Paper I represent the combined gravitational force from all of the mass components of the galaxy (dark matter, stars, gas). In order to determine the distribution of dark matter, the various contributions to the rotation curve need to be separated. To do this we use a Markov Chain Monte Carlo (MCMC) fitting procedure to model the rotation curves as a sum of the stellar, gaseous, and dark matter components: ${V}_{\mathrm{tot}}^{2}={V}_{\mathrm{stars}}^{2}+{V}_{\mathrm{gas}}^{2}+{V}_{\mathrm{DM}}^{2}$. We use the estimates of the stellar and gaseous rotation curves derived in Sections 2 and 3 respectively, and model the dark matter as a spherical halo following a generalized NFW (gNFW) profile:

$$\begin{eqnarray}&&\rho (r)=\displaystyle \frac{{\rho }_{0}}{{\left(\tfrac{r}{{r}_{s}}\right)}^{\beta }{\left(1+\tfrac{r}{{r}_{s}}\right)}^{3-\beta }}.\end{eqnarray} \tag{ 4 }$$

The gNFW profile includes the NFW scale radius r_s and characteristic density ρ₀, as well as an additional parameter β, that dictates the slope as $r\to 0$ (note that setting β = 1 returns the standard NFW profile). Typically, galaxies are described using a parameterization defined by the virial mass M₂₀₀ (given by the enclosed mass at r₂₀₀, the radius at which the mean enclosed density is equal to 200 times the critical density) and the concentration c₂₀₀ = r₂₀₀/r_s. In the case of the generalized NFW profile, we found that parameter covariance was reduced by recasting Equation (4) in terms of the parameter r₋₂ = (2 − β) r_s, which is the radius at which ρ ∝ r⁻² locally (see e.g., Newman et al. 2015; note that r₋₂ reduces to r_s for the NFW case where β = 1), and the concentration c₋₂ = r₂₀₀/r₋₂.

In addition to the three parameters describing the gNFW model, we also include the stellar mass-to-light ratio ϒ_* as a parameter with a prior based on the results of our population synthesis analysis from Section 2.2. Assuming Gaussian errors, the likelihood function L, including the full covariance matrix, is given by

$$\begin{eqnarray}&&\mathrm{log}L=-N\mathrm{log}j-\displaystyle \frac{1}{2}{\left({V}_{t}-m\right)}^{T}\bullet \ \left(\displaystyle \frac{{C}^{-1}}{{j}^{2}}\right)\bullet ({V}_{t}-m),\end{eqnarray} \tag{ 5 }$$

where • represents matrix multiplication, N is the number of points in the rotation curve data given by V_t, m represents the model of the total rotation curve (dark matter, stars, and gas), C is the N × N covariance matrix (obtained from the DiskFit output as described in Paper I), and j represents a jitter term. In order to control the errors, we included this "jitter" term, which is meant to aid in the cases where the model does not describe the data well by increasing the errors until a reasonable fit is obtained (see e.g., Price-Whelan et al. 2017). This means j ∼ 1 when the model adequately describes the data within the uncertainties, while j > 1 indicates that the model is inadequate or the uncertainties are underestimated. We also use the full covariance matrix from the bootstrapping analysis in order to take into account the correlated uncertainties among radial bins (see Paper I for details).

We use the MCMC ensemble sampler emcee (Foreman-Mackey et al. 2013), setting uniform priors of $0.0\lt \beta \lt 2.0,8\lt \mathrm{log}{M}_{200}\lt 12.0,1.0\lt {c}_{-2}\lt 40.0$ and Gaussian priors on $\mathrm{log}{{\rm{\Upsilon }}}_{* }$ and $\mathrm{log}j$. The prior for ϒ_* is centered on the SPS value with a width dictated by the uncertainty (see Table 1), while $\mathrm{log}j$ is centered on 0 with a dispersion of 0.7, corresponding to a factor of two in j. We first initialize the emcee walkers at positions drawn from these priors, then use a short initial run (500 steps) to estimate rough parameter values. These estimates are then used to initialize the final MCMC run of 50 walkers and 10,000 steps. We use only the second half of the chains in our analysis.

The MCMC procedure gives us a fiducial mass model for each galaxy that breaks down the contributions from the various mass components to the total rotation curve. Figure 4 shows the model of the PCWI rotation curves (comprising all components) as well as the decomposition into mass component contributions (dark matter, stars, atomic gas, molecular gas when available). While we generally find that our model of the total rotation curve is able to match the PCWI data, not all of our galaxies were fit well, as is to be expected.

In order to address the potential inadequacies of the model, we assign a grade from 1 to 3 to each galaxy. We take into account how well the total mass model fit the rotation curve, particularly the slope in the inner parts. We also note if there are strong radial motions, which are plotted in Figure 4 as well. If the radial velocities are a significant fraction of the tangential velocities, it could indicate the presence of irregular or non-circular motions, which may affect our ability to model such galaxies accurately. Finally, we consider the two-dimensional (2D) residuals from the DiskFit model.

As expected, there is a broad correspondence between our subjective grades and the quantitative value of the jitter term j. The galaxies with the largest jitter values received poor grades, while the grade 1 galaxies all have $\mathrm{log}j\leqslant 0.5$. This correspondence reinforces the validity of our subjective grades, which we prefer to use as the deciding criteria as the jitter term only quantifies the goodness of fit of the rotation curve and does not account for other factors, such as the fit quality of the velocity field or evidence of non-circular motions.

Galaxies assigned grades of 1 or 2 are fit sufficiently well to proceed with further analysis. The eight galaxies assigned a grade of 3 had mass models that we judge did not adequately describe the PCWI data and therefore cannot constrain the value of the inner dark matter density profile. We therefore exclude the grade 3 galaxies for the remainder of the paper. Moving forward only with grades 1 and 2 leaves us a final sample size of 18 galaxies. The grades of individual galaxies can be found in Table 3.

Table 3.  Inner Dark Matter Density Profile Slope

Galaxy	${\beta }^{* }\pm \sigma$	${\beta }^{* }\pm \sigma$	${\beta }^{* }\pm \sigma$	$\mathrm{log}j\pm \sigma$	Grade
		(Min. Disk)	(Max. Disk)
NGC 959	1.25 ± 0.12	1.21 ± 0.09	1.34 ± 0.24	0.28 ± 0.19	1
NGC 1035	0.58 ± 0.11	0.95 ± 0.21	0.51 ± 0.10	0.50 ± 0.16	1
NGC 2976	0.33 ± 0.05	0.42 ± 0.03	0.31 ± 0.08	−0.07 ± 0.14	1
NGC 4376	0.70 ± 0.10	0.70 ± 0.08	0.68 ± 0.15	0.37 ± 0.22	1
NGC 4396	0.84 ± 0.12	0.86 ± 0.10	0.81 ± 0.14	−0.15 ± 0.16	1
NGC 5303	0.54 ± 0.19	0.91 ± 0.11	0.50 ± 0.25	0.21 ± 0.26	1
NGC 5692	0.76 ± 0.24	1.09 ± 0.10	0.74 ± 0.46	0.35 ± 0.23	1
NGC 5949	0.84 ± 0.08	0.83 ± 0.04	1.11 ± 0.32	0.47 ± 0.21	1
NGC 6207	0.73 ± 0.13	0.79 ± 0.07	0.60 ± 0.14	0.35 ± 0.14	1
NGC 7320	1.03 ± 0.07	1.04 ± 0.16	1.05 ± 0.11	0.18 ± 0.17	1
UGC 4169	0.59 ± 0.16	0.58 ± 0.13	0.59 ± 0.17	0.36 ± 0.17	1
NGC 853	1.10 ± 0.25	0.95 ± 0.16	1.58 ± 0.45	0.94 ± 0.26	2
NGC 2644	0.70 ± 0.23	0.69 ± 0.16	1.14 ± 0.54	0.02 ± 0.22	2
NGC 6106	0.59 ± 0.14	0.68 ± 0.06	0.58 ± 0.15	0.57 ± 0.15	2
NGC 6503	0.72 ± 0.04	0.81 ± 0.03	0.70 ± 0.04	0.70 ± 0.12	2
UGC 1104	0.67 ± 0.35	1.21 ± 0.38	0.84 ± 0.46	0.41 ± 0.33	2
UGC 11891	1.00 ± 0.31	0.95 ± 0.28	1.01 ± 0.34	0.07 ± 0.20	2
UGC 12009	1.20 ± 0.54	1.26 ± 0.19	1.09 ± 0.50	0.19 ± 0.34	2
NGC 746	0.64 ± 0.32	0.67 ± 0.15	0.58 ± 0.23	0.24 ± 0.18	3
NGC 949	1.60 ± 0.42	1.41 ± 0.11	1.59 ± 0.46	0.73 ± 0.19	3
NGC 1012	1.09 ± 0.62	0.40 ± 0.09	0.30 ± 0.09	1.26 ± 0.18	3
NGC 3622	0.76 ± 0.38	0.70 ± 0.18	0.67 ± 0.21	0.44 ± 0.26	3
NGC 4451	1.48 ± 0.48	0.93 ± 0.22	1.52 ± 0.49	0.66 ± 0.25	3
NGC 4632	1.31 ± 0.15	1.33 ± 0.08	1.21 ± 0.47	0.92 ± 0.14	3
UGC 3371	0.47 ± 0.14	0.42 ± 0.13	0.50 ± 0.16	0.69 ± 0.22	3
UGC 8516	0.84 ± 0.56	0.38 ± 0.12	0.47 ± 0.28	1.10 ± 0.24	3

Note. See Section 4.3 for the definition of β* and Section 4.1 for details about $\mathrm{log}j$.

Download table as:  ASCIITypeset image

The mass models of grade 1 and 2 galaxies are presented in Figure 4, while the grade 3 galaxies are shown in Figure 12.

Although the value of β in the gNFW profile gives the slope of the dark matter density profile as $r\to 0$, this is not the best measured quantity as it is naturally asymptotic. Additionally, β is often covariant with one or more of the other parameters. Since many of our rotation curves do not fully reach the "flattened" part, we often cannot constrain all of the gNFW parameters well in our MCMC analysis. In particular, covariance with the scale radius r_s means the data can often be fit with an arbitrarily small β and small r_s, leading us only to an upper bound on β.

To remedy this, we decided to define a mean value of the inner slope over a fixed radial range, β*. Motivated by the resolution of our PCWI data and the CO data from Truong et al. (2017), as well as by the resolution achieved by simulations examining similar mass galaxies, we define β* over the range 300–800 pc. Taking the posteriors from the analysis described in Section 4.1, we calculate the average slope over this range as

$$\begin{eqnarray}&&{\beta }^{* }=-\displaystyle \frac{\mathrm{log}\tfrac{\rho (0.8\ \mathrm{kpc})}{\rho (0.3\ \mathrm{kpc})}}{\mathrm{log}\tfrac{0.8}{0.3}}.\end{eqnarray} \tag{ 6 }$$

We find that β* is much better constrained by the data than the asymptotic β from the gNFW profile. As seen in Figure 5, there is much less covariance between concentration c₋₂ and β* than with β, and we are not restricted to citing only an upper bound, instead deriving a tight constraint on β*.

Zoom In Zoom Out Reset image size

Table 3 gives the values of β* determined for each of our galaxies. We find a mean of β* = 0.74 ± 0.07 and an intrinsic scatter of ${0.22}_{-0.05}^{+0.06}$. We will discuss the interpretation of the β* distribution and compare it to other observations and theoretical expectations in Sections 6 and 7.

In order to compare this value to that of the NFW profile, we need to calculate the value of β* for a typical galaxy in our mass range if it were to follow NFW exactly. We choose a random sample of 10,000 masses drawn from our galaxies' stellar mass distribution and convert them to halo masses following Leauthaud et al. (2012). We then determine the concentrations following the relation and scatter given by Equation (8) of Dutton & Macciò (2014). These are used to construct NFW profiles for each of the 10,000 mock galaxies, which we then use to calculate β*. We find that a typical galaxy in our mass range that followed the NFW profile would have a value of ${\beta }_{\mathrm{NFW}}^{* }=1.05\pm 0.02$.

Our estimates of ϒ_* in Section 2.2 rely on the validity of the underlying stellar population synthesis models. In order to test the sensitivity of our dark matter profiles to the estimates of the stellar distribution, we want to explore the full range of kinematically plausible ϒ_*. As discussed in Section 2.3, we will consider a minimum and a maximum disk case, which bracket the range of kinematically permitted values of ϒ_*.

Figure 6 compares the value of β* from these two additional cases to our fiducial model for the 18 galaxies with grades 1 or 2. The values of β* for all cases can be found in Table 3. Typically the three values are close together and fall within one standard deviation of the fiducial value for each galaxy; however, there are a few cases where the minimum disk case is larger than the fiducial β* by more than one standard deviation. This is not surprising since this model represents a galaxy with no stellar or gaseous matter whatsoever (dark-matter-only) and is not expected to be a physically accurate description.

Zoom In Zoom Out Reset image size

While in every case, the maximum disk estimate of β* falls within one standard deviation of the fiducial case, it is somewhat surprising to note that the cases that disagree by the largest margin imply steeper profiles than the fiducial case. This is contrary to the usual expectations that the minimum and maximum assumptions bracket the kinematically permitted values of the inner slope.

We conclude that β* is indeed robust and not strongly influenced by our choice of ϒ_*, and thus any uncertainties in our estimate of ϒ_* do not have a significant effect on our final results. We will discuss this insensitivity of β* to ϒ_* further in Section 6.2.

Since the gas mass is sometimes neglected in rotation curve analysis and the relevant data are not always available, we chose to examine the effects of excluding gas from our mass models entirely. We proceed with the same MCMC analysis as before, but we instead model the rotation curves only as ${V}_{\mathrm{tot}}^{2}={V}_{\mathrm{stars}}^{2}+{V}_{\mathrm{DM}}^{2}$. We again keep ϒ_* as a free parameter in order to directly compare to our fiducial case.

We find that the value of β* effectively did not change when we removed the gas from the mass model, the largest difference being 0.09, with 14 of the 18 galaxies having differences of 0.02 or less.

This is not surprising, as we found the contribution to the rotation curves due to gas to be particularly insignificant in the inner portions, with atomic gas being essentially negligible in the 300–800 pc range in most cases. This leads us to conclude that excluding gas from mass models will not have a tangible effect on the shape of the estimated dark matter distribution, and that galaxies in our sample without molecular gas masses can be considered alongside the others.

In Section 4.5 of Paper I, we attempted to fit all of our galaxies with a bisymmetric kinematic model in order to determine if any of our sample was likely to contain a central bar. Only for four of the galaxies (NGC 949, NGC 2976, NGC 3622, NGC 4376) did the bisymmetric fit produce a rotation curve that was significantly different than the axisymmetric fit. For the remaining galaxies there was either no effective difference between the rotation curves produced by the two models, or the bisymmetric model was unable to converge. The median reduced χ² value was also lower for the axisymmetric fits, though only somewhat (${\chi }_{\mathrm{axisym}.}^{2}$/dof = 0.82, ${\chi }_{\mathrm{bisym}.}^{2}$/dof = 0.97).

Closer inspection of the residuals for NGC 3622 showed an unmistakable bar-like pattern, leading us to conclude the galaxy is barred. We thus used the bisymmetric rotation curve for the entirety of our analysis. We could not draw a similar conclusion for the remaining three, so we proceeded with the axisymmetric models for our analysis up until this point.

To examine what effect the presence of a bar might have, we repeat the mass modeling process using the bisymmetric rotation curves for the three candidate galaxies, in addition to the previously discarded axisymmetric fit for NGC 3622. Unfortunately, this new sample of rotation curves was not modeled well by the MCMC procedure (all received a grade of 3; mass models can be seen in Figure 13). The fiducial rotation curves for NGC 949 and NGC 3622 were also not adequately fit by our mass model. Thus we are not able to comment on the effects of the presence of a central bar, but we note that the change to β* was relatively small in all cases.

The presence of a bar in NGC 2976 was also considered by Spekkens & Sellwood (2007) and Adams et al. (2014) (for a comparison of rotation curves, see Paper I). Spekkens & Sellwood (2007) apply both radial and bisymmetric models to the velocity field of NGC 2976, taken from Simon et al. (2003). They found both fits to be adequate parameterizations of the kinematics, but favor the interpretation that NGC 2976 hosts a bar. Adams et al. (2014) also model their gas kinematics with both radial and bisymmetric fits, in addition to using stellar kinematics. They find the radial model shows a stronger core than the bisymmetric and stellar models, which are both somewhat steeper. This agreement with stellar kinematic data also leads them to favor the bisymmetric model. While the bisymmetric mass model did not result in a good fit to our data for NGC 2976, we do find the inferred dark matter slope to be slightly steeper than the axisymmetric model (${\beta }_{\mathrm{axisym}.}^{* }\,=0.33,{\beta }_{\mathrm{bisym}.}^{* }=0.42$), consistent with Adams et al. (2014).

Due to our inability to model the bisymmetric rotation curves well, we are unfortunately unable to draw any conclusions about the robustness of β* to the choice of kinematic model for the four galaxies where the kinematic models differed.

Figure 7 shows the distribution of β* for our sample, while the values are tabulated in Table 3. Our galaxies cover a range of β*, from as shallow as 0.33 to steeper than NFW at 1.25.

Zoom In Zoom Out Reset image size

The majority of our galaxies are found to have inner dark matter density profile slopes that are smaller than NFW (β* below 1.05), but many are only moderately shallower, deviating from the expected slope but not enough to be a nearly constant-density core. Around a third are consistent with NFW-like or even steeper profiles. The steeper profiles are likely due to baryonic effects such as adiabatic contraction (Blumenthal et al. 1986; Dutton et al. 2016), but it is notable because any theory put forth to explain the presence of shallow or cored dark matter profiles must also be able to account for steeper profiles in some galaxies.

We want to quantify the scatter in our distribution, which requires accounting for the different measurement uncertainties for each galaxy. We do this by hypothesizing that the true distribution is Gaussian, then we use the ensemble of posteriors on β* from each galaxy to derive constraints on the mean and dispersion of the true distribution. With this method, we find a mean $\langle {\beta }^{* }\rangle =0.74\pm 0.07$ and intrinsic scatter of $\sigma \,={0.22}_{-0.05}^{+0.06}$. This σ is clearly inconsistent with 0, demonstrating the significant diversity of dark matter profiles in our sample.

To investigate the stability of this mean, we also ran the above calculation on the minimum and maximum disk posteriors. We find $\langle {\beta }^{* }\rangle =0.85\pm 0.06$ and $\langle {\beta }^{* }\rangle ={0.70}_{-0.07}^{+0.08}$ for the minimum and maximum cases, respectively, both of which agree with the fiducial mean within their errors. We also considered only the galaxies that were assigned a grade of 1, finding $\langle {\beta }^{* }\rangle =0.74\pm 0.09$. The intrinsic scatter shows a similar consistency across the samples.

Overall, we find that the dark matter profiles in our galaxies are only moderately shallower than NFW, but they have a significant range of values of β*, showing that our sample of galaxies contains a diversity of dark matter profiles. This diversity is critical to looking for correlations between the inner slope and other galaxies properties, in addition to determining the cause of these different profiles, as any solution must be able to account for the range we find.

While our analysis so far has focused on the innermost region (300–800 pc), where the effects of baryonic feedback or new dark matter physics are expected to be strongest, it is also interesting to examine the shape of the rotation curves (and hence the dark matter distribution) as a whole. We would like to explore what, if any, large-scale dark matter properties can be derived from our sample.

To do this, we first compare the shapes of our baryon-subtracted rotation curves to NFW rotation curves. To determine the shape our rotation curves would take if our galaxies followed the NFW profile, we need to normalize to the maximum rotation velocity for each galaxy, as this value and the mass–concentration relation will determine the parameters in the NFW profile. We cannot always use the maximum velocity taken from our rotation curves as they are often still rising at the outermost measured radius, so we also consider maximum velocities derived from H i line widths (taken from HyperLeda) and use the larger of the two values. We again use the Dutton & Macciò (2014) mass–concentration relation to estimate the value of M₂₀₀ that reproduces the maximum velocity for each galaxy. We then use these NFW profiles to determine the corresponding dark matter contribution to the rotation curve, allowing us to compare directly to the dark matter components inferred from our mass models.

The left and center panels of Figure 8 show the dark matter contribution to the rotation curves plotted as a fraction of maximum velocity for the fiducial case and the maximum disk assumption, as well as the same for the minimum disk assumption (see Section 2.3). The curves are overlaid with the region spanned by NFW profiles following the Dutton & Macciò (2014) mass–concentration relation. The fiducial case typically has a smaller rotation velocity than NFW in the inner parts (r ≲ 1 kpc), which agrees with our expectation that our galaxies typically have less dark matter in their central regions, then approximately aligns with NFW at large radii, also as expected. The maximum disk case has a similar pattern, and both cases tend to show smaller velocities in the innermost region than NFW, though the fiducial case does so less. In contrast, the minimum disk case aligns well in the inner parts, deviating instead in the outer regions, where some galaxies overestimate the NFW velocities, again in line with expectations as this case assumes all of the galaxies' mass is due to dark matter.

Zoom In Zoom Out Reset image size

The right panels of Figure 8 give the ratio of the fiducial to maximum and minimum disk dark matter rotation curves. The flat lines in these panels indicate that our sample overall shows similar rotation curve shapes in the fiducial, maximum, and minimum disk cases, despite having different amplitudes. We see that while moving to a maximum or minimum disk model does systematically change the amount of dark matter in the inner kiloparsec, the slope in this region does not significantly change in most cases.

This seems to indicate that insensitivity of β* to the mass-to-light ratio is due to a combination of the relatively small differences between the maximum disk and SPS-predicted values of ϒ_* in many cases, and the similarity between the power-law slopes of the stars and dark matter in other cases.

Note that the only galaxy with a maximum disk ratio less than one is UGC 12009, for which the maximum disk estimate of ϒ_* is smaller than the fiducial value. The few galaxies with non-zero slopes in the ratio plot correspond to those with the largest deviations in β* (see Figure 6). In these few cases, the rotation curve permits a much larger ϒ_* than our SPS-based estimate.

This indicates that in our sample the dark matter density slope in the inner ∼1 kpc is more robustly measured than the dark matter fraction or concentration. Both depend on the amplitude of the rotation curves and their shape on larger scales, which are covariant with ϒ_*. Hence in this work, we focus on β* rather than the larger-scale dark matter properties.

The distribution in Figure 7 shows the variety of density profiles found in our sample of galaxies. This range enables us to look for correlations between the inner slope and other galaxy properties.

There are predictions from the Feedback in Realistic Environments (FIRE) simulations (Hopkins et al. 2014) that the inner slope of the dark matter profile should correlate with the stellar distribution and specific star formation rate. Strong gas outflows could reduce the gravitational potential, allowing both the dark matter and the stars to migrate outward. This process would be reflected by an increase in effective radius as well as a decrease in specific star formation rate (thus reddening the color; Oñorbe et al. 2015; El-Badry et al. 2016). Oñorbe et al. (2015) also find that in order to maintain a cored profile, a galaxy must have significant late-time star formation to keep the dark matter from re-accreting and forming a cusp. This indicates that core formation may be a continuous process, and the slope of the inner dark matter profile will evolve with time.

To look for signs of this feedback, we want to examine trends of the dark matter slope with specific star formation rate and effective radius. Figure 9 shows the correlations between β* and r_eff, as well as between β* and g−r color, which we use as a rough proxy for specific star formation rate, as well as the stellar surface density within 1 kpc, Σ_1kpc. In the case of g−r color, we find an intriguing suggestion of a correlation with the expected sign based on these physical considerations; however, we are unable to draw firm conclusions because of the low statistical significance of the correlation (1.6σ). For r_eff and Σ_1kpc, we find no evidence of a correlation. Note that we do not sample low specific star formation rate galaxies, for which tracers other than Hα would be required, so this result does not rule out the existence of such a correlation.

Zoom In Zoom Out Reset image size

If feedback is responsible for driving dark matter away through fluctuations of the gravitational potential, this mechanism would not be as efficient in more massive galaxies, as their deeper potentials weaken gas outflows. Discussion of a correlation with stellar mass in the context of existing literature as well as a more explicit comparison to the FIRE simulations can be found in Sections 7.2 and 7.3.

As discussed previously, Truong et al. (2017) used CO velocity fields to produce rotation curves for a sample of dwarf galaxies that overlaps with this work. The second paper in that series, P. N. Truong et al. (in preparation), uses the same methodology as this work to produce mass models and analyze the dark matter density profiles, allowing us to make a direct comparison of β* for the seven galaxies that were well fit by both samples.

We compare the slopes derived using Hα and CO tracers, as seen in Figure 10. Our results are consistent with those from P. N. Truong et al. (in preparation) and we do not find any systematic difference between the data sets. Note that although the CO value of β* for NGC 5303 appears to be more in agreement with NFW while the Hα shows a core, the CO measurement has a large-enough uncertainty that the results are consistent.

Zoom In Zoom Out Reset image size

Although one might expect CO to be a more faithful tracer, as it is kinematically colder than the ionized gas, the good agreement we find between CO and Hα shows that concerns about ionized gas as a tracer (e.g., outflows from H ii regions, thickness of the ionized gas disk—see Levy et al. 2019 for further discussion) do not affect the global kinematics enough to alter estimates of the dark matter profile.

Since we are tracing gas, our measurements are potentially susceptible to additional outward forces from gas pressure. If this effect is large enough, ignoring it could cause us to underestimate the enclosed mass at a given radius. While approximate corrections for pressure support exist, they depend on several estimations and assumptions. The limitations to the usefulness of these kind of approximations are significant, and are reviewed in detail in Dalcanton & Stilp (2010). However, Dalcanton & Stilp (2010) also indicate that pressure support is only likely to be relevant for galaxies with rotation speeds below 75 km s⁻¹. Of the galaxies in our final sample, only two (NGC 959 and UGC 12009) fall close to this limit, and in both of these cases we infer a cuspy central density profile. Since correcting for pressure support tends to produce steeper estimates of the density slope, it would not change our overall results for these galaxies. While we cannot draw a stronger conclusion about pressure support without resolved H i data, we can note that Adams et al. (2014) used both stellar and gas tracers to infer the dark matter profiles of a similar sample of dwarf galaxies and found general agreement between the two.

We conclude that the precise choice of tracer (ionized gas, molecular gas, stars) does not significantly affect the derived dark matter profile for late-type galaxies with M_* ≈ 10⁹ $\,{\text{}}{M}_{\odot }$.

There is much existing work studying the dark matter density profiles of low-mass galaxies, and we would like to place our results in the context of some of the larger data sets in order to compare and examine potential trends. The left panel of Figure 11 shows the inner dark matter slope versus stellar mass for a union of samples, described briefly here. We note that different groups use different definitions and methodologies for measuring the inner dark matter slope, so the comparison is somewhat imprecise, but we find it to be suitable for examining general trends.

Zoom In Zoom Out Reset image size

Simon et al. (2005) studied five galaxies using Hα velocity fields to derive rotation curves, which they modeled using a power-law dark matter model. They find a mean of $\beta =0.73$, in excellent agreement with ours. Two galaxies overlap with our sample, NGC 2976 and NGC 5949. For NGC 5949, Simon et al. (2005) determines an inner slope of 0.88, which again agrees well with our value of β* = 0.84. They found a very shallow inner slope of 0.01 ± 0.13 for NGC 2976, which is smaller than our value of β* = 0.33 ± 0.05; however, both studies find this galaxy to host the shallowest slope in their respective samples.

Adams et al. (2014) studied a sample of seven galaxies and measured both stellar and gas kinematics to derive dark matter density profiles. In order to directly compare to our results, we use their gNFW fits to convert the inner slopes to match our definition of β*, finding a mean of ${\beta }_{\mathrm{gas}}^{* }=0.74\pm 0.10$, in excellent agreement with ours. We use these converted values in our discussion and in Figure 11. The conversion typically increased the slope by around 0.15. We observe three of the same galaxies, allowing us to compare the individual slopes for those cases. NGC 2976 agrees well: Adams et al. (2014) finds a gas-traced β* of 0.44 ± 0.18, while we find 0.33 ± 0.05. We find a slope of 0.82 ± 0.08 for NGC 5949, which is consistent with the gas-traced slope of 0.65 ± 0.14 from Adams et al. (2014), as well as with their stellar-traced slope of 0.87 ± 0.11. The third galaxy, NGC 959, received a grade of 3 in our analysis and so no comparison of the inner slopes can be made, though we note the rotation curves were in excellent agreement (see Paper I).

The THINGS group (Walter et al. 2008) used H i observations to study galaxies over a wide range of masses, most of which are beyond the relevant scale for this work. The LITTLE THINGS group followed up THINGS and focused exclusively on 41 dwarf galaxies (Hunter et al. 2012). LITTLE THINGS re-observed four of the THINGS dwarf galaxies, so we include just those galaxies from THINGS that were not re-observed (Oh et al. 2011a). We include the subsample of 26 LITTLE THINGS galaxies studied by Oh et al. (2015), chosen from the total sample due to the regular rotation patterns in their velocity fields. They find a mean inner dark matter slope of 0.32 ± 0.24, indicating much shallower profiles than we find in our sample.

While Lelli et al. (2016) have constructed mass models for the SPARC data set, which includes dwarfs, we are not able to make a direct comparison as they do not quote logarithmic slopes. Di Cintio & Lelli (2016) note that there is evidence for cored profiles but better-quality rotation curves within 1 kpc would be needed to say more.

Collectively, the observational samples shown in the left panel of Figure 11 span a stellar mass range of $\mathrm{log}{M}_{* }/{\text{}}{M}_{\odot }\,\approx 7\mbox{--}9.5$ and over this range we find significant evidence for a trend with stellar mass (p-value of 0.001). The THINGS and LITTLE THINGS galaxies have shallow slopes, generally clustered around β = 0.5 and below, while several studies examining higher-mass galaxies have found steeper slopes. While there is significant scatter in the inner dark matter profiles at a given mass, the overall trend with mass is still apparent. This correlation provides an important quantitative constraint on dark matter models for future work.

An important question arises from the heterogeneity among the samples, relating to both the precise definition of the inner slope as well as the sample selection and method of analysis used. It is natural to be concerned that there might be systematic differences that explain some of the trends in Figure 11.

Regarding the method of measuring the inner slope, our data and the converted data from Adams et al. (2014) use a definition of the inner slope that is averaged over a particular radial range, namely 300–800 pc. The LITTLE THINGS data, however, fits a power law to a radial range that varies by galaxy, with the inner limit of this range typically falling between 100 and 400 pc. However, Oh et al. (2015) show that within their sample the inner slope does not correlate with the radial bounds of the fit. This implies that the mean of the LITTLE THINGS data is not sensitive to small differences in this radial range, and so the overall trend with mass cannot be explained by this mild inconsistency.

In our discussion of the various samples up to this point, we have focused on the difference in stellar mass, but there are other important considerations regarding the observations themselves. In particular, LITTLE THINGS use H i observations, which require careful treatment of beam smearing, and have velocity fields that are analyzed differently. For instance, they adopt tilted-ring models, in which the geometry of the gas distribution is refit at each radius, rather than fixing the geometric parameters throughout the disk as we do (see Paper I for details). There are also physical differences in the velocity fields, which are often more disordered than in higher-mass galaxies, and can sometimes lead to the inference of unphysical density profiles (e.g., density increasing with radius in the inner regions). The rotation curves are also analyzed differently. LITTLE THINGS fits a power law to only the innermost few data points (typically 2–3), which tends to minimize the value of the derived slope, whereas we fit a gNFW profile to the full rotation curve, then evaluate the slope over a fixed radial range. While it is possible that these differences in methodology could have an influence on the mass trend in Figure 11, it is beyond the scope of this paper to fully assess that possibility.

As mentioned previously, the trend with stellar mass seen in the observations in the literature is not entirely unexpected. Qualitatively, we might expect the effects of feedback to be lessened in more massive galaxies, which hints at the possibility of a mass limit above which core creation would be ineffective. It is interesting therefore to make a quantitative comparison with predictions from hydrodynamical simulations.

The right panel of Figure 11 shows the previously discussed results from observations in reduced opacity, overlaid with results from simulations. We include the results from Chan et al. (2015), who use the FIRE simulations to examine the influence of stellar feedback on core creation. They define the inner dark matter slope over the range 300–700 pc, very similar to our range of 300–800 pc used to calculate β*. We also include results from A. Lazar et al. (in preparation), which analyzes simulations run as part of the FIRE and FIRE-2 projects (El-Badry et al. 2017; Fitts et al. 2017; Chan et al. 2018; Hopkins et al. 2018; Wheeler et al. 2019; Garrison-Kimmel et al. 2019; Graus et al. 2019; Samuel et al. 2019). The slopes are measured from 600 to 800 pc, which is similar to our range but uses a larger lower limit to ensure convergence. The 600–800 pc slopes are 0.2 steeper, on average, than the slopes over 300–800 pc in the best resolved simulations of galaxies around M_* ∼ 10⁹. If we were able to resolve 300–800 pc in all the simulations, better matching the observations, it would likely produce even shallower slopes that only enhance the difference with the observations. The dashed line represents the slope from 300–700 pc (matching the range of Chan et al. 2015) for an NFW profile, calculated as in Section 4.3. The results of the simulations suggest that the dark matter slope should steepen with stellar mass over the range M_* ∼ 10⁷−10¹¹ $\,{\text{}}{M}_{\odot }$.

We see a similar trend among the overall literature samples. However, while the simulations generally agree with the inner slopes reported for lower-mass galaxies (M_* ∼ 10⁸ $\,{\text{}}{M}_{\odot }$), they produce slopes that are, on average, shallower for galaxies in our mass range. There are other properties of the simulated galaxies that differ from our sample. In particular, the simulated galaxies are more diffuse and dispersion supported than typical disk galaxies in this mass range (El-Badry et al. 2018, Chan et al. 2019).

Since strong outflow episodes will produce both shallower dark matter profiles and more diffuse galaxies (Di Cintio & Lelli 2016; El-Badry et al. 2016; Fitts et al. 2017), these could both be indications that the effects of feedback are too strong in these simulations around ${M}_{* }\sim {10}^{9}\,\,{\text{}}{M}_{\odot }$. As simulations improve, it will be interesting to compare the dark matter profiles to galaxies with more similar densities and morphologies.

The too-shallow slopes found by the FIRE simulations and also the NIHAO simulations (Tollet et al. 2016), which are consistent with FIRE (Bullock & Boylan-Kolchin 2017), suggest that the feedback in their models may be too effective around M_* ∼ 10⁹ $\,{\text{}}{M}_{\odot }$. Patel et al. (2018) also challenge the trend between specific star formation rate and effective radius predicted by FIRE for galaxies in the range M_* ∼ 10⁹−10^9.5 $\,{\text{}}{M}_{\odot }$, finding the opposite to be true in their observational sample.

It is also the case that some simulations that include baryonic feedback are unable to produce any shallow dark matter profiles. Benitez-Llambay et al. (2018), using the EAGLE simulations, were unable to create cores in dwarf galaxies, as not enough baryons were removed from the central region through feedback. They emphasize the importance of the density threshold for star formation in determining whether or not a core forms.

Beyond focusing on the inner dark matter density slope, Oman et al. (2015) emphasized that there is diversity in the shapes of observed rotation curves, even among galaxies with similar v_max, that is not reflected in the simulations. They emphasize that this dispersion cannot easily be explained through baryon-induced fluctuations to the gravitation potential, calling into question the effectiveness of current simulations.

The trends in dark matter slopes presented in this paper can be an important constraint (in addition to other observables) on feedback models, which must be able to account for a diversity in galaxy dark matter distributions, morphologies, and star formation rates in this mass range.

As mentioned in the introduction, SIDM could also explain the shallower profiles observed in dwarf galaxies. Simulations find that shallow dark matter profiles are a natural consequence of SIDM (Rocha et al. 2013; Elbert et al. 2015; Vogelsberger et al. 2014). Kamada et al. (2017) were able to use SIDM models to fit a variety of rotation curves with similar maximum velocities but very different shapes in the inner regions. Since our sample of galaxies span a fairly narrow range of stellar masses but show a wide range of inner slopes, it will provide a powerful test of whether SIDM models with a similar cross-section can naturally account for the diversity. We plan to examine this in our future work.