A low pre-infall mass for the Carina dwarf galaxy from disequilibrium modelling

Dark matter-only simulations of galaxy formation predict many more subhalos around a Milky Way-like galaxy than the number of observed satellites. Proposed solutions require the satellites to inhabit dark matter halos with masses 109–1010 Msun at the time they fell into the Milky Way. Here we use a modelling approach, independent of cosmological simulations, to obtain a pre-infall mass of Msun for one of the Milky Way's satellites: Carina. This determination of a low halo mass for Carina can be accommodated within the standard model only if galaxy formation becomes stochastic in halos below ∼1010 Msun. Otherwise Carina, the eighth most luminous Milky Way dwarf, would be expected to inhabit a significantly more massive halo. The implication of this is that a population of ‘dark dwarfs' should orbit the Milky Way: halos devoid of stars and yet more massive than many of their visible counterparts.

W hile the Cold Dark Matter paradigm for structure formation in the Universe has been very successful in reproducing observations on scales larger than B1 Mpc (refs 1-3), on galactic scales there have been long-standing puzzles. The discrepancy between the predictions of the Cold Dark Matter paradigm and the observed properties of the dwarf spheroidal galaxy satellites (dSphs) of the Milky Way has persisted for over a decade. Numerical models predict that thousands of dark matter subhalos should be found orbiting the Milky Way and Andromeda, yet only a few tens have been found to date 4,5 . This has become known as the 'Missing Satellites' problem. A popular resolution of this issue is to place stars only in the most massive satellite halos, implying a total or 'virial' mass for the Milky Way dSphs of B10 10 Msun (see Fig. 1). However, in such halos, the required central stellar velocity dispersions of the dSphs would be too high to be consistent with the Milky Way dSphs 6 . More recent work has shown further that more refined mappings between luminous and dark matter result in a population of massive satellites, which have inexplicably failed to form stars (the 'Too Big To Fail' problem) [7][8][9] . Several solutions have been proposed, including lowering the mass of the Milky Way halo 10 , or the central stellar velocity dispersion of the dSphs through the action of stellar feedback 11 . However, these still require the Milky Way dSphs to inhabit dark matter halos with pre-infall masses greater than B10 9 Msun 4 .
Here we use a new 'disequilibrium' model fitting algorithm to constrain the mass of a dwarf galaxy-Carina-that appears to be in the process of tidal disruption by the Milky Way 12 . Previous simulations of Carina, which assumed identical spatial distributions for the dark matter and stars, suggested that tides have played a role in the evolution of this system 13 . The extensive observed data for Carina, combined with our larger parameter space of non-equilibrium models, allows us to measure the mass of Carina over a far greater radial range than has been possible to date. Most significantly, we are able to 'wind the clock' back to estimate its mass before it fell into the Milky Way, without recourse to comparisons with cosmological simulations.

Results
Simulation procedure. Our method works by simulating the disruption of almost 19,000 N-body Carina models in a static Milky Way potential. To marginalize over the unknown model parameters, the N-body models are wrapped up inside a Markov Chain Monte Carlo (MCMC) pipeline (see Methods section). For the dark matter halo, we allowed both cusped and cored central density profiles. The former are found in cold dark matter simulations, while the latter give a better match to observations 14 . The main advantage of the method over previous studies is that by using full N-body simulations, we can model 'disequilibrium' systems like the tidally disrupting Carina dSph. A demonstration of the performance of our methodology on artificial data is shown in Fig. 2: we recover both the pre-infall and present-day mass of a mock dwarf.
Mass of the Carina dsph. In Fig. 3, as well as Table 1, we show our results for Carina. Figure 3 shows our best fit surface brightness (top) and projected velocity dispersion (bottom) profiles. Table 1 reports the median and 68% confidence intervals for all of our fitted parameters. To facilitate comparison with the predictions of cosmological simulations, we calculate the distribution of M 200 (the mass within the radius where the mean density of the dSph reaches 200 times the critical density of the Behroozi et al. 28 Moster et al. 29 Sawala et al. 16 Sawala et al. 16 Brook et al. 30 Carina  16 , using the halo mass before infall (BI: black) and at the present time (red). Unlike the other lines and data points, these two curves are not based on abundance-matching 16 . The purple points show the abundance-matching results between Local Group dwarfs and a 'constrained' simulation of our local volume 30 . The black and red circles are our pre-infall mass estimates (M 200 ) for Carina for cusped and cored dark matter halos, respectively. They are lower than any of the curves, only marginally consistent with the extrapolation to low luminosities of the relation found in ref. 28. Finally, the black (cusped) and red triangles (cored) are our present-day mass estimates within 1.5 kpc. We calculate the 1s error bars by including 68% of the good models around the median value as given in Table 1.

1,400
Cusped chains for mock dwarf The present-day mass is very well constrained (red histogram)-the red dashed line shows the actual value from the target model. The pre-infall mass (blue histogram) has more uncertainty but is still close to the actual value (solid blue line), despite the large range of masses explored by the chains. The values quoted in the legends of the histograms are the median values that the best models in the MCMC chains found. The likely reason for the larger uncertainty in the pre-infall mass is that the target model in this test was chosen to be on a very eccentric orbit and has undergone stronger tidal disturbance than Carina (see Table 1). In addition, the noise we add to its 'observations' (in order to make the uncertainties similar to those in Carina) admits models with a larger range of orbital eccentricities. Both histograms show models with w 2 o9.
Universe) values at the start of our simulations. We find pre-infall M 200 values of 3.9 Â 10 8 ( À 2.4; þ 3.9) and 3.37 Â 10 8 ( À 2.1; þ 3.8) Msun for Carina models with cusped and cored halos, respectively. Interestingly, this low mass estimate agrees with an earlier study, which found that the halos in the Aquarius cosmological simulations that reproduced the mass of Carina at the present time were those with pre-infall masses of o4 Â 10 8 Msun 4 . After infall, M 200 is no longer a meaningful quantity for a satellite galaxy, and we therefore use the mass within 1.5 kpc (the radial extent of our surface brightness and velocity dispersion data) as our present-day mass estimate. This quantity is very well constrained, with M(ro1.5 kpc) ¼ 7.1 Â 10 7 ( À 3.5; þ 2.8) and 9.7 Â 10 7 ( À 4.8; þ 4.9) Msun for the cusped and cored halos, respectively. We find that both cusped and cored models fit the data very well, suggesting that there is little power in the binned data to distinguish between the two. In general, protected by their higher central densities, cusped halos are able to withstand the stronger tidal forces experienced during closer perigalactic passages.
Our favoured models include both cases in which Carina inhabits a rather low-mass halo, showing significant tidal disruption in its outer parts and cases where it is protected from external tides by a massive halo. These latter models may be favoured by a recent study that found no evidence for tidal tails at a radius of B1 degree from the main body of Carina 15 . We explicitly tested which of our models show visible tidal tails when analysed similarly to that study and found that only the most tidally disrupted are inconsistent with their results. Table 1 presents results from the three sets of model chains, a cusped and a cored chain that use surface brightness, velocity dispersion and velocity gradient and an additional set of model chains with cusped haloes but which excluded the velocity gradient data. It is seen that ignoring the velocity gradient favours models with an even lower pre-infall mass (nevertheless, consistent with our other results within the errors). We expect our results for Carina to be relatively insensitive to the detailed properties of the Milky Way disk as the majority of successful models are not on diskcrossing orbits.

Discussion
The most striking aspect of our results is our upper bound on the mass of Carina both today and pre-infall. In Fig. 1, we compare this upper bound to predictions derived from 'abundance matching' schemes, where satellite luminosity is assumed to depend monotonically on the mass of the dark matter subhalos at infall, as well as to cosmological hydrodynamical simulations. Intriguingly, our data point for Carina lies to the low-mass side of all but one of the extrapolated relations. Even in this case, there is a tension at the B1s level. Our pre-infall mass estimate, being based on N-body simulations constrained by the observed data for the Carina dSph, does not rely on any assumed cosmological model and, unlike previous studies, it is entirely independent of  Figure 3 | A comparison between the best Carina models with cusped and cored dark matter halos and the observed data. The figure shows the surface brightness and projected velocity dispersion profiles (and the velocity gradient given in km s -1 in the legend) from two of our best models. The red data points are the observations with associated 1s errors. The blue and green curves come from our best fit cusped and cored models, respectively, with their 1s Poisson errors calculated in each bin.  ARTICLE any cosmological simulations. The low mass that we find thus provides a new and complementary insight into the mapping between low-mass haloes and low-luminosity galaxies and suggests that a simple monotonic mapping between light and dark in our standard cosmological model fails. One solution is to posit that for halo masses below B10 10 Msun, the physics of galaxy formation leads to a halo occupation function that is effectively stochastic, as suggested by some cosmological hydrodynamic simulations 16,17 .
Our findings can be extended in the near future: dwarf galaxy proper motions from the Gaia satellite will significantly decrease the orbital uncertainties for Carina and the other Local Group dwarfs 18 . Combined with deep photometric observations of dwarf outskirts, we will be able to obtain a pre-infall mass distribution for the whole population. This will address the important question of whether Carina is an outlier in the mass distribution. However, even taken in isolation, our determination of the mass of Carina before interacting with the Milky Way precludes any models that associate the most luminous dSphs with the most massive subhalos.

Methods
Markov chain Monte Carlo. We performed a large suite of N-body simulations that compared the final state of the disrupting dwarf spheroidal galaxies in the tidal field of the Milky Way to a host of observational data for the Carina dwarf galaxy. The simulations were performed within a MCMC framework 19 in order to sample the parameter space effectively and constrain the properties of the dSph. Our method is reminiscent of ref. 20 that uses a genetic algorithm combined with restricted N-body simulations to explore the tidal disruption of NGC 205 around the Andromeda galaxy. However, our method is more general, using full rather than restricted N-body simulations and using an MCMC algorithm that allows us to fully explore parameter degeneracies.
First, a two-component, non-rotating and spherical N-body model of the dSph was built with 2 Â 10 5 particles split equally between the dark matter and stellar components. The starting coordinates and velocities of the dSph were calculated by integrating the orbit of a point mass backwards in a static Milky Way potential using the present-day position and velocity. We then replaced the point mass with the live model and performed a full N-body simulation of its evolution around the Milky Way for 6 Gyr. At the end of the simulation, the radial profiles for the velocity dispersion (s) and the surface brightness (S), as well as the velocity gradient (rv) between the outermost bins along the major axis of the model dwarf were calculated, and a (reduced) w 2 -statistic was used to compare it with Carina data 12,21 . The likelihood ratio between the consecutive models was used to determine the more favourable region of the parameter space at each step, as the Markov Chain accepted either the new model or re-accepted the older one with the better likelihood. This is the learning process of the MCMC algorithm that chooses new initial conditions for the simulations at each step on the basis of the initial conditions of the last model that is accepted. In our pipeline, the proper motions are used as a prior on the basis of observational data and running N-body simulations that fit both the observational profiles and the proper motions, we are able to marginalize over orbit and halo properties simultaneously with the minimum number of assumptions for the dark matter mass and distribution of the satellite.
Supplementary Fig. 1 shows a schematic representation of the algorithm and the codes used at each step. MCMC chains that used models with cusped and cored halos were run separately, where each N-body simulation took B1 h on 32 processors on average. As a very large number of simulations needed to be performed (B19,000 in total for Carina), several chains were run in parallel. Before adding the chains together for the final analysis, the first 50-100 simulations of each of them were dismissed to account for the burn-in period of the algorithm.
Parameter space. The parameter space of initial conditions consisted of six free parameters (m a cos(d), m d, , M s , r s , M h and r h ) for which the allowed ranges are given in Table 1. Given the potential for systematic errors in the observed proper motions, the allowed range for the present-day proper motions (m a cos(d), m d ) used for the orbit integration was the 3s range of the observations (the values of the revised proper motions were obtained through private communication from the authors of ref. 22). The N-body models were generated with falcON 23 according to the split power profiles given in Equation 1, where a, b and g were fixed: 1,4,1 for the cusped halo; 0,4,1 for the cored halo; and 0.515, 4.45 and 0.287 for the stellar component.
The latter was chosen to have a functional form providing a good fit to the present-day surface density in the central regions, albeit with variable amplitude and scale radius that can take values between 0.2 and 2.5 times the original fit (r s ¼ 0.237 kpc). The initial stellar mass M s was allowed to be up to five times larger than the current one.
The initial mass of the dark matter halo M h could be as small as 10 6 Msun, making it only twice as massive as the stellar component, while the upper limit is high enough (10 10 Msun) to allow the large mass models predicted by the cosmological simulations to be tested. Similarly, while the scale radius of the dark matter halo, r h can be as large as 5 kpc, the lower limit is determined by that of the stellar component for each model.
Technical details of the N-body simulations and the Milky Way potential. The external Milky Way potential (pot 4a 24 ) implemented for the point mass orbit integrator 25 and the N-body code PkdGRAV-1 (ref. 26) was provided by the GalPot programme provided in the NEMO Stellar Dynamics Toolbox 27 . The mass of the Milky Way reaches 5 Â 10 11 Msun at 100 kpc, which is Carina's current distance. The simulation time of 6 Gyr was chosen to allow a large enough timescale for the dwarf's evolution while avoiding the complications because of the evolution of the Milky Way potential, and hence Carina's orbit 25 .
We keep the parameters describing the potential of the Milky Way fixed in our analysis. This is an additional systematic error that we do not currently marginalize over. However, much of the associated uncertainty is accounted for by the wide range of proper motions admitted by the large observed errors bars and so marginalizing over uncertainties in the Milky Way potential is not yet required. However, as Carina's proper motion errors shrink, this may become the leading error term. Such a situation would open up the possibility of actually constraining the Milky Way potential alongside the mass and orbit of Carina.
The parameters used in PkdGRAV were chosen on the basis of analysis of simulations with different number of particles, opening angle, time step and the softening length, taking into account the trade-off between CPU time and accuracy. Supplementary Table 1 shows a subset of these simulations where it is seen that the main factor that can affect the w 2 is the number of particles. Therefore, we chose N ¼ 2 Â 10 5 , the largest number of particles that were computationally affordablethe slow down of the simulation when run with 2 Â 10 6 particles instead of 2 Â 10 5 was a factor of 10. We underline that even for the very tidally disrupting model we present in Supplementary Fig. 2, the effects of changing these numerical parameters is not large enough to change the results of the MCMC.
Tests with mock data. To test the MCMC pipeline described above, we ran the chains first for a mock dSph instead of Carina. The mock dwarf had a cusped dark matter halo profile to start with and its final S, s and rv were similar to those of Carina. We modified the S and s profiles, adding Gaussian noise to each data bin as well as increasing the error bars to match the s.d. of the data that resulted in a Mock dwarf with similar uncertainties to Carina. The noise in the velocity gradient in Carina was calculated for B10 stars in the outermost data bin on the major axis. The error bars were obtained by using 100,000 random samples of 10 stars. Both the ideal N-body model and the noisy data are presented in Supplementary Fig. 3, where an increase in the velocity dispersion in the outer bins because of the eccentric orbit of the dwarf is observed.
Three sets of Markov Chains were run separately for the mock dwarf. Two chains where the likelihood comparisons between models were made using (S, s, rv) were run separately testing either cored or cusped dark matter halos, while the third chain with a cusped halo used only S and s in its 'observations'. Supplementary Table 2 shows the parameter values used in the mock data as well as those found by the three chains. The w 2 cuts were chosen to include as many as possible 'good' models for our statistical calculations to be meaningful. The numbers presented in the Supplementary Table 2 are derived from 2,469 (50 unique) models with w 2 (S, s,rv)o9 and 3,119 (335 unique) models with w 2 (S, s)o6 for the cusped chains.
It is seen from Supplementary Table 2 that the current mass M(ro1.5 kpc) is very well constrained by both the cored and the cusped chains within the radius probed by the data. It is also seen that when the velocity gradient was excluded from the w 2 calculations, the chains found models with larger and more massive dark matter halos and hence overestimated the both the initial and the final mass. Interestingly, although the Mock dwarf was cusped to start with, the cored chains recovered the initial mass better than those with cusps that found models with M200 twice as massive. Nevertheless, the results for these chains are very promising considering the wide ranges of masses the algorithm explored. Supplementary  Fig. 4 and Supplementary Fig. 5 show the distribution of the good models in the parameter space that was explored by the cusped (cored) chains, for all of the models as well as those within the w 2 o9 cut. The figure shows that the orbit is the worst constrained property of the models; hence, any improvement of the observed uncertainties in the proper motions will improve the results.
As seen in Fig. 2, the final mass is very well constrained by the chains within the radius probed by the data.
MCMC for Carina. The same three sets of MCMC chains were then run for Carina, which is ideally suited for our study because of the quality of available photometric and spectroscopic data extending to large radii, as well as the evidence of a velocity gradient along the major axis 12 . These new chains found more models with small w 2 than those for the mock dwarf. In order to test the validity of the w 2 cut used above, we re-analysed the Carina chains for different cuts. Supplementary  Figs 6 and 7 show the distribution of the models for a w 2 o6 cut, which included 445 models (22 unique) and the w 2 o9 cut that we use throughout this paper with 2,686 models (232 unique). As can be seen from the figures, these two options result in mass estimates consistent with each other. Therefore, as the overall power of the method we use is to provide a distribution of good models, we choose to use the same w 2 cut for Carina and the mock dwarf. As a final test, in Supplementary  Fig. 8, we compare two models with different w 2 for both cusped and cored halos and demonstrate that even at the highest end of the w 2 range the models used still provide a reasonable match for Carina.
Code availability. The simulation data, MCMC algorithm, pytipsy package necessary to read the binary data files and the codes analysing the simulation data are available on the project website: http://vo.aip.de/dwarfedmasses/.
The new version of PkdGRAV that was used for the N-body simulation can be obtained from https://hpcforge.org/projects/pkdgrav2/.
The falcON code used to generate the two-component N-body models is part of the NEMO package available at https://github.com/Milkyway-at-home/nemo/tree/ master/nemo_cvs/usr/dehnen/falcON.