Dark Matter Constraints from a Unified Analysis of Strong Gravitational Lenses and Milky Way Satellite Galaxies

Introduction

The ΛCDM cosmological paradigm assumes a cold, collisionless dark matter (CDM) particle and therefore predicts a plethora of dark matter structure and substructure on extremely small cosmic scales (e.g., Green et al. 2004;Diemand et al. 2005;Wang et al. 2020). It is often argued that small-scale structure measurements represent an outstanding test to this prediction (e.g., see Bullock & Boylan-Kolchin 2017 for a review); yet, our understanding of the distribution of dark matter structure on nonlinear scales is rapidly progressing. Recent analyses of Milky Way (MW) satellite galaxies using data over nearly the full sky-including the population of ultrafaint dwarf galaxies discovered by deep photometric surveys over the last decade-have only recently been performed (e.g., Drlica-Wagner et al. 2020;Nadler et al. 2020Nadler et al. , 2021)). Meanwhile, measurements of stellar streams from the Gaia mission are beginning to reach the requisite precision to infer the signatures of perturbations from nearby low-mass subhalos (Banik et al. 2021;Bonaca et al. 2019). On extragalactic scales, the number of compact-source strong gravitational lenses available for substructure analyses has drastically increased in recent years (e.g., Nierenberg et al. 2014Nierenberg et al. , 2017Nierenberg et al. , 2020)), and modeling efforts have advanced in step (e.g., Gilman et al. 2019Gilman et al. , 2020a;;Hsueh et al. 2020). Analyses of resolved distortion in extended strong lensing observations from adaptive optics and space-based imaging have also rapidly progressed (e.g., Hezaveh et al. 2016;Birrer et al. 2017;Vegetti et al. 2018).

All of the recent small-scale structure measurements outlined above are consistent with the CDM paradigm and have therefore been used to constrain microphysical properties of dark matter that would reduce its small-scale clustering (Banik et al. 2021;Gilman et al. 2020a;Nadler et al. 2021). Although analyses of different probes reach consistent dark matter constraints, to date they have been performed independently and with different modeling assumptions to address heterogeneous astrophysical systematics. Crucially, if evidence for a departure from the CDM paradigm arises, it must be confirmed across different redshifts and physical scales. It is therefore critical to jointly model and analyze small-scale structure probes. This effort will be particularly important to maximize small-scale structure measurements from next-generation surveys including the Rubin Observatory Legacy Survey of Space and Time (LSST), which will enable the discovery of vastly more strong gravitational lenses (e.g., Collett 2015) and revolutionize the search for dwarf galaxies and measurements of stellar streams in the local universe (e.g., Drlica-Wagner et al. 2019).

Jointly modeling the low-mass halo and subhalo populations relevant for various small-scale structure measurements requires precise theoretical predictions for the abundance and structure of these small systems-which probe highly nonlinear cosmological modes-as a function of redshift and environment. Even cosmological parameters play an important role given the precision of current data; for example, varying the running of the spectral index within Planck uncertainties significantly affects predictions for subhalo abundances (Stafford et al. 2020), while other cosmological parameters including Ω m and σ 8 have subleading effects that may become important to incorporate in models of next-generation smallscale structure data (Dooley et al. 2014). Moreover, a variety of other theoretical and numerical uncertainties must be marginalized over in joint likelihood analyses to robustly claim evidence for non-CDM physics. For example, specific systematics of interest for modeling the MW satellite galaxy population include the faint end of the galaxy-halo connection, the total mass of the MW halo, and the mass and accretion time of the Large Magellanic Cloud (LMC; Newton et al. 2020;Nadler et al. 2021). Meanwhile, the orbits of dark matter subhalos in the inner regions of the MW halo must be predicted precisely in a statistical sense while accurately modeling the effects of specific baryonic structures to infer dark matter properties from stellar stream measurements. For strong lensing, the mass-concentration relation in both CDM and alternative dark matter models is a key uncertainty that must be accounted for (Hezaveh et al. 2016;Gilman et al. 2020aGilman et al. , 2020b;;Minor et al. 2020), along with the host halo properties and selection functions of strong lenses, all while accurately modeling the differential signal contributed by substructure and small halos along the line of sight (Despali et al. 2018;Gilman et al. 2019).

Here, we perform a joint analysis of small-scale dark matter measurements by combining the results of recent strong gravitational lensing and MW satellite inferences in a multidimensional parameter space to break modeling degeneracies. 8 In particular, we combine these results in a parameter space that includes the mass scale describing the suppression of the subhalo mass function for thermal relic warm dark matter (WDM) and the amplitude of the projected subhalo mass function at the strong lensing host halo mass and redshift scale. In particular, we combine the constraints on these quantities derived from (i) the magnification and flux ratio data from quadruply imaged strong gravitational lenses presented in Nierenberg et al. (2014Nierenberg et al. ( , 2017Nierenberg et al. ( , 2020) ) and modeled in Gilman et al. (2020a), and (ii) the abundance and properties of MW satellite galaxies over ∼75% of the sky presented in Drlica-Wagner et al. (2020) and modeled in Nadler et al. (2020Nadler et al. ( , 2021)). We employ the semianalytic model Galacticus (Benson 2012;Pullen et al. 2014) to translate the subhalo population inferred from MW satellite measurements to the strong lensing host halo mass and redshift scale by calibrating its predictions using cosmological zoom-in simulations of MW-like halos Mao et al. (2015). Thus, our work lays the foundation for joint semianalytic models of small-scale structure that are benchmarked by high-resolution simulations at each halo mass and redshift scale.

Our joint analysis breaks degeneracies among the amplitude of the projected subhalo mass functions inferred from MW satellite and strong lensing observations, thereby improving limits on deviations from the CDM paradigm that have been derived independently from these data. Specifically, we show that our combined analysis improves the lower limit on the WDM particle mass by ∼30%. The framework we develop for combining MW satellite and strong lensing data is particularly important because strong lensing is potentially sensitive to the presence of halos with masses below the threshold for galaxy formation, a mass scale that dwarf galaxy observations constrain. We therefore quantify the observational and theoretical advances necessary to robustly infer the presence of such dark halos, showing that this outcome is within the reach of next-generation small-scale structure measurements.

This paper is organized as follows. In Section 2, we describe the analytic model of dark matter substructure that underlies our joint analysis. We then describe the MW satellite data and model that enters our analysis in Section 3 and the strong lensing data and model in Section 4. We combine these analyses in Section 5, present our results in Section 6, discuss key systematics and compare them to previous work in Section 7, and conclude in Section 8. Throughout, we adopt the following cosmological parameters, following both Gilman et al. (2020a) and Nadler et al. (2021): h = 0.7, Ω m = 0.286, Ω Λ = 0.714, σ 8 = 0.82, and n s = 0.96 (Hinshaw et al. 2013).

Projected Subhalo Mass Function

Strong lensing and MW satellites probe low-mass subhalos within host halos at different mass and redshift scales. Specifically, strong lensing probes both the projected dark matter substructure in the lens system and small-scale structure along the line of sight, while MW satellites probe the three-dimensional distribution of subhalos traced by luminous satellite galaxies within the MW. Because current strong lensing analyses are not highly sensitive to the radial distribution of subhalos within the host halo of the lens, we focus on the statistics of projected subhalo populations in this paper, although we will describe how observations of the radial distribution of MW satellites break model degeneracies.

To simultaneously predict the subhalo populations relevant for MW satellite and strong lensing studies, we construct an analytic model for projected subhalo abundances that depends on the host halo mass, M host , and redshift, z host . In particular, we express the projected SHMF-i.e., the differential number of subhalos within a host halo, in projection-by generalizing the form in Gilman et al. (2020a),

where M denotes subhalo mass, A denotes the unit area, M 0 = 10 8 M e , and α is the power-law slope of the SHMF. In Equation (1), Σ sub (M host , z host ) is the projected number density of subhalos within the virial radius of a host halo of mass M host at redshift z host , including the effects of baryonic physics, and

host host captures the dependence of the projected SHMF on M host and z host in CDM only (i.e., modulo the effects of baryonic physics) as described in Section 2.2. We discuss the impact of halo boundary definitions in Section 5.

Our model of the MW satellite population is based on the Nadler et al. (2021) analysis, which defines subhalo mass using the peak Bryan-Norman virial mass M peak (see Appendix A for details). Meanwhile, our strong lensing constraints are based on the Gilman et al. (2020a) analysis, which uses M 200 values relative to the critical density at z = 0, with subhalo masses evaluated at infall to compute the WDM SHMF and massconcentration relation. Here, we simply interpret the peak virial mass values from our MW satellite analysis as M 200 values at infall. The peak virial masses of subhalos relevant for our work are on average a factor of ∼2 larger than their M 200 values at infall in the cosmological zoom-in simulations our MW satellite analysis is based on, largely due to pre-infall tidal stripping (e.g., Behroozi et al. 2014;Wetzel et al. 2015). Thus, converting M peak to M 200 would further strengthen our joint WDM constraints; however, because changing this mass definition would nontrivially affect the abundance-matching model used in our satellite analysis, we leave a detailed investigation of this point to future work that combines satellite and lensing inferences at the likelihood level.

CDM Host Mass and Redshift Dependence

We model the dependence of the projected SHMF on host halo mass and redshift with the functions Σ sub (M host , z host ) and ( )  M z , CDM host host . Both of these terms play a key role in our joint analysis because they allow us to relate the subhalo populations corresponding to low-redshift, group-mass strong lens host halos (M lens ∼ 10 13 M e , z lens ∼ 0.5) to the regime of the MW halo today (M MW ∼ 10 12 M e , z MW = 0).

( )  M z , CDM host host captures the dependence of the SHMF on host halo mass and redshift in CDM only (i.e., without baryons), including the effects of tidal disruption by the dark matter host halo. This scaling depends on both the statistics of subhalo populations at infall, which can be predicted reasonably precisely using extensions of the Press-Schechter formalism (Press & Schechter 1974), and on the dynamical evolution of subhalos after infall into a host. Detailed semianalytic models calibrated to N-body simulations are necessary to model this evolution; we therefore follow Gilman et al. (2020a) in using the Galacticus model (Benson 2012;Pullen et al. 2014;Yang et al. 2020), which predicts

where k 1 = 0.88 and k 2 = 1.7. Section 2.3 describes our model for subhalo disruption due to baryons, which captures the leading-order corrections to ( )  M z , CDM host host . We do not model the impact of additional host halo, central galaxy, and environmental variables on the projected SHMF, noting that this is an important area for future work that ongoing observational efforts like the Satellites Around Galactic Analogs (SAGA) survey are informing at the MW-mass scale (Geha et al. 2017;Mao et al. 2021). However, we emphasize that the Nadler et al. (2021) MW satellite analysis our work is based on self-consistently uses simulations that are consistent with key secondary MW halo properties, including concentration, the existence of a realistic LMC analog system, and a formation history constrained by Gaia observations. Meanwhile, the Gilman et al. (2020a) lensing analysis is not sensitive to host-to-host variation in Σ sub beyond that modeled by ( )  M z , CDM host host given the current number of strong lenses studied and the information available per lens.

Subhalo Disruption Efficiency Due to Baryons

We model Σ sub (M host , z host ) with explicit host halo mass and redshift dependence to capture the impact of baryonic physics on the projected SHMF. This extra dependence relative to the CDM scaling is not captured by the ( )  M z , CDM host host term predicted by Galacticus, although baryonic effects can be modeled in future Galacticus implementations. Although Σ sub is not modeled with explicit host mass and redshift dependence in Gilman et al. (2020a), we include this dependence here because the subhalo populations probed by strong lensing and MW satellites are subject to baryonic effects that potentially impact the two regimes differently. Of these effects, the most important is tidal disruption due to the central galaxy, which suppresses the SHMF at the ∼50% level (e.g., Despali & Vegetti 2017;Garrison-Kimmel et al. 2017;Graus et al. 2018;Kelley et al. 2019;Richings et al. 2020;Webb & Bovy 2020). 9 Tidal disruption due to the central galaxy most strongly suppresses the abundance of subhalos in the inner regions of the host halo or (more precisely) subhalos that accrete early and have close pericentric passages (Garrison-Kimmel et al. 2017;Nadler et al. 2018).

The projected SHMF is largely driven by the plethora of subhalos in the outer regions of the host halo, which mitigates the impact of uncertainties in the strength and radial dependence of these baryonic effects on our probe combination. Nonetheless, our joint analysis is sensitive to both the amplitude of and differences in the efficiency of subhalo disruption due to baryonic physics as a function of host halo mass and redshift. We measure Σ sub (M host , z host ) in units of its value for strong lens host halos, and we define the differential subhalo disruption efficiency due to baryons as where Σ sub hereafter denotes the projected subhalo number density for strong lens host halos, following Gilman et al. (2020a), and Σ sub,MW denotes the same quantity for the present-day MW system. In Equation (3), q represents the efficiency of subhalo disruption due to baryonic physics in the MW at z = 0 in units of the efficiency of subhalo disruption due to baryonic physics in the group-mass halos and at the redshifts probed by strong lensing. Note that larger (smaller) values of q represent less efficient (more efficient) subhalo disruption in the MW relative to strong lenses and that differences in the radial dependence of subhalo disruption at these scales (which we do not model) do not affect our joint analysis of projected SHMFs. Motivated by the results of hydrodynamic simulations, we assume that subhalo disruption due to baryonic physics results in a mass-independent rescaling of the MW and strong lens projected SHMFs. This allows us to respectively express the projected SHMFs probed by strong lensing and MW satellite observations as As noted above, strong lenses typically have halo masses of M lens ≈ 10 13 M e and redshifts of z lens ≈ 0.5 and host massive elliptical galaxies (Gavazzi et al. 2007;Auger et al. 2010;Gilman et al. 2020a). In contrast, the MW has a halo mass of M MW ∼ 10 12 M e (e.g., Callingham et al. 2019;Cautun et al. 2020) at z MW = 0, and is largely typical for a spiral galaxy of its stellar mass, although it has a relatively quiescent formation history (e.g., Boardman et al. 2020;Evans et al. 2020). Subhalo disruption due to the central galaxy in hydrodynamic simulations of MW-mass systems reduces the amplitude of the SHMF by ∼50% relative to corresponding dark-matteronly simulations. This effect is roughly mass independent and is not a strong function of redshift at late times. Although hydrodynamic simulations of group-mass systems yield similar levels of SHMF suppression (Fiacconi et al. 2016;Graus et al. 2018;Richings et al. 2021), this regime is less well studied. We therefore adopt q = 1 in our fiducial analysis-i.e., equally efficient subhalo disruption due to baryons in the MW and strong lens host halos-and we also test values within a range of q ä [0.5, 2] when translating the projected SHMF amplitude inferred from MW satellites to the strong lens host halo regime. We emphasize that, in detail, subhalo disruption is expected to depend on the mass and formation history of the central galaxy along with host halo mass and redshift, and its efficiency will therefore differ among strong lenses. Although our phenomenological model for differences in subhalo disruption due to baryonic physics is very simple, we will demonstrate that the corresponding uncertainties do not significantly impact our joint dark matter constraints.

Warm Dark Matter Mass-Concentration Relation

The delay in the collapse of small-scale density perturbations in WDM suppresses the central densities of halos with masses near M hm , altering the mass-concentration relation for both field and subhalos. Because flux ratios in strong lenses are highly sensitive to the central densities of subhalos, the altered mass-concentration relation provides crucial information relevant for forwardmodeling strong lensing signals (Gilman et al. 2020b). We implement the WDM mass-concentration relation in a similar manner to the suppression of the SHMF (Gilman et al. 2020a),

where β(z) = 0.026z -0.04.

Halo concentrations are affected over an order of magnitude in mass above the turnover in the mass function set by M hm . Thus, the mass-concentration relation must be accounted for to self-consistently constrain WDM-like models using strong lensing data. Meanwhile, MW satellite abundances are relatively insensitive to the mass-concentration relation because subhalo disruption is mainly determined by subhalos' orbital properties. Moving beyond abundances, the internal dynamics of relatively bright MW satellite galaxies are often subject to baryonic effects that make it difficult to robustly infer halo concentration (e.g., Read et al. 2019). Meanwhile, it is difficult to obtain precise dynamical measurements given the limited number of spectroscopically confirmed stars associated with the faintest MW satellites; however, future spectroscopic measurements of these galaxies may reach the precision necessary to provide complementary constraints (Simon et al. 2019).

Galaxy-Halo Connection Model

To infer the present-day abundance of subhalos in the MW, Nadler et al. (2021) combine the simulations described above with an empirical model of the galaxy-halo connection (introduced in Nadler et al. 2019bNadler et al. , 2020)), which populates subhalos with satellite galaxies in a parametric fashion. By combining these predictions with observational selection functions derived from satellite searches in DES and PS1 data (Drlica-Wagner et al. 2020), the model is compared to observations assuming that satellites in each survey footprint populate the parameter space of surface brightness and heliocentric distance according to a Poisson process (see Nadler et al. 2020, Section 6). By marginalizing over the underlying Poisson rate in the calculation of the likelihood for each surface brightness bin, the galaxy-halo connection and dark matter model parameters are fit to data in a Markov Chain Monte Carlo (MCMC) framework.

The majority of the parameters in the Nadler et al. (2021) WDM analysis govern the relationship between satellite galaxies and the subhalos they inhabit. For example, these include the slope and scatter of the abundance-matching relation between galaxy luminosity and peak halo maximum circular velocity; the amplitude, scatter, and power-law slope of the relation between galaxy size and halo size; and parameters governing the fraction of low-mass dark matter halos that host observable galaxies. These parameters are not directly relevant for our strong lensing joint analysis because lensing is sensitive to the integrated amount of matter in the lens galaxy and along the line of sight, which is dominated by dark matter. However, they are crucial for robustly modeling the MW satellite population and are marginalized over in our probe combination.

Here, we highlight the aspects of the Nadler et al. (2021) model that are most relevant for our joint analysis:

satellite analysis is consistent with CDM predictions down to a characteristic halo mass scale referred to as the minimum halo mass ( min ). The minimum halo mass is defined as the peak virial mass of the smallest surviving subhalo inferred to host observed MW satellite galaxies and therefore represents the lowest mass down to which the SHMF is directly constrained by current MW satellite observations.  min is jointly inferred along with the fraction of halos that host observable galaxies, which is consistent with 100% down to  min . The upper limit on  min is calculated by marginalizing over the full posterior distribution, which yields

at 95% confidence in the Nadler et al. (2020) CDM fit (see Nadler et al. 2020, Sections 4.4, 7.4, and 7.5).

(ii) Baryonic disruption efficiency (): The efficiency of subhalo disruption due to the Galactic disk. Disruption probabilities due to baryonic physics for the subhalos in the dark-matter-only simulations described above are predicted using the To account for uncertainties resulting from the limited statistics of these training simulations, Nadler et al. (2021) parameterize the efficiency of subhalo disruption by assigning the following disruption probability to each subhalo in the MW-like zoom-in simulations: 6 M e (m WDM > 6.5 keV) at 95% confidence. We discuss the role of the MW halo mass in detail in Section 7.1.

Constraints from Milky Way Satellite Observations

Here, we rerun the Nadler et al. (2021) WDM MW satellite analysis using priors and a WDM SHMF parameterization chosen to match the Gilman et al. (2020a) lensing analysis, which allows us to self-consistently perform our multidimensional satellite-lensing probe combination. In particular, we rerun the MW satellite analysis adopting a uniform prior of ( )   0, 3 , which ensures that we match the shape of the Σ sub prior used in Gilman et al. (2020a) based on the linear relation between  and Σ sub we derive in Section 5.2. The use of a uniform (rather than lognormal) prior on  weakens the upper and lower limits of the marginalized posterior from Nadler et al. (2021) (i.e., < <  0.2 1.9) from 95% to 68% confidence constraints. In Appendix D, we show that the choice of this prior does not significantly impact our joint WDM constraints.

We also use the WDM SHMF and the ( ) ~ M log 5, 10 hm prior assumed in Gilman et al. (2020a). The resulting marginalized posterior distribution yields M hm < 10 7.4 M e (m WDM > 7.4 keV) at 95% confidence after MW host halo mass scaling, which is more constraining than the Nadler et al. (2021) result despite using a slightly less suppressed WDM SHMF. This is caused by the change in the lower limit of our M log hm prior, which is two orders of magnitude lower than that adopted in Nadler et al. (2021). As described in Section 6.1, the lower limit of the M hm prior is arbitrary unless we assume that WDM physics manifests at a particular halo mass scale. Thus, we also quote likelihood ratios for both our independent and combined constraints. We find that M hm = 10 7.9 M e (m WDM = 5.2 keV) is disfavored relative to the peak of the marginalized posterior at 10 5 M e with a 20:1 ratio, consistent with the Nadler et al. (2021) result.

Figure 1 shows the posterior from our updated WDM fit to MW satellite data in the two-dimensional parameter space of M hm versus , marginalized over seven other galaxy-halo connection parameters. In Figure 1 and subsequent plots, we do not scale parameters to account for MW host halo mass uncertainty unless explicitly noted. We reiterate that our MW satellite analysis only probes systems down to a peak halo mass threshold of ∼3 × 10 8 M e at 95% confidence and that M hm constraints below this mass scale are driven by the population statistics of halos near the minimum observable halo mass. This is demonstrated in Nadler et al. (2021, see their Figure 1), which shows that the WDM model ruled out by MW satellites at 95% confidence yields ∼25% suppression in subhalo abundances relative to CDM at the minimum halo mass, which is about one order of magnitude larger than M hm .

There is not a strong degeneracy between  and M hm in Figure 1 because  models the disruptive effects of the MW disk, which suppresses the inner radial distribution of MW satellites in an approximately mass-independent fashion, while M hm models the mass-dependent suppression of the projected SHMF caused by WDM free streaming. Figures 2 and7 illustrate the effects of  on the projected SHMF and radial distribution of our MW-like simulations.

Strong Lensing Model

The presence of small-scale structure in the lens and along the line of sight can perturb the magnified fluxes of unresolved quasar emission regions. The occurrence rate of the distribution of perturbed flux ratios between the multiple images is therefore sensitive to the underlying population of (potentially dark) subhalos within the host halo of the lens. Importantly, the strong lensing image position and flux ratio data described above are also sensitive to dark matter structure along the entire line of sight from the observer to the lensed quasars.

To perform the inference on the underlying subhalo and lineof-sight mass function population parameters, Gilman et al. (2020a) forward-model the quasar flux ratio with a large set of realizations of the small-scale lensing structure through a multiplane ray-tracing scheme, which accounts for the finite emitting source size and satisfies the astrometric constraints on the positions of the images. The likelihood for the individual lenses' population parameters was constructed with Approximate Bayesian Computation (ABC), and the joint posterior inference was performed by multiplying the individual likelihoods (see Gilman et al. 2020a Section 2).

Unlike in the case of subhalos, the line-of-sight dark matter structure is unaffected by tidal stripping and disruption. Thus, Gilman et al. (2020a) modeled the line-of-sight dark matter distribution using a Sheth-Tormen (Sheth et al. 2001) mass function with the same WDM SHMF and concentration suppression factors described above for subhalos, along with a contribution from the two-halo term near the main deflector's host halo and an overall scaling factor that allows for uncertainty in the halo mass function amplitude of 20% (see Gilman et al. 2020a, Section 5.3).

Constraints from Strong Lensing Observations

The Gilman et al. (2020a) strong lensing analysis is consistent with CDM predictions for the slope and amplitude of the halo and subhalo mass functions. In particular, Gilman et al. (2020a) derive constraints on the SHMF slope that are consistent with N-body simulations and find that the line-ofsight contribution is consistent with Sheth-Tormen mass function predictions. Moreover, Gilman et al. (2020b) demonstrate that these data are compatible with standard predictions for the CDM mass-concentration relation while self-consistently modeling the effects of tidal stripping on subhalos.

Here, we highlight the constraints from Gilman et al. (2020a) that enter our multidimensional MW satellite-lensing probe combination: 5)) using the host halo mass and redshift scaling predicted by Galacticus, evaluated at the average halo mass of our MW-like simulations with a slope of α = -1.92. Σ sub,MW is chosen such that the SHMF amplitude matches our simulations at the subhalo mass corresponding to the faintest observed MW satellites,  min (dashed vertical line). Dark (light) red contours show 68% (95%) confidence intervals from Galacticus for host halos with characteristics matched to our MW-like simulations. We impose the resolution cuts described in Appendix A on the simulation and Galacticus results.

confidence. This constraint results from the fact that warmer models suppress the abundance and concentrations of lowmass halos that contribute to the lensing signal (see Gilman et al. 2020a Sections 3.4 and 6.2).

Here, we reanalyze the Gilman et al. (2020a) marginalized M hm posterior using a slightly higher lower limit of the M log hm prior. We find M hm < 10 8 M e (m WDM > 4.9 keV) at 95% confidence, which is slightly less constraining than the Gilman et al. (2020a) result. Again, we also calculate likelihood ratios due to the ambiguity of the M hm prior and find that M hm = 10 8.7 M e (m WDM = 3.0 keV) is disfavored relative to the peak of the marginalized posterior at 10 6.4 M e with a 20:1 ratio, consistent with Gilman et al. (2020a). 10The right panel of Figure 3 shows the posterior distribution from the fit to strong lensing data in Gilman et al. (2020a) in the two-dimensional parameter space of M hm versus Σ sub , marginalized over the SHMF slope and line-of-sight mass function amplitude. There is a moderate degeneracy between Σ sub and M hm , particularly at high values of Σ sub ; in this regime, it is difficult to distinguish the coincident suppression of the WDM SHMF and mass-concentration relation relative to CDM from changes to the normalization of the CDM SHMF.

Inferring Σ sub from Milky Way Satellites

To connect the host halo mass and redshift regimes probed by strong lensing and MW satellites, we first construct a relation between the subhalo disruption efficiency  inferred from our MW satellite analysis and the projected subhalo number density Σ sub,MW predicted by evaluating Galacticus at the MW halo mass scale. We then translate Σ sub,MW to Σ sub at the strong lensing host halo mass and redshift scale 2020a) strong lensing analysis. The vertical band labeled "Σ sub Prior" shows the range of Σ sub inferred from the MW satellite posterior in our fiducial joint analysis (i.e., 0.015 kpc -2 Σ sub 0.03 kpc -2 ). In both panels, color maps show the probability density normalized to its maximum value in each parameter space, and solid (dashed) white lines indicate 1σ (2σ) contours for a two-dimensional Gaussian distribution.

using the dark matter substructure model described in Section 2.

Figure 2 shows the average z = 0 projected SHMF from the two realistic zoom-in simulations used in our MW satellite analysis for our fiducial disruption model calibrated to hydrodynamic simulations (i.e., =  1) and for bracketing values of the subhalo disruption efficiency (i.e., =  0 and =  2) that are ruled out at greater than 68% confidence by MW satellite data as discussed in Section 3. To construct the projected SHMF predicted by our analytic substructure model, we use Equation We construct a relation between Σ sub,MW and  by matching our analytic prediction from Equation (5) to the average projected SHMF predicted by our MW-like simulations as a function of , as illustrated in Figure 2. In particular, we match the amplitudes of the projected SHMFs inferred from our MW satellite analysis and predicted by Equation (5) at the minimum observable halo mass of 3.2 × 10 8 M e and within a fixed radius of 300 kpc (roughly corresponding to the virial radius of the MW host halo), chosen to match the Nadler et al. (2021) analysis. Our choice to match these SHMFs at the minimum halo mass is conservative because our MW satellite analysis is not sensitive to subhalos below this mass scale at 95% confidence. We then translate Σ sub,MW to the strong lensing regime using Equation 3, which yields

This relation allows us to transform the -M hm MW satellite posterior from Figure 1 into the Σ sub -M hm parameter space for a given value of q. Note that Σ sub scales linearly with MW halo mass because it measures the projected SHMF amplitude within a fixed physical radius in our model. The -Σ sub,MW relation constructed above is only based on the amplitude of the projected SHMF from our two realistic MW-like simulations, measured at the minimum halo mass scale of M peak = 3 × 10 8 M e . In general, the procedure to infer Σ sub,MW from an estimate of the MW SHMF should account for both host-to-host variation in the SHMF within the range of MW host halo properties allowed by observations and Poisson scatter in the SHMF given the range of subhalo masses probed by MW satellite measurements. Our SHMF matching procedure is intentionally simplistic because-as demonstrated in Appendix C-host-to-host and Poisson scatter in the projected SHMF near the minimum halo mass scale are both subleading uncertainties compared to the range of differential subhalo disruption efficiencies due to baryons that we explore. We therefore leave a statistically rigorous construction of the -Σ sub relation to future work that propagates such uncertainties into the joint analysis at the likelihood level.

The result of the transformation in Equation ( 14) is shown in the left panel of Figure 3 for our fiducial model of q = 1 (i.e., for equally efficient subhalo disruption in the MW and strong lensing host halos). The typical Σ sub values favored by the MW satellite posterior for this choice of q are significantly smaller than the largest values allowed by the Gilman et al. (2020a) lensing analysis; we return to this point below.

The lack of degeneracy observed between  and M hm in Figure 1, which results from the joint constraining power of the MW satellite radial distribution and luminosity function for subhalos near the minimum observable halo mass, persists in the Σ sub -M hm parameter space. On the other hand, strong lensing flux ratio statistics probe an integrated combination of subhalo masses and concentrations. The lensing analysis is currently less sensitive to subhalos in specific mass ranges than MW satellite population statistics and therefore exhibits a stronger Σ sub -M hm degeneracy in Figure 3. However, because lensing measurements do not depend on the connection between subhalos and luminous matter, they can probe subhalos below the minimum observable halo mass.

In Figure 2 and Appendix B, we also compare our simulation results to the SHMF in terms of both peak and present-day subhalo mass and to the radial subhalo distribution predicted by Galacticus for 14 halos selected to match the characteristics of our MW-like simulations. In all cases, we apply the same cuts on peak and present-day subhalo maximum circular velocity when comparing Galacticus to our zoom-in simulations; the details of these resolution cuts and our Galacticus runs are described in Appendix B. Note that these Galacticus predictions should be compared to our =  0 simulation results because the current Galacticus implementation does not model subhalo disruption due to central galaxies. The Galacticus-predicted SHMFs agree well with our simulations in terms of both peak and present-day subhalo mass, lending confidence to the choice of ( )  M z , CDM host host that enters our calibration of the projected subhalo number density at the MW scale.

Because we construct a Σ sub - relation based on our specific MW-like simulations, which have been shown to match the MW satellite population, we imposed several host halo selection criteria on the Galacticus runs used for the comparison above. These conditions include the existence of a realistic LMC analog system and a Gaia-Enceladus-like merger event. We emphasize that validating semianalytic halo and subhalo population predictions with self-consistent simulation suites of both MW-like and strong lens-like halos is an important avenue for future work. As discussed in Section 7.1, increasingly precise near-field observations and complementary data for strong lens systems will allow us to mitigate the impact of additional host halo properties including concentration (which is known to influence subhalo populations at fixed halo mass, e.g., Zentner et al. 2005;Zhu et al. 2006;Ishiyama et al. 2009;Mao et al. 2015;Fielder et al. 2019) and the characteristics of the local dark matter environment in future joint analyses.

Probe Combination Statistics

Having placed the lensing and MW satellite posteriors on the same footing, we now proceed to combine them to construct a joint Σ sub -M hm likelihood as follows. Formally, we write our joint MW satellite and strong lensing analysis as a combined Bayesian inference problem, where θ is the vector of parameters in both the lensing and satellite analyses including q (where the shared parameters Σ sub and M hm only appear once), θ MW (θ lensing ) are the parameters entering the MW satellite (strong lensing) inference, D = [D MW , D lensing ] is the joint datavector, and P(θ) is the prior distribution over all model parameters. Next, for a given value of q ä [0.5, 2], we marginalize over the independent parameters (i.e., the seven galaxy-halo connection parameters in the satellite analysis described in Section 3 and the SHMF slope and line-of-sight contribution in the lensing analysis described in Section 4) to arrive at a combined Σ sub -M hm posterior distribution, As described in Sections 3-4, we choose a lower limit of ( )  = M M log 5 hm because models with even lower M hm values are indistinguishable from CDM in both our MW satellite and strong lensing analyses. For simplicity, we repeat our analysis at several fixed values of q rather than marginalizing over this parameter. We note that our WDM limits marginalized over q are nearly identical to our fiducial (q = 1) result in the absence of a well-motivated, nonuniform prior for q based on hydrodynamic simulations (see Section 6.3).

Based on Equation (14), our MW satellite inference only samples Σ sub /kpc -2 ä [0.015q -1 , 0.03q -1 ] ä [0.0075, 0.06] given our prior of ( )   0, 3 and our assumed range of q ä [0.5, 2] (note that   0 by definition). Thus, our fiducial (i.e., q = 1) analysis is restricted to the range 0.015 Σ sub /kpc -2 0.03, labeled "Σ sub Prior" in Figures 3 and5. This range is set by combining the MW satellite posterior and zoom-in simulations with our analytic SHMF prediction and is narrower than the Σ sub range considered in Gilman et al. (2020a), which did not enforce priors based on cosmological simulations. This difference limits the range of Σ sub values from the lensing analysis relevant for our probe combination, but it does not formally affect our calculation of the marginal likelihood because the effective prior on Σ sub from the MW satellite analysis is nevertheless uniform.

Thus, exploiting the fact that our priors are uniform in both Σ sub and M log hm , we arrive at a joint marginal likelihood for these quantities in terms of the marginalized MW satellite and strong lensing posteriors, This joint marginal likelihood is illustrated in Figure 4 and analyzed in Section 6. Because the independently derived M hm distributions are consistent (Section 6.2) and the Σ sub distributions are only in mild tension (Section 6.4), we do not formally test for statistical consistency between the MW satellite and strong lensing analyses before constructing the joint likelihood.

Because we fix the slope of the projected SHMF in our analytic substructure model at α = -1.92, our joint analysis is effectively performed at a thin slice in α of the Gilman et al. (2020a) posterior (we remind the reader that α is defined in a CDM context well above the cutoff scale). We do not expect this choice to significantly affect our results because α is not highly degenerate with the parameters of interest (i.e., M hm and Σ sub ) in the Gilman et al. (2020a) analysis. However, in Section 7.1 we emphasize the importance of jointly inferring the SHMF slope in future work, and we discuss the role of the remaining line-of-sight dimension of the Gilman et al. (2020a) posterior, modeled by the amplitude of the Sheth-Tormen halo mass function, that is marginalized over in our analysis.

Conventions for WDM Constraints

To quantify the WDM constraints corresponding to the joint likelihood derived in Section 5.3, we marginalize over the Σ sub dimension and construct the following summary statistics from the marginal M hm likelihood:

(i) Confidence intervals: defined as the range of parameter values enclosing a particular fraction of the integrated marginal likelihood. Following common practice in the WDM literature, we quote upper limits on M hm and lower limits on m WDM at 95% confidence. (ii) Marginal likelihood ratios: defined as the parameter value at which the marginal likelihood probability density falls to a particular fraction of its peak value. Following Gilman et al. (2020a), we quote the M hm and m WDM values disfavored with 20:1 marginal likelihood ratios.

Although confidence intervals capture more information about the shape of the probability density and are commonly quoted in the WDM literature (e.g., Viel et al. 2013;Iršič et al. 2017), they depend on the arbitrary choice of a lower limit on the M hm prior (or equivalently, an upper limit on the m WDM prior) as noted above. In particular, small-scale structure data are currently consistent with CDM and therefore yield onesided limits on M hm or m WDM ; without assuming a preferred scale for a small-scale structure cutoff due to WDM (or other non-CDM) physics, this makes the lower limit of the M hm prior arbitrary. This situation motivated several authors (e.g., Enzi et al. 2020;Gilman et al. 2020a) to quote alternative summary statistics including marginal likelihood ratios that are less dependent on the choice of M hm prior, and we follow this practice here. Similarly, we follow both Gilman et al. (2020a) and Nadler et al. (2021) by adopting a logarithmic prior on M hm because any other choice would not be invariant to rescaling m WDM (e.g., see the discussion in Jethwa et al. 2018).

Fiducial WDM Constraints

We now present the results of our joint analysis for our fiducial subhalo disruption efficiency model of q = 1, which assumes equally efficient subhalo disruption due to baryons in the MW and in strong lens host halos, which is broadly compatible with the results of hydrodynamic simulations (see Section 2.3). The combined Σ sub -M hm marginal likelihood is shown in Figure 4 and the corresponding one-dimensional marginalized likelihoods for Σ sub and M hm are shown in Figure 5. The joint marginal likelihood retains the shape of the Σ sub -M hm distribution from the transformed MW satellite posterior and from the lensing analysis limited to the range of Σ sub inferred from our MW satellite analysis according to the procedure in Section 5.2. Moreover, the joint marginal likelihood visibly prefers lower values of M hm than either posterior alone, demonstrating the unique constraining power accessible when combining independent small-scale structure probes in a multidimensional parameter space. 4) for projected subhalo number density at the strong lensing scale (left panel) and WDM half-mode mass (right panel), assuming equally efficient subhalo disruption due to baryons in the MW and strong lens systems (q = 1). The marginalized MW satellite posterior is shown in blue, the marginalized strong lensing posterior is shown in red, and results obtained from our probe combination and marginalized over the remaining dimension are shown in purple. In the left panel, the vertical band labeled "Σ sub Prior" shows the range of Σ sub inferred from the MW satellite posterior in our fiducial joint analysis (i.e., 0.015 kpc -2 Σ sub 0.03 kpc -2 , slightly offset from the posteriors for visual clarity), and the dashed red line on the right panel shows the lensing half-mode mass posterior restricted to this range of Σ sub values. One-dimensional M hm Distributions q = 0.5 q = 1 q = 2 95% confidence level M hm (M e ) 10 7.2 10 7.1 10 7.0 10 6.9 95% confidence level m WDM (keV) 8.4 9.1 9.7 10.4 20:1 likelihood ratio M hm (M e ) 10 7.7 10 7.6 10 7.4 10 7.3 20:1 likelihood ratio m WDM (keV) 6.0 6.4 7.4 7.9

Note. q = 0.5 corresponds to twice as efficient subhalo disruption due to baryons in the MW relative to strong lenses, q = 1 (our fiducial model) corresponds to equally efficient subhalo disruption due to baryons, and q = 2 corresponds to twice as efficient subhalo disruption due to baryons in strong lenses.

Consistent with these qualitative aspects of the joint Σ sub -M hm likelihood, the upper limit of the marginal M hm likelihood shown in the right panel of Figure 5 is noticeably lower than either of the individual constraints from MW satellites or strong lensing. Quantitatively, the upper limit on M hm from our joint analysis improves upon those set by the MW satellite and strong lensing analyses individually by ∼60%, leading to a ∼30% increase in the strength of the lower limit on m WDM . Specifically, the 95% confidence limit of M hm < 10 7.4 M e (m WDM > 7.4 keV) from our MW satellite analysis improves to M hm < 10 7.0 M e (m WDM > 9.7 keV). We find a similar level improvement in terms of likelihood ratios, with M hm = 10 7.4 M e (m WDM = 7.4 keV) ruled out at 20:1 relative to the peak of the marginal likelihood at the lower limit of the prior at 10 5 M e .

To derive these limits, we conservatively increased the M hm values returned by our joint analysis by a factor of ∼25% to account for the maximum mass of the MW halo relative to the average host halo masses of our zoom-in simulations, following Nadler et al. (2021). As demonstrated in the following subsection, propagating the MW halo mass uncertainty into the Σ sub dimension would have a negligible impact on the results compared to uncertainties in the efficiency of subhalo disruption due to baryons, so we do not perform this scaling for simplicity.

Our fiducial constraint of m WDM > 9.7 keV at 95% confidence is one of the most stringent limits on the WDM particle mass set by small-scale structure observations to date. Moreover, it is set using only existing strong lensing and MW satellite measurements, underscoring the importance of unified, multidimensional small-scale structure analyses as the corresponding measurements continue to improve. Joint modelbuilding efforts that further incorporate Lyα forest (Viel et al. 2013;Iršič et al. 2017) and stellar stream (Banik et al. 2021) constraints while retaining the unique information provided by each probe will therefore be particularly fruitful.

Impact of the Differential Subhalo Disruption Efficiency Due to Baryons

We now explore the impact of the differential efficiency of subhalo disruption due to baryons on our WDM constraints. Table 1 lists the M hm and m WDM 95% confidence level and 20:1 likelihood ratio limits for q = 0.5, 1, and 2, and the right panel of Figure 6 shows the corresponding joint marginal likelihoods. In Table 1 and Figure 6, we also show the result of combining the fully marginalized one-dimensional M hm posteriors from our MW satellite and strong lensing analyses.

As demonstrated in the right panel of Figure 6, the joint marginal likelihoods for M hm become increasingly constraining as q increases. This is due to the fact that the transformed Σ sub -M hm posterior distribution from MW satellites (Figure 3 left panel) breaks the degeneracy between these parameters present in the strong lensing posterior. In particular, larger values of q correspond to more efficient subhalo disruption in strong lens host halos relative to the MW and yield lower inferred values of Σ sub at the strong lensing scale according to Equation (14). This shifts the region of two-dimensional parameter space in which we multiply the MW satellite and strong lensing posteriors toward lower values of Σ sub . Thus, because the low-Σ sub region of the lensing posterior does not allow for large values of M hm , larger values of q yield more stringent joint M hm constraints (and vice versa for smaller values of q). Indeed, as shown in the right panel of Figure 5, restricting the strong lensing posterior to the range of Σ sub inferred from our MW satellite analysis for q = 1 significantly strengthens the M hm constraint set by lensing alone.

Despite the qualitative effects of varying the differential subhalo disruption efficiency described above, varying q within a reasonably broad range only impacts the results of our probe combination at the ∼10% level in terms of m WDM . As discussed in Section 7.1, the differential efficiency of subhalo disruption due to baryons is one of several systematics that impact our probe combination at this level, all of which must be Figure 6. Left panel: the impact of systematics on the marginalized one-dimensional posterior distributions of projected subhalo number density at the strong lensing scale. The marginalized posterior distribution from our MW satellite analysis is shown in blue, the marginalized strong lensing posterior is shown in red, the dashed blue distributions indicate additional uncertainty in our MW satellite inference due to the mass of the MW halo, and the dotted-dashed green distribution illustrates the effects of systematic uncertainty in the differential efficiency of subhalo disruption due to baryons at the MW and strong lensing host halo scales. Right panel: joint marginal likelihood of WDM half-mode mass for our MW satellites plus strong lensing probe combination. Joint likelihoods are shown for equally efficient subhalo disruption in the MW and strong lens host halo mass and redshift regimes (q = 1, purple), twice as efficient disruption due to baryons in the MW relative to strong lens halos (q = 0.5, dotted-dashed green), and twice as efficient disruption in strong lens halos relative to the MW (q = 2, dashed green). The gray distribution shows the result of combining the fully marginalized one-dimensional M hm posteriors derived from strong lensing and MW satellites.

controlled in a joint modeling framework to claim a detection of non-CDM physics at the corresponding level of precision. Figure 6 and Table 1 demonstrate that combining the MW satellite and strong lensing posteriors with any value of q-i.e., performing the combination in multiple dimensions-is more constraining than combining the fully marginalized M hm posteriors, as expected.

Projected Subhalo Number Density Constraints

As demonstrated in the left panel of Figure 5, the marginalized posterior for Σ sub from strong lensing accommodates significantly larger values than we infer from our MW satellite analysis. In particular, our fiducial joint analysis yields a marginalized posterior distribution from MW satellites that peaks at Σ sub ≈ 0.025 kpc -2 ; moreover, Σ sub > 0.03 kpc -2 is not sampled because these projected SHMF amplitudes are larger than the average of our MW-like simulations.foot_3 Meanwhile, Σ sub values in this range are disfavored in the lensing posterior relative to its mild peak at Σ sub ≈ 0.067 kpc -2 by a ratio of ∼2:1. Although this is not a significant tension, it is worth exploring in future work that places Σ sub constraints at various host mass and redshift scales in the context of expectations from cosmological simulations. For example, Lazar et al. (2021) identify potentially significant contributions from backsplash halo populations near strong lenses beyond those captured by the two-halo term used in Gilman et al. (2020a), which (if modeled) may lower the inferred range of Σ sub and strengthen the corresponding WDM constraints. Furthermore, there are potential differences between the surviving subhalo populations inferred from our MW satellite and strong lensing analyses caused by tidal stripping, although heavily stripped halos do not dominate the signal in either case. Thus, although it is unlikely because the subhalos that contribute to strong lensing flux ratio statistics are usually tidally truncated well outside of their NFW scale radius (Gilman et al. 2020a;Minor et al. 2020), a careful analysis of whether these systems can be stripped severely enough such that their luminous content is affected warrants detailed investigation in future work.

In the left panel of Figure 6, we show how the Σ sub posterior from our MW satellite analysis shifts as a function of both the differential subhalo disruption efficiency due to baryons, q, and the MW halo mass, where we use the MW host halo mass uncertainties discussed in Section 3 and assume that Σ sub ∝ M MW . The Σ sub distribution from MW satellites is clearly sensitive to both of these systematic uncertainties, which we discuss further in Section 7. Because varying q changes the inferred Σ sub distribution in the strong lens halo mass and redshift regime, this quantity can potentially be constrained as the precision of Σ sub constraints from strong lensing increases. Although we do not attempt to constrain q here, we note that our results disfavor simultaneously high MW halo mass and low subhalo disruption efficiency due to baryons in the MW relative to strong lens host halos, which is physically reasonable.

Systematics

The analysis presented above casts MW satellite and strong lensing constraints in a shared, multidimensional parameter space for the first time. We emphasize that our WDM constraints (Section 6.2) are conservative due to our broad priors on key systematics and are robust to the modeling uncertainties directly addressed in the joint analysis at the ∼10% level (Section 6.3). Our work therefore provides important foundations for more detailed modeling frameworks that simultaneously constrain MW satellite and lensing observables at the likelihood level.

Nevertheless, our analysis makes several simplifying assumptions that circumvent a joint likelihood analysis. We regard these as crucial areas for future model-building work in preparation for next-generation facilities and surveys, both for the MW satellite-strong lensing probe combination presented here and to further combine these probes with analyses of stellar stream perturbations, the Lyα forest, and any other novel probes of small-scale structure. In general, joint small-scale structure constraints may be sensitive to additional "nuisance parameters" distinct from those governing non-CDM physics, which must be simultaneously measured to robustly claim evidence for a deviation from CDM. This underscores the importance of our multidimensional approach and of the following systematics, which we plan to build a joint model to simultaneously infer in future work.

SHMF slope. We assume a particular value of the SHMF slope α when constructing the -Σ sub relation in Section 5, thereby taking a thin slice through this dimension of the posterior from Gilman et al. (2020a). Although current MW satellite analyses do not strongly constrain the SHMF slope, future constraints from the MW satellite population probed by LSST may be sensitive to this quantity due to excellent observational sensitivity at the faint end of the satellite luminosity function throughout the MW virial radius (e.g., Ivezić et al. 2008;Hargis et al. 2014;Drlica-Wagner et al. 2019). Meanwhile, the Gilman et al. (2020a) strong lensing analysis already mildly constrains the SHMF slope, and this sensitivity will drastically increase with larger lens samples. Exploiting all of these data will require self-consistent suites of high-resolution simulations of both MW-like systems (including realistic LMC analogs) and strong-lens-like systems, which we are currently developing. Few such high-resolution zoom-in simulations at the group-mass scale have been performed, and these are particularly valuable to validate the predictions of semianalytic models like Galacticus used to inform strong lens substructure models. These studies must be coupled with detailed models for the impact of baryonic physics on smallscale dark matter structure because it is expected to significantly affect both the amplitude and slope of SHMF at low halo masses (Benson 2020).

Line-of-sight halo mass function. We marginalized over the amplitude of the line-of-sight halo mass function in our probe combination, noting that the Gilman et al. (2020a) lensing analysis our work is based on does not constrain this quantity within a broad prior range of ±20% relative to the mean Sheth-Tormen prediction. However, detailed zoom-in simulations of strong lens analogs coupled with realizations from cosmological simulations of the line-of-sight halo populations may provide more informative theoretical priors that-combined with upcoming strong lens discoveries and follow-up imaging and spectroscopy-will yield more decisive differential measurements of the line-of-sight and substructure contributions to the lensing signal (see Lazar et al. 2021 for a recent discussion). This will ultimately allow Σ sub to be measured more precisely, breaking degeneracies with WDM physics and facilitating a more direct combination with MW satellite data.

Subhalo disruption efficiency due to baryons. We combined MW satellite and strong lensing constraints at fixed values of the differential subhalo disruption efficiency due to baryons, q. Although q does not significantly affect the joint WDM limits presented here (Section 6), this quantity represents a key systematic that must be addressed in dedicated modeling work. In particular, it will be fruitful to analyze samples of hydrodynamic simulations at the MW and group-mass scales to refine subhalo disruption models that can be applied to larger simulation suites efficiently (e.g., Nadler et al. 2018). Constructing a physically motivated model for the differential efficiency of subhalo disruption due to baryons in strong lens systems and the MW will again allow for more informative theoretical priors in joint analyses, enabling robust constraints on deviations from CDM predictions.

Milky Way and strong lens host halo properties. The mass of the MW halo remains a key systematic for the interpretation of MW satellite measurements in terms of the underlying SHMF which then propagates into joint small-scale structure constraints. The MW halo mass is a particularly important nuisance parameter for setting non-CDM constraints because the (lack of a) turnover in low-mass subhalo abundances is inferred from the SHMF corresponding to MW satellite observations, while the SHMF amplitude scales linearly with host halo mass. In our analysis, uncertainty in the MW halo mass significantly affects our M hm and m WDM constraints, and we currently take a conservative approach to marginalize over this dependence. Forthcoming Gaia data releases will increase the precision of MW halo mass measurements, and combining detailed simulation suites of MW-like halos spanning the inferred mass range with next-generation observations of the MW satellite population will allow us to derive joint constraints on the MW halo mass and SHMF (e.g., see Newton et al. 2020).

Meanwhile, strong lensing measurements are less sensitive to host halo mass uncertainty because they probe both the SHMF, small-scale structure along the line of sight, and the concentrations of low-mass halos and subhalos. Nonetheless, the details of strong lens host halo selection functions are relatively unexplored (e.g., see Sonnenfeld et al. 2015) and will be better quantified using a variety of data including weak lensing and satellite velocity dispersion measurements. These efforts will lead to more precise constraints on the masses, secondary properties, and environments of strong lens host halos, further mitigating key theoretical uncertainties in forward models of strong lensing data. Enzi et al. (2020) recently presented a joint analysis of smallscale structure probes including MW satellite galaxies and gravitational imaging, with several distinct assumptions underlying the individual and joint modeling of these probes relative to our work. Here, we discuss the most important aspects of our individual models for MW satellites and strong lensing flux ratio statistics as well as our probe combination procedure relative to the Enzi et al. (2020) study.

Comparison to Recent Studies

For MW satellites, the Nadler et al. (2021) study upon which we base our analysis explicitly includes realistic LMC analog systems in the simulations used to perform the inference. This allows Nadler et al. (2021) to use the entire population of observed MW satellite galaxies-and particularly those within and near the DES footprint-without down-weighting systems based on the probability they are associated with the LMC, strengthening our dark matter constraints relative to the Newton et al. (2018Newton et al. ( , 2020) ) (2020). Importantly, these selection functions depend on satellite galaxy size, which is a crucial driver of satellite detectability that directly informs the translation from MW satellite observations to the underlying SHMF. This highlights the importance of including a model for the relationship between subhalo and satellite galaxy size like the one used in our analysis. As discussed in Section 3, we also marginalize over MW halo mass and the efficiency of subhalo disruption due to baryonic physics, which are both key systematics in the MW satellite inference. Our MW host halo mass marginalization procedure is analytic, unlike the simulation-based method employed in Newton et al. (2020), due to the limited statistics of MW-like simulations that include realistic LMC analogs. The significant improvements in sensitivity to non-CDM physics afforded by modeling the LMC satellite system further reinforce the importance of simulation suites of MWlike systems including realistic LMC analog systems.

On the strong lensing side, the Gilman et al. (2020a) study upon which we base our analysis uses flux ratio statistics that are significantly more constraining than the gravitational imaging data underlying the Enzi et al. (2020) joint analysis. This additional constraining power results from the fact that current gravitational imaging data probes ∼10 9 M e subhalos while flux ratio anomalies are sensitive to the presence of lower-mass subhalos. In terms of modeling, Gilman et al. (2020a) explicitly account for the host mass and redshift dependence of the SHMF using Galacticus-these are leading-order effects in predicting the SHMF for a given lens and its lens-to-lens variation-while the Vegetti et al. (2018) and Ritondale et al. (2019) analyses that the joint constraints in Enzi et al. (2020) are based on do not. In addition, Gilman et al. (2020a) self-consistently account for the reduction in halo concentration in WDM, which significantly increases the sensitivity of lensing observations to WDM effects and also models the effects of tidal stripping on subhalos after infall, which is again crucial to accurately forward-model flux ratio observations.

Finally, we emphasize the following key aspects of our probe combination relative to the procedure in Enzi et al. (2020), which combines fully marginalized one-dimensional M hm distributions from various small-scale structure probes including MW satellites and gravitational imaging to derive joint WDM constraints:

1. We cast the subhalo populations inferred from MW satellites and from the group-mass, z ∼ 0.5 host halos probed by strong lensing into a common, multidimensional parameter space of projected subhalo number density Σ sub versus WDM half-mode mass M hm ; 2. We combine these Σ sub -M hm distributions to construct a joint marginal likelihood that is strictly more constraining and informative than the joint M hm distribution resulting from fully marginalizing over all additional parameters (see the right panel of Figure 6), improving the precision of our joint analysis; and 3. We model the differential efficiency of subhalo disruption due to the central galaxies in the different host halo mass and redshift regimes probed by MW satellites and strong lensing, finding that our results are robust to uncertainties in these effects at the ∼10% level, which lends confidence to the robustness of our results.

The differences in the underlying data used in our inference -and particularly the inclusion of LMC-associated satellites in the MW satellite analysis and the use of strong lensing flux ratio statistics that probe lower-mass subhalos than current gravitational imaging data-therefore result in more precise joint constraints than those obtained in Enzi et al. (2020) and allow us to significantly improve upon their WDM limit. Moreover, the joint analysis choices described above lend to the robustness and accuracy of our results.

Conclusions

In this paper, we performed a multidimensional joint analysis of the distribution of small-scale dark matter structure inferred from MW satellite galaxies and strong gravitational lensing. In particular, we combined state-of-the-art dark matter substructure measurements derived from (i) the MW satellite galaxy population over ∼75% of the sky and (ii) the flux ratio statistics and image positions from eight quadruply imaged quasars. By combining constraints on the projected subhalo number density and the half-mode mass describing the suppression of the subhalo mass function in thermal relic WDM, we improved lower limits on the WDM particle mass derived independently by breaking degeneracies among the inferred subhalo distributions at each scale for the first time. Our m WDM constraint is more stringent than any limit set by independent analyses of small-scale structure probes to date.

Our key results are summarized below:

1. Our multidimensional joint analysis extracts information that was not accessed by MW satellite or strong lensing analyses independently, improving WDM constraints by ∼30%, with M hm < 10 7.0 M e (m WDM > 9.7 keV) at 95% confidence, or M hm = 10 7.4 M e (m WDM = 7.4 keV) disfavored with a 20:1 marginal likelihood ratio. (Figures 45); 2. Our joint WDM constraint is robust to uncertainties in the differential efficiency of subhalo disruption between the MW and strong lens host halo mass and redshift regimes at the ∼10% level; 3. Projected subhalo number density constraints from MW satellites and strong lensing flux ratio statistics are in mild tension but are sensitive to uncertainties in the efficiency of subhalo disruption in the corresponding host halo mass and redshift regimes; 4. We discuss key systematics that are conservatively marginalized over in the current analysis but which must be mitigated in future work to claim a detection of non-CDM physics from small-scale structure measurements. These systematics include the line-of-sight contribution to the strong lensing signal, the differential efficiency of subhalo disruption due to baryons at the MW and lensing host halo mass and redshift scales, and the properties of the MW and strong lens host halos (Figure 6); 5. Inferences of the small-scale dark matter structure from MW satellites and strong lensing are consistent despite the completely different nature of these probes and differences in their corresponding host halo mass and redshift regimes.

Recent studies have identified a variety of microphysical dark matter properties that suppress small-scale structure in a manner quantitatively similar to WDM, including the strength of velocity-independent interactions between dark matter and protons (Nadler et al. 2019a), the production mechanism of nonthermal dark matter in early matter-dominated cosmologies (Miller et al. 2019), and the dark matter formation redshift in models of "late-forming" dark matter (Das & Nadler 2021). Our jointly derived WDM constraints directly inform all of these properties. Dark matter models that feature qualitatively different suppression of small-scale structure compared to WDM can also be constrained by constructing a conservative mapping; for example, such mappings have been applied to constrain fuzzy dark matter (Schutz 2020), models with velocity-dependent self-and Standard Model dark matter interactions (Tulin & Yu 2018;Maamari et al. 2021), and models within the ETHOS framework (Bohr et al. 2020). Such dark matter physics may manifest differently in small-scale structure probes like MW satellites and strong lensing that are sensitive to halo abundances and concentrations in unique ways, and we regard this as a particularly compelling avenue for future work.

We expect the relative improvement offered by our probe combination to continue to increase as both techniques progress due to both additional data from existing instruments and nextgeneration observational facilities. Excitingly, the sample sizes of both nearby ultrafaint dwarf galaxies and quadruply lensed quasars are expected to drastically increase with LSST (Ivezić et al. 2008), Euclid Space Telescope (Laureijs et al. 2011), and Nancy Grace Roman Space Telescope (Spergel et al. 2015) observations. Forthcoming facilities including the Maunakea Spectroscopic Explorer (The MSE Science Team et al. 2019) will also help to confirm the nature of candidate MW satellites and faint dwarf galaxies throughout the Local Volume, while wide-aperture and extremely large telescopes (ELTs) will provide detailed information about the dynamical masses of these systems, which is key to refine galaxy-halo connection and WDM constraints (Simon et al. 2019). Meanwhile, the unprecedented sample of strong lenses expected to be discovered within 5-10 yr will yield precise measurements of the differential line-of-sight and substructure contributions to lensing signal and will allow the selection functions of strong lenses to be better quantified. Observations of extended source emission will also help constrain lens macromodels.

With sufficiently stringent limits on the minimum luminous halo mass from nearby dwarf galaxies and on the mass scale of a cutoff in the subhalo mass function, our procedure for combining satellite galaxy and strong lensing posteriors can potentially provide evidence for the existence of dark subhalos -i.e., subhalos devoid of observable baryonic componentswhich are a key, unverified prediction of many viable dark matter models. Nadler et al. (2020) estimate that the lowestmass halo expected to host a dwarf galaxy is more massive than ∼10 7 M e . Thus, the future observations discussed above, which are expected to constrain the subhalo mass function at and below these mass scales, will either yield evidence for a cutoff in galaxy (or halo) formation or evidence for halos devoid of observable baryonic matter. We plan to pursue these measurements by developing a multipronged theoretical framework to jointly infer the distribution of small-scale structure using heterogeneous data. zoom-in simulations for SHMFs evaluated using present-day subhalo virial mass. This indicates that the amount of stripping experienced by subhalos in our N-body simulations is well captured by the Galacticus model, on average. However, as shown in the right panel of Figure 7, our zoom-in simulations yield radial subhalo distributions that are slightly less concentrated than those predicted by Galacticus. This discrepancy is unchanged when comparing to the higher-resolution version of one of our simulations described above. We note that increasing the radial concentration of the subhalo distribution predicted by our simulations at fixed subhalo abundance would further strengthen our minimum halo mass and WDM constraints (Nadler et al. 2020). Because the radial distribution in our simulations and Galacticus are respectively subject to subtle numerical uncertainties including artificial subhalo disruption and semianalytic modeling of dynamical friction, we plan to explore this discrepancy systematically in future work.

Appendix C

Host-to-host, Poisson, and Projection Scatter in Σ sub,MW As we have emphasized, it is challenging to accurately estimate the host-to-host scatter in Σ sub,MW given the requirements we place on our realistic MW-like simulations. In this appendix, we take a very conservative approach by quantifying the scatter in Σ sub,MW for the entire suite of simulations from Mao et al. (2015) described in Appendix A. In particular, the left panel of Figure 8 shows the projected SHMF for all 45 of the Mao et al. (2015) simulations, shaded by their Σ sub,MW calculated according to the procedure in Section 5.2, and the right panel of Figure 8 shows the dependence of Σ sub,MW on host halo mass and concentration. These panels illustrate that the scatter in Σ sub,MW is at most ∼40% toward smaller values of Σ sub,MW than inferred from our realistic MW-like simulations, and at most ∼20% toward larger values of Σ sub,MW . The dependence of subhalo abundance on host halo properties among these zoom-in simulations is studied in detail by Mao et al. (2015) and Fielder et al. (2019).

Although the host-to-host uncertainty quoted above is not small in an absolute sense, the scatter in either direction is overshadowed by the factor of 2 uncertainty introduced by q in the translation from Σ sub,MW to Σ sub . Moreover, scatter toward lower values of Σ sub,MW (which is more common) would further strengthen our joint WDM constraints as described in Section 6.3. Furthermore, this estimate of the host-to-host scatter using the entire zoom-in simulation suite is an overestimate because it does not leverage additional information about the properties of the MW halo and because the Mao et al. (2015) hosts were chosen to span a cosmologically representative range of formation histories rather than being selected uniformly in host halo mass. We therefore regard our current analysis to be conservative because it accounts for the dominant uncertainties (i.e., M MW and q), and we plan to simultaneously infer q, M MW , and Σ sub,MW along with their associated uncertainties in future work.

The Poisson scatter in the projected SHMFs predicted by our simulations near the minimum halo mass is also moderate compared to the other systematic uncertainties discussed above. In particular, given our fiducial binning scheme, there are ∼50 subhalos per M peak bin near , corresponding to ∼15% Poisson scatter, which is again relatively minor compared to uncertainties in q and M MW . We refer the reader to Mao et al. (2015) for a detailed study of these subhalo populations that justifies the use of a Poisson distribution to describe their scatter.

Finally, we note that the scatter in the projected subhalo mass function induced by different orientations for the projection of the MW subhalo population is also small compared to the other sources of uncertainty we have discussed. For example, the subhalo population projected with half of the virial radius in our MW-like simulations varies at the percent level for different orientations. Dashed vertical lines approximately mark the radial range of observed MW satellite galaxies used in our analysis. In both panels, dark (light) red contours show 68% (95%) confidence intervals from Galacticus for a sample of halos with characteristics matched to our MW-like simulations (see Appendix B for details). To calculate the Galacticus radial distributions, we only consider halos with M peak > 10 8 M e in addition to the V peak and V max cuts described in Appendix A to facilitate a direct comparison to our simulation results.

Appendix D Subhalo Disruption Efficiency Prior

To formulate our probe combination in a statistically consistent way, we reran the Nadler et al. (2021) MW satellite analysis with a uniform prior on the baryonic subhalo disruption  as described in Section 3. However, the fiducial Nadler et al. (2021) model assumes a lognormal prior on this quantity centered around =  1 (i.e., the expectation for the efficiency of subhalo disruption from hydrodynamic simulations of MW-mass halos; also see Nadler et al. 2019bNadler et al. , 2020)). In this appendix, we explore the effects of performing the probe combination using this lognormal prior.

In particular, Figure 9 shows the posterior from the MW satellite analysis in the -M hm and Σ sub -M hm parameter spaces, translated according to Equation (3) with q = 1, assuming the fiducial Nadler et al. (2021) prior of ( ) m s ~= =   ln 1, 0.5 . It is visually evident that this prior favors a narrower range of Σ sub , as expected. Using this alternative prior and setting q = 1 does not change the results of our joint analysis, with M hm < 10 7.0 M e (m WDM > 9.7 keV) at 95% confidence and M hm = 10 7.5 M e (m WDM = 6.9 keV) disfavored with a 20:1 marginal likelihood ratio. This is due to a cancellation of effects: using a lognormal prior on  slightly strengthens our MW satellite constraint on M hm (compare the left panel of Figure 9 to 1), but also removes the low-Σ sub tail of the MW satellite posterior (compare the right panel of Figure 9 to the left panel of Figure 4). Because larger values of Σ sub lead to weaker joint constraints as described in Section 6.3, these effects push our joint WDM constraints in opposite directions and happen to be roughly equal in magnitude. ) anywhere within their virial radius, and circles show simulations that do not have an LMC analog. Colors indicate fractional differences relative to the average value of Σ sub,MW from the two MW-like simulations.