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

TL;DR

This study uses gravitational lensing of quasars to constrain dark matter properties, including free-streaming length and substructure abundance, aligning with Cold Dark Matter predictions and setting a lower limit on particle mass.

Contribution

First analysis to include line-of-sight haloes and complex galaxy structures in lens models, improving dark matter substructure constraints.

Findings

Substructure mass fraction f_sub = 0.012^{+0.007}_{-0.004}

Lower limit on thermal relic dark matter particle mass m_th > 5.58 keV

Results consistent with Cold Dark Matter simulations

Abstract

We present an analysis of seven strongly gravitationally lensed quasars and the corresponding constraints on the properties of dark matter. Our results are derived by modelling the lensed image positions and flux-ratios using a combination of smooth macro models and a population of low-mass haloes within the mass range 10^6 to 10^9 Msun. Our lens models explicitly include higher-order complexity in the form of stellar discs and luminous satellites, as well as low-mass haloes located along the observed lines of sight for the first time. Assuming a Cold Dark Matter (CDM) cosmology, we infer an average total mass fraction in substructure of f_sub = 0.012^{+0.007}_{-0.004} (68 per cent confidence limits), which is in agreement with the predictions from CDM hydrodynamical simulations to within 1 sigma. This result is closer to the predictions than those from previous studies that did not include line-of-sight haloes. Under the assumption of a thermal relic dark matter model, we derive a lower limit on the particle relic mass of m th > 5.58 keV (95 per cent confidence limits), which is consistent with a value of m_th > 5.3 keV from the recent analysis of the Ly-alpha forest. We also identify two main sources of possible systematic errors and conclude that deeper investigations in the complex structure of lens galaxies as well as the size of the background sources should be a priority for this field.

Tables

Table 1. Table 1: Summary of the different parameters that enter our analysis.

Parameter	Definition	Details
$\overrightarrow{\theta }$	mass function parameters	$f_{\mathrm{sub}},M_{\mathrm{hm}}$;
$f_{\mathrm{sub}}$	mass fraction in substructure
	w/i $2{\theta }_e$ projected area
$M_{\mathrm{hm}}$	half-mode mass	M⊙
$\overrightarrow{{\theta }_M}$	macro-model parameters	SIE+$\Gamma$, source position,
		mass of luminous satellite or stellar disc
$\overrightarrow{{\theta }_m}$	micro-model parameters	substructure: mass, position
		line-of-sight halo: redshift, mass, position
$𝐝$	observed data	lensed image position, flux ratio
$N$	number of low-mass perturbers	$N=N_{\mathrm{sub}}+N_{\mathrm{los}}$

Table 2. Table 2: Summary of the multiply imaged lensed quasars with radio or MIR flux measurements. The references for the lensed image positions and fluxes, and references for the evidence of more complex lens models (e.g., luminous satellites and stellar discs) are also listed. Additional information can also be found on the CASTLES lens database at https://www.cfa.harvard.edu/castles/. The symbol (†) indicates those systems that have radio flux-ratio measurements obtained from monitoring (Koopmans et al., 2003 ) .

Lens	Type	Radio flux	MIR flux	Satellite	Disc	References
CLASS B0128+437 †	fold	$\surd$	–	–	–	Phillips et al. (2000)
MG J0414+0534	fold	$\surd$	$\surd$	$\surd$	–	Falco et al. (1997); Minezaki et al. (2009)
HE 0435$-$1223	cross	$\surd$	–	$\surd$	–	Wisotzki et al. (2002); Jackson et al. (2015)
CLASS B0712+472 †	fold	$\surd$	–	–	$\surd$	Jackson et al. (1998); Hsueh et al. (2017)
HS 0810+2554	fold	$\surd$	–	–	–	Reimers et al. (2002); Jackson et al. (2015)
RX J0911+0551	cusp	$\surd$	–	$\surd$	–	Bade et al. (1997); Jackson et al. (2015)
PG 1115+080	fold	–	$\surd$	–	–	Weymann et al. (1980); Chiba et al. (2005)
CLASS B1359+154 †	fold	$\surd$	–	–	–	Myers et al. (1999); Rusin et al. (2001)
JVAS B1422+231 †	cusp	$\surd$	$\surd$	–	–	Patnaik et al. (1999)
CLASS B1555+375 †	fold	$\surd$	–	–	$\surd$	Marlow et al. (1999); Hsueh et al. (2016)
CLASS B1608+656	fold	$\surd$	–	–	–	Koopmans & Fassnacht (1999); Fassnacht et al. (2002)
CLASS B1933+503	fold	$\surd$	–	–	$\surd$	Sykes et al. (1998); Suyu et al. (2012)
CLASS B2045+265 †	cusp	$\surd$	–	–	–	Fassnacht et al. (1999); McKean et al. (2007)
Q2237+030	cross	–	$\surd$	–	–	Huchra et al. (1985); Minezaki et al. (2009)

Table 3. Table 3: Summary of the observational data for the seven gravitationally lensed quasars we used in our analysis. The lensed image positions are in units of arcsec. The values in the parentheses are the corresponding uncertainties. Note that MG J0414+0534 and PG 1115+080 do not have flux monitoring data from Koopmans et al. ( 2003 ) and their flux-ratio uncertainties are conservatively assigned to be 10 per cent.

Lens	Image	Positions		Flux ratio	Reference
		RA	DEC
CLASS B0128+437	A	$\equiv 0$	$\equiv 0$	$\equiv 1$
	B	+0.098 (0.003)	+0.094 (0.003)	0.584 (0.029)	Phillips et al. (2000)
	C	+0.520 (0.003)	$-$0.172 (0.003)	0.520 (0.029)
	D	+0.108 (0.003)	$-$0.250 (0.003)	0.506 (0.032)
MG J0414+0534	A1	+0.5876 (0.0003)	$-$1.9341 (0.0003)	$\equiv 1$
	A2	+0.7208 (0.0003)	$-$1.5298 (0.0003)	0.9027 (0.0903)	Katz et al. (1997)
	B	$\equiv 0$	$\equiv 0$	0.3890 (0.0389)
	C	$-$1.3608 (0.0007)	$-$1.6348 (0.0008)	0.1446 (0.0145)
CLASS B0712+472	A	$\equiv 0$	$\equiv 0$	$\equiv 1$
	B	+0.056 (0.003)	$-$0.156 (0.003)	0.843 (0.061)	Hsueh et al. (2017)
	C	+0.812 (0.003)	$-$0.663 (0.003)	0.418 (0.037)
	D	+1.174 (0.003)	+0.459 (0.003)	0.082 (0.035)
PG 1115+080	A1	+1.328 (0.003)	$-$2.034 (0.003)	$\equiv 1$
	A2	+1.477 (0.004)	$-$1.576 (0.004)	0.93 (0.093)	CASTLES,
	B	$-$0.341 (0.003)	$-$1.961 (0.003)	0.16 (0.016)	Chiba et al. (2005)
	C	$\equiv 0$	$\equiv 0$	0.21 (0.021)
JVAS B1422+231	A	+0.38925 (0.00005)	+0.31998 (0.00005)	$\equiv 1$
	B	$\equiv 0$	$\equiv 0$	1.062 (0.009)	Patnaik et al. (1999)
	C	$-$0.33388 (0.00005)	$-$0.74771 (0.00005)	0.551 (0.007)
	D	+0.95065 (0.00005)	$-$0.80215 (0.00005)	0.024 (0.006)
CLASS B1555+375	A	$\equiv 0$	$\equiv 0$	$\equiv 1$
	B	$-$0.0726 (0.001)	+0.0480 (0.001)	0.620 (0.039)	Hsueh et al. (2016)
	C	$-$0.4117 (0.001)	$-$0.0280 (0.001)	0.507 (0.030)
	D	$-$0.1619 (0.003)	$-$0.3680 (0.003)	0.086 (0.024)
CLASS B2045+265	A	$\equiv 0$	$\equiv 0$	$\equiv 1$
	B	$-$0.1338 (0.0001)	$-$0.2483 (0.0001)	0.578 (0.059)	McKean et al. (2007)
	C	$-$0.2877 (0.0001)	$-$0.7904 (0.0001)	0.739 (0.073)
	D	+1.6268 (0.0002)	$-$1.0064 (0.0002)	0.102 (0.025)

Table 4. Table 4: Each column represents (1) the lens name, (2) the lens redshift, (3) the source redshift, (4) the Einstein radius (in arcsec), (5) the opening angle, (6) to (8) the expectation value of substructures in the main halo for the idealized CDM model with f sub = 0.5 subscript 𝑓 sub 0.5 f_{\rm sub}=0.5  per cent, 1 per cent, 2 per cent, respectively, (9) the expectation value of line-of-sight haloes for the idealized CDM model within 2 ​ θ E 2 subscript 𝜃 𝐸 2\theta_{E} , (10) to (12) the expectation value of substructures for a thermal relic model with a particle mass of 8.0 8.0 8.0  keV with f sub = 0.5 subscript 𝑓 sub 0.5 f_{\rm sub}=0.5  per cent, 1 per cent, 2 per cent, respectively, and (12) the expectation value of line-of-sight haloes for the same thermal relic model. Redshift references: CLASS B0128+437; McKean et al. ( 2004 ) , MG J0414+0534; Tonry & Kochanek ( 1999 ) , CLASS B0712+472; Fassnacht & Cohen ( 1998 ) , PG 1115+080; Weymann et al. ( 1980 ); Tonry ( 1998 ) , JVAS B1422+231; Patnaik et al. ( 1999 ); Impey et al. ( 1996 ) , CLASS B1555+375; Fassnacht et al. in prep.; CLASS B2045+265; Fassnacht et al. ( 1999 ); McKean et al. ( 2007 ) .

					CDM				Thermal relic (8.0 keV)
					${\mu }_{\mathrm{sub}}$			${\mu }_{\mathrm{los}}$	${\mu }_{\mathrm{sub}}$			${\mu }_{\mathrm{los}}$
Lens	$z_l$	$z_s$	${\theta }_E$ (arcsec)	$\Delta \varphi$ (deg.)	0.5%	1%	2%		0.5%	1%	2%
CLASS B0128+437	1.145	3.12	0.21	123.3	7	15	29	240	2	4	6	29
MG J0414+0534	0.96	2.64	1.12	101.5	631	1262	2474	4666	128	255	510	552
CLASS B0712+472	0.41	1.34	0.69	76.9	104	209	417	206	21	42	84	25
PG 1115+080	0.3098	1.722	1.13	121.2	226	452	903	507	45	92	183	61
JVAS B1422+231	0.34	3.62	0.75	77.0	99	197	394	395	20	40	80	47
CLASS B1555+375	(1.0)	1.432	0.25	102.6	39	78	157	89	8	16	32	11
CLASS B2045+265	0.87	2.35	1.08	34.9	556	1113	2225	3497	112	224	450	414

Table 5. Table 5: Summary of our constraints on the dark matter mass fraction f sub subscript 𝑓 sub f_{\rm sub} (idealised CDM-only), the thermal relic particle mass m th subscript 𝑚 th m_{\rm th} and the half-mode mass M hm subscript 𝑀 hm M_{\rm hm} for different sub-samples of lensed quasars.

Sample	$f_{\mathrm{sub}}$ (CDM-only)		$m_{\mathrm{th}}$	$M_{\mathrm{hm}}$
7 quasar lenses	${0.012}_{-0.004}^{+0.007}$ (subs+LOS)	${0.023}_{-0.010}^{+0.018}$ (subs-only)	$>5.58$ keV (95% CL.)	M⊙ (95% CL.)
Exclude MG J0414+0534	${0.010}_{-0.004}^{+0.006}$	${0.019}_{-0.009}^{+0.008}$	$>5.30$ keV	M⊙
Exclude MG J0414+0534	${0.009}_{-0.004}^{+0.006}$	${0.018}_{-0.008}^{+0.013}$	$>4.77$ keV	M⊙
& CLASS B2045+265

Equations

$$n_{\mathrm{sub}}(m)=\frac{d^2{N}_{\mathrm{sub}}(m)}{dmdA}=n_0{\left(\frac{m}{m_0}\right)}^{-1.9}{\left(1+\frac{M_{\mathrm{hm}}}{m}\right)}^{-1.3},$$

$$M_{\rm{hm}}=\frac{4\pi}{3}\bar{\rho}\Big{(}6.97\lambda^{\rm{eff}}_{\rm{fs}}\Big{)}^{3}$$

$$\lambda^{\rm{eff}}_{\rm{fs}}=0.049~{}\Big{(}\frac{m_{\rm{th}}}{keV}\Big{)}^{-1.11}\Big{(}\frac{\Omega_{\rm{th}}}{0.25}\Big{)}^{0.11}\Big{(}\frac{h}{0.7}\Big{)}^{1.22}~{}{\rm Mpc~{}h^{-1}},$$

$$m_{\rm th}=2.32~{}\Big{(}\frac{M_{\rm hm}}{10^{9}M_{\odot}}\Big{)}^{-0.3}\Big{(}\frac{\Omega_{\rm th}}{0.25}\Big{)}^{0.4}\Big{(}\frac{h}{0.7}\Big{)}^{0.8}~{}{\rm keV}.$$

$$f_{\mathrm{sub}}=\frac{\Sigma m_{\mathrm{sub},\mathrm{proj}}}{M_{\mathrm{proj}}}=\frac{{\int }_{M_{\mathrm{low}}}^{M_{\mathrm{high}}}{n}_{\mathrm{sub}}(m)mdm\cdot A_{\mathrm{proj}}}{M_{\mathrm{proj}}},$$

$${\mu }_{\mathrm{sub}}={\int }_{M_{\mathrm{low}}}^{M_{\mathrm{high}}}{n}_{\mathrm{sub}}(m)dm\cdot A_{\mathrm{proj}}.$$

$$\begin{gathered} P(\vec{\theta }|\mathbf{d})\propto \sum _i^N\int P(\mathbf{d}|\vec{\theta },{\vec{\theta }}_{\mathrm{M}},{\vec{\theta }}_{\mathrm{m}},N)P({\vec{\theta }}_{\mathrm{m}}|\vec{\theta },{\vec{\theta }}_{\mathrm{M}})P(N|\vec{\theta },\overrightarrow{{\theta }_{\mathrm{M}}}) \\ \times P({\vec{\theta }}_{\mathrm{M}})P(\vec{\theta })d{\vec{\theta }}_{\mathrm{M}}d{\vec{\theta }}_{\mathrm{m}}. \end{gathered}$$

$$P(\mathbf{d}|\vec{\theta },{\vec{\theta }}_{\mathrm{M}},{\vec{\theta }}_{\mathrm{m}},N)\approx \mathcal{G}({\chi }_{\mathrm{tot}}^2(\vec{\theta },{\vec{\theta }}_{\mathrm{M}},{\vec{\theta }}_{\mathrm{m}},N)),$$

$$\mathcal{G}\propto \exp \left(-\frac{1}{2}\left({\chi }_{\mathrm{flux}}^2+{\chi }_{\mathrm{pos}}^2\right)\right).$$

$$P({\vec{\theta }}_{\mathrm{M}})⟶\frac{P({\vec{\theta }}_{\mathrm{M}})}{Q({\vec{\theta }}_{\mathrm{M}})}Q({\vec{\theta }}_{\mathrm{M}}),$$

$$\mathcal{P}(M_{\mathrm{hm}},f_{\mathrm{sub}})\propto M_{\mathrm{hm}}^{-1}{f}_{\mathrm{sub}}^{-1}\mu (M_{\mathrm{hm}},f_{\mathrm{sub}})\sqrt{\mu (M_{\mathrm{hm}},f_{\mathrm{sub}})+C},$$

Full text

††pubyear: 2019††pagerange: SHARP – VII. New constraints on the dark matter free-streaming properties and substructure abundance from gravitationally lensed quasars–A

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 et al., 2010a; Vegetti et al., 2010b, 2012; Nierenberg et al., 2014; Hezaveh et al., 2016; Gilman et al., 2018; Bayer 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 et al., 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^{8}$ and $10^{9}$ M⊙ with optical/infrared (IR) (Vegetti et al., 2010a; Vegetti et al., 2010b, 2012) and sub-millimetre 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^{7}$ 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^{8}$ and $10^{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, 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 Möller et al., 2003; Quadri et al., 2003; Chen et al., 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, 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 $\Lambda$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 $\Omega_{\rm M}=0.28$ and $H_{0}=70$ km s*-1* Mpc*-1*.

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

We assume the substructure population to be well described by spherical NFW haloes (Navarro, Frenk & White, 1997) with virial masses between $10^{6}$ and $10^{9}$ M⊙, and a concentration-mass relation from Duffy et al. (2008). We neglect the effect of tidal truncation and changes in the concentration-mass relation as a function of distance from the host halo centre, as both effects have been shown to be of secondary importance in terms of the lensing effects of these low-mass haloes (Despali et al., 2018). Following Schneider et al. (2012a) and Lovell et al. (2014), we parameterize the substructure mass function (i.e., the number density of substructures in the mass range $m,m+dm$ per unit area) as,

[TABLE]

where $M_{\rm{hm}}$ is the half-mode mass, that is, the mass scale at which the transfer function is suppressed by 50 per cent relative to CDM ($M_{\rm{hm}}=0$ for idealized CDM and $M_{\rm{hm}}\sim 10^{-6}M_{\odot}$ for WIMPs, Schneider et al., 2013). The normalization constants $n_{0}$ and $m_{0}$ can be related to the projected total mass fraction in subhaloes defined below. For thermal relic dark matter models, the half-mode mass is related to the mass of the dark matter particle by,

[TABLE]

(Viel et al., 2005; Schneider et al., 2012a), where $\bar{\rho}=\Omega_{M}\cdot\rho_{\rm crit}$ is the background density of the universe and $\lambda^{\rm{eff}}_{\rm{fs}}$ is the effective free-streaming length scale, which is given by,

[TABLE]

where $m_{\rm th}$ is the thermal relic mass of the dark matter particle.

The inverse conversion between $m_{\rm th}$ and $M_{\rm hm}$ can be expressed as (Nadler et al., 2019):

[TABLE]

In Figure 1, we show the substructure mass function for the idealized CDM and thermal relic models of different particle mass.

As lensing is sensitive to the total projected mass distribution within a cylinder, we define $f_{\rm sub}$ as the ratio between the total mass in substructure within a projected cylinder with a radius twice as large as the main lens Einstein radius, $\theta_{E}$, and the total mass of the main halo within the same projected cylinder. That is,

[TABLE]

where $A_{\rm proj}$ is the area within the aperture of the projected cylinder and $(M_{\rm{low}},M_{\rm{high}})=(10^{6},10^{9})~{}$M⊙. The expectation value of substructures, $\mu_{\rm sub}$, within the projected cylinder is expressed as,

[TABLE]

Following Xu et al. (2015), we assume the projected position of substructure on the plane of the lensed images to be uniform within $2\theta_{E}$.

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_{\rm 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^{13}$ M⊙ (e.g., Sonnenfeld et al., 2018). It should be also noted that we ignore the redshift dependence on $f_{\rm 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^{6}$ and $10^{9}$ 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 attentuation with $\left(1+M_{\rm hm}\right/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_{\rm los}>z_{l}+z_{\rm vir}$ or $z_{\rm los}<z_{l}-z_{\rm vir}$, with $z_{\rm vir}=10^{-3}$ (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).

2.3 The macro model

Each lensing galaxy is modelled as a singular isothermal ellipsoid (SIE), with the contribution of an external shear component $\Gamma$. 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 (SIS) 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, $\mathbf{x}$, of the lensed images, along with their uncertainties as the data, $\mathbf{d}$. The model parameters that we want to infer are the substructure mass fraction, $f_{\rm sub}$, and the half-mode mass, $M_{\rm 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.

Using Bayes theorem, we relate the posterior probability density of $\vec{\theta}$, marginalized over the nuisance parameters, to the likelihood function $P(\mathbf{d}|\vec{\theta},\vec{\theta}_{\rm M},\vec{\theta}_{\rm m},N)$ as,

[TABLE]

Under the assumption of Gaussian errors on the observed fluxes and positions, the log-likelihood function is related to $\chi_{\rm tot}^{2}=\chi_{\rm flux}^{2}+\chi_{\rm pos}^{2}$, such that,

[TABLE]

where

[TABLE]

For $P(N|\vec{\theta},\vec{\theta}_{\rm M})$, we adopt a Poisson probability function with an expectation number of substructures, $\mu_{\rm sub}$, and line-of-sight haloes, $\mu_{\rm los}$, respectively. The contribution from each population is given by the integral over the respective mass functions as defined in Sections 2.1 and 2.2. We apply a Monte-Carlo approach to compute the integral in eq. (7). To increase the computing efficiency, we introduce importance sampling of the macro-model prior, $P(\vec{\theta}_{\rm M})$, as,

[TABLE]

where $Q(\vec{\theta}_{\rm M})$ is obtained from an MCMC modelling of the data under the assumption of $N=0$. We then draw $\vec{\theta}_{\rm M}$ realisations from $Q(\vec{\theta}_{\rm M})$. The likelihood of each realisation is then weighted by $P(\vec{\theta}_{\rm M})/Q(\vec{\theta}_{\rm M})$. For the micro-model parameters, we adopt priors as discussed in Sections 2.1 and 2.2, respectively.

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~{}{\rm per~{}cent}<f_{\rm sub}<4~{}{\rm per~{}cent}$ and $10^{5}<M_{\rm hm}<10^{9}$ 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_{\rm hm}$ prior is due to the fact that the likelihood function flattens to a constant value below $\sim 10^{5}~{}\textrm{M}_{\odot}$, as the data lose sensitivity to the small difference in the number of detectable perturbers predicted by the different dark matter models (see Figure 4). Note that we also evaluate the CDM WIMP model at $M_{\rm hm}=10^{-6}$ M⊙* and interpolate the likelihood for values up to $M_{\rm hm}=10^{5}$ 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 Figure 4, it can bee seen that for dark matter models with $M_{\rm hm}>10^{8}~{}M_{\odot}$ the number of detectable line-of-sight perturbers quickly approaches zero, hence our choice of $M_{\rm hm}=10^{9}~{}M_{\odot}$ as an upper limit. We chose the prior that, without any knowledge of the Likelihood, is least informative on the scale of $M_{\rm hm}$, that is uniform in $\log$ space. We discuss this choice in Appendix A.

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}_{M}$ from $Q(\vec{\theta}_{\rm M})$ (importance sampling); 2. 2.

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_{\rm sub}$ via the projected mass of the main lens within an aperture of $2\theta_{E}$; 3. 3.

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})$; 4. 4.

we use pylens, a Python-based ray-tracing package that implements multi-plane 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 Figure 3.

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 LABEL:tab:lens for a summary of their general properties). Out of these fourteen 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 HS 0810+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). 2. 2.

CLASS B1359+154 and CLASS B1608+656 have multiple lensing galaxies within the Einstein radius and show strong flux-ratio anomalies (Rusin et al., 2001; Koopmans & Fassnacht, 1999; 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. 3. 3.

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. 4. 4.

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 a singular isothermal ellipsoid 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; 2. 2.

CLASS B0712+472 and CLASS B1555+375 both 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-fit models found by Hsueh et al. (2017, 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 prep). 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 elognated (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_{\rm sub}$, $M_{\rm hm}$ and the mass of the host galaxy, while the mean expected number of line-of-sight haloes is a function of $M_{\rm 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_{\rm 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 idealised CDM model (i.e. $M_{\rm 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}~{}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-sights 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_{\rm 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_{\rm 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 Figure 5, where the error bars represent the 1-$\sigma$ 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-$\sigma$ 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_{\rm sub}=0.008$; Xu et al., 2015, the value is recalculated to fit the definition of $f_{\rm 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).

Interestingly, from a sample of eleven 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 hydro-dynamical CDM simulations (Despali & Vegetti, 2017). Although Vegetti et al. (2014) infer a much smaller $f_{\rm 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 analyzed 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_{\rm sub}$ that had been inferred when fitting a single SIE model to these data.

From Figure 5, we notice that MG J0414+0534 is essentially an outlier, with a posterior probability that peaks at large values of $f_{\rm sub}$, outside our prior range. It is one of the most observed lensed quasars with a wide frequency coverage from radio, sub-millimetre 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-$\sigma$ 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 Figure 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-$\sigma$ level ($f_{\rm sub}$ = 0.008; Xu et al., 2015). Despali & Vegetti (2017) have quantified the suppression in the number density of substructure in hydro-dynamical 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-$\sigma$ 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 4 also shows that the expectation value of substructures and line-of-sight haloes are at the same level when $f_{\rm sub}=0.01$, for a CDM cosmology.

Figure 6 presents the comparison between the results from this paper and other studies that constrain $f_{\rm sub}$ from gravitational lensing. In this work, we define $f_{\rm sub}$ to be the substructure fraction within the aperture of $2\theta_{E}$. In Dalal & Kochanek (2002), $f_{\rm 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_{\rm 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\theta_{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_{\rm sub}$ scatter from $10^{-3}$ to $10^{-2}$, all of the estimates are marginally in agreement within the 1- to 2-$\sigma$ level.

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_{\rm 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_{\rm 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_{\rm 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_{\rm sub}=0.013\pm 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_{\rm sub}$ results with theoretical expectations. In the following, we present the marginalized constraints on $M_{\rm 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.

f

Figure 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-\sigma$ 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_{\rm hm}\leq 10^{5.5}~{}M_{\odot}$ 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_{\rm hm}=10^{9}~{}M_{\odot}$ as the number of detectable perturbers is equally consistent with zero for all models.

From Figure 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 toward 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_{\it th}>5.6$ keV at the same CL. In comparison with the latest 2-$\sigma$ constraints from the Ly$\alpha$ forest, $m_{\it th}>5.3$ keV or $m_{\it th}>3.5$ keV, depending on the assumption of the intergalactic medium thermal history (Iršič et al., 2017; Garzilli et al., 2018; Bolton et al., 2008), 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 et al., 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 over-estimate 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 hydro-dynamical 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, multi-wavelength observations at radio and sub-millimetre 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 multi-wavelength observations (e.g. Winn et al. 2004; Biggs et al. 2003; Mittal et al. 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 to 20 per cent for subhalo masses between $10^{6}$ and $10^{9}M_{\odot}$, 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; Green & van den Bosch, 2019; Errani & Peñarrubia, 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 realisations, and a detailed parametrisation as a function of the WDM model is still lacking.

Throughout the analysis, we have fixed the number density of all line-of-sight realisations to the average density of the Universe and thus neglected the fact that some lines of sight might be over- or under-dense. In a more stringent analysis, these realisation 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 parametrised the WDM mass function following the results of numerical simulations (Schneider et al., 2012b; Lovell et al., 2012). We extrapolated the fitting functions from these works below the resolution limit of the simulations ($M<10^{7}M_{\odot}$), 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 analyzed 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. 2. 2.

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_{\rm sub}$ into much closer agreement with the expectations from both CDM-only and hydro-dynamical 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). 3. 3.

The inclusion of extra complexity in the mass model of the lensing galaxies, although sub-dominant 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. 4. 4.

Under the assumption of a thermal relic dark matter model, we constrain the dark matter particle mass to be $m_{\rm th}>5.6$ keV at the 95 per cent confidence level. Our limits are in agreement with observations of the Ly$\alpha$ 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.

Appendix A Prior selection

Choosing a suitable prior is important in Bayesian inference. In this work we chose a prior that is uniform in $\log(f_{\rm sub})$, because it is shows the least information about the scale of a quantity. But in principle we could have chosen alternative forms for the prior.

A popular choice for priors is given by the Jeffreys prior (Jeffreys, 1946). Calculating the Jeffreys prior for the full Likelihood $\mathcal{P}(\vec{d}|M_{\rm hm},f_{\rm sub})$ requires a marginalization over data realizations $\vec{d}$, which implies a large number of model evaluations. As this is computationally prohibitive one might only consider the substructure part of the likelihood $\mathcal{P}(\vec{\theta}_{m},N|M_{\rm hm},f_{\rm sub})$. The Jeffreys prior in this case has the form:

[TABLE]

with $C\approx 2.3658\times 10^{10}$. Due to the large value of $C$ this can be approximated as: $\mathcal{P}(M_{\rm hm},f_{\rm sub})\approx M_{\rm hm}^{-1}f_{\rm sub}^{-1}~{}\mu(M_{\rm hm},f_{\rm sub})\propto\frac{1}{M_{hm}^{2}}$.

The advantage of the Jeffreys prior is that by construction every other parametrization of the mass function that is different to $(M_{\rm hm},f_{\rm sub})$ will have a Jeffreys prior that is related via a multiplication of the determinant of the Jacobian (while conjugate priors, for example, do not necessarily stay conjugate priors after a coordinate transformation). This property makes Jeffreys priors quite useful for parameters describing scales of some quantity.