SHARP – VII. New constraints on the dark matter free-streaming properties and substructure abundance from gravitationally lensed quasars

ABSTRACT

1 INTRODUCTION

Strong gravitational lensing has been shown to be a powerful tool to probe the presence of low-mass haloes in distant galactic haloes (Dalal & Kochanek 2002; Vegetti & Koopmans 2009; Vegetti, Czoske & Koopmans 2010a; Vegetti et al. 2010b, 2012; Nierenberg et al. 2014; Hezaveh et al. 2016; Bayer et al. 2018; Gilman et al. 2018), and along their lines of sight (Metcalf 2005; Despali et al. 2018; Gilman et al. 2019a). These low-mass haloes are dark matter dominated and therefore free from the uncertainty of baryonic process during structure formation. Hence, comparing to other approaches that focus on the local Universe, strong gravitational lensing provides an independent and promising approach to differentiate between alternative dark matter theories that modify the linear matter power spectrum and result in a different amount of low-mass haloes, such as cold (CDM; Springel et al. 2008), warm (WDM; Schneider et al. 2012a; Lovell et al. 2014), or fuzzy dark matter (Hui et al. 2017; Robles, Bullock & Boylan-Kolchin 2019). Mainly two approaches have been used to detect these low-mass haloes using strong gravitational lensing observations, which we now review.

The gravitational imaging technique constrains the projected position and effective mass of individual low-mass haloes via their effect on the surface brightness distribution of extended arcs (Koopmans 2005; Vegetti & Koopmans 2009). This technique has, so far, led to the detection of a few haloes in the mass regime between 10⁸ and 10⁹ M_⊙ with optical/infrared (IR; Vegetti et al. 2010a,b, 2012) and submillimetre imaging (Hezaveh et al. 2016), while future observations with long baseline interferometers are expected to lead to the discovery of haloes with masses lower than 10⁷ M_⊙ (McKean et al. 2015). Recently, Vegetti et al. (2014, 2018) and Ritondale et al. (2019) have used samples of 10–20 lenses from the SLOAN ACS Lens Survey (SLACS; Bolton et al. 2006) and the BOSS Emission Line Lens Survey (BELLS; Shu et al. 2016) to constrain the halo mass function in the regime between 10⁸ and 10¹⁰ M_⊙. They found that their results are consistent with the predictions from the CDM paradigm (e.g. Xu et al. 2015; Despali & Vegetti 2017), but more observations of higher quality are required to potentially rule out alternative WDM models.

The analysis of gravitationally lensed quasars uses the flux-ratio relation between the merging images to probe the total amount of low-mass haloes without inferring their individual positions and masses. Specifically, the presence of low-mass haloes is expected to change the relative fluxes of the multiple images compared to predictions from smooth lensing potentials. Currently, only those gravitational lens systems that produce four images of the background quasar can be used, as they provide enough constraints on the lens macro-model (typically, a smooth mass model plus an external shear). Mao & Schneider (1998) and Metcalf & Madau (2001) were the first to suggest that flux-ratio anomalies could be related to the presence of dark substructures contained within the dark matter halo of the foreground lensing galaxy, and therefore, these systems could be used to constrain the substructure fraction of distant galaxies. Follow-up studies, based on observations and numerical simulations, corroborated this idea and explicitly demonstrated the feasibility of flux-ratio anomalies as a means to detect low-mass substructure (Bradač et al. 2002; Metcalf & Zhao 2002; Dobler & Keeton 2006; Nierenberg et al. 2014).

Dalal & Kochanek (2002) presented the first homogeneous analysis of a small sample of seven lensed quasars, obtaining a result that was marginally consistent with the CDM paradigm. However, this result was contested first by Mao et al. (2004) and later by Xu et al. (2009), Xu et al. (2015), who showed that the level of the observed flux-ratio anomalies was significantly higher than expected from the subhalo population in high-resolution numerical simulations. Instead, they suggested that either line-of-sight haloes or complex mass distributions of the lensing galaxies were more likely the cause of the observed signal (see also Chen, Kravtsov & Keeton 2003; Möller, Hewett & Blain 2003; Quadri, Möller & Natarajan 2003; Metcalf 2005; Inoue & Takahashi 2012). Recently, significant progress has been made towards a better understanding of the origins of flux-ratio anomalies in these systems. In particular, Hsueh et al. (2016), Hsueh et al. (2017), Hsueh et al. (2018), and Gilman et al. (2018) have shown that baryonic structures, such as stellar discs, are a likely source of extreme flux-ratio anomalies, and it was demonstrated that deep-imaging observations were needed to break the degeneracy between low-mass haloes and other complexity in the lens mass distribution. Using mock observations, Gilman et al. (2019a) have shown the contribution from line-of-sight haloes to be significant and to provide extra constraining power on the halo mass function and the properties of dark matter (see also Metcalf 2005; Despali et al. 2018).

Since the properties of dark matter are inferred from deviations between the observed flux ratios and those predicted by the gravitational lensing mass model, reliable measurements of the flux ratios are needed. Historically, the analysis of gravitationally lensed quasars was restricted to systems with radio and mid-infrared (MIR) observations. This is because the radio jets produced from synchrotron emission and the thermal emission from the dusty torus of lensed quasars are expected to be free from dust extinction and stellar micro lensing (however, see Koopmans & de Bruyn 2000 for an exception). The small number of radio-loud lensed quasars and the difficulty in obtaining high-resolution MIR imaging from ground-based telescopes have limited the size of suitable samples. However, recent studies have shown a possible way forward to increase the sample size in the short term. In particular, Nierenberg et al. (2014) have demonstrated that narrow emission lines provide a new avenue for flux-ratio studies with near-infrared spectroscopy (see also Moustakas & Metcalf 2003), while Stacey & McKean (2018) have presented a new approach based on observations of cold dust and CO emission lines with the Atacama Large Millimetre/submillimetre Array (ALMA; see also Inoue et al. 2017). The launch of the James Webb Space Telescope (JWST) will also make the MIR flux measurements of lensed quasars faster and easier to obtain (Gardner 2012), while future large-scale surveys, such as with the Large Synoptic Survey Telescope (LSST; LSST Dark Energy Science Collaboration 2012), Euclid (Cimatti & Scaramella 2012), and the Square Kilometre Array (SKA; McKean et al. 2015), are expected to lead to the discovery of thousands of new gravitationally lensed quasars.

In this paper, we present a new analysis of the current sample of four-image gravitationally lensed quasars that have well-studied lens models and reliable radio or MIR flux measurements. For the first time, our analysis includes not only the contribution of substructure within the lensing galaxy but also that of stellar discs and line-of-sight haloes. Moreover, we use improved measurements of the observed flux ratios that have been obtained from monitoring campaigns to derive tighter constraints on the substructure and halo mass functions, and thereby the free-streaming properties of dark matter. In particular, our analysis focuses on thermally produced dark matter of which weakly interacting particles (WIMPs) are the best theoretically motivated CDM candidate. Our paper is organized as follows. In Section 2, we introduce our Bayesian modelling technique, while in Section 3, we describe the observational data used in the analysis. In Section 4, we present our new constraints on the dark matter mass function and discuss the implications for dark matter physics. In the first part of our analysis, we present the impact of line-of-sight haloes on the substructure mass function inference, assuming a concordant ΛCDM cosmology. The second part of our analysis provides constraints on the mass of a thermal relic dark matter particle. In Section 5, we summarize our results and discuss future extensions to this work.

Throughout out, we assume a flat cosmology with Ω_M = 0.28 and H₀ = 70 km s⁻¹ Mpc⁻¹.

2 METHODOLOGY

In this section, we describe the Bayesian methodology used to infer the low-mass-end of the halo mass function and the underlying dark matter properties. Specifically, we introduce the properties of the substructure and the line-of-sight halo populations in Sections 2.1 and 2.2, respectively. We discuss the specifics of the macro-models in Section 2.3. In Section 2.4, we present the posterior probability of the dark matter parameters, given the observed data, and in Section 2.5, we describe our analysis strategy.

2.1 The substructure population

Figure 1.

The subhalo mass functions of CDM and WDM models with different thermal relic particle mass, m_th. WDM models with a larger free streaming length, that is, smaller particle mass, lead to a suppression in the number of substructures at progressively larger substructure masses.

Open in new tabDownload slide

Gao et al. (2004) and Xu et al. (2015) have shown the substructure mass fraction to be a function of the host halo virial mass. Because we do not have virial masses for the lensing galaxies in our sample, we neglect this dependence. As a consequence, our constraints on f_sub should be interpreted as a mean value. While our sample is certainly not homogeneous in this respect, it is also unlikely to span a wide range of virial masses, as, statistically speaking, strong gravitational lens galaxies are more likely to reside in haloes of about 10¹³ M_⊙ (e.g. Sonnenfeld et al. 2018). It should be also noted that we ignore the redshift dependence on f_sub among our sample, since the systematic errors are considered to be larger than the effect of redshift evolution.

2.2 The line-of-sight halo population

Similar to what has been done for the substructure halo population, we include line-of-sight haloes as spherical NFW haloes with virial masses between 10⁶ and 10⁹ M_⊙, and again use the concentration–mass relation from Duffy et al. (2008). We apply the Sheth & Tormen (1999) halo mass function to calculate the number density of haloes per co-moving volume and within the mass range m, m + dm, and the best-fitting parameters are optimized for the Planck cosmology (Despali et al. 2016). We include the effect of the free-streaming properties of the dark matter particles with the same factor used for the substructure mass function; that is, an attenuation with |$(1+M_{\rm hm}/m)^{-1.3}$|⁠. We assume the normalization of the halo mass function to be constant; that is, we assume an average number density of haloes and neglect fluctuations amongst the different lines of sight. We discuss the implications of this choice in Section 4.

We only consider line-of-sight haloes inside two light-cones that share their base on the lens plane and have tips at the observer and at the redshift of the lensed quasar. The base of these cones is given by a circle of two times the Einstein radius, and is centred on the optical axis. To increase the computing efficiency, we consider multiple lens planes along the line of sight with an interval of dz = 0.05, and on each redshift plane, the haloes have projected positions drawn from a uniform prior. Unlike substructures, we assume line-of-sight haloes to be located outside the virial radius of the lensing galaxy; that is, z_los > z_l + z_vir or z_los < z_l − z_vir, with z_vir = 10⁻³ (see Fig. 2 for theoretical results and Sonnenfeld et al. 2018 for more details about the relation between virial radius and stellar mass of the host galaxy).

Figure 2.

The evolution of the virial radius in terms of redshift for different values of the halo virial mass (in units of redshift). The two vertical lines indicate the lowest and highest lens redshift in our sample.

Open in new tabDownload slide

2.3 The macro model

Each lensing galaxy is modelled as a singular isothermal ellipsoid (SIE), with the contribution of an external shear component Γ. Note that higher order terms are also introduced to the macro-model when a luminous satellite and/or a stellar disc is detected. For systems with a luminous satellite, we fix the centroid position of the satellite from optical/IR observations and assume a singular isothermal sphere mass density profile. The Einstein radius of the luminous satellite is then the only free parameter. Similarly, for systems with a detected stellar disc, we introduce an exponential disc component and assign the mass of the disc as a free parameter, while the other parameters are kept fixed at the values inferred from the corresponding light distribution.

2.4 Bayesian inference on dark matter

In the following, we refer to the observed flux ratios, f_i, and positions, |$\boldsymbol {x}$|⁠, of the lensed images, along with their uncertainties as the data, |$\boldsymbol {d}$|⁠. The model parameters that we want to infer are the substructure mass fraction, f_sub, and the half-mode mass, M_hm; these are collectively referred to as the target parameters, |$\vec{\theta }$|⁠. We consider the macro-model parameters; that is, the lensing galaxy mass distribution and the source position, |$\vec{\theta }_{\rm M}$|⁠, as nuisance parameters. Further nuisance parameters are the number of substructures and line-of-sight haloes, N, and the micro-model parameters, |$\vec{\theta }_{\rm m}$|⁠, which include their masses, projected positions, and redshifts. Table 1 summarizes the notations and definitions used in this work.

In the following subsection, we provide more details on the ray-tracing strategy that we have adopted to compute the likelihood and the posterior probability functions.

2.5 Analysis scheme

We assume a uniform prior in log  space for the target parameters in the ranges 0.2 per cent < f_sub < 4 per cent and 10⁵ < M_hm < 10⁹ M_⊙. We achieve this by defining a regular grid in logarithmic space, where all the points are equally weighted. The lower bound of the M_hm prior is due to the fact that the likelihood function flattens to a constant value below ∼10⁵ M_⊙, as the data lose sensitivity to the small difference in the number of detectable perturbers predicted by the different dark matter models (see Fig. 3). Note that we also evaluate the CDM WIMP model at M_hm = 10⁻⁶ M_⊙ and interpolate the likelihood for values up to M_hm = 10⁵ M_⊙. The loss of sensitivity is related to the size and structure of the background source and the opening angle of the merging images (we refer to Section 4.3.3 for a more detailed discussion). From Fig. 3, it can be seen that for dark matter models with M_hm > 10⁸ M_⊙ the number of detectable line-of-sight perturbers quickly approaches zero, hence our choice of M_hm = 10⁹ M_⊙ as an upper limit. We chose the prior that, without any knowledge of the likelihood, is least informative on the scale of M_hm that is uniform in log  space. We discuss this choice in the Appendix.

Figure 3.

The ratio between the expected number of detectable line-of-sight haloes in WDM models of different half-mode mass M_hm and idealized CDM. The expected number of perturbers are obtained by integrating equation (6) between 10⁶ and 10⁹ M_⊙, i.e. the range of halo mass that can be detected with the considered data.

Open in new tabDownload slide

To further increase computing efficiency, we set the quasars to be point sources. For each point on the grid, we calculate the likelihood with a Monte Carlo-based approach as follows:

• First, on each grid point of mass function parameters |$\vec{\theta }$|⁠, we draw a random set of macro-model parameters |$\vec{\theta }_{\rm M}$| from |$Q(\vec{ \theta }_{\rm M})$| (importance sampling).

• Then, we draw the corresponding set of line-of-sight haloes and subhaloes from |$P(N| \vec{ \theta },\vec{ \theta }_{\rm M})$|⁠, we stress here that only the number of substructures depends on f_sub via the projected mass of the main lens within an aperture of 2θ_E.

• For each subhalo and line-of-sight halo, we draw its mass, redshift, and projected position from |$P(\vec{ \theta }_{\rm m}| \vec{ \theta },\vec{ \theta }_{\rm M})$|⁠.

• We use pylens, a python-based ray-tracing package, that implements multiplane lensing with analytical mass profiles to derive the predicted image fluxes and positions, and calculate the relative likelihood.

In total, we generate 100 000 Monte Carlo realizations on each grid point of |$\vec{ \theta }$|⁠, and the posterior probability is then constructed from the summed likelihood. Since each lens is considered independent, we multiply the likelihood of each lens to obtain a joint inference on the model parameters. A schematic view of this strategy is provided in Fig. 4.

Figure 4.

A schematic view of our analysis strategy. On each grid point of target parameters (f_sub, M_hm), we run 100 000 simulations to collect statistical inference by comparing the predicted lens observable with the measured ones.

Open in new tabDownload slide

3 THE DATA

We have collected all radio or MIR observations that are available for the fourteen multiply-imaged quasars that have four lensed images (see Table 2 for a summary of their general properties). Out of these 14 lens systems, only seven are used in our full analysis, while the remaining seven are excluded for the following reasons:

• HE 0435−1223, HS 0810+2554, and RX J0911+080 are bright optical quasars where faint radio emission was detected by Jackson et al. (2015) from deep Very Large Array (VLA) and e-Multi-Element Remotely Linked Interferometer Network (e-MERLIN) imaging at cm wavelengths. Although they demonstrate the feasibility of detecting radio emission from radio-quiet quasars, the lensed images are partially resolved, and are, therefore, not suitable to be modelled as a point source or with a single Gaussian component. Very Long Baseline Interferometry (VLBI) observations of HS0810+2554 at mas-scale angular resolution has confirmed that the radio structure is indeed extended and composed of multiple compact components (Hartley et al. 2019).

• CLASS B1359+154 and CLASS B1608+656 have multiple lensing galaxies within the Einstein radius and show strong flux-ratio anomalies (Koopmans & Fassnacht 1999; Rusin et al. 2001; Suyu et al. 2009). Moreover, CLASS B1608+656 is a merging system, which may cause significant bias on the abundance of small-scale structures. Considering the strong coupling between the flux-ratio anomalies generated by multipole components in the macro-model and by substructure/line-of-sight haloes, we exclude these two systems from our analysis.

• CLASS B1933+503 is a 10-image system with a face-on spiral as the main lensing galaxy (Sykes et al. 1998; Suyu et al. 2012). We notice that the magnification of the lensed images close to the spiral arms has significantly strong deviations from the smooth model predictions. These strong anomalies are very likely due to the presence of the spiral arms. While we exclude this system from our current analysis, we plan to develop an algorithm that includes more complex baryonic structures from simulated disc galaxies into the lens modelling in the future.

• Q2237+030 is gravitationally lensed by the bulge of a low-redshift spiral galaxy (Irwin et al. 1989). The mass distribution of this system is, therefore, dominated by baryonic structures rather than a smooth dark matter distribution. Hence, as for CLASS B1933+503, we decided to exclude this system from our current analysis until we develop an appropriate description.

Our final sample includes the following lens systems: CLASS B0128+437, MG J0414+0534, CLASS B0712+472, PG 1115+080, JVAS B1422+231, CLASS B1555+375, and CLASS B2045+265. Table 3 summarizes the observational data we used in this work (positions and flux-ratios). Improved flux-ratio measurements are available from the MERLIN key programme (Koopmans et al. 2003) for CLASS B0128+437, CLASS B0712+472, JVAS B1422+231, CLASS B1555+375, and CLASS B2045+265. These measurements result in an improved flux-ratio uncertainty of less than 5 per cent. For the remaining systems, MG J0414+0534 and PG 1115+080, we adopt a flux uncertainty of 10 per cent. Each lens is modelled with an SIE plus external shear, except for

• MG J0414+0534 has a luminous satellite (object X) that is detected in optical imaging (Falco et al. 1997), which we include into the lens model and allow its mass to be a free parameter.

• Both CLASS B0712+472 and CLASS B1555+375 have an edge-on disc that lies across the merging images, where the flux-ratio anomalies are most significant (Jackson et al. 1998; Hsueh et al. 2016, 2017). We apply the best-fitting models found by Hsueh et al. (2017), Hsueh et al. (2016), and let the disc mass be the only free parameter. It should also be noted that for the lens system CLASS B1555+375, the redshift of the lensing galaxy is unknown. A recent detection of an emission line in the NIR spectrum suggests that the source redshift is z_s = 1.432 (Fassnacht et al. in preparation). Considering the red colour of the lensing galaxy, it is likely to be at high redshift and we assume the lens redshift to be z_l = 1.0. The disc mass fraction of CLASS B0712+472 and CLASS B1555+375 is about 15 per cent within the Einstein radius, which is consistent with the range of disc mass fraction in the SWELLS survey (Brewer et al. 2014).

The lens system CLASS B2045+265 shows a strong de-magnification on the central image of the cusp triplet (Fassnacht et al. 1999). This strong flux-ratio anomaly was thought to be due to the presence of a luminous satellite detected in NIR imaging, although the lens model was peculiar in that the satellite needed to be highly elongated (McKean et al. 2007). However, additional Keck adaptive optics imaging at a different epoch has shown proper motion of the luminous object, which indicates that it is very likely to be a star. Therefore, we do not include this additional component in the lens model.

For each lens system, Table 4 summarizes the redshift of the lens and source (where available), the Einstein radius, and the expected number of substructure and line-of-sight haloes for the idealized CDM model and for a 8.0 keV thermal relic model. As expected from Sections 2.1 and 2.2, the mean expected number of subhaloes is determined by f_sub, M_hm, and the mass of the host galaxy, while the mean expected number of line-of-sight haloes is a function of M_hm and the volume of the light-cone, that is, the redshift of the source and the lens. We notice that the expected number of line-of-sight haloes can be significantly larger than the number of substructures, depending on the length of the light-path. In the next section, we present how the contribution from these additional haloes affects our inference on f_sub.

4 RESULTS

We present the results of our analysis in two parts. In Section 4.1, we present our constraints on the mass fraction under the assumption of an idealized CDM model (i.e. M_hm = 0) and how the constraints change with the inclusion of line-of-sight haloes. In Section 4.2, we focus on thermal relic dark matter models, which also include the WIMP CDM model at |$M_{\rm hm} = 10^{-6}\, \mathrm{M}_\odot$|⁠. In Section 4.3, we discuss the effect of systematic uncertainties on our results.

4.1 Substructure mass fraction (CDM only)

In this section, we present the results of our analysis with and without line-of-sight haloes under the assumption of an idealized CDM cosmology. Our substructure-only analysis demonstrates the improvements on the precision of the constraints due to the more accurate flux-ratio measurements and lens macro-models that are now available. Importantly, we find that the inclusion of line-of-sight haloes now resolves the tension between the prediction from numerical simulations (Mao et al. 2004; Xu et al. 2015) and the large value of f_sub inferred by Dalal & Kochanek (2002).

4.1.1 CDM substructure only

We first derive constraints on the substructure mass fraction assuming a CDM model, that is, M_hm = 0, and excluding the contribution from line-of-sight haloes. The corresponding posterior probability distribution is presented as the curve marked with triangles in Fig. 5, where the error bars represent the 1σ uncertainties on the Monte Carlo integrals of the probability in each bin. From a Gaussian fit to the posterior curve of each lens, we derive a mean joint value of |$f_{\rm sub} = 0.023^{+0.018}_{-0.010}$| at the 68 per cent confidence level (CL). We find this value to be larger than, but within 2σ of, what is predicted by CDM-only numerical simulations for host galaxies with similar masses and redshifts, and substructure masses in the same range (f_sub = 0.008; Xu et al. 2015; the value is recalculated to fit the definition of f_sub in this work). We stress that due to galaxy formation processes, where baryonic effects may suppress the level of substructure, this discrepancy may be larger (Despali & Vegetti 2017).

Figure 5.

Joint constraints (triangles and thick curve) on the substructure mass fraction (f_sub) without (left) and including (right) line-of-sight haloes from seven gravitationally lensed quasars, under the assumption of an idealized CDM model. The shaded area shows the 1σ uncertainties of the probability distribution. The red solid vertical lines show the median constraint results from Dalal & Kochanek (2002), and the black dashed vertical lines show the upper limit from numerical simulations (Xu et al. 2015).

Open in new tabDownload slide

Interestingly, from a sample of 11 galaxies from SLACS, Vegetti et al. (2014) have inferred a fraction of |$f_{\rm sub} = 0.0064^{+0.0080}_{-0.0042}$|⁠, which is slightly smaller than, but is in much closer agreement with numerical predictions from both dark-matter-only and hydrodynamical CDM simulations (Despali & Vegetti 2017). Although Vegetti et al. (2014) infer a much smaller f_sub than our substructure-only result, we emphasize that the SLACS lenses are at relatively low redshift, and therefore, are less affected by line-of-sight haloes than our higher redshift sample of lensed quasars.

Our results are also consistent with those of Dalal & Kochanek (2002), who analysed a sample of lenses of comparable size and with some overlap in terms of lens systems. The improved precision of our results and an upper limit that is much closer to the theoretical predictions are mainly due to an improvement in the flux-ratio uncertainties from 20 to around 5 per cent. The improvement clearly demonstrates the importance of acquiring accurate flux measurements for more lensed quasars in the future. We also emphasize that, unlike Dalal & Kochanek (2002), we include in the macro-models of CLASS B0712+472 and CLASS B1555+375 an edge-on stellar disc that can explain most of the observed flux-ratio anomalies. The inclusion of this extra component brings down the upper limit on f_sub that had been inferred when fitting a single SIE model to these data.

From Fig. 5, we notice that MG J0414+0534 is essentially an outlier, with a posterior probability that peaks at large values of f_sub, outside our prior range. It is one of the most observed lensed quasars with a wide frequency coverage from radio, submillimetre molecular lines, far-infrared to the MIR. Recently, Stacey & McKean (2018) have shown that the flux ratios for this system do not change with frequency. As different wavelengths should be sensitive to the different mass scale of the perturbation (due to a change in the source size), we conclude that the large anomaly observed in MG J0414+0534 is more likely due to more massive structures, for example, a more complex macro-model. Further investigation and potentially deeper data for this system are required to conclusively understand the origin of the anomaly. When we exclude this system, we infer a mean mass fraction of |$f_{\rm sub} = 0.019^{+0.008}_{-0.009}$| at the 68 per cent CL.

As discussed in Section 3, we have not included the presence of a previously thought luminous satellite in the lens system CLASS B2045+265, given that it has now been identified as a foreground star from proper motion. If, together with MG J0414+0534, we also exclude this system, then we obtain |$f_{\rm sub}=0.018^{+0.013}_{-0.008}$|⁠, in agreement within 1σ of the expectations from CDM-only numerical simulations. We conclude, therefore, that complex macro-models can have a significant impact on the correct interpretation of flux-ratio anomalies, and that further investigations of the lens systems MG J0414+0534 and CLASS B2045+265 are required. A summary of the results from the different choices of sample sets is given in Table 5.

4.1.2 CDM substructure and line-of-sight haloes

Metcalf (2005) and Gilman et al. (2019a) have shown that low-mass haloes located along the lines of sight of the lens galaxies can have a dominant effect on the relative fluxes of multiply-imaged quasars (see Despali et al. 2018 for a similar result for extended lensed images). In this section, we discuss how our inference on the substructure mass fraction changes when the contribution from this population is taken into account, which is also shown in Fig. 5.

As expected, we find a significant drop of 50 per cent in the substructure mass fraction, to |$f_{\rm sub} = 0.012^{+0.007}_{-0.004}$|⁠, at the 68 per cent CL. The result is a reflection of the fact that once the dominant contribution from the line-of-sight haloes is included, a smaller abundance of substructure is needed to reproduce the observed flux ratios. In particular, our results imply that about half of the flux-ratio anomalies are produced by line-of-sight structures. This result may also explain the findings of Vegetti et al. (2014): if the degeneracy with the line-of-sight halo population has a significant impact on the substructure mass fraction inference, this effect is expected to be larger for the sample of high-redshift-lensed quasars considered here, than for the SLACS sample, because of the larger cosmological volume probed by the former.

Our new constraints on the total mass fraction in substructure are now consistent with CDM-only numerical predictions at the 1σ level (f_sub  = 0.008; Xu et al. 2015). Despali & Vegetti (2017) have quantified the suppression in the number density of substructure in hydrodynamical simulations relative to N-body-only simulations to be between 20 and 40 per cent, and the drop in the substructure mass fraction to be between 40 and 70 per cent. According to this correction, our constraints are also in agreement with CDM-hydrodynamical simulations at the 1.2σ level. When we exclude MG J0414+0534, we infer |$f_{\rm sub} = 0.010^{+0.005}_{-0.004}$|⁠, and when we exclude both MG J0414+0534 and CLASS B2045+265, we infer |$f_{\rm sub} = 0.009^{+0.006}_{-0.004}$|⁠. Table 3 also shows that the expectation value of substructures and line-of-sight haloes are at the same level when f_sub = 0.01 for a CDM cosmology.

Fig. 6 presents the comparison between the results from this paper and other studies that constrain f_sub from gravitational lensing. In this work, we define f_sub to be the substructure fraction within the aperture of 2θ_E. In Dalal & Kochanek (2002), f_sub is defined as one-half of the convergence at the critical radius, which is also the same definition used by Vegetti et al. (2014) and Hezaveh et al. (2016). This definition requires the perturbers to be close to the critical radius of the lens, which holds for the simulation designs of Dalal & Kochanek (2002) and the gravitational imaging technique. Following Xu et al. (2015), we recalculate the Dalal & Kochanek (2002) and Vegetti et al. (2014) results to match our definition of f_sub. The main difference between the definitions is that our work probes the gravitational effects of perturbers within a larger aperture, that is, to 2θ_E. However, these numbers are considered comparable within the region the convergence is close to the critical radius, as the projected positions of subhaloes are uniformly distributed. We also emphasize that, although the current estimates on f_sub scatter from 10⁻³ to 10⁻², all of the estimates are marginally in agreement within the 1 to 2σ level.

Figure 6.

Constraints on f_sub (CDM-only) derived in this paper and other strong lensing studies. DK02: Dalal & Kochanek (2002) analysis of seven lensed quasars; H16: Hezaveh et al. (2016) analysis of SDP.81 ALMA observations using the gravitational imaging technique; V14: Vegetti et al. (2014) analysis of 11 SLACS lenses using the gravitational imaging technique; Subs+LOS & Subs-only: results from this work with and without line-of-sight haloes, respectively. Xu15: Xu et al. (2015) analysis of re-scaled haloes in a CDM-only N-body simulation (Aquarius). All uncertainties are presented at the 68 per cent CL (solid) and 95 per cent CL (dashed), except for DK02, which is at the 90 per cent CL.

Open in new tabDownload slide

4.2 Inference on thermal relic dark matter

Schneider et al. (2012a) and Lovell et al. (2014) have shown that the effect of the free streaming of thermal relic dark matter particles can be well described by the half-mode mass M_hm. In terms of the subhalo mass function, this introduces an extra free parameter that is degenerate with the substructure mass fraction. Suppression in the number of low-mass haloes can either be related to a larger value of the half-mode mass or conversely to a lower dark matter fraction in the substructure. While the value of the half-mode mass solely depends on the dark matter physics, f_sub is related to the accretion history of the host halo and the efficiency of tidal disruption.

As we have assumed a mean line-of-sight halo mass function, we have ignored any degeneracy between the halo mass function normalization and M_hm. We expect that this degeneracy does not cause a significant effect on the results of this paper, as the systematic uncertainties in the macro-model for those galaxies that lack deep imaging information should be larger than the scatter introduced by varying lines of sight.

Under the assumption of a thermal relic dark matter model, we infer a mean dark matter fraction of f_sub = 0.013 ± 0.007, where the slightly increased value relative to our CDM-only results can be explained by the non-zero value of the half-mode mass. At present, we are not aware of any numerical simulation on cosmological scales of massive galaxies with a thermal relic cosmology that resolves small-scale haloes. Therefore, we cannot compare our f_sub results with theoretical expectations. In the following, we present the marginalized constraints on M_hm and discuss the implications for the free streaming properties of dark matter. In this part of the analysis, we have included the contribution of both substructure and line-of-sight haloes.

Fig. 7 presents our constraints on the free-streaming property of thermal relic dark matter from our sample of seven gravitationally lensed quasars. In both panels, the grey shaded area shows the 1 − σ uncertainty of the Monte Carlo integral on the probability distribution of the joint constraint. As anticipated, the likelihood function for all lenses flattens to a constant value for dark matter models colder than about M_hm ≤ 10^5.5 M_⊙ due to the small difference in the number of haloes that can be detected between the different dark matter models. Similarly, we expect the likelihood to reach a constant value above M_hm = 10⁹ M_⊙ as the number of detectable perturbers is equally consistent with zero for all models.

Figure 7.

The joint posterior probability distribution for the half-mode mass (left) and thermal relic particle mass (right) from our sample of seven gravitationally lensed quasars. The grey shaded area represents the 1σ uncertainty of the joint constraint. The black vertical line represents our upper limit on M_hm and lower limit on m_th at the 95 per cent CL. The black arrows show the direction of the allowed region at the 95 per cent CL from this work. The red dashed and solid lines represent the lower limits (95 per cent CL) from the latest Ly α forest constraints, assuming a smooth and non-smooth intergalactic medium temperature evolution, respectively (Iršič et al. 2017).

Open in new tabDownload slide

From Fig. 7, we also notice that the lens systems PG 1115+080 and CLASS B0128+437 show a clear preference for warmer models, especially towards the region where no detectable perturber is expected. This is due to the fact that these two systems show no flux-ratio anomaly, so that a smooth macro-lens model can successfully predict the image flux ratios. In contrast, the posteriors of other lenses increase towards the colder region until the limit of the data sensitivity. We interpret this bi-modal distribution as the result of small-sample statistics. Despite this, our results demonstrate that the analysis of gravitationally lensed quasars is a promising approach to further explore the parameter space of dark matter models.

Our joint analysis results in an upper limit at the 95 per cent CL for the half-mode mass of M_hm < 10^7.8 M_⊙, which corresponds to a lower limit on the thermal relic mass of m_th > 5.6 keV at the same CL. In comparison with the latest 2σ constraints from the Ly α forest, m_th > 5.3 keV or m_th > 3.5 keV, depending on the assumption of the intergalactic medium thermal histo (Bolton et al. 2008; Iršič et al. 2017; Garzilli et al. 2018), our current results provide a comparable level of constraints on the mass of the thermal relic dark matter particle, with the potential of significantly improving with upcoming large sample of lensed quasars (Gilman et al. 2019a). Our results are also in agreement with a recent analysis of gravitationally lensed quasars by Gilman et al. (2019b). However, due to a different choice of prior on the half-mode mass, at present, it is not clear how robustly the two results can be compared with each other.

4.3 Systematic uncertainties

Throughout our analysis, there are several factors that can introduce systematic uncertainties into our inference. We discuss these factors in this section and how they can be addressed in the future.

4.3.1 Stellar structures

Here, we have already included edge-on discs as a higher order component in the lens models for the cases of CLASS B0712+472 and CLASS B1555+375, since these two lenses show solid evidence of stellar discs in high-resolution IR imaging (Hsueh et al. 2016, 2017). However, Gilman et al. (2018) have shown that even early-type galaxies can have disc-like structures that can contribute to flux-ratio anomaly signals. An analysis of a sample of lens-like galaxies from the Illustris simulation also suggests that stellar structures in early-type lenses can increase the level of flux-ratio anomalies by around 10 per cent (Hsueh et al. 2018). In our sample, Keck adaptive-optics IR imaging of CLASS B0128+437 shows evidence of a face-on late-type galaxy (McKean et al. 2004; Lagattuta, Auger & Fassnacht 2010) and the lens model for CLASS B2045+265 requires a significant elliptical component to the mass model (McKean et al. 2007). By not accounting for possible stellar structures in the lens modelling, our analysis could potentially overestimate the abundance of low-mass haloes, and hence artificially favour colder dark matter models. Although evaluating the impact from undetected stellar structures is not possible with current observations, combining the stellar structures from the latest hydrodynamical simulations of lens-like galaxies, as for example, with the Illustris TNG suite of simulations, into the analysis or obtaining kinematic information on the lensing galaxies will help to reduce this source of systematic uncertainty. This will be the focus of a future theoretical work.

4.3.2 Source structures

To optimize the computing efficiency, we have assumed the source quasars to be with point like. There is evidence from VLBI observations that CLASS B0128+437 (Biggs et al. 2004), MG J0414+0534 (Ros et al. 2000), and CLASS B1555+375 (Hsueh et al. 2016) have extended background sources on mas-scales, but CLASS B0712+472 (Hsueh et al. 2017), JVAS B1422+231 (Patnaik et al. 1999), and CLASS B2045+265 (McKean et al. 2007) have compact structures; VLBI data for PG 1115+080 has been taken, but is yet unpublished. The size of the source affects the sensitivity of the data to small-scale structures, with larger sources being less sensitive to low-mass haloes. Whether this effect results in an artificial and significant preference for warmer dark matter models (potentially compensating for the effect of stellar structures) is unclear, as the size and internal structures of the quasar source can vary from one lens system to the other. In the near future, multiwavelength observations at radio and submillimetre wavelengths will help us to gain further information on the size and structure of the sources and correctly include them in our analysis.

4.3.3 Source variability

Any variation of the background radio source flux will be seen at different times in the multiple lensed images, meaning that flux measurements taken at a single epoch are sampling the intrinsic light curve of the quasar at different times for the different images. For this reason, in previous studies, such as Dalal & Kochanek (2002), the flux-ratio uncertainties were assigned to be 20 per cent for all of the systems with one-time flux measurements. The systematic uncertainties from quasar variability can be eliminated by averaging over a long period of monitoring. In this work, we quote the average flux-ratios from Koopmans et al. (2003), which bring the uncertainties in the intrinsic variation from 20 per cent to less than 5 per cent. Most of the radio-sources in our sample have also shown little variation throughout the monitoring. The clear improvement in the constraints produced by the analysis described in Section 4.1.1 emphasizes the importance of monitoring observations or some other technique for improving the precision of the flux-ratio measurements.

4.3.4 Propagation effects

Although propagation effects such as free–free absorption and scatter broadening can alter the properties of the different lensed images measured at radio-wavelengths, these effects have a strong wavelength dependence, and therefore, can be identified and corrected for with multiwavelength observations (e.g. Biggs et al. 2003; Winn, Rusin & Kochanek 2004; Mittal, Porcas & Wucknitz 2007). There is no clear evidence of propagation effects in our sample, except for CLASS B0128+437 (Biggs et al. 2004) where there is scatter broadening of the lensed images on VLBI scales. It is expected that the contribution from propagation effects to the systematic uncertainties in our analysis are small, but recently completed multi-frequency imaging campaigns with the VLA and VLBI will resolve this.

4.3.5 Properties of subhaloes and field haloes

We have assumed both field haloes and subhaloes to be well described by spherical NFW profiles with a concentration–mass relation from Duffy et al. (2008). While this is a reasonable assumption for small dark isolated field haloes, the profile of a subhalo may be affected by its interaction with the host halo and depends on its distance from the host centre and its properties at infall. In general, subhaloes are expected to be more concentrated than isolated haloes, especially those closer to the host centre due to tidal truncation (Hayashi et al. 2003; Springel et al. 2008; Moliné et al. 2017). However, Despali et al. (2018) found that, in terms of the lensing efficiency, neglecting the difference in the concentration–mass relation between subhaloes and field haloes leads to an uncertainty in the inferred mass of the order of 5–20 per cent for subhalo masses between 10⁶ and 10⁹ M_⊙, and an even smaller effect is caused by neglecting the dependence on distance. For this reason, we believe that neglecting the details of the subhalo concentration introduces only a second-order effect in our analysis. Moreover, recent works have shown that the properties of simulated subhaloes and the details of the disruption process, which influence the density profiles, are still not well constrained in numerical simulation at the low-mass end (van den Bosch & Ogiya 2018; van den Bosch et al. 2018; Errani & Peñarrubia 2019; Green & van den Bosch 2019).

Another relevant aspect is the effect of warm dark matter on the concentration–mass relation. Previous works (Ludlow et al. 2016) have shown that the concentration of low-mass (sub)haloes is depleted with respect to CDM, in a way that is similar to the effect on the mass function. This might result in a smaller lensing signal by low-mass (sub)haloes in WDM models, which we have not considered here. However, the (sub)halo mass is the main parameter that determines the strength of the lensing signal and the concentration has been shown to have a second-order effect (Despali et al. 2018) and thus we believe that our results are not significantly affected by not including a more complex description of the concentration–mass relation. Moreover, the concentration–mass relation has so far only been constrained for a limited number of specific WDM realizations, and a detailed parametrization as a function of the WDM model is still lacking.

Throughout the analysis, we have fixed the number density of all line-of-sight realizations to the average density of the Universe and thus neglected the fact that some lines of sight might be over- or underdense. In a more stringent analysis, these realization should reflect the environment of each lens and the characteristic variance of line-of-sight structures throughout the Universe. However, we have information on the environment and the field galaxies only for a few quasar lenses in our sample. Moreover, the link between the number density of luminous field galaxies and dark field haloes is ambiguous. In the future, we plan to use high-resolution large-scale simulations to quantify any variation between different lines of sights, and we will investigate the assumption that field haloes are isolated and do not contain significant subhaloes.

Finally, we parametrized the WDM mass function following the results of numerical simulations (Lovell et al. 2012; Schneider et al. 2012b). We extrapolated the fitting functions from these works below the resolution limit of the simulations (M < 10⁷M_⊙), where the functional forms result in a sharp decrease in the number of objects. It remains to be shown definitively whether this drop-off rate describes the WDM models accurately and higher resolution simulations would be required to confirm it.

5 CONCLUSIONS

We have analysed a sample of seven gravitationally lensed quasars and used the observed image positions and relative fluxes to probe the abundance of low-mass haloes within the potential of the lensing galaxies and along their lines of sight. Our results can be summarized as follows:

• We find that accurate flux ratio measurements are a key ingredient for the derivation of precise constraints on the (sub)halo mass function. By improving the flux-ratio uncertainties from 20 to better than 5 per cent, we substantially bring down both the upper limit and the uncertainty on the normalization of the subhalo mass function, when compared to a previous study by Dalal & Kochanek (2002), based on a sample of comparable size.

• Under the assumption of an idealized CDM model, we find that the degeneracy between the substructure and line-of-sight haloes has a significant effect on the inferred substructure mass fraction. In particular, the inclusion of a line-of-sight population brings our constraints on f_sub into much closer agreement with the expectations from both CDM-only and hydrodynamical simulations. This result also explains the long-standing discrepancy between the dark matter fraction inferred by Dalal & Kochanek (2002) and numerical simulations (Xu et al. 2015).

• The inclusion of extra complexity in the mass model of the lensing galaxies, although subdominant for the sample considered here, also plays an important role for a correct interpretation of flux-ratio anomalies. In particular, the inclusion of a stellar disc in the macro-model makes the edge-on disc lenses no longer outliers in terms of the strength of the observed anomalies. Deep imaging observations are therefore crucial to break the degeneracy between the stellar structures and small-scale dark matter perturbers. However, as the effect of other complex structures at a smaller scale, such as spiral arms and the intrinsic un-smoothness of elliptical galaxies, has not been properly quantified yet, we plan to address this issue in a follow-up paper.

• Under the assumption of a thermal relic dark matter model, we constrain the dark matter particle mass to be m_th > 5.6 keV at the 95 per cent confidence level. Our limits are in agreement with observations of the Ly α forest (Iršič et al. 2017), showing that the study of gravitationally lensed quasars can provide comparably strong constraints with the current sample size. Furthermore, compared to the uncertainties in the intergalactic medium thermal history, gravitationally lensed quasars provide a more direct and robust constraint on dark matter properties, which will be further improved with the analysis of increasingly larger sample of lens systems.

ACKNOWLEDGEMENTS

We thank Andrew Robertson, Alexey Boyarsky, Mark Lovell, Simon Birrer, Daniel Gilman, Chuck Keeton, and Dandan Xu for helpful discussions on this work. We also thank Ethan Nadler for important feedback on our preprint. CDF and JWH acknowledge support from the National Science Foundation under grant no. AST-1715611. JWH and JPM acknowledge support from the Netherlands Organization for Scientific Research (NWO) (project no. 629.001.023), Chinese Academy of Sciences (CAS) (project no. 114A11KYSB20170054), and the Dutch National Supercomputer Service. SV has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 758853). LVEK acknowledges support by a VICI grant (project no. 614.001.206) from NWO.

REFERENCES