Dancing in the Dark: Uncertainty in ultra-faint dwarf galaxy predictions from cosmological simulations

The galaxies with one star particles are, by definition, not resolved. In this section we set out to identify a set of halos that can be reliably compared against each other from each run. For every dark matter halo that contains a star particle at $z=0$, we trace back the most massive progenitor halo and identify the time step in which it has the maximum number of dark matter particles. In Figure 2 we show the halo mass and number of star particles inside the halo at that step. Solid lines connect matched galaxies across the MC and H₂ versions of the simulation. For halos above $10^{9}$ M_⊙, halos are well-matched, with a one-to-one correspondence between formed galaxies. Below $10^{9}$ M_⊙, however, it is clear that the MC run forms more galaxies than the H₂ run, and even in the matched cases the H₂ galaxies tend to form fewer stars at these low halo masses.

For halos with maximum halo mass above $10^{8}$ M_⊙ and N${}_{star}\geq 7$, there are MC halos that contain galaxies but with matched counterparts in the H₂ run that are completely dark. We have examined those dark halos in more detail to test whether their lack of star formation is expected for the H₂ model given the resolution. We have approached this in two ways. First, similar to Kuhlen et al. (2013), we compare the surface densities of our gas particles to the critical surface density $\Sigma_{crit}$ at which atomic H converts to molecular H. The metallicity of our non-star forming gas is $Z/Z_{\odot}=10^{-3}$ or lower. At $Z/Z_{\odot}=10^{-3}$, $\Sigma_{crit}=5700$ M_⊙ pc^-2. The non-star forming gas in the dark matched H₂ halos remains at surface densities $<10$ M_⊙ pc^-2. Thus, it would not be capable of forming H₂, and therefore should not form stars in the H₂ run (see also Sternberg et al., 2014).

Second, we can examine the timescale, $t_{H_{2}}$, for atomic H to convert to molecular H following Krumholz (2012). Following their equation 7:

where $t_{ff}$ is the free-fall time, $Z$ is the metallicity of the gas relative to solar, $C$ is a clumping factor equal to 10 in our simulations, and $n_{0}=\left\langle{n_{H}}\right\rangle/1$ cm^-3 where $\left\langle{n_{H}}\right\rangle$ is the mean hydrogen density. For the non-star forming gas in our matched H₂ halos, $\left\langle{n_{H}}\right\rangle$ is about 10 $m_{H}$ cm^-3. The free-fall time depends on $\left\langle{n_{H}}\right\rangle$ as well, and is 13.6 Myr at 10 $m_{H}$ cm^-3. The timescale for H₂ formation, $t_{H_{2}}$, is more than 10 Gyr for our gas with $Z/Z_{\odot}=10^{-3}$. Even at higher densities of $\left\langle{n_{H}}\right\rangle=100$ $m_{H}$ cm^-3, the timescale for H₂ formation is over 1 Gyr.

Thus, significant amounts of H₂ simply cannot form in these halos before reionization at this resolution (we discuss limitations of the resolution in Section 5). We conclude that our dark halos in the H₂ run are behaving as expected. We therefore restrict the following analysis to matched halos with M${}_{halo,max}>10^{8}$ M_⊙ and N${}_{stars}\geq 7$. It should be clear that we do not mean that these halos are “converged” across the two star formation recipes. The star formation in the two runs in the lowest mass halos is dramatically different, with the MC run forming stars while matched halos in the H₂ run remain dark. Even when the H₂ run forms stars there is a discrepancy in stellar mass with the matched halos in the MC run for halos with M${}_{halo,max}\lesssim 10^{9}$ M_⊙. This is exactly the difference we wish to explore in this work.

In Figure 3 we show the resulting stellar mass to halo mass relation for matched halos in both the MC (left panel) and H₂ (right panel) versions of the 40 Thieves. Each galaxy is color coded by its $V$-band magnitude at $z=0$. Filled circles are central (isolated) galaxies, while stars are satellite galaxies. Four of the lines show results from abundance matching in dwarf galaxies (Brook et al., 2014; Garrison-Kimmel et al., 2014; Read et al., 2017; Jethwa et al., 2018), and the fifth line shows the simulation results of Sawala et al. (2015). The stellar masses have been calculated following Munshi et al. (2013), based on photometric colors. To be consistent with abundance matching, the halo masses are the mass in the corresponding dark matter-only run. Note that we plot $z=0$ halo masses, while the abundance matching results use maximum halo mass. Thus, the satellite results should not be compared directly to the lines (though we note that all satellites must have peak halo masses above 10⁸ M_⊙ to be included in the sample shown here).

Again, it is immediately obvious from Figure 3 that there are galaxies residing in halos above M_halo $\sim$10⁹ M_⊙ that are produced in both runs. However, below $\sim$10⁹ M_⊙ the number of galaxies diverges, with twice as many halos containing galaxies in the MC run. We can now see from Figure 3 that many of these low mass halos are satellites (colored stars). There are five times as many satellites in the MC run than in the H₂ run.

All of the satellites in these runs are in the ultra-faint luminosity range. Like observed UFDs (Brown et al., 2012, 2014; Weisz et al., 2014a, b), they tend to form the bulk of their stars early, with their star formation trickling off soon after reionization (see Figure 5, discussed further below). In Figure 4 we compare the phase diagram for gas that forms stars within the first 1 Gyr of both simulations (i.e., before the end of reionization at $z\sim 6$). Although there is some overlap in the star formation temperatures and densities between the two runs, there is a clear offset in the regions where the majority of stars form. In the MC run (grey points), stars tend to form from gas that spans a range of temperatures (up to 10⁴ K) but is near the threshold density (100 $m_{H}$ cm^-3). In the H₂ run (contours and black points), star-forming gas forms from colder, denser gas ($<500$ K and $\gtrsim$ 1000 $m_{H}$ cm^-3). This difference was also discussed in CH12 (see their Figure 13).

The higher densities that stars form from in the H₂ model is tied directly to the subgrid H₂ formation model. At solar metallicities, the two models form star particles at similar densities, but the models diverge as metallicity decreases. As described in CH12, when metals are present the formation of H₂ is dominated by formation on dust grains. In the reionization epoch when metallicities are extremely low, formation may also proceed through gas-phase reactions. However, during the reionizaton epoch both of these processes are inefficient and slow. As a result, gas particles will frequently reach high densities through gravitational collapse prior to forming significant amounts of H₂. At the same time, the low metallicities reduce the amount of dust-shielding, which means that H₂ (and thus star formation) can only persist in high density gas. For most of the H₂ subhalos, significant amounts of H₂ can never be created before reionization removes the gas from the halos. In the MC run, because stars can form in lower density atomic gas, star formation can begin before reionization quenches it. This also explains why the MC run tends to form more stars in general for galaxies in halos with M_halo $<$ 10⁹ M_⊙ (see Figure 2).

Even at $z=0$, the most massive matched dwarfs in this work still show a difference in the effective star formation density threshold due to the impact of metallicity on dust-shielding (CH12; Christensen et al., 2014b). Despite this, dwarf galaxies in halos above $\sim$10⁹ M_⊙ generally still form similar numbers of star particles, due to their ability to self-regulate. This is consistent with the prior results discussed in the Introduction that found that the density threshold had little impact on the resulting SFR. However, a slight excess of stars seems to form in the H₂ model relative to the MC model at these masses (see Figure 2). This excess is likely due to the ability of the gas to shield in the H₂ model, protecting the gas from heating and leading to additional cold gas present in the H₂ simulations. The excess cold gas results in slightly higher stellar masses in the H₂ run in halos $\sim$10¹⁰ M_⊙.

Kuhlen et al. (2012) were the first to show that a model for H₂-based star formation could suppress star formation in dwarf galaxies primarily due to their lower metallicities. However, they found that their equilibrium H₂ model suppressed all star formation in halos below about M_halo $=$ 10¹⁰ M_⊙, while our non-equilibrium H₂ model forms stars in halos down to M_halo $=$ 10^8.5. Reionization played no role in Kuhlen et al. (2012), since their lowest mass halo that was able to form stars was well above the halo mass thought to be impacted by reionization. On the other hand, the changes in stellar mass that we see in halos below $\sim$10⁹ M_⊙ are explicitly due to the fact that these are halos in which reionization can strongly affect the evolution. The MC halos are able to start forming stars at lower densities than their H₂ counterparts, and thus are able to produce more stars before reionization removes their gas.

Early star formation is not the only reason why there are more satellites in the MC run than in the H₂ run: more of the satellites survive. Three of the surviving subhalos in the MC run are completely disrupted in the H₂ run after merging with a parent halo. Two of those halos had managed to form stars in the H₂ run before being fully destroyed. In fact, a close examination of Figure 5 shows that there is a surviving satellite in the MC run with a very extended star formation history (SFH; bolded red dashed line in the left panel). This satellite had a counterpart in the H₂ run with an extended SFH, but the subhalo is completely disrupted and part of the second largest halo in the H₂ run at $z=0$. All three of the surviving subhalos in the MC run are completely merged in the dark matter-only run as well.

Schewtschenko & Macciò (2011) showed that mergers occur later in simulations with baryonic feedback than in matched halos in an identical dark matter-only simulation. They hypothesized that the pressure of the hot halos found in baryonic simulations slows down accretion of subhalos. Our simulated dwarf galaxies do have hot halos of gas around them (Wright et al., 2018), with the mass in hot gas being up to an order of magnitude greater than the HI gas mass. Later infall of subhalos then leads to less tidal mass loss for subhalos simply due to the fact that there is less time for tidal stripping between infall and $z=0$ (Ahmed et al., in prep.). The end result should be that fewer satellites are fully destroyed at $z=0$ for a run with baryonic feedback. ¹¹1This assumes that additional tidal effects from the potential of the central galaxy are negligible. The disk potential is not negligible in Milky Way-mass galaxies (Peñarrubia et al., 2010; Zolotov et al., 2012; Arraki et al., 2014; Brooks & Zolotov, 2014; Garrison-Kimmel et al., 2017b), but is expected to be negligible in dwarf galaxies.

These results suggest that the feedback in the MC run is stronger, delaying mergers of subhalos. Indeed, another indication that feedback is stronger in the MC run is that the overall halo masses tend to be lower (see comparison of maximum M_halo for matched halos in Figure 2). As mentioned previously, feedback reduces the growth of halos (Munshi et al., 2013; Sawala et al., 2013). Despite the fact that the feedback recipe is the same in both the MC and H₂ runs, feedback is injected into lower density gas in the MC run because the MC run rarely reaches the higher densities found in the H₂ run (see, e.g., figure 7 of Christensen et al., 2014b). When comparing the nine most massive isolated dwarfs at $z=0$ (where the stellar masses are similar), the H₂ dwarf galaxies on average have 10% higher dark matter mass and more than twice the mass in hot gas within their virial radii (and more baryons generally) compared to their MC counterparts.²²2Since the hot gas mass of the MC galaxies is smaller than the H₂ runs, the delayed merging of subhalos does not seem to be a direct result of pressure from the hot halo (as was proposed in Schewtschenko & Macciò, 2011), but rather a response to the ejection of the gas itself. Thus, it seems that the MC galaxies have been able to remove more baryonic material from their halos, delaying their growth in dark matter as well.

However, this result is seemingly at odds with previous comparisons of these two models (Christensen et al., 2014a), which found that the H₂ model had more effective feedback and drove slightly more outflows. While the previous work used the Gasoline code, we use ChaNGa with an updated SPH treatment (Wadsley et al., 2017) and quantized feedback. It is not clear if these changes alter the feedback, making the MC run more effective at removing baryons. A full investigation is beyond the scope of this paper, but we find that all evidence points to stronger feedback in the current MC runs than in the H₂ runs.

In summary, these results suggest that feedback can lead to delayed mergers of the UFD satellites. Combined with the ability to form stars at early times in more halos in the MC run, this makes the difference in satellite numbers even more stark at $z=0$ between the two runs.

Consistent with the results presented in Section 3, we find that there is a delay in the onset of star formation in the H₂ run compared to matched halos in the MC run. We highlight a few cases of this in the right panel of Figure 5. In general (as can be seen in the left panel of Figure 5), most of our galaxies begin to form stars before $z=6$, i.e., in the first 1 Gyr after the Big Bang. However, there are two (central) galaxies in the H₂ run that only begin to form stars after 1 Gyr (shown in blue). Because their counterparts in the MC run can form stars without first generating significant amounts of H₂, the MC counterparts all begin star formation before the end of reionization (the counterparts are also shown in red in the right panel of Figure 5 – they are the two MC galaxies shown that begin star formation before 1 Gyr).

However, post-reionization onset of star formation cannot be attributed to only a single star formation model. There are four galaxies (two centrals and two satellites, shown in the right panel of Figure 5) in the MC run that do not begin to form their stars until after 1 Gyr. In each of the four cases, their counterparts in the H₂ run are dark halos that never formed stars. In that case, we might say that the star formation is so delayed in the H₂ run that it does not occur before $z=0$.

In general, all observed galaxies have early star formation, though the error bars on the time of onset can be 1 Gyr or more depending on how far down the main-sequence resolved stars are detected (Brown et al., 2012, 2014; Weisz et al., 2014a, b). It is not clear if our galaxies that delay star formation until after reionization are consistent with observations. However, the dwarf galaxies in these simulated volumes are much further away from a massive galaxy than any observed galaxies with resolved star formation histories that have pushed below the oldest main-sequence turnoff. Our dwarf galaxies are $\sim$5 Mpc away from a Milky Way-mass galaxy. It remains to be seen if environment plays any role in onset of star formation.

In summary, while the H₂ run has consistently later star formation (or none at all by $z=0$) compared to matched galaxies in the MC run, there is no obvious trend in onset time that could be used to discriminate models based on observations. Rather, it is the number of galaxies and stellar mass of those galaxies that discriminates the models (discussed in the next section).

The factor of two difference in the number of faint dwarf galaxies that form between our two models will lead to different faint ends of the stellar mass function (SMF). Additionally, the fact that the MC run forms higher stellar masses for matched galaxies with M${}_{halo,max}<10^{9}$ M_⊙ (see Figure 2) will also alter the SMF.

First, we demonstrate the difference in the SMF that arises due to the different number of low mass galaxies in each model alone. For both the MC and H₂ models, we populate a SMF following the method outlined in Garrison-Kimmel et al. (2017a) in which the scatter, $\sigma$, in the stellar mass at a given halo mass grows as

where $\gamma$ is the rate as which the scatter grows, and M₁ is a characteristic halo mass above which the scatter remains constant. We adopt the scatter model for field galaxies from Garrison-Kimmel et al. (2017a), where M${}_{1}=10^{11.5}$ M_⊙ and $\gamma=-0.25$. We populate the $z=0$ ELVIS catalogs (Garrison-Kimmel et al., 2014) with this SMHM relation, and then assign galaxies as “dark” based on the fraction of luminous to non-luminous halos shown in Figure 1 for each of our models. The resulting SMFs are shown in Figure 6 as the blue (MC) and green (H₂) lines. We do this $1000$ times in order to estimate our errors (dashed lines), and normalize the results so that each has 12 galaxies with stellar mass comparable to Fornax (this is roughly the number of Fornax-mass galaxies within 1 Mpc McConnachie, 2012). Figure 6 demonstrates that there is a difference in slope below M${}_{star}\sim 10^{5}$ M_⊙, with the H₂ simulation predicting the shallower faint-end slope.

The blue and green lines adopt the same slope and scatter of the SMHM relation, and thus only reflect the change in the fraction of “occupied” dark matter halos shown in Figure 1. However, the H₂ model tends to form fewer stars in matched halos at the low mass end, which will yield a steeper SMHM slope that will impact the SMF. As can be seen in Figure 3, the Garrison-Kimmel et al. (2017a) SMHM slope appears to be a decent fit to the central galaxies in the MC run (left panel). However, the Brook et al. (2014) SMHM slope appears to be a better match to the central galaxies in the H₂ run. Our simulations do not have enough galaxies to independently define the slope and scatter of the SMHM relation for each prescription, so we adopt the slope and scatter of Garrison-Kimmel et al. (2017a) to describe the MC run, and Brook et al. (2014) to describe the H₂ run, in order to show the additional affect that the steeper SMHM relation will have on the SMF. We assume the same growing scatter as in the Garrison-Kimmel et al. (2017a) relation but applied to the Brook et al. (2014) slope, and again adopt the galaxy occupation fraction for the H₂ simulation. The result is shown as the orange line in Figure 6. It can now be seen that there is a dramatic difference in the number of expected UFD galaxies in the two models.

The SMF for the two runs is indistinguishable above M${}_{star}\sim 10^{5}$ M_⊙. Since this is the mass range that is currently best constrained, current observations cannot constrain the two models. However, galaxies at lower stellar masses are being discovered by surveys like DES and HSC-SSP, and potentially hundreds will be found near the Milky Way with the Large Synoptic Survey Telescope (LSST; Tollerud et al., 2008; Walsh et al., 2009; Newton et al., 2018), and out to greater distances using integrated light surveys (Danieli et al., 2018). The masses of the UFDs discovered in DES is highlighted by the grey region in Figure 6. It can be seen that pure number counts of faint dwarfs in LSST (i.e., UFDs found out to $\sim$1 Mpc) can help us to constrain how star formation proceeded at high redshift.

As discussed above, the MC simulation contains five times more satellites than the H₂ run at $z=0$. Here we use this result to estimate the predicted frequency of satellites around dwarfs in the Local Group in the two models.

We follow the methodology in Wheeler et al. (2015) and use the ELVIS suite of collisionless zoom-in simulations of Local Group-like environments (Garrison-Kimmel et al., 2014) combined with the results of our two simulations to predict the frequency of a satellite around a dwarf galaxy with M${}_{vir}\sim 10^{10}$ M_⊙. To summarize the procedure: (1) we select isolated dwarf galaxies from the ELVIS suite from both the isolated and paired Milky Way-mass halo simulations and (2) we estimate the frequency of subhalos with peak $M_{vir}\geq 10^{8}$ M_⊙ within 50 kpc of the dwarf host. Figure 7 shows the results.

As expected, the MC model predicts that UFD satellites of dwarf galaxies are more abundant than the H₂ model. The MC run produces non-negligible frequencies for finding at least one satellite around a dwarf host, 46%. The H₂ model predictions are roughly a factor of 4 lower, at 16%. The H₂ model also suppresses galaxy formation to the point that no Local Group dwarf with M${}_{vir}\sim 10^{10}$ M_⊙ is expected to host more than one luminous satellite companion, while multiple UFD satellites are possible in the MC run.

The Wheeler et al. (2015) results adopt the FIRE 1 star formation and feedback prescription (Hopkins et al., 2014). FIRE 1 adopted a threshold for star formation of $>$ 100 $m_{H}$ cm^-3, while the updated FIRE 2 adopts a higher density threshold of $>$ 1000 $m_{H}$ cm^-3 (Hopkins et al., 2018). Both FIRE 1 and FIRE 2 use the equilibrium model from Krumholz & Gnedin (2011) to determine the H₂ fraction in a given particle, which is used to calculate the self-shielding of the gas particle and the cooling rate. The presence of H₂ for star formation is an explicit requirement, which is usually ensured by the simultaneous requirement of a high density threshold. However, because FIRE adopts the equilibrium model of H₂, there is no delay in star formation due to the H₂ formation time.

Figure 7 demonstrates that the predictions in Wheeler et al. (2015) are closer to our MC results than H₂ results. This is in line with our expectations given the similarity in star formation threshold (and resolution) between the FIRE 1 model and the MC model. Of course, both the feedback and reionization models in FIRE are different than those used here, which may also play some role in the results. We note, though, that we have adopted a reionization model (Haardt & Madau, 2012) that is known to lead to earlier reionization than that used in FIRE (Faucher-Giguère et al., 2009), as was shown in Oñorbe et al. (2017). Despite the stronger, earlier heating in our MC model, we predict slightly more satellites than Wheeler et al. (2015).

The results presented in Figure 7 serve to underscore that current predictions for the number of UFD satellites expected around Local Group dwarfs should be approached with caution. This work highlights the range of values we expect to see in state-of-the-art zoom simulations that reach the UFD mass/luminosity range.³³3 Though we note that these results are all currently based off of simulations of classical field dwarfs, and no one has yet simulated UFDs around a Milky Way-mass galaxy in a cosmological context. However, we do not currently know what star formation physics model is “correct.” Until future observations better constrain the models, we must find independent constraints from existing observations.

It is unlikely that either of the explored models, MC or H₂, is “correct.” The MC model does not take into account whether the gas can shield when forming stars, and shielding of the gas is likely to be a necessary ingredient for star formation to occur. Similarly, it has been argued that H₂ is not necessary for star formation, but that its presence is correlated with the ability of the gas to shield (Glover & Clark, 2012; Mac Low & Glover, 2012; Krumholz et al., 2011; Krumholz, 2012; Clark & Glover, 2014). Hence, it is possible that neither model accurately reflects the physics of star formation in the first halos. Instead, a model that links star formation to shielding may be more appropriate (Byrne et al., 2019).

We emphasize that not all H₂ models for star formation will behave in the same way as our adopted non-equilibrium H₂ model. As discussed throughout this paper, some H₂ models assume equilibrium between the formation and destruction of H₂ in the interstellar medium. Because H₂ formation is not explicitly followed, there will be no dependence on the formation time of H₂ in equilibrium models. Thus, equilibrium H₂ models may behave more like the MC model, though this remains to be tested.

Recently, the Dark Energy Survey (DES), has nearly doubled the number of known UFDs in the Milky Way (Bechtol et al., 2015; Drlica-Wagner et al., 2015; Kim & Jerjen, 2015; Kim et al., 2015; Koposov et al., 2015; Luque et al., 2017). Many of the DES dwarfs are thought to potentially be satellites of the Magellanic Cloud system (Deason et al., 2015; Yozin & Bekki, 2015; Jethwa et al., 2016; Sales et al., 2017; Kallivayalil et al., 2018; Li et al., 2018). We note that we do not have an LMC-mass galaxy in this simulation volume. The halo mass of the LMC is estimated to be at least M${}_{halo}>10^{11}$ M_⊙ (Besla, 2015; Dooley et al., 2017a; Peñarrubia et al., 2016; Laporte et al., 2018; Erkal et al., 2018). A lot of work has been done recently to determine if the DES dwarfs are satellites of the LMC or Magellanic System (e.g., Deason et al., 2015; Jethwa et al., 2016; Sales et al., 2017). The majority of studies find that 25-50% of the newly discovered DES dwarfs are likely to have come in with the Magellanic Clouds. Kallivayalil et al. (2018) and Pace & Li (2018) derive proper motions from Gaia data to test whether the UFD satellites have kinematics consistent with falling in with the Clouds. Between the two studies, seven UFDs are thought to be associated. It is not immediately clear if these high numbers are compatible with our predictions, given our lack of LMC-mass galaxies.

The models of Dooley et al. (2017a) predict $\sim$5-15 UFDs with M${}_{star}>3000$ M_⊙ (the lower limit for at least seven star particles in our well-matched galaxy sample) around an LMC-mass galaxy, and $\sim$2-9 around an SMC-mass galaxy. These results are consistent with the number of UFDs that are potentially associated with the Magellanic Cloud system. However, the abundance matching results of Dooley et al. (2017b) put the SMC in a halo with M${}_{halo}>10^{11}$ M_⊙, more massive than our most massive simulated galaxy. Based on their abundance matching results, a more direct comparison of our most massive simulated galaxies is to WLM, for which Dooley et al. (2017b) expect to find roughly one UFD companion. This suggests that our results are generally consistent with the estimates in Dooley et al. (2017b, a), as they should be as long as our simulation results are reasonably described by one of the SMHM relations that they explore. However, it is unlikely that even our MC model would predict seven satellite companions around an LMC-mass galaxy. To reach such high numbers, the satellites must be associated with the whole Magellanic system (SMC+LMC), rather than just the LMC.