The Imprint of Cosmic Reionization on the Luminosity Function of Galaxies

In the ΛCDM model of structure formation, small dark matter halos are already forming profusely during the dark ages—the period following the (re)combination of hydrogen at redshift z ∼ 1100, when the gas becomes neutral. Neutral gas is able to cool into these dark matter halos and form stars and galaxies, bringing the dark ages to an end. UV radiation from the first sources of light reionizes the hydrogen, heating it up to a temperature of $\sim {10}^{4}$ K, raising its entropy, and preventing it from cooling into halos with effective temperature T_vir ≲ 10⁴ K (Doroshkevich et al. 1967; Couchman & Rees 1986). Thus, the reionization process temporarily halts the formation of galaxies in low-mass halos (Rees 1986; Efstathiou 1992; Loeb & Barkana 2001). Galaxy formation resumes some time later when sufficiently massive halos, with virial temperatures well above 10⁴ K, begin to form.

The temporary suppression of galaxy formation as a result of reionization is reflected in the abundance of dwarf galaxies today. No galaxies form below a present-day halo mass of a few times 10⁷M_⊙, and only a fraction of the halos with masses between this value and ∼10¹⁰M_⊙ form a galaxy (Sawala et al. 2013, 2016b; Fitts et al. 2017). In halos that do form a galaxy, the growth of stellar mass is further limited by supernovae feedback (Larson 1974; White & Rees 1978; White & Frenk 1991). One consequence of these processes is that the number of these "chosen few" galaxies is much smaller than the number of dark matter subhalos predicted to be orbiting around the Milky Way in cosmological N-body simulations (e.g., Kauffmann et al. 1993; Bullock et al. 2000; Benson et al. 2002a, 2002b; Somerville 2002; Font et al. 2011; Sawala et al. 2016a), thus readily explaining away the so-called "missing satellites problem" often deemed to afflict the ΛCDM model (Klypin et al. 1999; Moore et al. 1999).

The two critical features of reionization that impact on the abundance of dwarf galaxies are:

• 1.  
  The time when reionization happened.
• 2.  
  The characteristic scale below which gas could no longer cool in dark matter halos.

These two features are linked because the epoch of reionization determines how long low-mass halos are able to continue forming stars before reionization inhibits further gas cooling in them. Understanding these features is therefore crucial to an understanding of galaxy formation.

An important constraint on the epoch of reionization can be derived from the polarization of the cosmic microwave background radiation. The polarization data directly constrain the electron scattering optical depth to recombination, τ, which can be converted to an equivalent redshift of reionization by assuming a model for the redshift evolution of the ionization fraction. Recent estimates from the Planck satellite data imply that the universe was 50% ionized by ${z}_{\mathrm{re}}={8.8}_{-1.4}^{+1.7}$ (Planck Collaboration et al. 2016). Theoretical work suggests that reionization proceeded relatively quickly, with the ionization fraction increasing from 20% to 90% over ∼400 Myr between 9 ≲ z ≲ 6 (Robertson et al. 2015; Sharma et al. 2018). An alternative probe of the ionization state of the IGM comes from the spectra of QSOs: the absence of a Gunn–Peterson trough in the absorption spectra of quasars at z ≲ 6 (e.g., Fan et al. 2000), and its presence in spectra at z ≳ 6 (e.g., Becker et al. 2001; Bolton et al. 2011), suggest that the universe completed the transition from neutral to ionized at around that time. A third source of evidence is the decline in Lyα emission observed from galaxies at z > 6, attributed to absorption by intervening H i gas (Stark & Ellis 2006; Fontana et al. 2010; Caruana et al. 2012, 2014; Treu et al. 2013; Pentericci et al. 2014; Schenker et al. 2014; Tilvi et al. 2014; Mason et al. 2018).

There is still considerable uncertainty regarding the characteristic scale below which galaxies are significantly affected by the photoionizing background. Rees (1986) suggested that halos with circular velocities ∼30 km s⁻¹ would be able to confine photoheated gas in stable equilibrium (i.e., with photoheating balanced by radiative cooling), an idea recently corroborated by Benítez-Llambay et al. (2017) in the Apostle hydrodynamical simulations (Fattahi et al. 2016; Sawala et al. 2016a). Gnedin (2000) expressed the characteristic scale in terms of a filtering mass (corresponding to a circular velocity of ∼50 km s⁻¹) that sets the scale over which baryonic perturbations are smoothed over in linear perturbation theory (see also Thoul & Weinberg 1996, who reached a similar conclusion using 1D hydrodynamical simulations). Okamoto et al. (2008) used high-resolution hydrodynamical simulations to estimate the loss of baryons from low-mass halos resulting from photoionization and revised the filtering mass scale down to ∼25 km s⁻¹ (corresponding to a halo mass of ∼6 × 10⁹M_⊙). Recent radiation-hydrodynamic simulations of reionization give a halo mass of $\sim 2\times {10}^{9}\,{M}_{\odot }$, below which the effects of photoionization become important (Ocvirk et al. 2016).

In this paper we propose a new probe of the physics of reionization: the shape of the (differential) dwarf galaxy luminosity function. We show explicitly that this function encodes both when reionization happened and the characteristic scale below which it had a significant impact. We focus specifically on the luminosity function of satellites of both Milky Way and LMC-mass galaxies because these are easier to measure observationally than the luminosity function of the dwarf galaxy population as a whole. However, all the features of the satellite luminosity functions that we highlight here are also present in the general dwarf galaxy luminosity function. Current observational surveys like the SDSS (Adelman-McCarthy et al. 2007; Alam et al. 2015), the Dark Energy Survey (DES, Bechtol et al. 2015; Drlica-Wagner et al. 2015; Kim et al. 2015; Koposov et al. 2015), and Pan-STARRS1 (Laevens et al. 2015; Chambers et al. 2016) are rapidly improving the census of faint satellites in the Milky Way. The total luminosity function of Milky Way satellites can be readily inferred from a partial census (Koposov et al. 2008; Tollerud et al. 2008; Hargis et al. 2014; Newton et al. 2018). Future surveys like DESI and LSST will measure properties for large samples that may enable estimates of the luminosity function of the dwarf galaxy population as a whole and of satellites of LMC-mass galaxies, which will test the ideas developed in this paper.

This paper is organized as follows. In Section 2, we introduce the theoretical aspects of this work, including a systematic investigation of how reionization shapes the luminosity function of satellites (Section 2.2.2). In Section 3, we combine the satellite populations of the Milky Way and M31 to test if the imprint of reionization can be detected in the observed luminosity function of dwarf galaxies. In Section 4, we present the cumulative and differential luminosity functions of the Milky Way satellites predicted by our models and compare them to the data. We also provide the predictions of our models for the satellite luminosity function of LMC-mass halos (Section 4.2). Finally, our conclusions are summarized in Section 5.

In this section we provide an overview of the N-body simulations used in this work. We describe the semi-analytic model of galaxy formation, Galform, used to populate dark matter halos in the simulation with galaxies. We also explore how reionization shapes the luminosity function of dwarf satellite galaxies.

2.1. The Copernicus Complexio Simulations

2.2. The Galform Semi-analytic Model of Galaxy Formation

2.2.1. The Semi-analytic Philosophy

Semi-analytic models of galaxy formation are very instructive theoretical tools for understanding the physics of galaxy formation. Although, assuming spherical symmetry, semi-analytic models are unable to follow the evolution of gas in galaxies in full generality, as is done in hydrodynamical simulations, they are much cheaper computationally. A great advantage of this is that, in addition to generating a large, statistical sample of galaxies, it is possible to explore rapidly the parameter space describing the physical processes implemented in the model. This makes it straightforward to examine how model predictions are affected by turning on or off particular mechanisms; this is a feature we exploit in Section 2.2.2. For further discussion of the methodology and philosophy behind semi-analytic modeling, we refer the reader to the review by Baugh (2006).

The Durham semi-analytic model of galaxy formation, Galform, was first presented in Cole et al. (1994, 2000); it incorporates the various physical processes thought to be important for galaxy formation, such as the cooling of gas in halos; star formation in galactic disks and central starbursts; metal enrichment of the interstellar medium (ISM); chemical evolution of stellar populations; feedback from stellar winds, supernovae, and active galactic nuclei (AGNs). As the demands on galaxy formation models have increased due to the availability of better observational data, the Galform model has undergone several upgrades. For example, Baugh et al. (2005) introduced a top-heavy IMF in starbursts in order to reproduce the observed number counts of submillimeter galaxies. To explain the exponential tail at the bright end of the galaxy luminosity function, Bower et al. (2006) introduced AGN feedback as a means to suppress star formation in bright galaxies. Motivated by the improved observational understanding of the link between star formation rates and the gas content of galaxies, Lagos et al. (2011) introduced a star formation law that depends on the molecular gas content of the ISM. Galform employs the Maraston (2005) stellar population synthesis model to compute broadband luminosities and magnitudes from the stellar SEDs of galaxies.

In Galform (and, to some extent, in some of the other semi-analytic models currently in use, e.g., Menci et al. 2002; Monaco et al. 2007; Somerville et al. 2008; Guo et al. 2011; Benson 2012; Henriques et al. 2015), the free parameters of the model are set by requiring a good match to a small selection of properties of the local galaxy population, in particular for Galform: (1) the optical and near-IR luminosity functions at z = 0; (2) the H i mass function at z = 0; (3) galaxy morphological fractions at z = 0; (4) the normalization of the black hole–bulge mass relation at z = 0. (A comprehensive list may be found in Section 4.2 of Lacey et al. 2016, L16 hereafter). In this sense, while semi-analytic models do contain free parameters, the degree to which they can be "tuned" is limited by demanding that a small set of observational data be always reproduced.

2.2.2. Reionization in Galform

Reionization in Galform is implemented using a simple two parameter model: to describe the effect of reionization from a global UV background, the cooling of gas in a halo with circular velocity, ${V}_{c}$, is turned off if V_c < ${V}_{\mathrm{cut}}$ at z < ${z}_{\mathrm{cut}}$, where ${V}_{\mathrm{cut}}$ and ${z}_{\mathrm{cut}}$ are input parameters. In this scheme, ${z}_{\mathrm{cut}}$ controls when reionization happens and ${V}_{\mathrm{cut}}$ determines which halos are affected by reionization. While this treatment may appear oversimplified at first, the ${V}_{\mathrm{cut}}$–${z}_{\mathrm{cut}}$ approach is in fact a good approximation to a comprehensive, self-consistent calculation of reionization in Galform performed by Benson et al. (2002b). In fact, it was shown by Font et al. (2011) that the ${V}_{\mathrm{cut}}$–${z}_{\mathrm{cut}}$ method remains a good approximation even when local ionizing sources are included in addition to the global ionizing background. This approach has the added advantage that investigating the effect of changing the small number of parameters controlling reionization on the predicted satellite luminosity function is relatively simple.

We now explore in detail the effects of changing ${z}_{\mathrm{cut}}$ and ${V}_{\mathrm{cut}}$ on the shape of the luminosity function of satellites. In what follows, we will treat the L16 model (in which ${z}_{\mathrm{cut}}$ = 10 and ${V}_{\mathrm{cut}}$ = 30 km s⁻¹) as the "fiducial" model against which all qualitative changes will be compared. All parameters of the galaxy formation model, apart from ${V}_{\mathrm{cut}}$ and ${z}_{\mathrm{cut}}$, are kept fixed at their L16 values. Throughout this paper, satellites are defined as Galform galaxies that are located within r₂₀₀ of their host halo center. Since we follow the merger trees obtained from the coco N-body simulation, the effects of tidal stripping and dynamical friction on infalling subhalos are automatically taken into account. Galform keeps track of "orphan" galaxies—those that have "lost" their subhalo after infall due to limited numerical resolution. The orbits of these galaxies are followed by tracking the most bound particle of the subhalo from the last snapshot in which this particle is associated with a resolved object (Simha & Cole 2017, see also Appendix C in Newton et al. 2018).

Figure 1 shows the effect of changing ${z}_{\mathrm{cut}}$ (the redshift at which reionization happens) on the differential luminosity function. Each of these models assumes ${V}_{\mathrm{cut}}$ = 30 km s⁻¹. The general shape of the curves is similar: the abundance of satellites slowly increases faintward of M_V = −16, peaking at ${M}_{V}=-10$. Fainter than this, all models exhibit a "valley," the location of which depends only weakly on the value of ${z}_{\mathrm{cut}}$. The depth of this valley (and the number of satellites fainter than this) differs significantly as ${z}_{\mathrm{cut}}$ varies. In particular, the earlier the redshift of reionization, the lower the abundance of galaxies fainter than M_V = −5 (M_⋆ ≈ 10⁴ M_⊙). This figure shows that for a fixed value of ${V}_{\mathrm{cut}}$, the location of the peaks of the two populations carved out by the reionization valley is largely independent of the choice of ${z}_{\mathrm{cut}}$.

Zoom In Zoom Out Reset image size

The interpretation of this dependence of the number of faint satellites on ${z}_{\mathrm{cut}}$ is straightforward: when reionization occurs very early (say at z = 12), very few halos with circular velocity exceeding ${V}_{\mathrm{cut}}$ = 30 km s⁻¹ have formed. As a result, cooling is suppressed in a significant fraction of halos, preventing the formation of new satellites (although halos in which gas has already cooled can continue to form stars and become brighter). A later redshift of reionization allows many more faint galaxies to form before cooling is prevented. The parameter ${z}_{\mathrm{cut}}$ therefore affects the amplitude of the differential luminosity function fainter than the location of the "valley." It is also important to note that changing ${z}_{\mathrm{cut}}$ has no effect on the bright end of the luminosity function, as these galaxies primarily assemble much later, long after reionization has ended.

The effect of ${z}_{\mathrm{cut}}$ on the assembly history and abundance of faint galaxies is demonstrated in Figure 2, which shows the fraction of stellar mass assembled before z = 6 by satellites identified at z = 0 in the L16 model (${z}_{\mathrm{cut}}$ = 10), and a variant of L16 where we have set ${z}_{\mathrm{cut}}$ = 6 (to which we refer hereafter as the L16-z6 model). It is evident from this figure that the smallest galaxies are the ones that form earlier (as expected from the hierarchical build-up of structure in CDM). The faintest population of satellites (M_⋆ ≤ 10⁵ M_⊙), on average, assembles $\sim 60 \% \mbox{--}70$% of their stellar mass prior to reionization in the L16-z6 model. By contrast, in the L16 model, in which reionization occurs at z = 10, the faintest satellites, on average, have assembled 100% of their present-day mass by z = 6. These differences explain the larger abundance of faint satellites in the L16-z6 model compared to L16, and by extension, the systematic effect of changing ${z}_{\mathrm{cut}}$ on the amplitude of the luminosity function fainter than the valley.

Zoom In Zoom Out Reset image size

The location of the valley itself is instead primarily controlled by ${V}_{\mathrm{cut}}$, the threshold that determines which halos are affected by reionization. This is demonstrated in Figure 3. Here, we have fixed the value of ${z}_{\mathrm{cut}}$ = 10. The L16 model, by default, assumes ${V}_{\mathrm{cut}}$ = 30 km s⁻¹. The effect of changing ${V}_{\mathrm{cut}}$ is dramatic across the entire range of luminosities. Increasing ${V}_{\mathrm{cut}}$ shifts the location of the valley to brighter luminosities. As gas is prevented from cooling in larger and larger halos, fewer and fewer bright galaxies form.

Zoom In Zoom Out Reset image size

Reducing the value of ${V}_{\mathrm{cut}}$ to 15 km s⁻¹, below the fiducial L16 value (the red line in Figure 3) leaves the abundance of galaxies brighter than M_V = −10 unchanged, but increases the number of faint satellites. Reionization now only affects very small halos, allowing galaxies in halos in the range ${V}_{c}=[15\mbox{--}30]\,\mathrm{km}\,{{\rm{s}}}^{-1}$ to grow in stellar mass and become brighter than in the fiducial case. The bottom of the "valley" shifts to much fainter magnitudes (M_V ≈ 0 in the example in Figure 3).

A summary of the numerical experiments performed in this section is provided in Figure 4, which shows a schematic illustration of the effects of changing ${z}_{\mathrm{cut}}$ and ${V}_{\mathrm{cut}}$ on the shape of the differential luminosity function of satellites. In short, ${z}_{\mathrm{cut}}$ (when reionization takes places) determines the abundance of satellites fainter than the "valley" in the luminosity function, leaving the abundance of bright galaxies unaffected. ${V}_{\mathrm{cut}}$, on the other hand, determines where exactly the "valley" is formed, and can influence both the faint and bright ends of the luminosity function. As the abundance of bright satellites of the Milky Way is well known, the range of allowed values for ${V}_{\mathrm{cut}}$ is better constrained than the value of ${z}_{\mathrm{cut}}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While ${z}_{\mathrm{cut}}$ and ${V}_{\mathrm{cut}}$ are input parameters specific to Galform, they are parameterizations of very general properties of the physics of reionization: the time when reionization happens and the mass scale of halos that are affected by it. In this sense, the effects described in this section are generic, and not specific to Galform or semi-analytic models in general. Indeed, using a formalism similar to that in Galform to calculate the properties of galactic subhalos and a simple prescription for assigning a stellar content to subhalos that crudely models the sort of processes that we have considered here. Koposov et al. (2009) also identified two populations of satellites. However, their models make rather different predictions to ours for the properties of the two populations. In Section 3, we investigate whether the general features described in this section are present in the observed satellite luminosity functions of the Milky Way and M31.

2.2.3. The Lacey et al. (2016) and Hou et al. (2016) Models of Galform

The most recently published version of Galform, presented in Lacey et al. (2016), includes all of the revisions mentioned in Section 2.2.1 in a single, unified model. This model has been shown to reproduce a wide range of observational relations at various redshifts, such as the fraction of early-type galaxies, the Tully–Fisher relation, the far-IR number counts, the evolution of the K-band luminosity function to z ∼ 4, and the far-UV luminosity functions at z ∼ 3–10. In this paper we treat the L16 version of Galform as the fiducial model against which we compare variations of Galform. In the published version, L16 assumes ${z}_{\mathrm{cut}}$ = 10 and ${V}_{\mathrm{cut}}$ = 30 km s⁻¹. This value of ${V}_{\mathrm{cut}}$ is consistent with the hydrodynamical simulations of Okamoto et al. (2008).

A shortcoming of the L16 model is that the choice of ${z}_{\mathrm{cut}}$ = 10 is not self-consistent. The condition for reionization may be defined as the redshift at which ∼6 ionizing photons are produced per hydrogen nucleus.³ Counting the total number of ionizing photons produced by galaxies as a function of redshift in the L16 model implies that the universe is reionized at z ≈ 6, later than the redshift of reionization inferred from Planck data ($z={8.8}_{-1.4}^{+1.7}$, Planck Collaboration et al. 2016). This discrepancy can be traced back to the strong supernova feedback implemented in the L16 model, which was required to reproduce the faint end of the z = 0 galaxy luminosity function. Strong feedback not only suppresses the number of ionizing photons produced at high redshift, but also results in the metallicities of Milky Way satellites being too low when compared with observations (Hou et al. 2016, H16 hereafter).

To remedy these problems, H16 proposed a feedback prescription in which the feedback strength varies not only as a function of halo circular velocity (as in all models of Galform) but also with redshift. In particular, feedback becomes weaker at high redshift (z > 4), resulting in the production of more ionizing photons, and therefore an earlier redshift of reionization. The redshift dependence was chosen to emulate the dynamical supernova feedback model of Lagos et al. (2013), which attempts to capture the relationship between the efficiency of feedback and properties of the ISM, including gas density, metallicity, and molecular gas fractions. Reionization in the H16 model occurs at z = ${z}_{\mathrm{cut}}$ = 7.9; it also sets ${V}_{\mathrm{cut}}=30\,\mathrm{km}\,{{\rm{s}}}^{-1}$. This model also produces an acceptable luminosity function (M_V ≲ −4) and metallicity–luminosity relation for Milky Way satellites. We will present detailed predictions of the L16 and H16 models of Galform in Section 4.

The (relatively) large number of satellites now known to orbit around the Milky Way and M31 invites us to investigate if the features in the luminosity function of satellite galaxies that our models predict (see Section 2.2.2) are present in the data.

A total of 54 satellites around the Milky Way have now been detected in the SDSS and DES. This census is incomplete because the combined sky coverage of SDSS and DES is only ∼47%, and both surveys are subject to detection limits that depend on the satellite luminosity and distance.

Extrapolations based on N-body simulations (Koposov et al. 2008; Tollerud et al. 2008; Hargis et al. 2014) have suggested that the total population count is at least ${70}_{-30}^{+65}$ (at 98% confidence) for satellites brighter than M_V = −2.7. Recently, Newton et al. (2018) applied a new Bayesian method to a sample that includes the newly detected satellites in SDSS DR9 and DES. They estimated a total of ${124}_{-28}^{+40}$ (at 68% confidence) satellites brighter than M_V = 0. The Newton et al. estimate is particularly important for our test because the faint end of the luminosity function is especially sensitive to the physics of reionization.

To test the predictions of Galform against this data set, we assume that the luminosity function of the pre- and post-reionization satellite populations can each be approximated by a Gaussian. Each Gaussian has three free parameters controlling the location of the peak, the width (standard deviation) and height; in total, a general double-Gaussian model has six free parameters. However, our a priori standard theoretical model, which has ${V}_{\mathrm{cut}}=30\,\mathrm{km}\,{{\rm{s}}}^{-1}$, has two fewer degrees of freedom because the locations of the two peaks (at M_V ≈ −10 and ${M}_{V}\approx -3$) are predicted by the model (approximately independently of halo mass and the value of ${z}_{\mathrm{cut}}$; see Figure 1). We refer to this as the constrained double-Gaussian model. Finally, we compare the double-Gaussian models to the simplest possible model of the satellite luminosity function, a power law.

To determine which model is preferred by the data, we apply the Akaike Information Criterion (AIC; Akaike 1974). The AIC gives a measure of the relative quality of different models given the data and is therefore very useful for model selection. It penalizes models with a larger number of parameters; the model with the lowest AIC value is preferred. For two models, A and B, with corresponding AIC values AIC_A and AIC_B, the quantity $\exp \left[\tfrac{1}{2}({\mathrm{AIC}}_{A}-{\mathrm{AIC}}_{B})\right]$ may be interpreted as the relative likelihood of model A over model B. For our analysis we consider a variant of the AIC that corrects for small sample size (Burnham & Anderson 2003).

Figure 5 shows fits of our three models to the Milky Way satellite luminosity function and the associated AIC values. While the goodness-of-fit is best for the general double-Gaussian, the AIC penalizes that model for having six free parameters (compared to only two for the power law). The constrained double-Gaussian, where two parameters are fixed according to predictions of our galaxy formation model, strongly improves the AIC value. In the Milky Way data alone, therefore, there is evidence for the presence of a bimodal population of satellites just as our model predicts.

To maximize the statistical power of the test, we combine the satellite populations of the Milky Way and M31 using the strategy that we now describe. We combine the satellite luminosity function for the Milky Way estimated by Newton et al. with the satellite luminosity function for M31 compiled from the Pan-Andromeda Archaeological Survey (PAndAS; McConnachie et al. 2009; McConnachie 2012; Ibata et al. 2014; Martin et al. 2016). PAndAS has surveyed the region within ∼150 kpc (in projection) of the center of M31, but the census of satellites is by no means complete. Limiting the sample to satellites brighter than ${M}_{V}=-9.13$, Ferrarese et al. (2016) find a total of 19 satellites in M31. Here, we make the assumption that the PAndAS sample of satellites is complete to about M_V = −8 (c.f. McConnachie 2012), which extends the sample of M31 satellites to 23.

To combine the satellite luminosity functions of the Milky Way and M31 we assume that the abundance of satellites scales with the mass of the central galaxy's halo (Wang et al. 2012). We consider four values for the ratio of the masses of the two halos: ${M}_{200}^{{\rm{M}}31}/{M}_{200}^{\mathrm{MW}}=[0.5,1.0,1.5,2.0]$. For a given value of ${M}_{200}^{\mathrm{MW}}$ we can then derive a corresponding value of r₂₀₀ for M31 and extrapolate the PAndAS luminosity function to r₂₀₀ by assuming a radial profile for the distribution of satellites. Following Newton et al., we assume that the radial number density of satellites, n(r), follows an Einasto (1965) profile:

$$\begin{eqnarray}&&\displaystyle \frac{n(\chi )}{\langle n\rangle }=\displaystyle \frac{\alpha {c}_{200}^{\ 3}}{3{\left(\tfrac{\alpha }{2}\right)}^{\tfrac{3}{\alpha }}\gamma \left(\tfrac{3}{\alpha },\tfrac{2}{\alpha }{c}_{200}^{\ \alpha }\right)}\exp \left[-\displaystyle \frac{2}{\alpha }{({c}_{200}\chi )}^{\alpha }\right],\end{eqnarray} \tag{ 1 }$$

where χ = r/r₂₀₀, $\langle n\rangle$ is the mean number density within r₂₀₀, c₂₀₀ = 4.9, α = 0.24, and γ is the lower incomplete Gamma function.

For satellite magnitudes sampled in both galaxies and for a given value of ${M}_{200}^{\mathrm{MW}}$, the combined luminosity function is the average of the Milky Way and M31 estimates extrapolated to r₂₀₀, with the latter rescaled by the ratio of the M31 and Milky Way halo masses. As the estimate for M31 satellites is limited to galaxies brighter than M_V = −8, only the Newton et al. estimate for the Milky contributes to bins fainter than this.

Figure 6 shows fits of our three models to the combined Milky Way and M31 satellite luminosity function for different values of ${M}_{200}^{{\rm{M}}31}/{M}_{200}^{\mathrm{MW}}$. Here, we have fixed ${M}_{200}^{\mathrm{MW}}={10}^{12}{M}_{\odot }$. In each panel we list the corresponding AIC values for each of the best-fitting power-law, double-Gaussian, and constrained double-Gaussian models. It is clear that the constrained double-Gaussian model is preferred by the data. As we have already seen from Figure 5, the smaller number of free parameters in this model results in a significant improvement in its AIC value over the general double-Gaussian.

Zoom In Zoom Out Reset image size

The results of this test for other choices of ${M}_{200}^{\mathrm{MW}}$ are summarized in Table 1. In all cases, the presence of two populations in the observed luminosity function is significantly preferred over a single power law. The value of the AIC for the power law varies very little with either the mass of the Milky Way halo or the ratio of the masses of the M31 and Milky Way halos. The values for the double-Gaussian and constrained double-Gaussian, on the other hand, decrease significantly with ${M}_{200}^{{\rm{M}}31}/{M}_{200}^{\mathrm{MW}}$ and the difference relative to the power law is largest when the M31 halo is assumed to be twice as massive as the Milky Way halo.

Table 1.  Summary of AIC Statistics Comparing the Power Law (PL), Double-Gaussian (DG), and Constrained Double-Gaussian (CDG) Models of the Combined Milky Way and M31 Satellite Luminosity Functions

${M}_{200}^{\mathrm{MW}}$	Ratio $[{M}_{200}^{{\rm{M}}31}/{M}_{200}^{\mathrm{MW}}$]
[M⊙]	0.5	1.0	1.5	2.0
5 × 1011	PL = 4.15	PL = −1.05	PL = −3.07	PL = −4.23
	DG = 7.38	DG = −7.23	DG = −11.45	DG = −16.87
	CDG = 0.30	CDG = −15.97	CDG = −20.43	CDG = −22.26
1 × 1012	PL = 4.66	PL = 0.23	PL = −1.91	PL = −3.19
	DG = 8.30	DG = −3.93	DG = −9.18	DG = −11.66
	CDG=1.61	CDG = −12.35	CDG = −18.06	CDG = −20.65
1.5 × 1012	PL = 5.57	PL = 1.11	PL = −1.09	PL = −2.44
	DG = 9.83	DG = −1.47	DG = −7.33	DG = −10.28
	CDG = 2.98	CDG = −9.64	CDG = −16.07	CDG = −19.22
2 × 1012	PL = 6.27	PL = 1.79	PL = −0.45	PL = −1.84
	DG = 10.98	DG = 0.52	DG = −5.73	DG = −9.04
	CDG = 4.10	CDG = −7.44	CDG = −14.33	CDG = −17.91

Note. Each column corresponds to a different choice for the mass of the M31 halo relative to that of the milky way halo; each row corresponds to a different choice for the mass of the Milky Way halo. The preferred model in each case (highlighted in Bold) is the model with the lowest AIC value.

Download table as:  ASCIITypeset image

4.1. The Luminosity Function of the Milky Way-M31 Satellites

We have seen in Section 3, that both the Milky Way and the combined Milky Way-M31 satellite luminosity functions are best described as the sum of two distinct populations, each characterized by a Gaussian. We now consider how the predictions of specific Galform models compare with the data. For definitiveness, we will take ${M}_{200}^{{\rm{M}}31}/{M}_{200}^{\mathrm{MW}}=1$, but our conclusions are not affected by this choice: the inclusion of the M31 satellites affects only the bright end of the satellite luminosity function, not the faint population that includes most of the satellites.

Figure 7 compares the cumulative satellite luminosity function for a variety of Galform models with the combined Milky Way-M31 estimate (gray data points). In addition to the L16 and H16 models described in Section 2.2.3, we have also included the L16-z6 model, which, we recall, is identical to L16 except that ${z}_{\mathrm{cut}}$ = 6, which is the self-consistent value for the redshift of reionization in the L16 model. It should be noted that according to our approximate method to determine when reionization actually occurs in the model (based on counting the total number of ionizing photons produced per hydrogen nucleus), the L16-z6 model does, indeed, reionize the universe at z = 6, as in the default L16 model. The reason for this is that changing the value of ${z}_{\mathrm{cut}}$ only affects the abundance of galaxies with M_⋆ ≲ 10⁵M_⊙ (see Section 2.2.2 and Figure 2), far below the scale of the dominant sources of ionizing photons at z ≥ 6 (typically ${M}_{\star }\geqslant {10}^{7}{M}_{\odot }$, Hou et al. 2016; Sharma et al. 2016). Therefore, reasonable changes to the value of ${z}_{\mathrm{cut}}$ do not affect the time at which the model satisfies the condition for reionization.

Zoom In Zoom Out Reset image size

Each panel in Figure 7 shows the luminosity function for different ranges of the assumed mass of the Milky Way host halo; the shaded regions mark the 10th–90th percentile spread in the predicted number counts for that mass bin. The dotted and dashed curves in magenta, respectively, represent the mean predictions from the Auriga (Simpson et al. 2018) and Apostle (Fattahi et al. 2016; Sawala et al. 2016a) hydrodynamical simulations in which the host halo mass is ∼10¹²M_⊙. These curves are only plotted down to stellar masses where the simulations are well-resolved. Neither Apostle nor Auriga have sufficient resolution to follow the luminosity function to fainter magnitudes.

It is interesting to note that all three variants of Galform seem to prefer relatively low masses for the Milky Way (and M31) halo (M₂₀₀ = [0.7 − 1.4] × 10¹² M_⊙). While the L16, H16, and L16-z6 models roughly match the cumulative number counts for satellites brighter than M_V = −10, significant differences can be seen at fainter magnitudes.

The L16 model, for example, is in good agreement with the data down to M_V ∼ −5 (M_⋆ ∼ 10⁴ M_⊙), but underpredicts the number of satellites fainter than this magnitude. While this may at first appear to be a consequence of the strong feedback employed in the L16 model, it is, in fact, a result of the choice of ${z}_{\mathrm{cut}}$ = 10. This can be seen by comparing the prediction of L16 to L16-z6, which agrees very well with the data down to the faintest magnitudes. In this model, the strength of supernova feedback as a function of halo mass is identical to that in L16; the only difference is that reionization now occurs later, at z = 6, rather than at z = 10. Since reionization is delayed, gas can now cool in halos with V_c < ${V}_{\mathrm{cut}}$ = 30 km s⁻¹ for a longer period of time, allowing the abundance of the faint galaxy population to build up.

While the total number of satellites brighter than the M_V = 0 predicted by the H16 model is consistent (within ∼2σ) with the observed total number, the overall shape of the predicted luminosity function is not consistent with the data. This difference is clearly seen in Figure 8, which shows the differential satellite luminosity function. This figure shows that the H16 model overpredicts the abundance of galaxies in the range $-9\leqslant {M}_{V}\leqslant -5$ and vastly underpredicts it at fainter magnitudes. This behavior can be attributed to the weaker feedback at high redshift in H16, which allows faint galaxies to build up their stellar mass, shifting their occupancy from fainter to brighter magnitudes in the luminosity function. This results in a shape that is inconsistent with the Milky Way-M31 data. Considering this evidence, the H16 model is ruled out by the Milky Way data. A similar case can be made for the default L16 model, which, as shown in Figure 8, greatly underpredicts the abundance of satellites fainter than M_V = −7.

Zoom In Zoom Out Reset image size

As we have already seen in Figure 7, the L16-z6 model provides an excellent match to the data for the Milky Way. This is demonstrated in greater detail in the lower panels of Figure 8, where we can see that even the shape of the observed differential luminosity function is captured almost perfectly by the L16-z6 model, particularly for the lowest mass bins. We recall that ${z}_{\mathrm{cut}}$ = 6 was not chosen to provide a good match to the luminosity function; it is the self-consistent value for the redshift of reionization appropriate to the L16 model. Although the remarkable level of agreement between this model and the data may well be fortuitous given the noisy data, it is interesting, nevertheless, that the main features of the shape of the observed luminosity function are reproduced by the model. The discovery of new, ultra-faint satellites will help to confirm or exclude the predictions of this model. Importantly, other predictions of the L16-z6 model (such as the field luminosity functions, Tully–Fisher relation, etc.) that the L16 model reproduces are unaffected by this change in ${z}_{\mathrm{cut}}$. This is because the bright galaxies that these observations probe are not sensitive to when exactly reionization happens, as we demonstrated in Section 2.2.2. Finally, it is worth noting that, as Figure 8 shows, the general shape of the differential luminosity function predicted by a given Galform model is independent of the host halo mass.

4.2. Predictions for the Luminosity Function of LMCs

To conclude this section, we provide predictions for the luminosity function of satellites of LMC-mass galaxies, in anticipation of future surveys like the LSST (Ivezic et al. 2008) and WFIRST (Spergel et al. 2015) that may detect the satellites of such systems. The luminosity functions of lower-mass hosts are especially interesting, as the abundance of their faint satellites is particularly sensitive to the details of reionization.

Figure 9 presents the cumulative satellite luminosity function of LMC-mass systems predicted by the L16, H16, and L16-z6 models of Galform. As in Figure 7, each panel shows the luminosity function in different bins of mass for the LMC-mass host. The systematic difference between the three models is consistent with the trends seen for Milky Way-mass hosts, with the L16-z6 model predicting the largest number of satellites brighter than M_V = 0 in each mass bin. For example, the L16 model, on average, predicts fewer satellites in the highest bin of host halo mass (12) than the L16-z6 model does in the lowest mass bin (18). The halo mass for the LMC itself is close to $\sim 2.5\times {10}^{11}\,{M}_{\odot }$ (e.g., Peñarrubia et al. 2016), for which the L16-z6 model predicts 26 ± 10 (68% confidence) satellites brighter than M_V = 0.

Zoom In Zoom Out Reset image size

Differences in the predictions of the Galform models are revealed explicitly in Figure 10, which shows the PDF of the LMC satellite luminosity function (i.e., the probability that a satellite occupies a particular magnitude bin). This figure is simply the differential luminosity function of LMC satellites normalized by the total number of satellites. The PDFs in each of the L16, H16, and L16-z6 models are qualitatively similar, exhibiting a bimodal population in all cases. All three models peak at ${M}_{V}\approx -10$, before displaying the characteristic "reionization valley" we have previously seen in Section 2.2.2. The distributions then peak once more at magnitudes fainter than the location of this valley, whereas the fraction of galaxies in both peaks is comparable in the L16 model, the L16-z6 model predicts ∼10 times as many faint satellites (M_V ≥ −5) as bright ones.

Zoom In Zoom Out Reset image size

Recently, Dooley et al. (2017) made use of abundance matching to infer the total satellite population around LMC-mass hosts. As shown by these authors, the different abundance matching models available in the literature predict very different numbers of satellites, particularly at the faint end of the luminosity function. In their work, one of the abundance matching models tested is the one calibrated by Garrison-Kimmel et al. (2017), and predicts ∼16 satellites more massive than 10³M_⊙, comparable to what is predicted by the L16-z6 model (∼18) for the lowest LMC-mass bin. However, Sawala et al. (2015) have shown, using the Apostle hydrodynamical simulations of the Local Group, that standard abundance matching prescriptions, such as those on which these numbers are based, are invalid for galaxies with stellar masses less than $\sim {10}^{6}{M}_{\odot }$ or halo masses less than ∼3 × 10⁸M_⊙, because only a decreasing fraction of halos below this mass host a visible galaxy.

The luminosity function of dwarf satellites is one of the most informative statistics of the galaxy population that can be measured from local observations. The total number of satellites is sensitive to the physics of reionization, the strength of supernova feedback, the host halo mass, and the nature of the dark matter. In this paper, we have explored the ways in which reionization influences both the amplitude and the shape of the satellite luminosity function.

To obtain a well-resolved sample of Milky Way and LMC-mass halos, we made use of the high-resolution Copernicus Complexio (coco) suite of simulations (Bose et al. 2016; Hellwing et al. 2016). The merger trees of the dark matter halos in coco were populated with galaxies using the Durham semi-analytic model of galaxy formation, Galform (Cole et al. 2000). Galform is a flexible tool that allows the parameter space of galaxy formation models to be explored efficiently. In this model, reionization is characterized by two parameters: ${z}_{\mathrm{cut}}$, which determines the redshift at which reionization is complete, and ${V}_{\mathrm{cut}}$, which controls the mass scale of halos that are affected by reionization. To emulate the net effect of an ionizing background, gas cooling in halos with circular velocity, V_c < ${V}_{\mathrm{cut}}$, is suppressed at redshifts z < ${z}_{\mathrm{cut}}$. Benson et al. (2002b) showed that this simple prescription agrees remarkably well with the results of a detailed, self-consistent model for the coupled evolution of the global properties of the intergalactic medium and the formation of galaxies in the presence of a photoionizing background due to stars and quasars. Thus, while ${z}_{\mathrm{cut}}$ and ${V}_{\mathrm{cut}}$ are parameters specific to Galform, they quantify general features of the effects of reionization on galaxy formation.

In this paper we have considered two recent versions of Galform: the fiducial Lacey et al. (2016) (L16) model, which assumes ${z}_{\mathrm{cut}}$ = 10 and ${V}_{\mathrm{cut}}$ = 30 km s⁻¹, and the Hou et al. (2016) (H16) model, which assumes ${z}_{\mathrm{cut}}$ = 7.9 and ${V}_{\mathrm{cut}}$ = 30 km s⁻¹. The two models differ only in their treatment of supernovae feedback: whereas the strength of feedback in L16 is a function of halo circular velocity only, supernovae feedback in H16 varies as a function of circular velocity and redshift, becoming weaker at z > 4 for the reasons explained in Section 2.2.3. To understand the effects of reionization on the satellite luminosity function, we additionally considered departures from the L16 model, varying ${V}_{\mathrm{cut}}$ and ${z}_{\mathrm{cut}}$ about their fiducial values. Figure 4 illustrates the effects of varying ${z}_{\mathrm{cut}}$ and ${V}_{\mathrm{cut}}$ on the amplitude of the faint end and the overall shape of the satellite luminosity function. The general picture that emerges is as follows:

• 1.  
  The general shape of the differential satellite luminosity function exhibits two peaks: one corresponding to a population of faint galaxies that were mostly assembled before reionization and one corresponding to a population of bright galaxies that were mostly assembled after reionization. These features are generic and do not depend on the details of the Galform model.
• 2.  
  Between these peaks there is a "valley" whose location depends on the mass scale at which reionization affects the cooling of gas in halos (${V}_{\mathrm{cut}}$ in our parameterization).
• 3.  
  The abundance of satellites fainter than the position of the dip is determined by when reionization occurred (${z}_{\mathrm{cut}}$ in our parameterization); earlier reionization inhibits the build-up of a significant population of faint satellites and vice versa.
• 4.  
  The abundance of satellites brighter than the position of the dip is unaffected by the redshift of reionization, as these galaxies typically assemble the bulk of their stellar mass long after reionization.

In principle, the signatures of reionization described in (1)–(4) are measurable. Remarkably, our general prediction that the satellite luminosity function is made up of two distinct components seems to be validated by the observed satellite luminosity function: by combining a recent estimate for the Milky Way (Newton et al. 2018) with an estimate for M31 based on the PAndAS survey (McConnachie et al. 2009), we find that the presence of a bimodal distribution is preferred over a simple power law (Figure 6). Although with somewhat larger uncertainties, we also find that the existence of two populations can be inferred from the Milky Way data alone (Figure 5).

The observed number of Milky Way satellites brighter than ${M}_{V}\approx -8$ is well reproduced in both the L16 and H16 Galform models (Figure 7), but both vastly underpredict the abundance of galaxies fainter than this magnitude (Figure 8). This is because in these models, ${z}_{\mathrm{cut}}$ is large, 10 in L16 and 7.9 in H16. With such large values of ${z}_{\mathrm{cut}}$, gas cooling shuts off too early, preventing the formation of faint galaxies after this time. The large values of ${z}_{\mathrm{cut}}$ adopted in these Galform models were chosen by reference to the value of the redshift of reionization inferred from early Planck data, z_re ∼ 11 (Planck Collaboration et al. 2014). In the latest Planck data analysis, this value has come down to ${z}_{\mathrm{re}}={8.8}_{-1.4}^{+1.7}$.

Perhaps not coincidentally, a variant of the L16 model (L16-z6) where ${z}_{\mathrm{cut}}$ = 6 provides an excellent match to both the cumulative and differential versions of the observed luminosity function down to the faintest magnitudes; in particular, it produces many more satellites fainter than M_V = −7 than L16, with the faintest satellites assembling the bulk of their stellar mass prior to reionization (Figure 2). The choice of ${z}_{\mathrm{cut}}$ = 6 is also appealing, as it is the self-consistent value for the redshift of reionization in the L16 model (see Section 2.2.3); in addition, L16-z6 retains all the successes on large scales of the L16 model. While ${z}_{\mathrm{cut}}$ = 6 may appear "too late" compared to the latest Planck value of z_re, it should be noted that in the Planck analysis z_re is defined as the time when the universe is 50% ionized. By contrast ${z}_{\mathrm{cut}}$ is more readily interpreted as the time when reionization is complete. Given the large quoted uncertainty in Planck's z_re, these two values are compatible. Furthermore, ${z}_{\mathrm{cut}}$ = 6 is also consistent with the inference from the absorption spectra of QSOs that reionization should have been completed by z ∼ 6.

The epoch by which a region of the universe is completely reionized depends on its environment (see, e.g., Alvarez et al. 2009; Font et al. 2011; Dawoodbhoy et al. 2018). Our assumed value of ${z}_{\mathrm{cut}}$ = 6, which results in a good match to the luminosity functions of the Milky Way and M31, could differ for galaxies located in regions of higher or lower overdensities due to the presence of a larger or smaller population of local ionizing sources. While this would affect the number of galaxies predicted to have formed prior to reionization, even relatively large changes to ${z}_{\mathrm{cut}}$ have very little impact on the scale of the transition between the pre- and post-reionization population of satellites (c.f. Figure 1).

Finally, we have predicted the number of satellites expected to be present around galaxies similar in mass to the LMC (Figure 9). The L16-z6 model, which provides the best match to the combined Milky Way and M31 satellite data, predicts 26 ± 10 satellites (68% confidence) brighter than ${M}_{V}=0$. As shown in Figure 10, the majority of the contribution to this population is from galaxies with ${M}_{V}\sim -3$, or M_⋆ ≈ 10³ M_⊙.

With the continuing investment in observational efforts to compile a census of satellites around galaxies other than our own (see, e.g., the recent results from the SAGA survey; Geha et al. 2017), the statistical significance of the features detected in the satellite luminosity function of the Milky Way and M31 may be confirmed. The prospect of also detecting satellites around less massive galaxies, such as the LMC, offers the possibility of a further test of current ideas about some of the most fundamental physical processes involved in galaxy formation.

We thank the anonymous referee for useful suggestions that improved our paper and Marius Cautun and Andrew Cooper for valuable comments on an early draft. We are also grateful to Alan McConnachie for useful comments on the completeness of the PAndAS catalog, and for making it publicly available in the first place. S.B, is supported by Harvard University through the ITC Fellowship. A.J.D. is supported by a Royal Society University Research Fellowship, and she and C.S.F. are supported by the STFC Consolidated Grant for Astronomy at Durham (ST/L00075X/1). This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grants ST/H008519/1 and ST/K00087X/1, STFC DiRAC Operations grant ST/K003267/1, and Durham University. DiRAC is part of the National E-Infrastructure. This research was carried out with the support of the HPC Infrastructure for Grand Challenges of Science and Engineering Project, co-financed by the European Regional Development Fund under the Innovative Economy Operational Programme. The data analyzed in this paper can be made available upon request to the author.