Lensing substructure quantification in RXJ1131-1231: A 2 keV lower bound on dark matter thermal relic mass

TL;DR

This study uses realistic simulations and Bayesian analysis of gravitational lens RXJ1131-1231 to statistically constrain dark matter particle mass, ruling out warm dark matter below 2 keV at 2σ confidence.

Contribution

It introduces a forward modelling approach combining semi-analytic substructure generation with ABC to constrain dark matter properties from lensing data.

Findings

Ruling out warm dark matter below 2 keV at 2σ confidence.

Developing a simulation-based framework for dark matter substructure analysis.

Applying the method to HST data of RXJ1131-1231.

Abstract

We study the substructure content of the strong gravitational lens RXJ1131-1231 through a forward modelling approach that relies on generating an extensive suite of realistic simulations. We use a semi-analytic merger tree prescription that allows us to stochastically generate substructure populations whose properties depend on the dark matter particle mass. These synthetic halos are then used as lenses to produce realistic mock images that have the same features, e.g. luminous arcs, quasar positions, instrumental noise and PSF, as the data. We then analyse the data and the simulations in the same way with summary statistics that are sensitive to the signal being targeted and are able to constrain models of dark matter statistically using Approximate Bayesian Computing (ABC) techniques. In this work, we focus on the thermal relic mass estimate and fix the semi-analytic descriptions of the substructure evolution based on recent literature. We are able, based on the HST data for RXJ1131-1231, to rule out a warm dark matter thermal relic mass below 2 keV at the 2$\sigma$ confidence level.

Equations

$$r_T(k)={\left[P(k{)}_{\text{ΛWDM}}/P(k{)}_{\text{ΛCDM}}\right]}^{1/2}.$$

$$r_T(k)={\left[1+{\left(\alpha k\right)}^{2\nu }\right]}^{-5/\nu },$$

$$\alpha =0.049{\left(\frac{m_{\chi }}{1\text{keV}}\right)}^{-1.11}{\left(\frac{{\Omega }_{\chi }}{0.25}\right)}^{0.11}{\left(\frac{h}{0.7}\right)}^{1.22}{h}^{-1}\text{Mpc}.$$

$${\sigma }^2(M)=\frac{1}{2{\pi }^2}{\int }_0^{\infty }4\pi k^{2}P(k)W^2(k|M)dk,$$

$$W(k|M)=\{\begin{array}{ll}1& \text{if }k\le k_s(M)\\ 0& \text{if }k>k_s(M).\end{array}$$

$$M=\frac{4\pi }{3}\bar{\rho }{R}^3=\frac{4\pi }{3}\bar{\rho }{\left(\frac{a}{k_s}\right)}^3.$$

$$c=A_c{\left(\frac{M}{M_{\text{pivot}}}\right)}^{B_c}(1+z{)}^{C_c}.$$

$$\dot{m}=-A_{ts}\frac{m}{{\tau }_{\text{dyn}}}{\left(\frac{m}{M}\right)}^{\zeta },$$

$$n_{\text{dyn}}\equiv \frac{{\tau }_{\text{df}}}{{\tau }_{\text{dyn}}}=A_{df}\frac{(M/m{)}^{b_{df}}}{\ln (1+M/m)}\exp [c_{df}\eta ]$$

$$P(\eta )\propto {\eta }^{1.22}(1-\eta {)}^{1.22}.$$

$$a_{\text{df}}=a_{\text{infall}}\exp \left[n_{\text{dyn}}/10\right].$$

$$y(x)=\frac{2^{\alpha }{x}^{\beta }}{(1+x{)}^{\alpha }}$$

$${\vec{\alpha }}_{dm}^s={\vec{\alpha }}_{dm}\ast g({\varphi }_{\text{min}}).$$

$$\Delta {\vec{\alpha }}_{dm}={\vec{\alpha }}_{dm}-{\vec{\alpha }}_{dm}^s.$$

$${\vec{\alpha }}_{\text{mock}}={\vec{\alpha }}_{\text{smooth}}+\Delta {\vec{\alpha }}_{dm}.$$

$${\vec{\alpha }}_{\text{pert}}(x_i,y_i)={\vec{\alpha }}_{\text{smooth}}+{\vec{\alpha }}_{\text{clump}}(x_i,y_i).$$

$${\alpha }_{\text{clump}}(r)=\{\begin{array}{ll}{\theta }_{\text{E,clump}}& \text{if }r\le r_{\text{trunc}}\\ {\theta }_{\text{E,clump}}\left(2-\frac{r}{r_{\text{trunc}}}\right)& \text{if }{r}_{\text{trunc}}<r\le 2r_{\text{trunc}}\\ 0& \text{if }r>2r_{\text{trunc}}\end{array}$$

$$(x_i^{\prime },y_i^{\prime })=(x_i,y_i)-{\vec{\alpha }}_{\text{smooth}}(x_i,y_i),$$

$$R_i^{\text{pert}}=\sum _{\text{pixel}}\frac{1}{{\sigma }^2}{\left(I-I_i^{\text{pert}}\right)}^2,$$

$$R_i^{\text{smooth}}=\sum _{\text{pixel}}\frac{1}{{\sigma }^2}{\left(I-I_i^{\text{smooth}}\right)}^2,$$

$$R_i^{\text{sens}}=\sum _{\text{pixel}}\frac{1}{{\sigma }^2}{\left(I_i^{\text{pert}}-I_i^{\text{smooth}}\right)}^2,$$

$$\Delta R_i=R_i^{\text{pert}}-R_i^{\text{smooth}}.$$

$$D(C_1,C_2)={\int }_{0"}^{2.5"}{\left(C_1(dr)-C_2(dr)\right)}^{2}dr,$$

Full text

1 Introduction

Dark matter is one of the core ingredients in cosmology and evidence for its existence can be found in many areas including the cosmic microwave background (CMB) temperature anisotropies [1], galaxy clustering [2], weak gravitational lensing [3, 4], supernovae of type Ia [e.g. 5, 6] and clusters [7]. The smallest scale phenomenological tests come from the Lyman-alpha forest [see e.g., 8, 9]. The standard paradigm is that dark matter consists of cold particles that interact only through gravity. On large cosmological scales this model makes predictions that remain consistent with data. However on small, sub-galactic, scales tensions between theoretical expectation and observational measurements seem to appear. Discrepancies have been reported predominantly in the number, phase space densities and density profiles when comparing simulations of dark matter substructure with observations of luminous satellite galaxies in our Milky Way (MW) see e.g. [10, 11, 12, 13, 14] and a recent review [15]. There also remain open questions about the inner slopes in galaxy clusters [e.g. 16]. These tensions attract attention due to the fact that some of these features could be explained by introducing dark matter particle properties and interactions beyond those of cold dark matter (CDM) of the standard cosmological model.

Warm dark matter (WDM), a thermal relic with a free streaming length, has been proposed to solve the discrepancies in the abundance of small scale structure [17, 18]. The empirical constraints on this model, from the Lyman-$\alpha$ forest [19], lead to a lower bound on the thermal relic mass $m_{\text{TH}}~{}=3.3$keV at the 2$\sigma$ confidence level. The latest constraints from [20] results to $m_{\text{TH}}~{}=5.3$keV based on their chosen priors on the thermal evolution. Allowing for a non-smooth evolution of the temperature of the IGM with sudden temperature changes reduces their lower limit to $m_{\text{TH}}~{}=3.5$keV. The limits obtained from dwarf galaxy counts also disfavors particle masses below $m_{\text{TH}}~{}=2.3$keV at the 2$\sigma$ confidence level [21, 22]. Another important probe of substructure properties is strong gravitational lensing [23, 24, 25]. For example, strong lensing has been used to detect luminous and dark substructure in strong lens systems [26, 27]. This has led to the detection of individual clumps down to masses of roughly $2\times 10^{8}M_{\odot}$ [28, 29, 30]. In these studies the approach is to detect substructure on an object by object basis. However, to make precision measurements of WDM matter properties, we need to move to statistical methods that are able to capture information from the collective action of the population of subhalos within a parent halo. Statistical approaches to quantify substructure have been done in simple lens configurations, for instance by analyzing flux ratios in multiple imaged lensed quasars [see e.g., 31, 32, 33, 34, 35]. Anomalous flux ratios have been reported relative to a prediction based on simple smooth lens models, known as flux-ratio anomalies. Constraints on warm dark matter from weak lensing in anomalous quadruple lenses [36] lead to $m_{\text{TH}}~{}=1.3$keV under the assumption that the lens profile is perfectly known. These approaches requires multiple strong lens systems as the statistics is rather weak for a single lens system as well as strong assumptions about the lens mass profile.

In this paper, we take a different statistical approach to quantifying substructure in strong gravitational lenses. Our method is based on an extensive forward modelling scheme that relies on simulating lensing imaging data using a suite of halos, where the statistics of the subhalo population are set by the dark matter properties. These generated mocks are tuned to be similar to the data [see e.g. 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, for other application of forward modelling in cosmology]. A simple feature of a forward modelling approach is that the same analysis tools are run on both simulations and the data. This then allows us to use an approach known as Approximate Bayesian Computing (ABC) [50, 51, 52, 53, 54] to perform Bayesian inference calculation even for cases where the likelihood cannot be computed. This is because the variability in the mocks allows us to propagate errors to the uncertainties on the parameters being estimated. A key step to this approach is to find diagnostic measures, known as summary statistics, that are sensitive to the signal we wish to target, which is substructure in our case. A recent ABC framework for cosmological purpose has been implemented by [55]. We apply our method to the quadrupole lens system RXJ1131-1231 , that was discovered by [56], to quantify the substructure content and to differentiate between different dark matter models.

This paper is structured as follows: In section 2, we describe the smooth lens model and source reconstruction of the lens RXJ1131-1231 . In section 3, we describe our substructure modeling and how we create realistic mock lenses for different dark matter models. We then describe how we perform a substructure analysis and how we compare the statistical features in section 4. In section 5, we present our dark matter model constraints based on the analysis of the lens RXJ1131-1231 based on hundreds of different simulations. We discuss the results in section 6, compare them to the literature and discuss possible extensions to this work. Further technical details about the analysis are provided in the appendices. Throughout this work, we assume a flat $\Lambda\text{CDM}$ cosmological model with $H_{0}=69.31$ km s*-1Mpc-1*, $\Omega_{b}=0.049$, $\Omega_{m}=0.315$, $\sigma_{8}=0.829$ and $n_{s}=0.968$.

2 Smooth lens model of RXJ1131-1231

In this work, we focus on substructure inferences based on the strong lens system RXJ1131-1231 . The redshift of the lens is $z_{l}=0.295$ and of the background quasar source $z_{s}=0.658$. These were determined spectroscopically by [56]. Throughout this work, we make use of the MultiDrizzle product from the HST archive in bands ACS WFC F814W and F555W. We use a 1602 pixel image centered at the lens position with pixel scale 0.05”. This corresponds to a FOV of 8”$\times$ 8”.

In targeting the signatures of substructure we separate our analysis into global features and small scale structure, where we focus on the latter to measure dark matter properties. Furthermore, two components need to be set for the global model, the lens and the source. For the smooth lens model, we use the same model parameters as in [57]. This model is composed of an ellipsoidal power law mass distribution, an SIS profile centered on the visible lens substructure and external shear parameters. Additionally, we use shapelet potential perturbations, as introduced in [58] based on the shapelet basis set [59]. We set the scale parameter of the lens model shapelets to be the Einstein radius of the lens. In total, we use 21 shapelet parameters (corresponding to $n_{\text{max}}=5$), which enables us to model the lens mass distribution down to 0.4” resolution.

For the global luminous source structure, we use shapelets up to order $n_{\text{max}}=50$ with scale parameter $\beta=0.18"$ centered at the position of the quasar, as described in [58]. In a second step, we increase the source resolution with adaptive shapelets. We do this by identifying multiple distinct features in the residual map that correspond to source structure on smaller scales. We traced back the residual features to the source plane and identify 15 distinct regions corresponding to flux coming from small scale structure. For each of these specific regions, we place additional shapelets with order $n_{\text{max, clump}}=3$ and scale $\beta_{\text{clump}}$ set to the magnified resolution limit of the data. The source reconstruction is performed independently for the two imaging bands, but the common configuration and order of the basis sets is used for the source.

We reconstruct the two HST ACS filter bands F814W and F555W simultaneously and combine the likelihood of the smooth lens model given the imaging data. We use the framework presented in [58] to fit the high dimensional non-linear parameter space using a Particle-Swarm Optimization (PSO) method [60]. The best fit lens model result has a reduced $\chi^{2}$ value of $\chi^{2}_{\text{red}}=1.03$.

In modeling the global feature that do not come from substructure, care needs to be taken in two respects. The first is that the smooth lens model needs to be sufficiently flexible to capture the scales that are larger than the scales where dark matter substructure has an impact. In appendix A.1, we show how the choice of the smooth model can affect the substructure analysis. The second consideration is that the substructure signal depends on the source surface brightness variation. This means that we need to ensure that the simulations have the right level of small scale variation in the source to give reliable statistics in the mocks. We also need to allow the model to have sufficient small scale features as present in the data so as to not lead to biases. In appendix A.2, we present illustrative examples and a discussion of how these model biases affect substructure signatures.

3 Simulations

We anchor our comparison of the statistical features in the image residuals with the features of different dark matter models on simulations. For our analysis, we target deflection angle perturbations arising from substructure within the lens relative to a smooth lens model. To investigate the statistical significance of a signal, we use multiple realizations of the same physical model (WDM, CDM) and different halo masses. We use semi-analytic descriptions to generate a large number of realizations of a given model. The details on how we produce the set of simulations and how we incorporate WDM in the semi-analytic description is given in this section.

3.1 Dark matter substructure model

To compute a sample of expected deflection perturbations, we use a model based on the extended Press-Schechter [61, 62] formalism. In particular, we use a merger tree and a subhalo evolution and disruption prescriptions, tuned to N-body CDM simulations [see e.g. 63, 64, 65, 66, 67].

3.1.1 Power spectrum

We follow [8] in describing a power spectrum for WDM. The relative transfer function $r_{T}(k)$, defined as the ratio of WDM to CDM model can be defined as

[TABLE]

In the case of WDM, the relative transfer function is given by [17]

[TABLE]

where $\alpha$ is the scale of the free streaming break of the WDM particle and $\nu$ was fixed by [8] to $\nu=1.12$ and

[TABLE]

$m_{\chi}$ is the DM thermal relic mass and $\Omega_{\chi}$ is the normalized dark matter density. The mass variance $\sigma^{2}(M)$ can be computed as

[TABLE]

which depends on the window function $W$ used.

A sharp-$k$ filter can capture the impact of a truncated power spectrum (such as WDM) since the suppressed high $k$ modes leads to a flattening of $\sigma^{2}(M)$. We use the sharp-$k$ filter

[TABLE]

A definition of the mass $M$ according to the filter scale $k_{s}$ and a clear physical relation on the $k_{s}(M)$ relation is missing. We follow [68] and relate $k_{s}$ to the mass radius $R$ of a real space top-hat filter through a scale parameter $a$ with $k_{s}=a/R$ as

[TABLE]

We set $a=2.5$ as stated by [68] to predict the turnover in the halo mass function.

3.1.2 Merger tree

To generate stochastic merger trees based on the power spectrum and mass variance, we follow [68], which is based on the original merger tree for CDM structure of [64]. Using a sharp-$k$ filter to compute $\sigma^{2}(M)$ will change the halo mass function at high mass. To avoid over-predicting high mass halos, the barrier for collapse has to be increased by a factor of 1.197 [see 68]. Furthermore, WDM changes the collapse threshold below a characteristic mass [69]. Based on work of [69, 68], a threshold of collapse, $\delta_{c}(M,t)$, from WDM and CDM can be fit. For this, we adopt equation 7-10 of [69]. We do not, for simplicity, incorporate the non-Gaussian walks of a moving barrier in WDM, as proposed by [68]. However, we set the resolution limit of the merger tree to the characteristic mass. For the subhalo population, we do not expect a significant impact of the simplified implementation. This description allows us to produce many realizations of halos and the in-fall histories of the accreted sub-halos.

3.1.3 Mass-concentration relation

The concentration of the dark matter halos depends on the mass assembly history. Halos that assemble earlier are more concentrated [e.g. 70, 71, 72]. Taking a mass-concentration-redshift relation based on CDM simulations may not reflect the fact that WDM halos form later and are therefore expected to have lower concentrations. We use the $c(M,z)$ relation of [73] for CDM cosmologies, namely

[TABLE]

Throughout this work, we set $A_{c}=5.22$, $B_{c}=-0.072$, $C_{c}=-0.42$ and $M_{\text{pivot}}=2\cdot 10^{12}M_{\odot}$, taken from [73]. In the WDM case, we set the $c(M,z)$ relation by a mapping of the same formation histories of CDM halos as described in [74]. Halos of different WDM models with the same formation history, but different masses, have the same concentration. This leads to a lower concentration in a WDM case compared to CDM as a result of the later formation of the WDM halo.

3.1.4 Substructure evolution

To describe the evolution of the subhalos from the time they are assigned to a more massive halo, the information in a merger tree is insufficient since subhalos get disrupted by tidal stripping and dynamical friction.

We use the description of [75], which is based on [76], to track the average mass-loss rate, $\dot{m}$, of a dark matter subhalo expected by tidal stripping. The parameterized form of the subhalo decay $\dot{m}$ depends on the instantaneous mass ratio $m/M$ of the subhalo $m$ and the parent halo $M$ and the dynamical time $\tau_{\text{dyn}}$ of the halo as

[TABLE]

where $A_{ts}$ and $\zeta$ describe the normalization and mass dependents of the subhalo decay. Throughout this work, we use $\tau_{\text{dyn}}=0.1H^{-1}$. We take the parameter values derived from [77], $A_{ts}=1.54$ and $\zeta=0.07$. We do not introduce scatter in this relation [as done by 75] as we expect that the scatter introduced by the merger tree will have this effect.

Dynamical friction has also a significant impact on the subhalo disruption. We use the description of [78], which provides a fitting formula for this process based on simulations. The number of dynamical time scales to fully merge the subhalo into the parent $n_{\text{dyn}}$ is given by [78] as

[TABLE]

where $\eta=j/j_{c}(E)$ is the orbital circularity. We ignore the dependence on orbital energy. The fitting parameters are tuned to simulations and given by [78] as $A_{df}=0.216$, $b_{df}=1.3$ and $c_{df}=1.9$. [79] showed that the circularity $\eta$ of subhalos at infall can be well fit by

[TABLE]

We draw from this distribution for every subhalo at in-fall to compute the dynamical friction time according to equation 3.9. The scale factor of disruption due to dynamical friction $a_{\text{df}}$ is then given by

[TABLE]

Subhalos whose dynamical friction time $a_{\text{df}}$ is shorter than that needed to survive to the redshift of interest (in our case to the redshift of the lens $a_{\text{lens}}$) are assumed to be fully disrupted.

For the evolution of the subhalo internal structure, we use the description of [80, 81]. The ratio of $y=r_{\text{max}}/r_{\text{max,infall}}$ and $y=v_{\text{max}}/v_{\text{max,infall}}$ are both parameterized as a function of mass loss $x=m/m_{\text{infall}}$ as

[TABLE]

where $\alpha$ and $\beta$ are fitting coefficients. The coefficients from [80, 81] for $y=r_{\text{max}}/r_{\text{max,infall}}$ are $\alpha=-0.3$, $\beta=0.4$ and for $y=v_{\text{max}}/v_{\text{max,infall}}$ they are $\alpha=0.4$, $\beta=0.3$.

The semi-analytic prescription used to compute the disruption of halos does not have a specific concentration dependence. In the case of WDM, where the concentration of halos is expected to be lower than for CDM, ignoring the concentration as a parameter in the disruption process might slightly under-predict the strength of the disruption in the WDM case. For all of our conclusions, the bounds on the thermal relic mass of a WDM particle are conservative with respect to the subhalo disruption effect. In appendix B, we show the subhalo mass function and maximum circular velocity function for different WDM models predicted with the stated descriptions.

The spacial position of the subhalos is computed based on the orbit at infall. The orbit is computed based on the orbital circularity distribution (see equation 3.10) and provides a radial distribution function from which we draw from. Our tests show that the radial distribution of substructure mass follows the NFW profile of the parent halo.

The prescriptions and assumptions stated here allow us to generate mock halos of different masses with substructure based on different WDM thermal relic masses. Figure 1 shows a set of realizations of halos in the range 12.0-13.5 log$(M/M_{\odot})$ and thermal relic masses 1.0-10.0 keV. The enhanced amount of substructure of higher WDM mass is visible on the figure.

3.2 Mock image generation

In our approach, we aim to separate the effects of global features from those coming from substructure. In generating mock images, we have therefore decided to use a hybrid approach. For the global features, we use the deflection angles calculated by fitting the smooth lens to the data, $\vec{\alpha}_{\text{smooth}}$. This allows us to generate mock strong lens systems that are in the same regime as that of the data. The small scale features, which is where dark matter properties have the most significant impact, are modeled using the halos generated by the merger tree procedure described in section 3. To do this, we first calculate the full deflection angles for the mock halos $\vec{\alpha}_{dm}$. We then generate a smoothed field $\vec{\alpha}^{s}_{dm}$ through a convolution with a Gaussian kernel $g$ of width $\phi_{\text{min}}$,

[TABLE]

We then approximate the perturbations to the deflection angles $\Delta\vec{\alpha}_{dm}$ as

[TABLE]

We chose the kernel width $\phi_{\text{min}}$ to match the smooth lens model shapelet deflection scale, $\theta_{\text{min}}=0.4"$. The result is a deflection angle map that has no net mass beyond scales of $\theta_{\text{min}}$111There is no unique way to filter scales below $\phi_{\text{min}}$. We also performed a median wavelet transform up to the level matching $\phi_{\text{min}}$ on a subset of the halo realizations. We find no significant impact on the summary statistics we apply. We conclude that the details of the filter does not significantly impact our analysis.. As an example, Figure 2 shows the deflection angle perturbation map of one of the components of $\Delta\alpha_{dm}$ for the cases shown in Figure 1.

The deflections used to generate the mock images are therefore

[TABLE]

We also use the same source model, PSF model and weight map as inferred from the data and we add Poisson and Gaussian noise on the mock image.

4 Substructure model comparison

The tools introduced in the previous sections 3 and 2 allow us to generate mock data that is very similar to the actual imaging data. Our aim is to construct a statistical test that will allow us to select from different dark matter models. The Approximate Bayesian Computing (ABC) [53, 54, 55] framework allows us to construct a posterior distribution even for cases where the likelihood, probability of data given a model, can not be calculated. Instead, ABC relies on detailed forward modeling to capture the impact of errors. It requires distance measures (metric distances) between the output of different simulations, based on summary statistics. In general, the summary statistic can take many forms but needs to map to one or several distance metrics.

Our goal is to construct distance metrics that we can measured from the data and to the mock data and that are sensitive to dark matter properties. This means that comparisons of different realizations of the same dark matter model should result in small distances. On the other hand, comparing realizations of different dark matter models should result in large distances.

The features (substructure in the lens galaxy) that we are targeting emerge due to astrometric anomalies of surface brightness in the reconstructed image compared with a smooth lens model. As a complication, the intrinsic source surface brightness distribution is not known in the actual data. Therefore, anomalies can only be quantified by simultaneously reconstructing the lens (with and without small scale perturbations) and the source.

Instead of trying to identify individual structures within the lens, we use a scanning strategy that introduce local perturbations in the lens model as described below in 4.1. These scans provide us with information that we further compress into a summary statistics (4.2) to discriminate between different models of dark matter with ABC.

4.1 Substructure scanning procedure

Our scanning procedure works by fitting a lens model that includes a single subclump perturber, with corresponding deflection $\vec{\alpha}_{\text{clump}}$, at a position ($x_{i},y_{i}$) as well as the smooth global model $\vec{\alpha}_{\text{smooth}}$. The combined deflection angle is then

[TABLE]

We have chosen the perturber to have an SIS profile up to a radius $r_{\text{trunc}}$, followed by a linearly decreasing deflection such that beyond $2r_{\text{trunc}}$, any perturbative deflection disappears. The functional form of the perturber is

[TABLE]

where $\theta_{\text{E,clump}}$ is the Einstein radius of the clump and $r_{\text{trunc}}$ its truncation radius. For our analysis, we set $\theta_{\text{E,clump}}=0.01"$ and $r_{\text{trunc}}=0.1"$. The total mass enclosed in the truncation radius is $M_{\text{clump}}=10^{8.35}M_{\odot}$ for the lens in consideration.

We then move the position of the subclump perturber in increments of 0.05” over the image in both the horizontal and vertical directions. At each position ($x_{i},y_{i}$), we make a best fit reconstruction of the lensed image $I$ (this can be a fit to data or to simulated data) using the perturbed model ($\vec{\alpha}_{\text{pert}}$). We call this reconstructed image $I^{\text{pert}}_{i}$. We also perform fits using the smooth model $\vec{\alpha}_{\text{smooth}}$ leading to the reconstructed images $I^{\text{smooth}}_{i}$. When performing a fit we again use the source basis sets as described in section 2. In addition, we add source surface brightness shapelets of order $n_{\text{max}}=5$ at the source position

[TABLE]

which corresponds to the source plane position being mapped through the subclump perturber. We set the scale $\beta$ of these additional shapelets to the average magnification within the truncated SIS profile at ($x_{i},y_{i}$).

From the set of images that we generate, we can define three different residual measures ($R$), related to the relative agreement of the images $I$ (the data or mock generated data), $I^{\text{pert}}_{i}$ (the best fit reconstruction of $I$ with a smooth model and one additional perturber) and $I^{\text{smooth}}_{i}$ (the best fit reconstruction of $I$ with a smooth model ), which are effectively a $\chi^{2}$ measures, as

[TABLE]

[TABLE]

and

[TABLE]

where the sum is over the pixels of the image and $\sigma$ is the noise associated with each pixel. We do not expect any of those measures to agree (i.e reduced $\chi^{2}\approx 1$), i.e neither a smooth model nor a smooth model with one added truncated perturber can describe the data (or simulated data) sufficiently. We expect the sensitivity measure $R^{\text{sens}}_{i}$ to be the closest to agreement. The two models compared in $R^{\text{sens}}_{i}$ only differ in the addition of one truncated perturber in one model and the relative residual comes from a very targeted presence of one perturber. This measure provides us with an estimate of the maximal additional residual measure to be excepted in the presence of a subclump. In this regard, to quantify the deflection anomalies at particular positions, we further introduce a relative excess distance $\Delta R_{i}$ defined as

[TABLE]

Negative $\Delta R_{i}$ show that an addition of a subclump perturber improves the reconstruction of the data (or simulated data) while a positive $\Delta R_{i}$ shows that the smooth model performs better in reconstructing the original image. Figure 3 shows the scanning results, in terms of $\Delta R$ values, for the HST data (left column), a model with a halo mass of $10^{13.5}M_{\odot}$ (second column), a model with a halo mass of $10^{13.0}M_{\odot}$ (third column). The last column shows the sensitivity map, $R^{\text{sens}}_{i}$. The different rows indicate the results for different observing bands. Band F814W results are shown in the top row, F555W in the middle and the combined results F814W+F555W at the bottom.

4.2 Distance metric to compare simulations

Given the scan maps that we are able to generate for both data and mocks, the challenge is to construct a distance metric that allows us to compare two different scan maps and encapsulate the statistical features imprinted in them. In principle, any metric can be applied to the scann maps to compare two distributions. The challenge is to contract the information as much as possible not to find our self to reject too many simulations while keeping most of the contained information about the quantities of interest. The imprint of dark matter properties is effectively mapped to the abundance and density profile of subhaloes at different masses. These primary quantities of interest in the deflection pattern result in the abundances of deflection anomalies and their spacial patterns in the substructure scans.

To probe the abundances and masses of the subhalo population imprinted in the deflection anomalies, our primary statistics is the abundance and spacial pattern of negative excess distance $\Delta R_{i}<0$ (equation 4.7 described in section 4.1). In this work, we restrict ourself to the pattern emerging from $\Delta R_{i}<-4$, which is more stable to artifacts in the noise modelling of the simulations compared to the data.

We find that feasible distance measures for our aim can be constructed based the spherical averaged two-point correlation function $C(dr)=\left<\Delta R_{i}(r)\Delta R_{i}(r+dr)\right>_{r}$ of the negative residuals $\Delta R_{i}<0$. This function contains information about the number of anomalies (normalization) and their spatial patterns (slope). We chose the metric to comparing two correlation functions of two realisations of data $C_{1}(dr)$ and $C_{2}(dr)$ as

[TABLE]

where the integral goes over scales between 0” and 2.5” with equal weight on all scales. We emphasise that the expressions in equations 4.8 is not a likelihood. The ABC framework will turn a rejection based on that metric in a likelihood.

Figure 4 shows the one-dimensional correlation function $C(dr)$ of several realisations for six different dark matter models, where the panels are for different parent halo masses. The thick black line shows the statistics extracted from the data.

5 Dark Matter model constraints

In this work, we focus on thermal relic mass constraints for dark matter. From the diagnostics of section 4 and Figure 4, we see that the mass of the parent halo also has a significant impact on our distance metrics and therefore needs to be included as a free parameter in our analysis. Hence, we focus on the two parameters thermal relic mass $m_{\text{TH}}$ and parent halo mass $M_{h}$. The goal is to calculate the posterior distribution in $m_{\text{TH}}$ , marginalized over halo mass $M_{h}$. We set the prior on the halo mass as flat in log-space in the range [$10^{12}M_{\odot}$, $10^{13.5}M_{\odot}$] and for $m_{\text{TH}}$ uniform in [$1$keV, $5$keV]. For all practical purposes in our simulations and the quality of the data, $m_{\text{TH}}~{}>5$keV corresponds to CDM. To explore the parameter space in an ABC analysis, a large number of simulations are required.

After an initial grid based study with a total of 480 realisations in the range $M_{h}\in[10^{12}M_{\odot},10^{13.5}M_{\odot}]$ and $m_{\text{TH}}$ $\in[1,10]$keV, we are able to exclude any halo mass below $10^{13}M_{\odot}$ by more than $2\sigma$ confidence222This statement is sample limited and is confirmed in the more dense sampling of a narrower region in parameter space. and found that a good prior for our full ABC calculations to be $M_{h}\in[10^{13}M_{\odot},10^{13.5}M_{\odot}]$ and $m_{\text{TH}}$ $\in[1,5]$keV. We then randomly sample from this prior to produce a total of 2483 mock data sets.

In the ABC process, we accept the simulations, whose summary statistics in the cumulative distribution and the correlation function reside in the best 5% of the sample. Setting the threshold in this way does not guarantee that we have reached convergence, however, since lowering the threshold on the distance measures will give tighter constraints if our calculations have not converged, the results we present are conservative estimates. Figure 5 (left) shows the two-dimensional distribution of the accepted samples. We see that high halo masses and high thermal relic masses are preferentially selected through our cut in the summary statistics of ABC. Figure 5 (right) shows the $m_{\text{TH}}$ marginal distribution. We can state that thermal relic mass models below 2keV are disfavored by 2$\sigma$, in agreement with a CDM scenario. This result is comparable in strength and in agreement with the other probes of small scale structure. The latest constraints from the Lyman-$\alpha$ forest [19] results in a lower bound on the thermal relic mass $m_{\text{TH}}~{}=3.3$keV at the 2$\sigma$ confidence level. The limits obtained from dwarf galaxy counts dis-favors particle masses below $m_{\text{TH}}~{}=2.3$keV at the 2$\sigma$ confidence level[21, 22].

6 Discussion

We have performed a statistical analysis of the substructure content in the strong lens RXJ1131-1231 and report a lower limit for the mass of thermal relic dark matter of $m_{\text{TH}}~{}=2$keV at 2$\sigma$ confidence level. This is comparable to earlier limits derived from Lyman-$\alpha$ clustering and MW dwarf counts. Our statistical method significantly improves upon clump-by-clump strong lensing analyses and can easily be extended to multiple lens systems. Our forward modelling analysis of RXJ1131-1231 uses an ABC method and targets statistics that are sensitive measures of dark matter properties.

The model shows that multiple substructures within the same mass range have a significant impact on the lensing statistics. This means that probing substructure on a one-by-one basis at this mass limit may not be feasible. However, when pushing the limits to lower substructure masses, a statistical approach like ours is able to account for the effects of hundreds of subclumps simultaneously. For instance, this can be important for ALMA. ALMA can provide higher resolution data than HST images, which leads to greater potential for measuring smaller lens perturbations [see e.g. 30, 82, 83] if the information from these higher resolution data can be properly tapped.

Our current constraints are mainly limited by the statistics of the single lens and the moderate sample size of our simulations. The simulations and there analysis requires large computational resources. A larger sample of simulations would allow us to apply stricter cuts in the ABC framework which may lead to tighter constraints. The results we get are therefore conservative estimates of the likelihood. Better data (more lenses, better quality data) can discriminate models with higher significance. Tackling the partial degeneracy between parent halo mass and dark matter thermal relic mass can also be done by incorporating additional, independent, priors on the halo mass. For instance coming from abundance matching or galaxy-galaxy lensing. What we do have from the strong lensing measures is an accurate measure of the total mass (dark matter + baryonic mass) within the Einstein radius, which for RXJ1131-1231 is $M_{<\theta_{E}}=10^{11.9}M_{\odot}$. The host galaxy of the lens is a very massive elliptical early-type galaxy. A significant fraction of the mass within $\theta_{E}$ comes from baryonic matter. These observations allow us to set a conservative lower limit on the expected halo mass of the lens RXJ1131-1231 to be $10^{13}M_{\odot}$ from abundance matching and forward modeling of the galaxy population through cosmic time [e.g. 84, 85].

Systematics in the lens modeling could be another limiting factor. In this work, we focused on the source reconstruction scale and the intermediate lens model scale descriptions. Further effects that we do not include could mimic lensing substructure effects, such as micro-lensing by stars of the lensing galaxy, luminous structure of dwarf galaxies, substructure displaced along the line-of-sight or dust extinction. For the current constraints, we do not expect a major impact of such effects on the thermal relic mass. For lower mass halos to infer higher thermal relic mass constraints, substructure along the line-of-sight can contribute significantly to the observed substructure signatures [86, 87]. Its first order effect is a renormalization of the subclump mass function, which is effectively a renormalization of the parent halo mass in our description, or a weak lensing effect that we capture in the highly parameterized smooth lens model. Our current constraints are marginalized over the parent halo mass within our prior.

Our method can be extended to probe other statistics of lensing substructure, such as mass-concentration relations and disruption efficiencies and might provide a window to probe suggested solutions of the cusp-core discrepancies. In this work, we have fixed the semi-analytic description of the subhalo evolution based on recent literature. Covariances between the thermal relic constraints and other physics, such as tidal disruption of subhaloes have not been investigated in the current work. Systematics in the semi-analytic description might translate to systematics in the thermal relic mass constraints. We do not expect these systematics to dominate our results at the current level of precision.

Such studies are in reach with the current HST data and with a parallel advance in predictive power of the physical models. With the increasing number of discovered strong lens systems and their follow-up campaigns, significant advances can be expected from strong lens studies in the near future to shed light on the nature of dark matter and the disruption processes of satellite halos.

Acknowledgments

SB thanks Aurel Schneider for useful discussions. We acknowledge the import, partial use or inspiration of the following python packages: CosmoHammer [88], FASTELL [89], numpy 333www.numpy.org, scipy 444www.scipy.org, astropy 555www.astropy.org, matplotlib [90], FASTELL [89], triangle 666https://github.com/dfm/triangle.py.

Appendix A Systematics in the modeling

Spurious artifacts in the reconstruction modeling may lead to signatures visible in the summary statistics (see section 4) and may therefore be interpreted as caused by substructure. In this appendix, we investigate to what extent the resolution of the smooth lens model and on the source surface brightness reconstruction can impact our statistics.

A.1 Substructure deflection perturbations

We use a smooth lens model described in section 2 with multiple shapelet potentials. This basis set enables us to model the intermediate scales of the lens model. Intermediate smooth mass distribution may give rise to signatures in our summary statistics similar to what small clumps may give rise to as they give rise to astrometric and magnification anomalies as well.

To test the impact of intermediate scales on the summary statistics of section 4, we model the data with different number of lens model shapelets in addition to the elliptical power-law profile. We fit the data with 0 (no additional intermediate parameterization) 10 ($n_{\text{max}}=3$) and 21 ($n_{\text{max}}=5$) additional shapelets. Figure 6 shows the $\left(\Delta R_{i},R^{\text{sens}}_{i}\right)$ (section 4.1) distributions of those models. We see that the statistics changes significantly when adding the first 10 shapelets. This means that the most of the “signal” captured in our summary statistics in the smooth model can be attributed to features arising from scales captured by the first 10 lens shapelets. Increasing the shapelet number to 21 does not lead to a significant change in the statistics presented. This means that there are intermediate scales we have to explicitly model before applying our summary statistics. Our underlying assumption is that the features present in the scanning of the model with 21 shapelets arise from scales smaller than being captured by the lens shapelets.

A.2 Source surface brightness

The specific source reconstruction model may also have an effect on the summary statistics of section 4. Potential signal of substructure relies on the flux and its variation of the source surface brightness. Simulations require to have a source surface brightness very accurately reflecting those of the real data to give rise to the same substructure sensitivities.

To investigate this dependence, we take the lens model with 21 additional (our default model) lens shapelets and perform the scanning based on three different source reconstructions. First, we only reconstruct the source with the global shapelets with $n_{\text{max}}=50$. Second, we add the nested shapelet description at the position of the lens perturber (see section 2). Third, additionally to the nested shapelets, we add 15 fixed high resolution clumps to the source model (also described in section 2). Figure 7 shows the scanning statistics of those three reconstruction models. Substructure perturbation can change the magnification locally significantly. When the source reconstruction can not resolve the existing scales, substructure can significantly help in reconstructing the image because it can demagnify those regions such that the source reconstruction description better match the scales involved. This effect is an artifact. A substructure detection method requires to be able to describe the smallest scales involved in the source surface brightness. We tested our method by further enhancing the nested shapelets and increasing the additionally modeled source clumps. Neither of those pushes to smaller scales did significantly change the substructure scanning statistics. The underlying assumption in our inference is that we are able to match the smallest scales in the source reconstruction description relevant to match the features present in the HST images.

Appendix B Subhalo statistics

Our semi-analytic description to render a subhalo population described in section 3.1 provides the statistics of the subhalo population. Figure 8 shows the cumulative subhalo mass function and $v_{max}$ function for different WDM models for a parent halo $M_{h}=10^{12.5}M_{\odot}$ at $z=0$. The distributions are computed as the average of five realizations of each WDM model. Our results are in good agreement with the N-body simulations of [91].