Too Big To Fail in Light of Gaia

TL;DR

This paper reveals an anti-correlation between dark matter density and orbital proximity in Milky Way dwarf galaxies, suggesting new insights into the too-big-to-fail problem and galaxy formation.

Contribution

It identifies a correlation between dark matter density and orbital history of dwarf galaxies, challenging previous assumptions and prompting a reevaluation of the too-big-to-fail issue.

Findings

Less dense dark matter in distant dSphs like Fornax and Carina.

Closer dSphs like Draco have higher dark matter densities.

Extended stellar distributions correlate with lower dark matter densities.

Abstract

We point out an anti-correlation between the central dark matter (DM) densities of the bright Milky Way dwarf spheroidal galaxies (dSphs) and their orbital pericenter distances inferred from Gaia data. The dSphs that have not come close to the Milky Way center (like Fornax, Carina and Sextans) are less dense in DM than those that have come closer (like Draco and Ursa Minor). The same anti-correlation cannot be inferred for the ultra-faint dSphs due to large scatter. Including ultra-faints, a trend that dSphs with more extended stellar distributions tend to have lower DM densities emerges. A fresh look at solutions to the too-big-to-fail problem is warranted in light of these observations.

Figures

Tables

Table 1. Table 1: Results from the Bayesian inference of the DM content in the nine bright MW dSphs with the refined Jeans analysis described in the text. From left to right we report in order: the name of the dSph, the data-based estimator for the fourth-order virial shape parameter adopted in our analysis (16- t ​ h 𝑡 ℎ th ,50- t ​ h 𝑡 ℎ th ,84- t ​ h 𝑡 ℎ th percentile), the inferred V max subscript 𝑉 max V_{\rm max} and R max subscript 𝑅 max R_{\rm max} pair (16- t ​ h 𝑡 ℎ th ,50- t ​ h 𝑡 ℎ th ,84- t ​ h 𝑡 ℎ th percentile), χ 2 superscript 𝜒 2 \chi^{2} /d.o.f. related to equation ( 4 ) (yielding a rough estimate on the goodness of the residuals in the fit to stellar line-of-sight velocities), and the difference in the Akaike Information Criterion and Bayesian deviance Information Criterion (see text) between NFW and cISO scenarios. Positive differences of more than few units  (Kass & Raftery, 1995 ) a preference for cISO and negative ones point to a preference for the NFW profile.

MW dSph	$\widehat{{𝝂}_{𝒔\mathrm{𝟏}}}\left[{\mathrm{𝟏𝟎}}^{\mathrm{𝟑}}{\text{km}}^{\mathrm{𝟒}}{\text{/s}}^{\mathrm{𝟒}}\right]$		${\mathrm{𝐥𝐨𝐠}}_{\mathrm{𝟏𝟎}}\left({𝑽}_{\mathrm{𝐦𝐚𝐱}}/\left[\text{km/s}\right]\right)$	${\mathrm{𝐥𝐨𝐠}}_{\mathrm{𝟏𝟎}}\left({𝑹}_{\mathrm{𝐦𝐚𝐱}}/\left[\text{kpc}\right]\right)$	${𝝌}^{\mathrm{𝟐}}/\text{d.o.f.}$	$𝚫\text{AIC}$	$𝚫\text{BIC}$
Ursa Minor (Umi)	${3.73}_{-0.62}^{+0.75}$	NFW	${1.26}_{-0.03}^{+0.07}$	$-{0.19}_{-0.46}^{+0.46}$	0.95	0	0
		cISO	${1.24}_{-0.02}^{+0.03}$	$-{0.75}_{-0.77}^{+0.41}$	0.94
Draco	${3.38}_{-0.46}^{+0.54}$	NFW	${1.42}_{-0.07}^{+0.08}$	${0.50}_{-0.24}^{+0.23}$	1.0	-3	-3
		cISO	${1.36}_{-0.04}^{+0.05}$	${0.03}_{-0.13}^{+0.12}$	0.99
Sculptor	${3.94}_{-0.37}^{+0.43}$	NFW	${1.43}_{-0.08}^{+0.10}$	${0.62}_{-0.25}^{+0.29}$	0.95	8	4
		cISO	${1.35}_{-0.03}^{+0.04}$	${0.06}_{-0.08}^{+0.11}$	0.99
Sextans	${1.84}_{-0.29}^{+0.37}$	NFW	${1.19}_{-0.06}^{+0.07}$	${0.31}_{-0.34}^{+0.27}$	0.95	-1	-1
		cISO	${1.18}_{-0.04}^{+0.07}$	${0.05}_{-0.28}^{+0.21}$	0.95
Fornax	${7.12}_{-0.48}^{+0.52}$	NFW	${1.36}_{-0.05}^{+0.04}$	${0.80}_{-0.19}^{+0.12}$	0.96	10	4
		cISO	${1.33}_{-0.04}^{+0.06}$	${0.45}_{-0.14}^{+0.17}$	0.99
Carina	${1.79}_{-0.25}^{+0.33}$	NFW	${1.28}_{-0.08}^{+0.06}$	${0.57}_{-0.26}^{+0.16}$	1.0	-1	-2
		cISO	${1.14}_{-0.05}^{+0.07}$	$-{0.09}_{-0.32}^{+0.18}$	1.0
Leo II	${1.46}_{-0.32}^{+0.41}$	NFW	${1.16}_{-0.03}^{+0.06}$	$-{0.58}_{-0.49}^{+0.47}$	1.0	1	2
		cISO	${1.17}_{-0.03}^{+0.04}$	$-{1.30}_{-0.46}^{+0.68}$	1.0
Leo I	${3.42}_{-0.56}^{+0.67}$	NFW	${1.42}_{-0.12}^{+0.13}$	${0.70}_{-0.38}^{+0.35}$	1.0	-2	-2
		cISO	${1.33}_{-0.06}^{+0.07}$	${0.12}_{-0.19}^{+0.17}$	0.99
Canes Venatici I (CVnI)	${3.13}_{-0.77}^{+1.33}$	NFW	${1.24}_{-0.10}^{+0.15}$	${0.33}_{-0.60}^{+0.44}$	0.98	0	-1
		cISO	${1.19}_{-0.06}^{+0.15}$	$-{0.09}_{-0.41}^{+0.38}$	0.99

Equations

$$\left(\frac{d\log \left(\nu {\sigma }_r^2\right)}{d\log r}+2\beta \right){\sigma }_r^2=-G_N\frac{M}{r}.$$

$$\beta (r)=\frac{{\beta }_0+{\beta }_{\infty }(r/r_{\beta }{)}^{\eta }}{1+(r/r_{\beta }{)}^{\eta }}.$$

$$\Sigma {\sigma }_{los}^2={\int }_{R^2}^{\infty }\frac{d{r}^2}{\sqrt{r^2-R^2}}\left(1-\beta \frac{R^2}{r^2}\right)\nu {\sigma }_r^2.$$

$${\mathcal{L}}_{v_{los}}=\prod _{i=1}^{N_{\star }}\frac{1}{\sqrt{2\pi \left(\delta v_{los,i}^2+{\sigma }_{los}^2(R_i)\right)}}\exp \left(-\frac{1}{2}\frac{v_{los,i}^2}{\delta v_{los,i}^2+{\sigma }_{los}^2(R_i)}\right),$$

$$v_{1s}\equiv \frac{1}{2}{\int }_0^{\infty }d{R}^2\Sigma ⟨v_{los}^4⟩$$

$$v_{2s}\equiv \frac{1}{2}{\int }_0^{\infty }d{R}^2{R}^2\Sigma ⟨v_{los}^4⟩$$

$$\widehat{v_{s1}}\equiv \frac{1}{2\pi N_{\star }}\sum _{i=1}^{N_{\star }}{v}_{los,i}^4\simeq \frac{G_N}{5}{\int }_0^{\infty }d{r}^{2}M\left(5-2\beta \right)\nu {\sigma }_r^2.$$

$${\mathcal{L}}_{\text{corr}}=\prod _{i=1}^9\frac{\exp \left(-\frac{1}{2}{\chi }_{\text{corr},i}^2\right)}{\sqrt{2\pi \left(\Delta x_i^2+\delta y_i^2\right)}},{\chi }_{\text{corr},i}^2=\frac{{\left(y_i-10^q{x}_i^m\right)}^2}{\Delta x_i^2+\delta y_i^2},\Delta x_i=m10^q{x}_i^{m-1}\delta x_i,$$

Full text

††pubyear: 2019††pagerange: Too Big To Fail in Light of Gaia–C

1 Introduction

We are in an era where the remarkable success of the standard cosmological model with cold and collisionless dark matter (CDM) (Ade et al., 2016) can be tested by observations of dwarf galaxies. The densities probed by dwarfs are sensitive to a wide range of particle physics, including the mass of fermionic DM (Tremaine & Gunn, 1979), DM self-interactions (Spergel & Steinhardt, 2000; Kaplinghat et al., 2016), ultra-light bosonic DM (Hu et al., 2000; Hui et al., 2017) and DM as a superfluid  (Berezhiani & Khoury, 2015). Milky Way (MW) dwarf spheroidal galaxies (dSphs) are a crucial part of this test (Bullock & Boylan-Kolchin, 2017). Photometry and spectroscopy of these galaxies (McConnachie, 2012; Walker, 2013; Battaglia et al., 2013; Simon, 2019) have led to the formulation of the too-big-to-fail (TBTF) problem (Boylan-Kolchin et al., 2011, 2012), based on a key progress in high-resolution simulations of the MW (Madau et al., 2008; Springel et al., 2008) and mass estimators of the dSphs (Walker et al., 2009b; Wolf et al., 2010; Ullio & Valli, 2016; Errani et al., 2018) tested against hydro simulations (Campbell et al., 2016; González-Samaniego et al., 2017).

The dominant interpretation of the TBTF of MW dwarf satellites has been that the CDM N-body subhalos are overly dense. This can be alleviated if the MW halo mass is smaller than $10^{12}\rm M_{\odot}$, but the prevalence of this problem in M31 and in the field (Boylan-Kolchin et al., 2012; Garrison-Kimmel et al., 2014b; Papastergis et al., 2015) makes this an incomplete solution. Another possibility is that the tidal effect of the disk can reduce the inner densities of most of the bright MW satellites (Brooks & Zolotov, 2014; Zolotov et al., 2012). If satellites have a constant density core in DM, then the impact is even larger (Penarrubia et al., 2010). Note that this resolution implicitly assumes dSphs have orbits that take them close to the disk.

Cosmological simulations of the Local Group based on the CDM model, with and without strong feedback effects, have claimed to solve the TBTF problem satisfactorily (Dutton et al., 2016; Sawala et al., 2016; Wetzel et al., 2016). Solutions in the context of non-trivial DM physics including warm DM (Lovell et al., 2017; Bozek et al., 2018), DM-radiation couplings  (Vogelsberger et al., 2016; Schewtschenko et al., 2016), fuzzy DM (Schive et al., 2014) and self-interacting dark matter (SIDM) (Vogelsberger et al., 2012; Zavala et al., 2013) have also been investigated. The SIDM models that alleviate TBTF require $\sigma/m\gtrsim 1~{}\rm cm^{2}/g$ (Zavala et al., 2013), which can also explain the DM densities inferred in field galaxies (Rocha et al., 2013; Kamada et al., 2017; Ren et al., 2018). As with CDM, a complete resolution of the TBTF within SIDM (or other models) will require a deeper understanding of environmental effects (Dooley et al., 2016; Valli & Yu, 2018).

In this Letter, we provide a new observational handle on this issue. We first infer the DM content in the bright MW dSphs through a spherical Jeans analysis (Battaglia et al., 2013) supplemented with fourth-order projected virial shape parameters (Merrifield & Kent, 1990; Richardson & Fairbairn, 2014; Read & Steger, 2017) and up-to-date photometric information (Munoz et al., 2018). We find that both cuspy Navarro-Frenk-White (NFW) profile (Navarro et al., 1997) and cored isothermal profile (cISO) can fit the data well and the results are summarized in figure 1. In neither case, do we have to impose cosmological priors to infer the DM densities. For both profiles, we recover the well-known TBTF problem.

We show in figure 2 that the inferred inner density of the bright MW dwarfs is tightly anti-correlated with the MW pericenter distances (Fritz et al., 2018), estimated from the most recent observational data provided by the Gaia collaboration (Brown et al., 2018; Helmi et al., 2018). We do not find the same strong anti-correlation for the ultra-faint dSphs; see figure 3. However, all the dSphs taken together show a clear trend of decreasing DM densities at larger half-light radii in figure 4. Our results provide a new perspective for understanding the formation and evolution of the MW dSphs, as well as a constraint on solutions to the TBTF problem. We discuss possibilities that could reproduce the observed trends in both CDM and SIDM scenarios.

2 Jeans analysis and fourth-order virial theorem.

We adopt the standard spherical Jeans analysis of dSph stellar kinematics (Binney & Tremaine, 2008), and refine it with the inclusion of estimators for fourth-order projected virial theorems (Merrifield & Kent, 1990; Richardson & Fairbairn, 2014). We apply this to the bright dSphs (V-band luminosity $>10^{5}\ \rm L_{\odot}$ as in Boylan-Kolchin et al. (2011)), which include the eight classical dSphs and Canes Venatici I (Zucker et al., 2006). The number of available spectroscopic members greatly exceeds a hundred stars for these nine dSphs, which allows for a detailed Jeans analysis. We analyze the corresponding samples in refs. (Simon & Geha, 2007a; Mateo et al., 2008; Walker et al., 2009a, 2015; Spencer et al., 2017), following the methodology in (Walker et al., 2006; Strigari et al., 2008b, a; Martinez et al., 2009) and in the more recent works (Pace & Strigari, 2019; Petac et al., 2018). We perform a Bayesian fit (Foreman-Mackey et al., 2013; Goodman & Weare, 2010) of dSph line-of-sight velocities, varying a total of six parameters: the pair of $V_{\rm max}$ and $R_{\rm max}$, i.e., the maximal circular velocity of the halo and its associated radius, fully determining the DM mass profile, and four parameters to model the stellar orbital anisotropy as in Baes & Van Hese (2007). We describe the stellar density profile by the Plummer model (Plummer, 1911) with half-light radii from the recent photometric collection in Munoz et al. (2018). To impose projected virial theorems, we only use the “$v_{s1}$" estimator, proportional to $\langle v_{los}^{4}\rangle$, developed in refs. (Richardson & Fairbairn, 2014; Read & Steger, 2017). More details specific to our analysis are provided in appendix A.

The mass of MW dSphs is inferred under two separate assumptions for the DM density profile: i) cuspy NFW, $\rho=\rho_{s}(r_{s}/r)(1+r/r_{s})^{-2}$, which encapsulates the predictions of N-body CDM simulations, and ii) cored isothermal (cISO) profile, which is the prediction of the SIDM model for moderate cross sections when baryons are dynamically irrelevant (Kaplinghat et al., 2016). Both profiles have two free parameters. For the cISO profile, they are the central density $\rho_{0}$ and the constant velocity dispersion $\sigma_{0}$. These two determine the DM density profile $\rho(r)$ through the Poisson equation: $\sigma^{2}_{0}\mathop{}\!\mathbin{\bigtriangleup}\ln\rho=-4\pi G_{N}\rho$ with the “cored" boundary conditions at $r=0$: $\rho=\rho_{0}$ and vanishing spatial derivative. Note that both NFW and cISO parameters can be related to ($V_{\rm max}$,$R_{\rm max}$) pair by: i) $\rho_{s}\simeq 1.721{V_{\rm max}^{2}}/({G_{N}R_{\rm max}^{2}}),~{}r_{s}\simeq 0.462R_{\rm max}$; ii) $\rho_{0}\simeq 2.556{V_{\rm max}^{2}}/({G_{N}R_{\rm max}^{2}}),~{}\sigma_{0}\simeq 0.630V_{\rm max}$.

3 The too-big-to-fail problem revisited

In figure 1, we show preferred $V_{\rm max}$ and $R_{\rm max}$ values for the classical MW satellites from our NFW (upper) and cISO (lower) fits, compared to prediction for the 50 most massive subhalos (up-pointing triangles for the most massive ones, i.e. $V_{\rm peak}>30$ km/s, and down-pointing triangles for the others) in the high-resolution ELVIS CDM simulation named Kauket, which has the lowest host-halo mass in the suite (Garrison-Kimmel et al., 2014a). While the choice of the latter goes in the direction of minimizing the TBTF problem, it does not solve it. As figure 1 shows, the host halo profiles inferred from both NFW and cISO fits are systematically less massive than those predicted in CDM simulations. It should be noted that both NFW and cISO scenarios provide equally good fits to dSph line-of-sight velocity datasets, while respecting the global constraint from the fourth order projected virial theorem. Details about our fits are in appendices A and B, where we also exploit two notions of information criterion (Akaike, 1974; Gelman et al., 2013) for model comparison (Kass & Raftery, 1995).

The TBTF problem for the cISO profile is particularly interesting in connection with the study in Valli & Yu (2018), where the classical MW dSphs were analyzed using the SIDM halo model of Kaplinghat et al. (2014, 2016). The inferred $\sigma/m$ spanned a large range from $0.1\textup{--}0.3~{}{\rm cm^{2}/g}$ (Ursa Minor and Draco) to values greater than $10~{}{\rm cm^{2}/g}$ (Fornax and Sextans), with a CDM concentration-mass prior from Vogelsberger et al. (2016). This hierarchy of $\sigma/m$ values reflects the diversity of the DM content inferred for the host subhalos, which is not evident in the field halos (Kamada et al., 2017; Ren et al., 2018). Thus, to fully address the MW TBTF problem in CDM or SIDM models, we must understand the physics which is unique to MW dSphs, but is not present for field galaxies. As recently investigated in Robles et al. (2019), the interactions of dwarf satellites with the MW disk could play a crucial role in reducing the number of satellites, removing mass from the outskirts of the subhalos and reducing its inner densities both in CDM and SIDM context. Then, it would be natural to expect the dynamical properties of MW dSphs to be correlated with MW pericenter distances, as we discuss next.

Pericenters from Gaia: a novel diagnostics. The tidal field of the MW could create a correlation between the central density of a satellite and its pericenter passages (Gnedin et al., 1999; McConnachie & Irwin, 2006; Read et al., 2006; Dooley et al., 2016). In figure 2, we show the inferred DM densities at $150$ pc, $\rho_{150}$, from our fits versus the pericenter distances (Fritz et al., 2018), $r_{\rm P}$, estimated from the recent Gaia data (Brown et al., 2018; Helmi et al., 2018) for a MW model with mass of $0.8\times 10^{12}$ M*⊙*. For both NFW and cISO fits, there is a tight correlation between the central densities and the pericenter distances of the satellites, in such a way that the denser dSphs have small pericenter distances. The correlation persists for a MW model that is twice as heavy or if we adopt density measurements from Read et al. (2018),while it is no longer present if one trades the pericenter distance with the dSph heliocentric distance, which has been discussed recently, see Hammer et al. (2018). These checks and more details on the $\rho_{150}\textup{--}r_{\rm P}$ correlation are available in appendix C, where we also show the lack of correlation of $\rho_{150}$ with stellar mass and star formation shutoff time.

The existence of this correlation indicates that the galactic tides play a significant role in shaping the host subhalos of the bright MW bright dSphs, which is to be expected (Brooks & Zolotov, 2014). It is, however, surprising that the central DM density is anti-correlated with the pericenter distance because the tidal effects become more significant as the pericenter distance decreases. One might expect a “survivor" bias to higher densities in subhalos that venture closer to the MW, but this is not evident in the N-body simulations for the most massive subhalos (upward-pointing triangles in figure 2). Adding a disk will preferentially reduce the densities of the smaller $r_{\rm P}$ subhalos and change the orbits. This does not seem to lead to the required correlation (Robles et al., 2019), but more investigation is needed to firmly establish this point.

Note that $150~{}\rm pc$ is chosen as a compromise between the expectations driven by the dynamics (which constrains the mass within the half-light radius best (Walker et al., 2009b; Wolf et al., 2010)) and the need for a common scale at where dSph densities can be compared. From figure 2, we see that Fornax, Carina and Sextans have lower $\rho_{150}$ values in cISO fits than NFW ones, because their half-light radii are considerably larger than $150$ pc. For the other bright dSphs half-light radii are closer to $150$ pc and consequently $\rho_{150}$ is less sensitive to the choice of the density profile.

4 Including ultra-faint dSphs

A natural question is whether the MW ultra-faint dSphs (Simon & Geha, 2007b; Martin et al., 2007; Koposov et al., 2015) follow the trend set by the bright dSphs. Here we focus on ultra-faint dSphs that have a pericenter distance measurement (Fritz et al., 2018; Simon, 2018), $>2\sigma$ evidence for non-zero average velocity dispersion, and no proposed association with the Large Magellanic Cloud (Walker et al., 2016; Kallivayalil et al., 2018). From this list, we discard four dSphs: Aquarius II and Pisces II due to a small number of spectroscopically-confirmed members; Willman I due to the large uncertainty on the dynamical-equilibrium condition (Willman et al., 2011); Eridanus II given the chance that it may not be a MW-bound satellite (Li et al., 2017; Fritz et al., 2018). Finally, we also include the recently discovered dSph Antlia II (AntII), with dispersion and pericenter distance estimates from Torrealba et al. (2018). To infer $\rho_{150}$ conservatively we fit only to the mass estimator of Wolf et al. (2010), using measurements and uncertainties tabulated in Munoz et al. (2018), except for Reticulum II, for which we follow Minor et al. (2018), that includes corrections due to binary systems. Note that Boötes I may show some evidence for multiple kinematically distinct populations (Koposov et al., 2011), not accounted for in the adopted average velocity dispersion. For the NFW profile, we additionally impose the concentration-mass relation in Moliné et al. (2017) by adding a statistical weight: \chi^{2}_{\rm cosmo}=\big{(}0.905-\{\log_{10}(c_{200})+[0.146-0.101\,\log_{10}(M_{200}\,10^{-12}/h)]\}/0.3\big{)}^{2}.

Figure 3 shows the inferred $\rho_{150}$ versus pericenter for the selected MW ultra-faint dwarfs. Some of them have significantly different $\rho_{150}$ between NFW and cISO fits because the half-light radii quite off from $150~{}{\rm pc}$. Our key result with regard to the ultra-faint sample is the larger scatter in the $\rho_{150}\textup{--}r_{\rm P}$ plane, compared to the anti-correlation for the bright MW dSphs. Such a scatter could be due to increased scatter in the stellar-to-halo mass relation at the ultra-faint end and/or reflective of underestimated errors for the ultra-faint dSphs, given their low velocity dispersions and smaller spectroscopic datasets. It may be also connected to DM physics, in particular to DM self-interactions, as discussed later.

The large scatter is evident even without Crater II (CraII) and AntII, characterized by their extremely low DM density and large stellar extent (Caldwell et al., 2017; Torrealba et al., 2016; Torrealba et al., 2018). This begs the question of why CraII and AntII are so diffuse and underdense (in both DM and stars) compared, for instance, to Draco, whose luminosity and pericenter are similar to those of the two ultra-faint dSphs.

To get a global view of the stellar sizes and mass densities of the dSphs, we show in figure 4 the mean density at the half-light radius, $\langle\rho_{1/2}\rangle\equiv M_{1/2}/(4/3\pi r_{1/2}^{3})$, versus the stellar half-light radius for all MW dSphs considered in this work. Note that the inference of $\langle\rho_{1/2}\rangle$ is very robust against the assumed density profile. Consequently, for the bright dSphs we show only NFW fits. For the ultra-faint dSphs we directly derive $\langle\rho_{1/2}\rangle$ from the estimator in Wolf et al. (2010) using what tabulated in Simon (2019).

In figure 4 dashed gray lines represent the outcome from CDM ELVIS simulation Kauket (Garrison-Kimmel et al., 2014a). From the figure one can clearly notice how more extended dSphs favor lower-mass subhalos compared to more compact objects. This interesting trend – seemingly independent from the pericenter information – constitutes another unexplained facet of the MW satellites. The TBTF problem, captured in the upper panel of figure 1, also stands out in the $\langle\rho_{1/2}\rangle\textup{--}r_{1/2}$ plane of figure 4. In what follows, we discuss the TBTF problem in light of the $\rho_{150}\textup{--}r_{\rm P}$ correlation inferred for the bright MW dSphs.

5 The TBTF problem in CDM

To explain TBTF for the bright MW dSphs in light of the Gaia data using collisionless DM models, one must host Draco and Ursa Minor in the densest subhalos with pericenter distances close to 30–40 kpc, while explaining why Carina, Fornax and Sextans are hosted by lower-density subhalos that did not come close to the MW disk. If the TBTF is solved by reducing the MW mass, then the impact of the disk on the subhalos with smaller $r_{\rm P}$ may render them inconsistent with Draco or Ursa Minor stellar kinematics. Indeed, the effect of a disk in N-body simulations is to reduce the inner densities of the small $r_{\rm P}$ subhalos, but this would leave the other subhalos still too dense compared to data, see e.g. Robles et al. (2019). Therefore, a full resolution of the TBTF may likely involve physics beyond the tidal impact of the disk.

Feedback from supernovae is one possibility, but its role in setting or diluting the $\rho_{150}-r_{\rm P}$ anti-correlation is unclear. Most of the cosmological simulations indicate that below about $10^{7}~{}\rm M_{\odot}$ in stellar mass, feedback is not efficient in changing the density distribution at $1\textup{--}2\%$ of virial radius (Bullock & Boylan-Kolchin, 2017). For subhalo masses around $\rm M_{200}=10^{9}~{}\rm M_{\odot}$, the relevant radius is about $300~{}\rm pc$. So it may be still possible for feedback to change the inner density at $150$ pc. However, we have verified that a correlation with $r_{\rm P}$ can still be assessed for densities inferred at $300~{}\rm pc$. In addition, Draco, Ursa Minor, Carina and Sextans have similar V-band luminosities (below $10^{6}\ \rm L_{\odot}$) but they bracket the range of inferred densities. This would indicate that there is another variable that controls the impact of feedback on DM halos, for example, the “burstiness" of star formation. In this regard, we observe only a weak trend of our inferred $\rho_{150}$ values with the star formation shutoff times collected in Read et al. (2019), see the appendices. Interestingly, the difference in the $\rho_{150}$ seen in the sample of dwarfs in Read et al. (2019) may reflect the differences in the satellite versus field populations, with the field population having lower $\rho_{150}$. Clearly, more work along these lines is warranted.

6 The TBTF problem in SIDM

SIDM models with cross section over mass $\sigma/m\gtrsim 1~{}\rm cm^{2}/g$ have been investigated as promising solutions of the TBTF problem (Vogelsberger et al., 2012; Zavala et al., 2013). In this regime, a simple model has emerged for the SIDM halo profile (Kaplinghat et al., 2014, 2016; Elbert et al., 2018; Sameie et al., 2018) that explains the diversity of galaxy rotation curves across the full range of galaxy masses (Kamada et al., 2017; Ren et al., 2018; Creasey et al., 2017). However, this halo model does not single out a consistent SIDM cross section when addressing the TBTF problem in the MW (Valli & Yu, 2018).

Recently, Nishikawa et al. (2019) suggested that for $\sigma/m\gtrsim 5\ \rm cm^{2}/g$ some of the MW satellites may be in the core-collapse phase, i.e., the central density would be increasing with time (Balberg et al., 2002; Elbert et al., 2015). Interestingly, this process is naturally correlated with how close the satellite comes to the center of the MW. In fact, according to this scenario, Draco would be in the core-collapse phase and hence dense. On the other hand, Fornax, Carina and possibly Sextans would be in the core-expansion phase, similar to the field halos in their evolutionary phase but the central density would be higher due to tidal-mass loss (Nishikawa et al., 2019). We point out that for a given pericenter, a high-concentration subhalo is more likely to experience core collapse and develop a higher central density than a low-concentration one. This is due to the fact that the core-collapse time scale, $t_{c}$, depends on concentration, $c_{200}$, as $t_{c}\propto(\sigma/m)^{-1}M^{-1/3}_{200}c^{-7/2}_{200}$ (Essig et al., 2018; Nishikawa et al., 2019). For $\sigma/m\simeq 5~{}\rm cm^{2}/g$, the relevant time scale can be much larger than the age of the Universe or comparable to it depending on the concentration of the MW subhalo. The low-concentration subhalo, which is not undergoing core collapse, is more vulnerable to tidal effects, further reducing the core density (Penarrubia et al., 2010). For higher concentration subhalos, we expect the density to be larger due to core collapse. This interplay between tidal effects and thermalization of the inner halo due to DM self-interactions is a generic prediction of SIDM models with fairly large cross sections.

7 Conclusions

We have fitted stellar kinematics of the nine most luminous MW dSphs and reassessed the TBTF problem for both cupsy and cored DM density profiles. We found a strong anti-correlation between the inner densities of the bright dSphs and their pericenter distances. The ultra-faint dSphs show a much larger scatter and the presence of an anti-correlation (if any) is muted. These results show that the TBTF problem in light of the Gaia data is still a challenge. We argued that the tidal field of a disk is important but bursty star formation and consequent feedback on DM subhalos may be required in solutions based on CDM models. In SIDM models, MW dSphs could be in either core-expansion and core-collapse phases, depending on their concentrations and pericenter distances. Whether this predicted diversity in the DM content of SIDM subhalos is consistent with data remains to be seen. A more precise knowledge of the central DM densities, orbital motions, stellar sizes and star-formation histories for the MW dSph satellites promises to provide an incisive test of the DM nature.

Acknowledgements

We are very grateful to Andrew Pace and Matthew Walker for their help with the data and useful discussions. We thank Mike Boylan-Kolchin, James Bullock, Shea Garrison-Kimmel and Joshua Simon for interesting comments. M.K. and M.V. are supported by the NSF Grant No. PHY-1620638. H.B.Y. acknowledges support from U.S. Department of Energy under Grant No. de-sc0008541 and UCR Regents’ Faculty Development Award.

Appendix A Jeans analysis refined by fourth-order virial shape parameters

The spherical Jeans analysis has been extensively used in literature to analyze the dynamics in MW dSphs via the measured stellar kinematics; see refs. (Bonnivard et al., 2015b; Petac et al., 2018) for recent analyses reviewing in detail this approach. Under spherical symmetry approximation,111This is well-motivated by MW dSph ellipticities inferred from photometric data, see for instance Munoz et al. (2018). Small systematics related to departure from spherical symmetry have been assessed via numerical simulations (Bonnivard et al., 2015a; Campbell et al., 2016; González-Samaniego et al., 2017; Errani et al., 2018), and Schwarzschild modeling (Jardel et al., 2013; Breddels & Helmi, 2013; Kowalczyk et al., 2017, 2018, 2019). Axi-symmetric Jeans analysis (Cappellari, 2008) has been also performed under simple specific ansatz on the stellar orbital anisotropy (Hayashi & Chiba, 2012; Hayashi et al., 2016). Jeans equations (Binney & Tremaine, 2008) reduce to:

[TABLE]

with $\nu$ being the stellar density profile of the system, $\beta$ the stellar velocity dispersion anisotropy and $M$ is the total mass profile, which is well-approximated by the DM component for dSphs (Walker, 2013; Battaglia et al., 2013). The parameterization of $\beta$ adopted from Baes & Van Hese (2007) is:

[TABLE]

Solving equation (1) for $\sigma_{r}^{2}$, one can make contact with observations through the line-of-sight projection:

[TABLE]

In the above, $\Sigma$ is the surface brightness of the system, related to $\nu$ via an Abel transform (Wolf et al., 2010; Ullio & Valli, 2016) and constrained by available photometric data, while $\sigma_{los}$ is the line-of-sight velocity dispersion profile of stars inferred from spectroscopy of individual stars in the dSph. In this work, we model the surface brightness with a Plummer model, i.e. $\Sigma\propto(1+R/R_{1/2})^{-2}$. The structural parameters of the Plummer model, projected half-light radii $R_{1/2}$ and ellipticities, are taken from the recent work of Munoz et al. (2018) (see also Simon (2019)). We modify the semi-major half-light radius by a factor of $\sqrt{1-\epsilon}$ to be consistent with our assumption of spherical symmetry. In the present Jeans analysis for the brightest MW dSphs we also include the contribution of the stars to the potential well of the system for the objects where such an effect turns out to be more important than a few percent level for the inference of $\rho_{150}$. This is actually the case for Fornax, Leo I and Sculptor, for which we adopt the total stellar mass reported in Read et al. (2019), describing the stellar contribution to the internal dynamics of the system with the Plummer model.

The sample of line-of-sight velocities for the nine bright MW dSphs is publicly available from refs. (Simon & Geha, 2007a; Mateo et al., 2008; Walker et al., 2009a, 2015; Spencer et al., 2017) (for Ursa Minor data, M.G. Walker, private communication) and using stellar-membership cuts obtained in Pace & Strigari (2019) (A.Pace, private communication). We adopt dSph RA-Dec coordinates and heliocentric distance from Munoz et al. (2018) and infer the average bulk velocity $\overline{V}$ and dispersion $\overline{\sigma}_{los}$ using the Gaussian estimator proposed in Walker et al. (2006). We then subtract the estimated $\overline{V}$ from the line-of-sight velocities in order to remove the bulk motion of the dSph.

The test statistic to model the probability of obtaining the dataset $\{v_{los}\}_{i=1,\dots,N_{\star}}$ given the model parameters involved in our analysis corresponds here to the Gaussian likelihood (Strigari et al., 2008b; Martinez et al., 2009):

[TABLE]

where $\delta v_{los,i}$ and $R_{i}$ are the velocity measurement error and position of the $i$-th star, and $\sigma_{los}^{2}(R_{i})$ corresponds to the intrinsic dispersion predicted by equation (3). We use the public package emcee (Foreman-Mackey et al., 2013) in order to perform Monte Carlo Markov Chain (MCMC) sampling. Our theoretical model comprises of 6 parameters: the four characterizing $\beta$ in equation (2), and ($V_{\rm max}$, $R_{\rm max}$) pair determining the DM halo profiles considered in this work. For the stellar anisotropy $\beta$, we adopt the priors of Valli & Yu (2018), which take into account the conditions stemming from the requirement of non-negative phase-space distribution functions. For the DM halo, we impose the broad priors: $-2\leq\log_{10}(V_{\rm max}/[\textrm{km}/\textrm{s}])\leq 2$ and $-2\leq\log_{10}(R_{\rm max}/\textrm{kpc})\leq\log_{10}(r_{\textrm{J}}/\textrm{kpc})$. We demand that $R_{\rm max}$ be smaller than the instantaneous Jacobi radius $r_{\textrm{J}}$, which is a reasonable description of $R_{\rm max}$ for subhalos in CDM simulations (Springel et al., 2008). This restriction is mostly relevant for the NFW profile and it removes some of the large $R_{\rm max}$ solutions, which is evident in the contour shape at high $R_{\rm max}$ in figure 1.

We set the Jacobi radius $r_{\textrm{J}}\equiv(G_{N}M_{\textrm{dSph}}D^{2}/(2\sigma_{\textrm{MW}}^{2}))^{1/3}$ similarly to what was done previously (e.g., ref (Strigari et al., 2008a)). We approximate the distance $D$ of the satellite from the MW by the dSph heliocentric distance, and take 200 km/s as an estimate for the MW velocity dispersion $\sigma_{\textrm{MW}}$. Since the restriction on $R_{\rm max}$ does not have a significant impact on the physics discussion, we use estimates of the total mass of the dSphs $M_{\textrm{dSph}}$ from Errani et al. (2018) for the cored and cuspy cases. Finally, we note that for the cISO fit of CVn I, a second mode is found at extremely low $R_{\rm max}$ values: we have cut this second solution restricting for this case the prior on $R_{\rm max}$ to values $\log_{10}R_{\rm max}\geq-1$.

We refine the standard Jeans analysis with the constraints from the fourth-order projected virial theorems (Merrifield & Kent, 1990) for approximately spherical systems (Richardson & Fairbairn, 2014):

[TABLE]

where $\langle v_{los}^{4}\rangle$ denotes the fourth moment of the line-of-sight velocity distribution. The left-hand side of equation (5) involves only observational quantities and is determined by the following means $\{v_{los,i}^{4}\}_{i,\dots,N_{\star}}$ and $\{R_{i}^{2}v_{los,i}^{4}\}_{i,\dots,N_{\star}}$ (Read & Steger, 2017).

We compute mean and variance for the estimators of $v_{1s}$ and $v_{2s}$ by generating $10^{6}$ sets of the samples $\{v_{los,i}^{4}\}_{i,\dots,N_{\star}}$ and $\{R_{i}^{2}\,v_{los,i}^{4}\}_{i,\dots,N_{\star}}$ via the distribution encoded in equation (4) for $v_{los}$, namely a Gaussian with mean 0 and variance approximated by $\overline{\sigma}_{los}^{2}$. Another possibility suggested by Read & Steger (2017) would be to directly compute the surface brightness integral in equation (4) within a given functional parameterization for $\Sigma$ and a data-based interpolation for the radial profile of $\langle v_{los}^{4}\rangle$. We tested this alternative procedure and found the estimate of $v_{1s}$ to be quite stable against the method adopted. The estimate of $v_{2s}$, instead, turns out to be strongly sensitive to the assumed outer behavior for $\Sigma$ and $\langle v_{los}^{4}\rangle$ profiles at large radii (e.g.: for a Plummer surface brightness profile and a mildly varying $\langle v_{los}^{4}\rangle$ profile, $\nu_{2s}$ would formally diverge), implying there is a large uncertainty on the corresponding data-based estimator. Therefore, we only use constraints from $v_{1s}$:

[TABLE]

From our Monte Carlo sampling of $v_{los}$, the data-based estimator $\widehat{v_{s1}}$ can be reasonably approximated by a Gaussian distribution. We report the 16-$th$, 50-$th$ and 84-$th$ percentiles of the inferred distribution in table 1. We then use the estimated mean and variance for $\widehat{v_{s1}}$ to build a Gaussian weight, $\mathcal{L}_{v_{s1}}$, in our fit that supplements the observational information from the likelihood in equation (4) and refines our Jeans analysis according to equation (6).

In the same table 1 we show the results from our Bayesian inference of $(V_{\rm max},R_{\rm max})$ pair for both cuspy and cored DM profiles considered in this work, reporting again 16-$th$, 50-$th$ and 84-$th$ percentiles of the related posterior distribution. For a rough estimate of the goodness of the fit, we compute $\chi^{2}/\textrm{d.o.f.}\equiv-2\log\widetilde{\mathcal{L}}_{v_{los}}/(N_{\star}-(6-1))$, where $\widetilde{\mathcal{L}}_{v_{los}}$ corresponds to the product of the Gaussian distributions in equation (4) without the normalization factor. As shown in table 1, the constraint from $\eqref{eq:fourth_order}$ does not forbid an overall good fit to dSph stellar kinematics in any of the cases considered. To compare the cored and cuspy density profile fits, we estimate the information criteria: (i) the frequentist Akaike Information Criterion (AIC) (Akaike, 1974), related to the log-likelihood $\log\mathcal{L}_{v_{los}+v_{1s}}$ at the global mode of the fit, and (ii) the Bayesian deviance Information Criterion (BIC) (Gelman et al., 2013), obtained from the mean and variance of the posterior distribution of $\log\mathcal{L}_{v_{los}+v_{1s}}$. Note that, differently from the$\chi^{2}/\textrm{d.o.f.}$ defined above, both AIC and BIC are sensitive to the goodness of fit to the line-of-sight velocity dispersion and the global constraint from the estimated fourth-order virial-shape parameter $\nu_{s1}$. For both AIC and BIC, the standard scale of evidence reported e.g. in Kass & Raftery (1995) suggests that differences (in absolute value) of $\mathcal{O}(10)$ or greater in the information criteria point to data-based preference for the model with the lower value of AIC or BIC.

Interestingly, we observe that $\Delta\textrm{AIC}\equiv\textrm{AIC}_{\rm NFW}-\textrm{AIC}_{\rm cISO}$ seems to highlight some mild preference for the cored halos in Sculptor and Fornax. We verified that this would no longer be the case if we would remove the information stemming from $v_{s1}$ in the corresponding fits. At the same time, $\Delta\textrm{BIC}\equiv\textrm{BIC}_{\rm NFW}-\textrm{BIC}_{\rm cISO}$, while confirming the trend highlighted by the AIC, never exceeds $\mathcal{O}(1)$ values, pointing to the fact that the assumed cored and cuspy halo profiles provide roughly the same level of goodness in the description of dSph stellar kinematics. We conclude that the inclusion of the information from $\nu_{s1}$ in the spherical Jeans analysis does not effectively break the well-known $M$ - $\beta$ degeneracy for MW dSphs.

Finally – in order to highlight the effects of fourth-order projected virial theorems in our analysis – we show in figure 5 the 68% and 95% highest probability density region for $(V_{\rm max},R_{\rm max})$ with and without the information form the fourth-order virial shape parameters. The results for both cored and cuspy halo profiles are shown in figure 5.

Appendix B Too-big-to-fail for cuspy and cored halo profiles

We have seen how cuspy and cored DM halos in our detailed investigation lead to the very same quality of the fit on dSph stellar kinematics in the absence of any cosmological prior imposed on ($V_{\rm max}$,$R_{\rm max}$) pair. Here we reinforce the picture already drawn in the main body of this manuscript about how both cuspy and cored DM scenarios are both subject to the too-big-to-fail problem. We present in figure 6 the $68\%$ region of the posterior distribution of the circular velocity profile obtained for both cuspy and cored fits. In the same figure, we show the band encompassing the prediction for the circular velocity profiles of the most massive subhalos ($V_{peak}>30$ km/s) from the CDM ELVIS simulation Kauket (Garrison-Kimmel et al., 2014a). We also present the mass estimator data point $V_{1/2}$ at the deprojected half-light radius (Wolf et al., 2010). In order to compare with the literature in a quite agnostic manner, we update the values reported in Wolf et al. (2010) using the average line-of-sight velocity dispersions collected in Simon (2019) (our own $\overline{\sigma}^{2}_{los}$ results compatible within errors) and, most importantly, including the new photometric dataset from (Munoz et al., 2018). Figure 6 shows the tension between the CDM N-body prediction and the outcome of our fits with cored and cuspy DM halos for the classical dSphs. It also shows that the mass estimator works well for both profiles. This is a non-trivial check since the mass estimator was derived with only information from the dispersion and not the fourth-order projected virial theorem.

Appendix C Details and cross-checks on the $\rho_{150}\textup{--}r_{\rm P}$ correlation

We provide here further details and validation of the correlation discussed in the main text for the nine bright MW dSphs, involving a relation between DM densities at 150 pc, $\rho_{150}$, and their pericenter distance, $r_{\rm P}$, estimated from the Gaia data (Fritz et al., 2018). To formally establish such correlation, we perform a power-law fit for the nine bright MW dSphs based on the basis of the following test statistic:

[TABLE]

where for each of the objects: $x=r_{\rm P}/\textrm{kpc}$, $y=\rho_{150}/(10^{7}M_{\odot}\textrm{kpc}^{-3})$ and $\delta x$, $\delta y$ are the estimated errors. In order to infer $\rho_{150}$ and $\delta\rho_{150}$ we compute the 16-$th$, 50-$th$ and 84-$th$ percentile of the posterior distribution of the DM density at 150 pc, inferred from the Bayesian analysis of classical dSphs described in detail in appendix A. The quoted values for $r_{\rm P}$ and $\delta r_{\rm P}$ are taken from Fritz et al. (2018), who provide this information for both a low-mass and high-mass MW model. For both density and pericenter distance we symmetrize errors and shift median values accordingly. We perform a MCMC analysis using again the affine-invariant sampling algorithm of Goodman & Weare (2010), implemented in the package emcee (Foreman-Mackey et al., 2013).

In figure 7 we show the outcome of our inference on $q$ and $m$ parameters appearing in equation (7). For both NFW and cISO cases, there is a tight $\rho_{150}$ - $r_{\rm P}$ anti-correlation: no correlation case is well beyond the 99.99$\%$ highest probability density region (the largest contour drawn for the joint posterior probability distribution in the same figure). In both the low-mass and high-mass MW cases, we observe that the analysis with the NFW profile yields a more precise determination of two parameters. We also show the results of adopting the pericenter distances corresponding to the high-mass MW model from Fritz et al. (2018), as a point of comparison for the results in the main text. The same information that is in figures 2 - 3 is plotted in the left and central panel of figure 7 for the MW model with mass $1.6\times 10^{12}$ M*⊙*. We note the inference of this correlation does not critically depend on the assumed MW model.

In the same figure we provide a direct comparison of the inferred densities at 150 pc with the ones obtained in Read et al. (2019). CVn I has not been included in that study and therefore we restrict to the other eight bright dSphs. Figure 8 shows that the inferred $\rho_{150}$ for the classical dSphs in this work are in good agreement with the inferences of Read et al. (2019).

Finally, in figure 9 we provide further tests of the correlation inferred for the nine bright dSphs. In the left panel of the figure, we show our $\rho_{150}$ estimate versus the star formation shutoff time presented in Read et al. (2019); in the central panel we inspect the possible correlation of dSph central densities with their stellar mass counterpart; in the right panel of figure 9 we replace the pericenter distance with the heliocentric distance. While the first two panels may be informative with respect to what investigated in Read et al. (2018), the last one may be useful in light of the recent claims provided in Hammer et al. (2018). As mentioned in the main text, we highlight that no strong correlation emerges from the comparison of the central density of the bright dSphs against these three different quantities.