Constraints on the properties of warm dark matter using the satellite galaxies of the Milky Way

The satellite galaxies of the Milky Way (MW) are effective probes of the underlying dark matter (DM) substructure, which is sensitive to the nature of the DM particle. In particular, a class of DM models have a power spectrum cut-off on the mass scale of dwarf galaxies and thus predict only small numbers of substructures below the cut-off mass. This makes the MW satellite system appealing to constrain the DM properties: feasible models must produce enough substructure to host the number of observed Galactic satellites. Here, we compare theoretical predictions of the abundance of DM substructure in thermal relic warm DM (WDM) models with estimates of the total satellite population of the MW. This produces conservative robust lower limits on the allowed mass, $m_\mathrm{th}$, of the thermal relic WDM particle. As the abundance of satellite galaxies depends on the MW halo mass, we marginalize over the corresponding uncertainties and rule out $m_\mathrm{th} \leq 2.02\, \mathrm{keV}$ at 95 per cent confidence independently of assumptions about galaxy formation processes. Modelling some of these - in particular, the effect of reionization, which suppresses the formation of dwarf galaxies - strengthens our constraints on the DM properties and excludes models with $m_\mathrm{th} \leq 3.99\, \mathrm{keV}$ in our fiducial model. We also find that thermal relic models cannot produce enough satellites if the MW halo mass is $M_{200}\leq 0.6\times 10^{12}\, \mathrm{M_\odot}$, which imposes a lower limit on the MW halo mass in CDM. We address several observational and theoretical uncertainties and discuss how improvements in these will strengthen the DM mass constraints.


Introduction
Recent astrophysical observations have provided tentative indirect evidence for a candidate dark matter (DM) particle with mass in the keV range, e.g. [1,2]. Such a particle would be incompatible with the mass range proposed for candidate cold DM (CDM) particles and could have very different clustering properties on small scales [1][2][3][4]. This, together with a lack of any experimental detection of a CDM particle despite considerable advances in particle detector technology [5,6], has motivated a renewed interest in possible alternatives to CDM [7][8][9]. These seek to replicate the success of CDM on large scales and to explain the observed small-scale features of ΛCDM [10] with less reliance on 'baryonic processes'. One family of these alternative DM models posits a 'warm' DM (WDM) particle that would have a much higher thermal velocity than its CDM counterpart at early times in the evolution of the Universe. These 'thermal relics' are formed in equilibrium with the primordial plasma with masses such that they are relativistic at decoupling but non-relativistic by matter-radiation equality [11,12]. Such particles would free-stream out of small-scale primordial density perturbations, preventing their condensation into small haloes and producing a cut-off in the linear matter power spectrum on astrophysically relevant scales. Detecting this suppression of structure relative to CDM predictions would provide a means of discriminating between the prevailing cosmological paradigm and viable WDM models. The goal of this paper is to use visible tracers of the DM substructure to rule out thermal relic WDM models that do not produce enough subhaloes to host the observed number of low-mass satellites of the Milky Way (MW).
Low mass, DM-dominated galaxies provide an excellent probe of the 'small-scale' DM structure [13][14][15][16]. The smallest and faintest of these can be observed best in the environs of the MW; however, the current census of ∼60 satellite galaxies is highly incomplete as extant surveys do not cover the entire sky to sufficient depth and large parts of it are partially or totally obscured by the MW itself [17][18][19]. Simple volume corrections to the observed complement of satellite galaxies have been used already to constrain the viable parameter space of thermal relic WDM models by comparing the number of DM substructures in MW-mass haloes with the number of observed satellites [20,21]. Such approaches make assumptions about the completeness of the surveys, which could lead to a misestimation of the real satellite population. More recent estimates of the satellite galaxy luminosity function that account for the stochasticity of observational data and uncertainties arising from the variability of host haloes at fixed halo mass suggest that the size of the total complement of MW satellites could be several times larger than previously assumed [19,22,23].
This paper improves on previous work and strengthens the methodology used to constrain the properties of candidate WDM particles in several important ways, which we demonstrate using the thermal relic class of WDM models. First, we use one of the most recent estimates of the total satellite population of the MW, which takes advantage of recent observational data to infer a population of 124 +40 −27 satellites brighter than V =0 within 300 kpc of the Sun [22]. This properly accounts for the incompleteness of current surveys; the method used to obtain this estimate has been tested robustly using mock observations. Secondly, our results account for resolution effects in N-body simulations that prevent the identification of DM subhaloes that survive to the present day but fall below the resolution limit of subhalo finders or are destroyed by numerical effects that enhance tidal stripping [24][25][26]. These significant effects have been overlooked in previous studies which, as a result, produce constraints on the viable parameter space of WDM models that are too restrictive. Finally, we incorporate the uncertainty in the total number of satellite galaxies, which has not been included in previous analyses.
We organize this paper as follows. In Section 2, we describe the method to constrain the properties of WDM models by comparing their predictions of the abundance of subhaloes in MWmass haloes with estimates of the total number of MW satellite galaxies from observations. We apply this methodology to thermal relic WDM and present our main results in Section 3. We investigate further the effect of reionization on the constraints that we obtain in Section 4. In Section 5, we discuss the implications of our results and consider some of the limitations of our method; we present concluding remarks in Section 6.

Methodology
Our goal is to use the satellite luminosity function of our Galaxy to constrain the properties of WDM models using a minimal set of assumptions. In DM cosmologies, galaxies of all masses form almost exclusively within DM haloes. The abundance of these can be probed readily with numerical simulations [31] which provide a useful tool to investigate the predictions of different models; we introduce these in Section 2.1. A DM model is viable only if it forms enough subhaloes to host each MW satellite galaxy. To test for this condition we need two ingredients. First, we need an accurate estimate of the MW satellite galaxy luminosity function, which we discuss in Section 2.2. Secondly, Dwarf galaxies can also form during the collision of gas-rich massive galaxies and these are known as 'tidal dwarf galaxies', e.g. [27][28][29][30]. These are low mass and possess negligible DM content; consequently, they are thought to be short-lived. As our Galaxy has not experienced any recent major mergers, the MW is unlikely to contain a significant population of tidal dwarf galaxies. we need a model to predict the number of substructures given the properties of the WDM particle and the mass of the host DM halo, which we describe in Section 2.3.

N-body simulations
We calibrate our predictions for the number of substructures in a WDM model using high-resolution DM-only N-body simulations of cosmological volumes. The Copernicus Complexio (COCO) suite consists of two zoom-in simulations: one of ΛCDM that we refer to as COCO-COLD [32], and the other of 3.3 keV thermal relic WDM, hereafter COCO-WARM [33]. These two versions differ only in the matter power spectra used to perturb the simulation particles in the initial conditions. Both COCO-COLD and COCO-WARM are simulated in periodic cubes of side 70.4 ℎ −1 Mpc using the 3 code that was developed for the Aquarius Project [24]. The high-resolution regions correspond approximately to spherical volumes of radii ∼18 ℎ −1 Mpc that each contain ∼1.3 × 10 10 DM particles of mass, =1.135 × 10 5 ℎ −1 M . Haloes at the edges of these regions can become contaminated with high-mass simulation particles that disrupt their evolution. We identify these contaminated haloes as having a low-resolution DM particle within 3 × 200 of the halo centre at =0. The cleaned catalogues provide large samples of haloes in both cosmological models and both simulations resolve the subhalo mass functions of DM haloes down to masses ∼10 7 M . The cosmological parameters assumed for this suite of simulations are derived from the WMAP seventhyear data release [34] In N-body cosmological simulations the discreteness of the simulation particles can give rise to gravitational instabilities that produce artificial structures. Models such as WDM that impose a cut-off in the primordial matter power spectrum are especially susceptible to these effects [21,35,36]. The instabilities are resolution-dependent and lead to the artificial fragmentation of filaments, giving rise to small 'spurious' haloes that create an upturn at the low-mass end of the WDM halo mass function. Ref. [21] developed a method to identify and prune these objects from the halo merger trees using their mass and particle content. The onset of numerical gravitational instabilities translates into a resolution-dependent mass threshold. Haloes that do not exceed this during their formation and subsequent evolution are likely to be spurious. This coarse requirement is refined further by a second criterion on the particles that compose the halo when its mass is half that of its maximum value, max / 2. In the initial conditions of the simulation, the Lagrangian regions formed by the particles in spurious haloes are highly aspherical. Their shapes are parametrized by half-max = / , where a and c are the major and minor axes of the diagonalized moment of inertia tensor of the DM particles in Lagrangian coordinates. These considerations were applied to the COCO-WARM simulation by ref. [33] who find that almost all spurious haloes can be removed by applying the criteria: max <3.1 × 10 7 ℎ −1 M and half-max <0.165. The details of the calculation of these threshold values can be found in section 2.3 of ref. [33]. Applying such simple criteria means that some genuine haloes can be removed while some spurious haloes remain; however, the numbers of each are extremely small and do not affect our results. Therefore, we follow this prescription to 'clean' the COCO-WARM catalogues of spurious haloes for use throughout the rest of this paper.
The resolution of a simulation also affects the identification of subhaloes in the inner regions of simulated haloes, e.g. [24,25]. Subhaloes that fall below the resolution limit at any time are discarded by some substructure finders, and some other subhaloes are disrupted artificially by numerical effects [26,[37][38][39]. Consequently, these objects do not appear in the subhalo catalogue, even though they may still exist at the present day. We correct for this by identifying such subhaloes in COCO-COLD and COCO-WARM before they are accreted and tracking them to =0, and restoring them to the subhalo catalogues that we use to calibrate our methodology. In Appendix A, we discuss in more detail the procedure we use to recover these objects, and the effect that excluding them has on the halo mass function.

Model-independent radial density profile of the MW satellites
To obtain the best constraints on the WDM particle mass we need a complete census of the Galactic satellites. The satellite population is dominated by ultra-and hyperfaint galaxies with absolute magnitudes fainter than V = −8 (see e.g. [19,22,40]), which can be detected only in deep surveys. This means that large areas of the sky remain unexplored and that currently we have only a partial census of the MW satellites. However, there are several methods that use the current observations to infer the total satellite count of our Galaxy (see e.g. [17,40]). Here, we use the estimates from ref. [22] that are based on a Bayesian formalism that has been tested robustly using mock observations. These results were obtained by combining the observations of the Sloan Digital Sky Survey (SDSS) [41] and the Dark Energy Survey (DES) [42,43], which together cover nearly half the sky, and estimating the MW satellite luminosity function down to a magnitude, V =0. This roughly corresponds to galaxies with stellar mass higher than 10 2 M [44].
The method of ref. [22] (code implementing this is available from [45]) takes two input components. First, it uses the sky coverage of a given survey and the distance from the Sun within which a satellite galaxy of a given magnitude can be detected. This depends on the depth of the survey and the satellite detection algorithm. Secondly, the ref. [22] method requires the radial probability distribution function of satellite galaxies. Simulations of DM-only CDM haloes show that subhaloes selected by peak , the highest maximum circular velocity achieved in their evolutionary histories, have the same radial number density profile as that of the observed satellites (see ref. [22], and discussion therein). Furthermore, CDM simulations (e.g. [24,32]) have shown that the radial distribution of satellites is largely independent of their mass as well as of the host mass when expressed in terms of the rescaled distance, / 200 , where and 200 denote the radial distance and the host halo radius, respectively. This is studied further in Figure 1, where we compare the normalized radial number density profiles of stacked populations of subhaloes in the COCO-WARM and COCO-COLD simulations. The fiducial populations were obtained by selecting subhaloes with peak ≥ 20 km s −1 and identifying and including subhaloes that would exist at =0 if they had not been prematurely destroyed or missed by substructure finders (for details see Appendix A). We apply this correction after pruning the spurious haloes from the merger trees (see Section 2.1) to ensure that they are not inadvertently restored. Figure 1 illustrates that both CDM and WDM predict the same radial distribution of satellites, which means that we can use the satellite distribution inferred from CDM to make predictions for WDM models. This is beneficial as CDM simulations sample better the inner radial profile, to which the ref. [22] result is particularly sensitive.
To summarize, in this paper we infer the satellite galaxy luminosity function of the MW within 200 for assumed host halo masses in the range, 200 = [0.5, 2.0] × 10 12 M , using the Bayesian methodology presented in ref. [22]. As we mentioned above, this requires two components: 1. a tracer population of DM subhaloes with a radial profile that matches that of the observed satellites; and, 2. a set of satellite galaxies detected in surveys for which the completeness is characterized well.
For the former, we use the same peak -selected peak ≥10 km s −1 fiducial CDM subhalo populations as used in ref. [22]. These are obtained from five high-resolution ΛCDM DM-only N-body simulations of isolated MW-like host haloes from the Aquarius suite of simulations [24]. For the latter, we use the observations of nearby dwarf galaxies from the SDSS and DES supplied in appendix A of ref. [22] (compiled from [43,[46][47][48][49][50][51][52][53][54]). Later work to infer the luminosity function using more recent observational data and a better characterization of the DES completeness function is in good agreement with the ref. [22] results [23,55].

Estimating the amount of halo substructure
Estimates of the average number of subhaloes in MW-like DM haloes can be obtained using the Extended Press-Schechter (EPS) formalism [56][57][58][59][60]. In this approach, the linear matter density field is filtered with a window function to identify regions that are sufficiently dense to collapse to form virialized DM haloes. In CDM models the filter employed takes the form of a top-hat in real space. However, applying this to models such as WDM in which power is suppressed at small scales leads to an over-prediction of the number of low-mass haloes [61]. This occurs because the variance of the smoothed density field on small scales becomes independent of the shape of the linear matter power spectrum if the latter decreases faster than −3 . Consequently, the halo mass function continues to increase at small masses rather than turning over [62, 63 section 3.1], making the top-hat filter an inappropriate choice. Using a sharp k-space filter seemed to address this by accounting for the shape of damped power spectra at all radii [61,64]; however, subsequent work by ref. [63] demonstrates that this over-suppresses the production of small haloes. They find that using a smoothed version of the sharp k-space filter produces halo mass functions in best agreement with N-body simulations. Throughout this paper, we use the ref. [63] smooth k-space filter for the WDM models that we consider. To obtain our estimates of the number of substructures, sub , within 200 of MW-like haloes we follow the approach described by ref. [65] that was subsequently modified in section 4.4 of ref. [66] for use with sharp k-space filters. Using the ref. [63] filter, a conditional halo mass function, SK , is generated from the primordial linear matter power spectrum. Ref. [12] showed that WDM power spectra, WDM ( ) , are related to the CDM power spectrum, CDM ( ) , by WDM ( ) = 2 ( ) CDM ( ) , where ( ) is the transfer function given by (2.1) Here, = 1.12 and is described by ref. [67] as being a function of the WDM particle mass, th , given by Ref. [66] showed that integrating the conditional halo mass function over the redshift-dependent spherical collapse threshold of a given progenitor, ( ), gives the subhalo mass function where is the filter mass and norm is a normalization constant. The latter term, which is a free parameter, corrects the total count for progenitor subhaloes that exist at multiple redshifts which are counted more than once. Using the ref. [63] filter introduces two other free parameters,ˆandˆ, that control the 'smoothness' and the mass-radius relationship of the filter function.
We calibrate the free parameters of the EPS formalism by comparing its predictions of DM substructure with the fiducial subhalo populations of COCO haloes in the mass bin 1. the EPS estimate of the mean number of CDM subhaloes with mass ≥ 10 9 M must equal the mean number of objects with peak ≥ 10 9 M in COCO-COLD haloes; and, 2. the EPS prediction of the mean number of WDM subhaloes with ≥ 10 6 M must equal the mean number of objects with peak ≥ 10 6 M in COCO-WARM haloes (i.e. all subhaloes).
Here, peak is determined using the definition of halo mass [68,69] and represents the highest mass achieved by the subhaloes at any time during their evolutionary histories. Typically, haloes achieve peak just before infall into a more massive halo. In the second calibration criterion, we compare the mass functions at 10 6 M as this is below the turnover in the WDM power spectrum used in COCO-WARM. We obtain excellent agreement between the mean EPS estimates and the COCO simulation results by setting norm = 2.59,ˆ= 4.6, andˆ= 3.9. This is shown in Figure 2, which is discussed below.
The EPS formalism predicts only the mean number of subhaloes in DM haloes of a given mass, and not the host-to-host scatter in the subhalo count. As we will discuss later, including this scatter is very important to obtain unbiased results and thus needs to be accounted for. We do this using the results of cosmological N-body simulations that have shown that the scatter in the subhalo mass function is modelled well by a negative binomial distribution [70,71]. This takes the form where N is the number of subhaloes and Γ( ) = ( − 1)! is the Gamma function. The variable, = / σ 2 , where and σ 2 are, respectively, the mean and the dispersion of the distribution. This scatter in the subhalo count can be described best as the convolution of a Poisson distribution with a second distribution that describes the additional intrinsic variability of the subhalo count within haloes of fixed mass, such that σ 2 = σ 2 Poisson + σ 2 . The parameter r then describes the relative contribution of each of these two terms: = σ 2 Poisson / σ 2 I . We find that the scatter in the subhalo count of haloes in the COCO suite is modelled well by σ =0. 12 , as depicted in Figure 2. We use this approach to characterize the scatter associated with the EPS predictions throughout the remainder of this paper. In both cases, unfilled symbols represent objects from a subhalo catalogue where the 'prematurely destroyed' subhaloes have not been recovered, and filled symbols indicate the same haloes using the subhalo catalogue after restoration of the 'prematurely destroyed' subhaloes.
In Figure 2, we compare the EPS predictions for haloes in the mass range [0.5, 2.0] × 10 12 M to the number of subhaloes in individual COCO haloes of the same mass. We obtain excellent agreement with the N-body results across the entire halo mass range of interest for this study. In particular, our approach reproduces very well both the mean number of subhaloes and its halo-to-halo scatter, which are represented by the grey shaded region and the vertical error bars, respectively.

Calculating model acceptance probability
We rule out sections of the viable thermal relic WDM parameter space by calculating the fraction, v , of WDM systems that have at least as many subhaloes as the total number of MW satellites. We denote with EPS the probability density function of the number of DM subhaloes predicted by the EPS formalism. Then, the fraction of haloes with MW sat or more subhaloes is given by However, as we discussed in Section 2.2, the inferred total number of MW satellite galaxies is affected by uncertainties. We can account for these by marginalizing over the distribution of MW satellite counts, MW MW sat . Combining everything, we find that the fraction of WDM haloes with at least as many subhaloes as the MW satellite count is given by (2.6) Figure 3. The fraction, v , of WDM systems with at least as many DM subhaloes, sub , as the inferred total number of MW satellites, sat , for a DM halo with 200 = 1 × 10 12 M . Thermal relic masses for which v ≤ 0.05 are ruled out with 95 per cent confidence. Earlier works that do not account for the uncertainty in sat or the scatter in sub at fixed halo mass (thin lines) artificially exclude too many thermal relic particle mass values. In this work (thick line) we include both sources of uncertainty in our calculation. The horizontal dotted line indicates the 5 per cent rejection threshold that we use to rule out parts of the WDM parameter space.
While not explicitly stated, both the number of MW satellites and the number of subhaloes (e.g. see Figure 2) depend on the assumed MW halo mass [72], which is still uncertain at the 20 per cent level (e.g. [73]). This means that the fraction of valid WDM haloes depends strongly on the assumed mass of the Galactic halo. Note that the inferred total number of MW satellites depends weakly on the MW halo mass when calculated within a fixed physical distance, e.g. within 300 kpc from the Galactic Centre (see fig. 10 in ref. [22]); however, here we calculate the expected number of satellites within 200 for each MW halo mass.
This approach to calculating the fraction of viable WDM systems for the first time incorporates the scatter in sub at fixed halo mass and the uncertainty in the inferred total MW satellite population. This is important, as excluding one, or both, of these sources of uncertainty produces constraints on th that are too strict. We demonstrate this in Figure 3 where, for each WDM particle mass, we plot the fraction of haloes with mass 200 =10 12 M that contain enough DM substructure to host the inferred population of MW satellite galaxies. We derive our constraints on th from the intersection of these cumulative distributions with the 5 per cent rejection threshold indicated by the horizontal dotted line. In this example, neglecting both sources of uncertainty excludes thermal relic DM with particle masses th ≤2.4 keV, which is ∼15 per cent more restrictive than our reported value of th 2.1 keV (thickest solid line). Some previous analyses (e.g. [21,74]) account for some of the uncertainty by modelling the scatter in the number of DM subhaloes at fixed halo mass. This weakens the constraint; however, the results are still artificially ∼5 per cent more stringent than they should be with our more complete treatment of the uncertainties. In addition to these complications, earlier works also suffer from incompleteness in the =0 subhalo catalogues due to numerical resolution effects. This contributes to a much more significant overestimation of the constraints and we discuss this in detail in the next section.

Constraints on the thermal relic mass
Here we present the results of our analysis obtained using the EPS formalism calibrated to fiducial subhalo populations from the COCO-COLD and COCO-WARM simulations. Our most robust result assumes that all DM subhaloes that form host a galaxy, thereby making no assumptions at all about galaxy formation processes.

Thermal relic particle mass constraints
We calculate the model acceptance distributions of DM haloes in the mass range 200 =[0.5, 2.0] × 10 12 M for several thermal relic WDM models. We rule out with 95 per cent confidence all combinations of 200 and th with v ≤ 0.05. Problems arising from resolution effects persist even when using high-resolution simulations, and these effects are amplified as the resolution decreases. In addition to incorporating the scatter in sub and the uncertainty in sat , we account for resolution effects in the N-body simulations with which we calibrate the EPS formalism by including subhaloes that have been lost below the resolution limit at higher redshifts or destroyed artificially by tracking the most bound particle of these objects to =0 (for details see Appendix A).
The results that we obtain using this approach are displayed in Figure 4. The shaded region represents the parameter combinations that we rule out with 95 per cent confidence. Independently of MW halo mass, we find that all thermal relic models with particle mass th ≤1.80 keV are inconsistent with observations of the MW satellite population. The exact constraints vary with the MW halo mass, such that for lower halo masses we exclude heavier DM particle masses.
Recent studies, especially using Gaia mission data, have provided more precise measurements of the MW halo mass (for a recent review, see [73]). We can take advantage of these results to marginalise over the uncertainties in the MW halo mass. For this, we use the ref. [75] estimate of the MW mass, which we illustrate in Figure 4 with two horizontal dashed lines indicating their 68 per cent confidence interval. This estimate is in good agreement with other MW mass measurements, such as estimates based on the rotation curve or on stellar halo dynamics [73,76]. Marginalising over the MW halo mass, we rule out all models with th ≤2.02 keV. These constraints do not depend on uncertain galaxy formation physics and therefore they are the most robust constraints to be placed on the thermal relic particle mass to date. A more realistic treatment of galaxy formation processesthe effect of which would be to render a large number of low-mass subhaloes invisible -would allow us to rule out more of this parameter space as fewer WDM models would produce a sufficient number of satellites to be consistent with the inferred total population. We consider this possibility in more detail in Section 4.1.
In Figure 4, we include for comparison the constraints obtained by refs. [21] and [74] who use similar analysis techniques. These constraints suffer from resolution effects that suppress the identification of some substructures that survive to the present day. The dotted line demarcates the exclusion region that we would obtain in our analysis if we did not account for these prematurely destroyed subhaloes. Such issues are not revealed by numerical convergence tests that are typically used to assess the reliability of particular simulations. For example, even the 'level 2' simulations of Aquarius haloes, which are some of the highest resolution DM-only haloes available, are not fully converged. We explore this in more detail in Appendix A.

The effects of galaxy formation processes 4.1 Modelling galaxy formation
In the preceding sections we described an approach that gives a highly robust, albeit conservative, lower limit on the allowed mass of the WDM thermal relic particle. This ignores the effects of galaxy formation processes on the satellite complement of the MW. These mechanisms play an important role in the evolution of the satellite galaxy luminosity function but still are not fully understood. Semi-analytic models of galaxy formation enable the fast and efficient exploration of the parameter space of such processes and thus help us to understand how they affect the WDM constraints. [77,78] is one of the most advanced semi-analytic models that is currently available and is tuned to reproduce a selection of properties of the local galaxy population. A complete summary of the observational constraints used to calibrate the model parameters is provided in section 4.2 of ref. [79]; hereafter L16. Of particular interest to our study is the reionization of the Universe, which is the main process that affects the evolution of the faint end of the galaxy luminosity function. The UV radiation that permeates the Universe (and that is responsible for reionization) heats the intergalactic medium and prevents it from cooling into low-mass haloes, impeding the replenishment of the cold gas reservoir from which stars would form.
In , the effect of reionization on haloes is modelled using two parameters: a circular velocity cooling threshold, cut , and the redshift of reionization, reion . The intergalactic medium is taken to be fully ionized at a redshift, = reion , whereafter the cooling of gas into haloes with circular velocities, vir < cut , is prevented. This simple scheme has been verified against more sophisticated calculations of reionization, with which it has been shown to produce a good agreement [80,81]. Recent studies by e.g. ref. [82] have characterized the sensitivity of the satellite galaxy luminosity function to changes in these parameters: a later epoch of reionization allows more faint satellites to form, and a smaller circular velocity cooling threshold permits those faint satellites to become brighter.  [94], and 3.8 keV [95] envelopes of the most robust constraints on the thermal relic particle mass obtained from modelling of the Ly α forest. The shaded region shows the 68 per cent confidence interval on the mass of the MW halo from ref. [75].
We use to explore the effect of different parametrizations of reionization on the number of substructures containing a luminous component around the MW. Several previous works that have adopted a similar approach [20,83] used the L16 model, which has reion =10 and cut =30 km s −1 ; however, this combination of parameters is now disfavoured by more recent theoretical calculations and the analysis of recent observational data, e.g. [84,85]. Additionally, others have noted that using reion =10 is not self-consistent and that a modified L16 model with reion =6 is a more appropriate choice [82]. In light of these theoretical and observational developments, for this study we consider parametrizations of reionization in the ranges 6 ≤ reion ≤ 8 and 25 km s −1 ≤ cut ≤ 35 km s −1 (see [81,[84][85][86][87][88][89]).

Constraints using GALFORM models
Our exploration of different prescriptions for reionization assumes the L16 model as a reasonable description of various feedback and evolutionary processes in galaxy formation. We vary the reionization parameters in the ranges described in Section 4.1 and apply to Monte Carlo merger trees calibrated as closely as possible to the COCO suite. The Monte Carlo algorithm used in cannot be calibrated to match exactly the N-body results as it lacks sufficient free parameters to match both the high-and low-mass ends of galaxy formation. Where a discrepancy exists between the Monte Carlo and N-body luminosity functions, we remap the V values of Monte Carlo satellite galaxies to new values such that the resulting luminosity function is consistent with the N-body results. Using these, we obtain predictions for the dwarf galaxy luminosity function for 500 realizations of each MW halo mass, allowing us to compute the model acceptance distributions in the same manner as before (see Section 2.4). Details of the merger tree algorithm and the functions to remap the Monte Carlo satellite galaxy −band magnitudes are provided in Appendix B. here, the cooling threshold is fixed at cut =30 km s −1 and the dotted, solid, and dashed lines represent constraints obtained assuming reion =6, 7, and 8, respectively. High values of reion produce more stringent constraints on the thermal relic mass at fixed MW halo mass. Right panel: here, reionization is assumed to have ceased by reion =7, and the dotted, solid, and dashed lines represent the constraints obtained assuming cooling thresholds of cut =25, 30, and 35 km s −1 , respectively. Higher cooling thresholds produce more stringent constraints on the thermal relic mass.
In Figure 5, we plot our constraints on thermal relic WDM models assuming a fiducial model of reionization with reion =7 and cut =30 km s −1 . This is a viable parametrization that is consistent with observations and resides in the centre of the parameter ranges that we explore. In this model, we rule out all thermal relic WDM particle masses with th ≤2.95 keV independently of the MW halo mass. When marginalising over the uncertainties in the estimate of the MW halo mass from ref. [75], our constraints strengthen and we exclude with 95 percent confidence all models with th ≤3.99 keV. Our fiducial constraints are considerably stronger than our model-independent result and produce more stringent constraints in different MW halo mass regimes compared with work by refs. [20,83], who also model the effects of galaxy formation processes. More recently, ref. [90] carried out a similar analysis and obtained tighter constraints on the WDM particle mass than we find in this work. We discuss the reasons behind this and its implications in Section 5. In Figure 5, we have also included for comparison the most conservative constraints derived from the Ly α forest by refs. [91][92][93][94][95], which our results complement.
In Figure 6, we explore the effect on the constraints of varying cut or reion while holding the other parameter constant. The left panel shows the effect of varying the redshift at which reionization concludes while fixing cut =30 km s −1 . An epoch of reionization that finishes later, characterized by a lower value of reion , allows more faint galaxies to form in low-mass DM haloes, which weakens the constraints that can be placed on thermal relic WDM models. The right panel shows the effect of curtailing further star formation in low-mass haloes after reionization finishes at reion =7. As the cut cooling threshold increases, a larger fraction of the low-mass galaxy population is prevented from accreting new cold gas from the intergalactic medium after the end of reionization. Consequently, the reservoir of cold gas available for further star formation in these galaxies depletes over time, limiting how bright these objects become by =0. When the cooling threshold is large, fewer faint galaxies evolve to become brighter than V =0 and populate the MW satellite galaxy luminosity function, leading to stronger constraints on the thermal relic mass. For completeness, in Appendix C we provide the constraints obtained for the three values of cut assuming two scenarios with reion =6 and 8, respectively.

Discussion
We have placed new conservative and highly robust constraints on the mass of the thermal relic WDM particle by comparing EPS predictions of the DM subhalo content of WDM haloes with the total number of MW satellite galaxies inferred from observations. We obtain estimates of the total satellite complement using the ref. [22] approach including recent observations of satellites from the SDSS and DES. To calibrate the EPS formalism, we use DM haloes from the COCO simulation suite with masses in the likely MW halo mass range 200 =[0.5, 2.0] × 10 12 M . We improve upon previous constraints by incorporating for the first time the uncertainty in the size of the total MW satellite population and by accounting for unresolved or numerically disrupted subhaloes in N-body simulations (see Appendix A). In a separate analysis we also explore the effect of various assumptions about galaxy formation processes on the constraints that we can place on the WDM particle mass.
We find that, when marginalizing over uncertainties in estimates of the MW halo mass, thermal relic models with th ≤2.02 keV are ruled out with 95 per cent confidence (see Figure 4). This result is independent of assumptions about galaxy formation physics, as for our purposes we treat all DM subhaloes as hosts of visible galaxies. This ensures that the constraints provide a robust lower limit on the mass of the thermal relic WDM particle, improving on the results reported in ref. [21] across the entire MW halo mass range considered (see Figure 4). Our results are competitive with but slightly less restrictive than the constraints obtained by ref. [74] because we account for subhaloes that exist but are missing for numerical reasons from the =0 halo catalogues.
The resolution of a simulation can affect the population of haloes at =0 in two major ways. First, haloes close to the resolution limit of a simulation experience stronger tidal disruption due to numerical effects that can destroy the halo. Secondly, some structure finders stop tracking haloes that fall below a mass threshold at any time during their evolution. Haloes composed of few particles can occasionally fall below this and recover later, with the result that the object is permanently excluded from the final catalogue even if it survives to the present day. Omitting these objects significantly affects the constraints on the WDM parameter space, strengthening them artificially (see Figure 4). This effect worsens as simulation resolution decreases, so constraints that are obtained using lowerresolution simulations and methods that do not account for 'prematurely destroyed' subhaloes will be significant overestimates.
The processes responsible for the formation of galaxies are complex and are yet to be understood fully; nevertheless, they play an important role in shaping the luminosity function of the dwarf galaxies of the MW. Incorporating the effects of these mechanisms into our approach allows us to refine the constraints on the properties of the DM and rule out many more WDM models. In a modified version of the L16 model with reion =7 and cut =30 km s −1 (our fiducial model) we rule out, with 95 per cent confidence, thermal relic models with th ≤3.99 keV when marginalizing over uncertainties in the MW halo mass (see Figure 5). Furthermore, we rule out all thermal relic WDM particle masses with th ≤2.95 keV independently of MW halo mass. These improve on our modelindependent results and are consistent with the constraints obtained in previous works that adopted similar approaches. This result also compares favourably with complementary constraints derived from the Ly α forest by refs. [91][92][93][94][95].
Recently, ref. [90] conducted a similar analysis to constrain the particle mass of thermal relic WDM using the inferred luminosity function of MW satellite galaxies from ref. [55]. Their constraints on the DM particle mass are stricter than all of our results spanning the parametrizations of reionization considered in this work (see Figures 6 and 10, and Table 1). Two factors contribute to this discrepancy. First, ref. [90] use an abundance matching technique extrapolated to very faint magnitudes to populate substructure with galaxies. Such techniques adopt a model to describe the relationship between the DM structure and the luminous component; however, they may not capture the full complexity of galaxy formation physics at the faint end [14]. Semi-analytic models like the one used in this work are physically motivated and fare better at modelling the baryonic processes taking place on small scales, encapsulating more of the complexities of galaxy formation in this regime; however, they are not entirely free of simplifying assumptions. Secondly, the ref. [90] results are based on a combination of Pan-STARRS and DES data whereas our analysis uses satellite galaxy data from SDSS and DES. The Pan-STARRS survey data are not as deep as those from SDSS, particularly at the faint end of the satellite galaxy luminosity function. Consequently, the size of the satellite population inferred from the Pan-STARRS data, and hence the WDM constraint derived from this, is more sensitive to modelling uncertainties in the inner halo. The discrepancy in the DM particle mass constraints between these two approaches demonstrates the role that uncertainties in galaxy formation physics play in analyses of this type and motivates continued efforts to further our understanding of these processes. It also shows that the incompleteness of existing surveys of the MW satellite galaxy population contributes to analysis uncertainties. Future surveys such as the Legacy Survey of Space and Time (LSST) to be carried out by the Vera C. Rubin Observatory will improve the sky coverage and depth of extant surveys of the MW halo and help to tighten the uncertainties on DM particle mass constraints.
Two important aspects of the reionization of the Universe affect the formation of the low-mass galaxy population. The timing of the end of reionization influences how many low-mass DM haloes are able to accrete cold gas for use in star formation prior to reionization. The later reionization finishes, the more time is afforded for faint galaxies to form in such haloes. After this, further cold gas accretion is limited to those haloes that are massive enough that the gas can condense out of the intergalactic medium and onto the galaxy. The star formation that this facilitates enables the faintest galaxies to become brighter, changing the shape of the faint end of the satellite galaxy luminosity function [82]. These processes are reflected in our constraints (see Figures 6 and 10), where we find that an epoch of reionization that finishes earlier (i.e. at higher values of reion ) and a larger cooling threshold ( cut ) produce the most stringent constraints on the thermal relic particle mass. At high MW halo mass well away from the lower limit of the constraint envelope, the value chosen for cut has the largest effect on the number of substructures with a luminous component, in agreement with previous work, e.g. refs. [20,83]. However, close to the MW halo mass favoured by ref. [75], we find that the choice of reion has a significant effect on the constraints that can be placed on thermal relic models.
Our key results (see Figure 4) assume MW halo masses in the most likely range 200 =[0.5, 2.0] × 10 12 M . The constraints have only a moderate dependence on host halo mass because the number of MW satellite galaxies within a fixed radius inferred from observations scales much less strongly with halo mass than the number of subhaloes predicted by DM models (see Section 2.3). Better measurements of the mass of the MW halo will improve the constraining power of this approach; in the most extreme case, a MW halo with mass at the lowest end of the likely range would rule out thermal relic models with th ≤2.4 keV independently of galaxy formation physics. This estimate does not account for the effect of the central baryonic disc of the host halo that destroys subhaloes [96][97][98][99][100][101][102], which would exclude more of the WDM parameter space. For our fiducial galaxy formation model (see Figure 5), we also find that all DM particle masses are excluded for MW halo masses, 200 ≤ 0.6 × 10 12 M . This arises from the failure of the models to produce enough faint galaxies to be consistent with observations of the MW satellites, even in very cold thermal relic models where the number of low-mass subhaloes does not differ significantly from CDM predictions. Therefore, this threshold can be interpreted as a lower-mass limit for our Galaxy within the CDM model (see also [103,104]).
Recently, the EDGES collaboration announced the detection of a global 21 cm absorption line in measurements of the cosmic microwave background radiation [105]. This shows promise as a potential complementary probe of WDM models at high redshift because its shape and location ( =17.2) depend partly on the abundance of low-mass structures that act as sites of early star formation [106]. Currently, this epoch is inaccessible to other observational techniques [107]. Unfortunately, the 21 cm signal is very sensitive to uncertainties in the modelling of the Galactic foreground and in our understanding of the physics of star formation at early times. Therefore, the current data cannot constrain the properties of the DM [108][109][110]. Future studies of the statistics of the spatial distribution of the 21 cm signal and further work to understand stellar evolution at high redshift will overcome these difficulties [109,110].
The size of the satellite population inferred by the ref. [22] method is a lower limit to the true population as it cannot account for spatially-extended dwarf galaxies that fall below the surface brightness threshold of the surveys. Additionally, it does not encompass the contribution of the former satellites of the Large Magellanic Cloud that lie outside the DES footprint that could increase the size of the satellite complement still further. Taken together, these caveats strengthen the robustness of our lower limits on the thermal relic particle mass as a larger inferred satellite complement would rule out an even larger region of WDM parameter space.

Conclusions
In the continued absence of the direct detection of a DM particle or the observation of an astrophysical phenomenon that unambiguously constrains its properties, the debate about its exact nature and the acceptability of the current cosmological paradigm will continue. The discussion of 'small-scale' challenges to ΛCDM -perceived discrepancies between the observations of low-mass galaxies and predictions of DM substructure -has renewed impetus in this regard and has encouraged further exploration of alternative DM models. One class of these, which are broadly termed WDM models, produces a cut-off in the linear matter power spectrum that leads to a suppression in the formation of DM haloes on the scale of (and smaller than) those that would usually host dwarf galaxies in ΛCDM. The location and nature of this suppression depends sensitively on the properties of the DM particle. One method to constrain the parameter space of these models is the use of sophisticated hydrodynamic simulations to simulate self-consistently the formation and evolution of dwarf galaxies in the Local Group, and around MW-like hosts in particular. However, the resolution that would be required to achieve this in a volume that is large enough to attain high statistical power is, at present, computationally challenging. The development of other approaches to explore efficiently the viability of different cosmological models on these scales is, therefore, important.
In this work, we improve a method to constrain the properties of WDM models by comparing Extended Press-Schechter (EPS) predictions of the amount of substructure within MW-mass WDM haloes with the most recent estimates of the size of the satellite population of the MW (see Sections 2.3 and 2.4). This approach is complementary to previous work and for the first time accounts fully for limitations in the resolution of N-body cosmological simulations, incorporates the scatter in the number of substructures inside haloes at fixed DM halo mass, and includes the uncertainty associated with estimates of the number of satellite galaxies in the MW. The constraints that can be produced by this method are efficient at ruling out WDM models independently of any particular choice of galaxy formation physics, making the results highly robust.
We demonstrate the utility of this approach by applying it to thermal relic WDM models to constrain the DM particle mass (see Section 3.1). Our most robust constraint rules out, with 95 per cent confidence, thermal relic WDM particles with masses th ≤2.02 keV when marginalizing over uncertainties in estimates of the MW halo mass. This is competitive with existing limits that also use the abundance of MW satellite galaxies to constrain the WDM parameter space with minimal assumptions; however, our approach accounts for small subhaloes in N-body simulations that are not identified by substructure finders for numerical reasons, even though some of them actually survive to =0. Excluding them from the subhalo catalogue reduces the number of subhaloes that are available to host dwarf galaxies, artificially strengthening restrictions on the viable thermal relic model parameter space (see Figure 4). This effect worsens as the simulation resolution becomes poorer, so constraints that are obtained using lower-resolution simulations without accounting for the 'prematurely destroyed' subhaloes are significant overestimates.
All methods that seek to constrain the properties of DM models using visible tracers of the underlying substructure must make assumptions about galaxy formation processes that affect the satellite complement of the MW. Here, to obtain our highly robust constraints on the allowed properties of candidate WDM particles independently of galaxy formation physics, we have made the minimal and conservative assumption that a galaxy forms in all DM haloes. This allows us to place stringent lower bounds on the parameter space of thermal relic WDM models. In reality, baryonic physics mechanisms are important to determine the fraction of DM haloes that go on to host visible galaxies at late times, leaving many small subhaloes 'dark' [111]. While the details of these processes are still not understood fully, they are now constrained quite well. Accounting for these physical processes in models reduces the effective size of the satellite complement and in our analysis this improves significantly the constraints on the WDM particle properties.
Of particular interest to this study, the reionization of hydrogen in the early Universe, and the size of DM haloes in which it suppresses galaxy formation, dominates the formation and evolution of low-mass galaxies and imprints a characteristic signature on the luminosity function of MW satellite galaxies. We use the Durham semi-analytic model to explore several possible descriptions of this process and examine how different parametrizations affect the constraints on thermal relic WDM (see Section 4.1). By assuming that reionization is complete by reion =7 and that galaxy formation is suppressed in DM haloes with circular velocity vir <30 km s −1 , we rule out with 95 per cent confidence thermal relic DM with mass th ≤3.99 keV, when marginalizing over uncertainties in estimates of the MW halo mass (see Figure 5). We also find that a MW halo mass below 200 =0.6 × 10 12 M would not permit any thermal relic models that are warmer than CDM. This improves on the ref. [20] result and is competitive with conservative astrophysical limits from recent analyses using the Ly α forest. Furthermore, we find that the redshift at which reionization is assumed to cease has a significant effect on the constraints near to the most likely MW halo mass; however, for large MW halo masses the value chosen for the cooling threshold is more important (see Figure 6). Continued efforts to constrain further the probable ranges of the reionization parameters and the mass of the MW are therefore crucial if we wish to place ever more stringent constraints on the viability of alternative models to the ΛCDM paradigm.
While a DM particle candidate remains undetected, WDM models remain a feasible alternative to CDM. The satellite galaxy system of the MW provides a powerful means of probing structure formation on small scales and can help to discriminate between different cosmological models. However, the MW may not be typical of most DM haloes of similar mass. Hydrodynamic simulations that self-consistently model star formation and gas physics on the scale of dwarf galaxies will facilitate more robust astrophysical tests of this; however, achieving sufficient resolution is computationally challenging at present. A complementary means of testing the predictions of structure formation from different cosmological models is to consider their predictions of the evolution of structure across a range of mass scales and in a variety of environments, and to compare these with observations. Currently, this is challenging as it is difficult to identify the faintest and most extended objects at vast distances against observational backgrounds. Future improvements in observational capability will offer the prospect of further constraining the parameter space of viable WDM models.

A Resolution effects in numerical convergence studies
Numerical simulations are a useful tool to study the physical behaviour of cosmological models in the non-linear regime, where analytical approaches are unable to account fully for the complexity at these scales. While the dynamic range of such simulations is vast, spanning many orders of magnitude, N-body simulations are limited by the resolution at which their smallest objects can be self-consistently modelled. It is important to understand whether the phenomena that are observed in the simulations occur for physical reasons, or whether they arise because of this limitation.
The traditional approach to identify the onset of resolution effects has been to conduct convergence studies, e.g. [112,113]. These entail re-running the same simulation at different resolution levels and comparing the results: those that are unaffected by an increase in the resolution are deemed to be converged. A number of studies using several different N-body simulations support this conclusion and suggest that the subhalo present-day mass function of DM haloes is converged down to approximately 100 simulation particles per object, e.g. [24,25,114]. Some of these low-mass subhaloes could be disrupted by numerical effects from the limited resolution of the simulation [26,[37][38][39]. Ref. [25] also show that configuration space structure finders are ineffective at identifying all substructure near the centre of simulated haloes. This resolution-dependent deficiency of the halo finding algorithms implies that some substructures may be missed. These effects complicate attempts to understand and characterize the convergence of the subhalo peak mass function, which is of interest for this study as peak mass correlates more strongly with the formation of a luminous component than the present-day halo mass. It also affects directly the calibration of the EPS formalism that we use to estimate the amount of substructure in MW-mass haloes.
In the peak mass function, resolution limitations can also affect the high-mass end as even haloes with large peak mass can be excluded from the =0 halo catalogue if they fall below the resolution limit. This could occur after many orbits of the host during which the subhalo experiences continuous tidal stripping of mass. It is important to correct for missing and 'prematurely disrupted' subhaloes as these can bias our results: as we discuss in the main text, under-predicting the true subhalo count produces overly stringent constraints on the WDM particle mass. We are also careful to distinguish these from the spurious haloes found in N-body WDM simulations, which are produced by artificial fragmentation of filaments due to numerical effects and should be removed from the halo catalogues.
The 'prematurely destroyed' subhaloes may be recovered relatively easily by tracing their constituent particles through the simulations and identifying whether they survive to the present day. Details may be found in appendix B of ref. [22]. Briefly, we use the ref. [115] merging scheme implemented in to carry out this procedure. This tracks the most bound particle of objects that fall below the resolution limit from the last epoch at which they were associated with a resolved subhalo. From this, a population of substructures is recovered that contains the 'prematurely destroyed' subhaloes and other objects that are disrupted by physical processes. We remove the latter from the recovered population if they satisfy one of the following criteria:

1.
A time has elapsed after the subhalo fell below the resolution limit, which is equal to or greater than the dynamical friction timescale.
2. At any time, the subhalo passes within the halo tidal disruption radius.
In both cases, the effects of tidal stripping and of interactions between orbiting subhaloes are ignored. The size of this correction to the COCO suite is not easy to ascertain as COCO does not have counterpart simulations with different resolution levels with which to conduct a similar convergence study. Instead, we use the Aquarius suite [24], the constituent simulations of which span a range of resolution levels that encompass that of COCO, to estimate the size of the effect of excluding the prematurely destroyed subhaloes. In Figure 7, we compare the subhalo peak mass functions of the Aquarius A halo simulated at four different resolution levels: 2, 3, 4, and 5. Aq Level 5 is simulated coarsely, with a DM particle mass, p = 3.14 × 10 6 M . The simulation resolution improves with decreasing level number, such that Aq Level 2 is simulated with a DM particle mass, p = 1.37 × 10 4 M (i.e. a factor of ∼200 times better mass resolution). The figure shows the subhalo count before and after recovering the population of missing and prematurely destroyed subhaloes. At high halo mass, the original and 'corrected' curves are consistent with the highest resolution simulation. As the resolution degrades, the lower-resolution simulations peel away from the Level 2 curves, with the lowest-resolution simulation turning off at the highest value of peak . This demonstrates the major consequence of limited resolution, which is particularly acute for low-mass objects: in the cases considered here for haloes with mass 200 ≥ 1.5 × 10 12 M , restoring the missing population increases the total subhalo abundance by an order of magnitude. However, as we discussed earlier, resolution effects are not confined to the low-mass regime and can also affect higher masses. Massive haloes can experience considerable tidal stripping after being accreted by a host, which can lead to their exclusion from the =0 halo catalogue. The resulting discrepancy between the original and corrected mass functions at high masses indicates that this population of 'missing' objects composes a non-negligible fraction of the subhaloes even in the high-mass regime. Therefore, 'traditional' convergence studies that do not account for missing and prematurely destroyed subhaloes cannot properly characterize these numerical effects on the peak mass function.
In Figure 7, we also plot for comparison the average mass function of COCO-COLD haloes with masses similar to the Aquarius A halo. The COCO-COLD and COCO-WARM simulations have a DM particle mass resolution that lies between that of the Aq Level 3 and Level 4 runs. Comparing the subhalo mass functions of the incomplete subhalo catalogues of COCO-COLD and Aq Level 3 suggests that subhaloes with peak 3 × 10 8 M are resolved well. However, after recovering the prematurely destroyed subhaloes, a comparison of the mass functions implies consistency at masses peak 5 × 10 6 M , approximately two orders of magnitude better than before. This is consistent with the correction to the Aq Level 4 simulation, which suggests that the same correction for prematurely disrupted subhaloes that we have shown to work well for the Aquarius Level 2 to 5 runs is also applicable to the two COCO simulations.

B Calibrating the Galform merger tree algorithm
Monte Carlo merger trees are generated within using an implementation of the ref.
[60] merger tree algorithm, which iteratively splits the present-day halo mass into different progenitor haloes as it progresses to higher redshifts. The algorithm depends on three free parameters: 0 =0.57, a normalization constant; 1 =0.38, which controls the mass distribution of the progenitor haloes; and 2 = − 0.01, which controls the halo-splitting rate. Ref. [60] calibrated these parameters by comparing the Monte Carlo progenitor halo mass functions at several redshifts with those from the Millennium simulation [116]. This follows the evolution of 2160 3 particles with mass, = 8.6 × 10 8 ℎ −1 M , resolving the halo mass function to ∼1.7 × 10 10 ℎ −1 M , which is three orders of magnitude larger than the regime of interest for this study. The merger trees produced from the best-fitting free parameter values derived from the Millennium simulation predict a factor of two times more galaxies at the faint end of the cumulative luminosity function in MW-mass haloes than is obtained by applying to the COCO suite.
To attempt to address this overestimate, we performed the ref. [60] calibration procedure using the COCO simulations and found best-fitting values of 0 =0.75, 1 =0.1 and 2 =−0.12. The resulting Monte Carlo merger trees produce a better match with the COCO merger trees; however, they remain discrepant across the range in satellite brightness. Consequently, when applying on the new Monte Carlo merger trees, this produces an overestimate of the faint end of the cumulative satellite galaxy luminosity function by a factor of ∼1.6 compared with that obtained by applying on the COCO merger trees directly (cf. the COCO+ and Monte Carlo luminosity functions in Figure 8). This discrepancy can be improved self-consistently only by altering the ref. [60] algorithm, which would require more thorough investigation and possibly the introduction of one or more additional free parameters; this is beyond the scope of this work.
Instead, to obtain a satellite luminosity function for the Monte Carlo merger trees that is in agreement with the cosmological predictions, we map the satellite magnitude, V , predicted in the 'Monte Carlo merger trees + ' case to that of the 'COCO + ' case by matching objects at fixed abundance, i.e. fixed sat per host. By carrying out this procedure, we construct a remapping relationship between the 'old' V and new values that are consistent with the N-body results. In Figure 9, we plot these relationships calculated for the CDM and 3.3 keV thermal relic WDM models in three bins in halo mass. For clarity, we plot only the remapping functions obtained for our fiducial model with reion =7 and cut =30 km s −1 . The error bars (CDM) and shaded region (WDM) indicate the bootstrapped 68 per cent confidence intervals on the remapping relationships in the halo The remapping functions are in excellent agreement across the range in halo mass in both DM models, and across almost the entire range in satellite brightness, although there is a small discrepancy between the CDM and WDM functions at the faint end. This corresponds to low-mass subhaloes near the cut-off scale in the WDM power spectrum, whose properties differ the most from their equal mass CDM counterparts. The differences in the formation histories of such objects in WDM and CDM models are modest [21], which explains the similarly modest discrepancy between the remapping functions of these models calculated here. We find similar results for the other parametrizations of reionization that we consider (not shown). Therefore, when calculating the results presented in Section 4.1 we use the CDM remapping relationships appropriate for each parametrization of reionization to adjust the -produced absolute magnitudes of dwarf satellite galaxies.
To check the remapping technique, we plot the corrected 'Monte Carlo merger trees + ' satellite luminosity function as a green curve in Figure 8. By construction, the mean satellite count matches the COCO predictions. More importantly, the 68 per cent scatter (represented by the green shaded region) is also in good agreement with the N-body results despite this not having been calibrated.

C Thermal relic mass constraints for different Galform results
Reionization plays an important role in the formation of low-mass dwarf galaxies and shapes the star formation history of the Universe more widely. In , reionization is described in terms of two key variables: the redshift by which reionization has ceased, reion , and the circular velocity cooling threshold, cut , below which galaxies and DM haloes are prevented after reionization from accreting cool gas from the intergalactic medium with which they might form more stars. To understand better how reionization affects the constraints on the thermal relic particle mass, we considered nine parameter combinations that span the allowed parameter range given our current observational constraints on reionization and galaxy formation models: cut =[25, 30, 35] km s −1 and reion = [6,7,8]. In the main paper, we showed how the DM particle mass constraints change when varying cut assuming reion =7, and when varying reion assuming cut =30 km s −1 , the results of which are presented in Figures 5 and 6. Here, we provide the constraints for parameter combinations assuming reion =6 and reion =8 (see Figure 10, left and right panels, respectively). In both cases, we also plot our fiducial constraint as a thicker solid line to facilitate easier comparison with these results.
The dependence of the constraints on reion and cut demonstrated in Figure 6 also holds for the parameter choices shown here. If reionization finishes later (Figure 10 left panel), the strength of the constraints weakens considerably and the choice of cut becomes significantly more important. In Table 1, we provide the particle masses at and below which thermal relic WDM models are excluded at 95 per cent confidence for each combination of reionization parameters that we consider in this study. Figure 10. Constraints on the thermal relic particle mass obtained assuming three values of cut for reion =6 (left panel) and reion =8 (right panel) within the galaxy formation model. As in Figure 5, parameter combinations to the left of and beneath the envelopes are ruled out with 95 per cent confidence. The thicker solid lines indicate the constraint envelope of our fiducial model with reion =7 and cut =30 km s −1 . The shaded regions indicate the 68 per cent confidence interval on the mass of the MW halo from [75].