Exploiting flux ratio anomalies to probe warm dark matter in future large scale surveys

The Cold Dark Matter paradigm (CDM) is one the of the most successful models in cosmology. The existence of a massive, non-relativistic particle that interacts only via gravity allows us to explain a wide variety of astronomical observations, including the distribution of galaxies over scales that span many orders of magnitude (de la Torre et al., 2013; Anderson et al., 2012, e.g.). Despite the success of this model, the lack of any detection of new particles in terrestrial experiments around the weak scale means that, astronomically, our ability to gain further information on the particle nature of CDM further is limited (LUX Collaboration et al., 2013; Aprile et al., 2012). As such we are diversifying our search, looking for new signatures that might lead us away from the CDM paradigm and give us telling insights in to its nature.

Extensions to the CDM paradigm are becoming increasingly commonplace. Models of dark matter that invoke a self-interaction (Harvey et al., 2018b, a; Spergel & Steinhardt, 2000, e.g.), or assume an ultra-light state (Hu et al., 2000; Hui et al., 2017, e.g) or models that do not assume that dark matter is generated at non-relativistic velocities (Bode et al., 2000; Boyarsky et al., 2019; Kusenko, 2009) have become popular predicting discriminate signatures. In this letter we will study the imprint of two popular Warm Dark Matter (WDM) models on the expected flux ratios observed in strongly lensed quasars, identifying what key observation will we require in order to significantly detect or rule out these particular models.

In this paper we will concentrate on WDM. Interest in WDM has grown since the detection of an unidentified X-ray emission line in clusters of galaxies (Bulbul et al., 2014; Boyarsky et al., 2014; Franse et al., 2016; Urban et al., 2015) that is consistent with a 7 keV sterile neutrino (Boyarsky et al., 2009a; Shi & Fuller, 1999a; Asaka et al., 2005; Laine & Shaposhnikov, 2008; Boyarsky et al., 2009b). Unlike CDM, WDM is produced relativistically at radiation-matter equality, generating very different growth of structure at scales with wavenumber $k>1$ $\mathrm{Mpc}/h$. With higher particle velocities, WDM is able to free-stream out of the smallest scale perturbations, suppressing structure. The mass-scale of this suppression depends on the mass of the dark matter particle, $m_{\chi}$ (Shi & Fuller, 1999b; Abazajian et al., 2001; Asaka et al., 2007; Ghiglieri & Laine, 2015; Venumadhav et al., 2016). Thus if we are able to probe the dark matter haloes down to $\sim 10^{8}$, we will be able to constrain WDM as a plausible dark matter candidate.

Our goal in this paper is to generate an observable that will allow us to distinguish clearly between CDM and WDM. We will focus on two specific models of WDM; however, the results have much wider implications than just these models. Regardless of whether the dark matter is a thermal relic or is generated in some more complicated process, as long as it shares as a similar primordial transfer function to those that we use here the results will be applicable.

When the geodesics from distant sources intersect compact objects they are bent and deformed. In the rare event that the intervening compact object is dense enough it can cause the geodesics to split, resulting in multiple images of the same source. Strong gravitational lensing occurs at many different scales, for example distant galaxies can be lensed by massive foreground clusters of galaxies, which allows for the precise measurement of the total mass in their cores (e.g. Jauzac et al., 2018) and the study of the nature of dark matter (e.g. Harvey et al., 2018b). It is also observed on galaxy scales where background light from a distant galaxy is lensed by a foreground one, distorting the galaxy into a giant arc surrounding the foreground lens. As this light from the distant source is lensed, it is slightly perturbed by the small-scale structure along the line of sight and therefore can allow for the sensitive measurement of galaxy structure (e.g. Nightingale et al., 2019; Gilman et al., 2017), including any small substructures and hence where it may be sensitive to the dark matter model (Vegetti et al., 2018; Metcalf & Zhao, 2002; Metcalf & Madau, 2001; Metcalf & Silk, 1999; Gavazzi et al., 2007; Nierenberg et al., 2017). This manifests itself as small deviations in the arcs and can be directly observed.

In a similar way, distant quasars can experience the same distortion, where the light from a point source is split into discrete lines by a foreground galaxy, resulting in two or four images of the same quasar. In this situation perturbations along the line of sight cause the flux of each image to deviate from what is expected of a smooth foreground halo model. Given that the flux ratios between different images can be predicted very precisely, any anomalies in these ratios can be attributed to small-scale structure (Metcalf & Madau, 2001; Metcalf & Silk, 1999).

Measurements of flux anomalies have been already been observed in quadruply imaged quasars (Metcalf & Zhao, 2002; Morgan et al., 2006; Hsueh et al., 2019), predicting that in order to account for the data there must be $\sim 5-10$ per cent of the lensing mass in a substructure. However, recently it was shown that assuming an oversimplified initial lens model for the foreground galaxy can result in a biased estimate of the amount of substructure required to account for the data (Xu et al., 2015); similarly a baryonic disk could mimic flux anomalies (Hsueh et al., 2018). As a result, in this study we adopt a method that takes a very different approach to the modelling of the foreground lens, which assumes only that the cosmological simulations we use are a good representation of the population of observed lensing galaxies. Integrating over all possible lenses, lens and source redshifts, we would be able to produce a prediction of what a large scale survey would observe should it observe a complete sample of multiply imaged quasars. This way observations of flux ratios can be directly compared to the simulations with no modelling required.

We solve for strong lensing assuming the projected densities from simulations are at a single redshift, and hence we use the thin-lens approximation (Bartelmann, 2010). We initially obtain a solution to the lens equation:

which is the deflection map $\alpha_{i}(\theta_{i})$ that maps lens-plane positions $\vec{\theta}$ to source plane positions $\vec{\beta}$. The solution is obtained using a Fast Fourier Transform (FFT) as a poisson solver. The deflection angle is then given by the gradient of an obtained lens potential:

where the lens potential is

where the $D$ variables are the angular diameter distances between the observer, the lens (L) and the source (S), and $\Phi$ the Newtonian 3D potential. Finally the magnification is given by the laplacian of the potential, being equal to the jacobian of the coordinate transformation: as a consequence of the FFT choice, we have taken care that there is sufficient empty space around the lens, in order not to get unphysical results owing to the periodic boundary conditions.

This deflection map then needs to be numerically inverted (using simply bi-linear interpolation), to obtain a regularly spaced sampling of the source plane positions. Specifically, we invert the vertices of the lens plane and find those cells in the source plane that lie within each image plane cell. In this way we are implicitly simulating a finite-sized source, which is by default the pixel size of the image plane ($100$pc). Although this is much larger than the size of a quasar, we are limited by the mass resolution of the simulation, and any higher resolution source size would provide no extra information in the results. As such we assume that there are no structures smaller than this $100$pc scale that would impact the flux ratios. We therefore have at each position in the source plane, one or more lens-plane positions at which the source will be observed by the observer behind the lens.

From the inverted map, we collect all the pixels that have exactly two or three (with the central image) images on the lens plane. We discard the central images, and we discard quadruply or more lensed source pixels. The remaining data are hence our bare samples of doubly imaged source plane pixels, with their associated magnifications. In other words, we take a uniform prior on the position of a source in the source plane, by sampling all positions and assigning them the same probabilistic weight.

For a given lens at a fixed redshift $z_{\rm lens}$, we are interested in its probability distribution of flux ratios for all possible source positions. We therefore populate the source plane and then nominally integrate this plane to infinity. In practice this is impossible and therefore choose to integrate to some maximum source redshift. However, at high redshift $dV/dz$ (for volume $V$) tends to zero and hence choosing a high-redshift cut-off will not affect the results. In the presence of a actual survey, it would be trivial to integrate over the redshift distribution of the quasar sample; however, here we are agnostic about the survey choice and therefore integrate out to a redshift of $z_{\rm source}=8.0$. Similarly, when we combine the statistics from lenses across multiple lens redshifts, we weight the distribution from each redshift according to the analytic integral over comoving volume, such that the number of flux ratios for a given lens represents the total observable volume at the redshift of that lens.

We will thus focus this paper on doubly imaged quasars. This is motivated two-fold. First, with the limited number of simulated haloes, concentrating on doubly imaged quasars means we can construct a simple observable and garner large statistics. Second, and similarly, in any large scale surveys the number of doubles will outnumber quads by a factor of $\sim 20$. We found in this study that there were 12, 22 and 26 times more doubly imaged source positions than quadruply imaged ones for CDM, $L_{6}=8$ and $L_{6}=11$ respectively. The premise of this method is in measuring the flux ratios of many quasars, and therefore relies on large numbers of quasars in any future survey (e.g. Laureijs et al., 2011). Throughout this paper we will measure the ratio between the fainter and brighter quasar image,

It is known that quasar flux varies on timescales of $\sim 10$ days (MacLeod et al., 2012). When a variable quasar is strongly lensed and the light beams are split the differential length path of each geodesic results in a time delay in each photon arrival. As a result, at any given time there may exist an anomalous flux ratio that is due to the delayed arrival of quasar variability, which may mimic small structures. In order to incorporate these into our predictions we first state that the magnitude change due to quasar variability, $\Delta m$ is

where $f^{\prime}$ is the varying flux and $f$ is the mean flux. Hence the impact on a flux ratio, $F$, would be:

Following this, we can work out the impact on the probability distribution of the final flux ratio, $P_{F}$,

where $P_{F}(F)$ is the flux ratio PDF in an unvarying quasar and $P_{\Delta m}$ is the PDF for magnitude change between two images. Given that:

it follows that:

We find that the resulting probability distribution is just a convolution of the intrinsic flux ratio PDF and the PDF of the varying quasar magnitude. To calculate the final PDF we assume that the quasar variability PDF, $P_{\Delta M}$, can be modelled by Gaussian (in magnitude space), with zero mean. We calculate the final PDF for three standard deviations of $\Delta m$=0, (unvarying) $\Delta m=0.5$ (expected) and $\Delta m=1$ (Millon, 2019). Since the quasar variability will smear out any signal, a larger $\Delta m$ constitutes as a more conservative estimate of our ability to distinguish between CDM and WDM. In the left-hand panel of Figure 1 we show the predicted CDFs for the flux ratios of CDM for these three cases of varying quasar flux, with $\Delta m=0$ (green line, no variability), $\Delta m=0.5$ (orange line, realistic variability) and $\Delta m=1.0$ (black line, strong variability).

In addition to the non-linear distortion by the lens galaxy, the light rays from a lensed quasar will also be perturbed by intervening structures that happen to lie along the line-of-sight. Unlike the distortions from the lens, these perturbations are linear and considered in the “weak” regime.

Moreover, even though the lensing is strong, we are still in the regime where a point source remains a point source (i.e.infinitesimally small) after lensing. As a consequence, the various sources of magnification are multiplicative, and the order of multiplication does not matter. For a single image, its observed flux $\mathcal{F}_{\rm obs}$ can be decomposed into three components,

with $\mathcal{F}_{\rm intrinsic}$ the intrinsic flux of the source, $\mu_{\rm lens}$ the magnification due to the strong lens (computed as per Section 3.1), and $\mu_{\rm LoS}$ the magnification due to the line-of-sight structure.

For the computation of the line-of-sight contribution, we will assume that there is a negligible probability of having a structure in the line of sight that is of a comparable mass to the main lens such that compound strong lensing occurs. Indeed, compound lenses are almost absent from current samples of strong lenses. In this case, all structures in the line of sight will have a smaller effect, which justifies the weak lensing approximation. Moreover we assume that since the two images of the quasar are separated by the approximately the size of a galaxy, any halo that is larger than this will affect the two fluxes equally and hence have no affect on the flux ratio. As such we therefore consider a mass function between $10^{13}M_{\odot}<M_{\rm los}<10^{4}M_{\odot}$. Although this mass will be redshift dependent, by taking a large mass it will smooth out any differences and cause us to overestimate the number of lenses required to discriminate. As such in this sense it is a conservative limit.

In order to calculate the contribution of lensing by haloes along the line-of-sight we first derive a probability density distribution of magnification by large scale structure, $P_{G}(z_{\rm source})$. To do this we use TurboGL (Kainulainen & Marra, 2009a, b, 2011), a stochastic approach to cumulative weak lensing that models the line of structure through the halo model. We modify the input mass function of TurboGL in order to simulate the effect of line of sight structures in a WDM cosmology. In order to do this we adopt the WDM correction,

where we fit the two free parameters $\beta$ and a characteristic mass $M_{\rm s}$ (Schneider et al., 2012; Lovell et al., 2014) (note that this is not the equivalent half-mode mass since this shape of these two functions are very different). We find that L8 $(z=0.5)$ $M_{\rm S}=1.6\times 10^{9}$ and $\beta=0.35$, L8 $(z=0.2)$ $M_{\rm S}=1.5\times 10^{9}$ and $\beta=0.37$, L11 $(z=0.5)$ $M_{\rm S}=4.1\times 10^{9}$ and $\beta=0.55$, L11 $(z=0.2)$ $M_{\rm S}=1.5\times 10^{9}$ and $\beta=0.52$. We convert the CDM mass functions in TurboGL into L8 and L11 equivalents using the fits found here and calculate the weak lensing probability distribution by line-of-sight structure. The central panel of Figure 1 shows examples of $P_{G}$ at redshift $z=8$ with the ratio to CDM in the bottom panel; the example of $z=8$ is chosen to show a maximal difference between each PDF. We note that given that the estimate of the mass function is carried out over the simulation box, which itself is an overdense region, means that we may bias our estimate of $M_{\rm s}$ and $\beta$. We therefore test the sensitivity of our final distributions to the two parameters and find that a difference of 2 orders of magnitude in $M_{\rm s}$ results in a difference of $\sim 0.1$% in the PDF.

In this study we have estimated the PDF of flux ratios in a way that is independent of any survey specific parameters. It is subsequently trivial to take the results and apply them, for example, to a given source redshift distribution of observed quasars. With our distributions we are able to estimate the number of flux ratios that would need to be observed in a survey of specified volume and depth in order to rule out these particularly models of WDM.

The first step is to draw $N$ flux ratios randomly from the WDM PDF, draw a further $N$ flux ratios from the CDM PDF. We then calculate the likelihood of each sample, given the mean predicted by the CDM PDF, assuming a Poisson Probability Mass Function. i.e.

where $k_{i}$ is the observed number of flux ratios from a given dark matter model in the $i$th flux ratio bin, $F_{\rm n}$ is the number of flux ratio bins and $\lambda_{i}$ is the expected number from CDM for the same $i$th flux ratio bin. Since we randomly draw from the WDM and CDM PDFs, we Monte Carlo this test 1000 times. Following this step we compare the likelihoods by estimating the mean change in the Bayesian Information Criterion, BIC (with the same number of degrees of freedom):

In other words we are using Bayesian statistics to work out how large $N$ needs to be for the WDM and CDM distributions to be unambiguously different. Figure 3 shows the predicted $\Delta$BIC for each WDM-CDM model pair. It is generally accepted that a $\Delta$BIC$>10$ is very strong preference for a model. Thus we estimate the number of flux ratios whereby the $\Delta$BIC$>10$ for three different values of quasar variability $\delta m=0$ (solid line), $\delta m=0$ (dashed line) and $\delta m=1.0$ (dotted line). We find that in the case of no quasar variability we would require $\sim 300-600$ flux ratios, in a realistic quasar variability scenario we would require almost double this, and then for a sample of strongly varying quasars we would require $\sim 1000$ flux ratios. These numbers are definitely within reach of forthcoming surveys (e.g. Laureijs et al., 2011; Ivezic et al., 2008) and hence should be considered a promising test of dark matter.

In this letter we have measured the expected flux ratio PDF from double-imaged strong lensed quasars. We identified quasar variability and line of sight structures as two of the main sources of systematic uncertainty associated with this technique. From this we have estimated the number of flux ratios required to rule out the two WDM models in question here. However, there exist further systematics that still need to be addressed.

The primary theoretical systematic error will be associated with the “modelling” of the halos in the N-body simulations. Sub-grid modelling of the baryonic processes will affect the inner distribution of matter, altering the distribution of observed flux ratios. Hsueh et al. (2018) showed that what was once interpreted as subhalos causing anomalous flux ratios, were more likely to be over-simplification of the lens models used. Here we model the ensemble distribution of lenses using state-of-the-art simulations, not individual ones. However, in the event that these lenses do not reproduce the true lenses, the comparison between observed and simulated distribution will be biased. We therefore clearly state here that in order for this statistical method to be profitable we require simulations that can produce all the known features of a true lens (e.g. baryonic disks). Hence, it will be important when comparing to observations that the nature of these uncertainties are understood.

Furthermore, it will be important to ensure the selection function of observed lenses matches those that are simulated. It is not clear that the selection function of observed lensed quasars is indiscriminate, and we hypothesise that a particularly type of galaxy will be more efficient at lensing than another. Therefore it may not be sufficient to integrate over all galaxies as we have done here, and thus future work will need to identify which type of galaxy, with which set of parameters, is more likely to lens a quasar and incorporate those in to these findings.

Finally, we have not attempted to model the contribution of dust emission and/or absorption to the lensing maps (Trayford et al., 2015). In the presence of dust within the lens, a single image may appear de-magnified whilst another that has no intervening dust may appear brighter, mimicking a flux ratio. Although important to be considered in observations, modelling it here is beyond the scope of this paper.