The lensing properties of subhaloes in massive elliptical galaxies in sterile neutrino cosmologies

TL;DR

This study compares the lensing signals of subhaloes in massive elliptical galaxies under Cold and Sterile Neutrino dark matter models using high-resolution simulations, highlighting differences in subhalo properties and detectability.

Contribution

It provides the first detailed parametrization of subhalo properties and lensing signals in sterile neutrino cosmologies relative to CDM, aiding future observational tests.

Findings

Subhalo convergence follows a log-normal distribution.

Lensing signal in SN models is 10-80% lower than in CDM.

Approximately 30 lenses are needed to distinguish between models.

Abstract

We use high-resolution hydrodynamical simulations run with the EAGLE model of galaxy formation to study the differences between the properties of - and subsequently the lensing signal from - subhaloes of massive elliptical galaxies at redshift 0.2, in Cold and Sterile Neutrino (SN) Dark matter models. We focus on the two 7 keV SN models that bracket the range of matter power spectra compatible with resonantly-produced SN as the source of the observed 3.5 keV line. We derive an accurate parametrisation for the subhalo mass function in these two SN models relative to CDM, as well as the subhalo spatial distribution, density profile, and projected number density and the dark matter fraction in subhaloes. We create mock lensing maps from the simulated haloes to study the differences in the lensing signal in the framework of subhalo detection. We find that subhalo convergence is well described by a log-normal distribution and that signal of subhaloes in the power spectrum is lower in SN models with respect to CDM, at a level of 10 to 80 per cent, depending on the scale. However, the scatter between different projections is large and might make the use of power-spectrum studies on the typical scales of current lensing images very difficult. Moreover, in the framework of individual detections through gravitational imaging a sample of ~30 lenses with an average sensitivity of M_sub=5x10^7M_sun would be required to discriminate between CDM and the considered sterile neutrino models.

Tables

Table 1. Table 1: Summary of halo properties in the twelve hydro simulations: ID, halo mass M 200 ​ c subscript 𝑀 200 𝑐 M_{200c} and radius r 200 ​ c subscript 𝑟 200 𝑐 r_{200c} , stellar mass of the central galaxy M ∗ subscript 𝑀 M_{*} . In the next four columns we list the total number of subhaloes with subfind mass M > M min ≃ 3 × 10 7 ​ M ⊙ ​ h − 1 𝑀 subscript 𝑀 min similar-to-or-equals 3 superscript 10 7 subscript M direct-product superscript h 1 M>M_{\mathrm{min}}\simeq 3\times 10^{7}\mathrm{M_{\odot}}\mathrm{h^{-1}} , those with M > 10 8 ​ M ⊙ ​ h − 1 𝑀 superscript 10 8 subscript M direct-product superscript h 1 M>10^{8}\mathrm{M_{\odot}}\mathrm{h^{-1}} and the total number of ‘luminous’ satellites, given two different thresholds for stellar mass. For the L8 and L11 runs, we also show in brackets the subhalo abundances in the DMO runs. Finally, we have the stellar effective radius r ∗ , e subscript 𝑟 𝑒 r_{*,e} , calculated both from the 3D distribution of stars and averaging the projected distribution over different orientations – these last two quantities are calculated with the stellar particles belonging to the main galaxy and located within 300 (100) kpc from the halo centre. The masses are in units of M ⊙ subscript M direct-product \mathrm{M_{\odot}} , r 200 ​ c subscript 𝑟 200 𝑐 r_{200c} in units of comoving kpc, while the effective radii are expressed in physical kpc, similarly to the observational data in general. ID’s 1, 2, 3, and 4 were identified in Oppenheimer et al. ( 2016 ) as Grp008, Grp009, Grp005, and Grp002 in their Table 1.

ID	$M_{\text{200c}}$	$r_{\text{200c}}$	$M_{\ast }$	$N_{\text{sub}}$	$N_{\text{sub}}\left(M>{10}^8\right)$	$N_{\text{sub}}\left(M_{\ast }>{10}^6\right)$	$N_{\text{sub}}\left(M_{\ast }>{10}^7\right)$	$r_{\ast ,e}$(3D)	<$r_{\ast ,e}$(2D)>
	[$M_{\odot }$]	[kpc]	[$M_{\odot }$]					[kpc]	[kpc]
CDM
1	1.06 ${10}^{13}$	401	9.45 ${10}^{10}$	2209	482	84	25	12.76 (8.45)	9.31 (8.12)
2	1.05 ${10}^{13}$	401	1.00 ${10}^{11}$	2577	596	112	36	22.31 (13.31)	27.11 (13.93)
3	6.40 ${10}^{12}$	340	5.07 ${10}^{10}$	1436	333	61	30	6.24 (4.64)	5.80 (4.24)
4	3.99 ${10}^{12}$	290	5.48 ${10}^{10}$	829	189	30	16	6.46 (5.09)	5.15 (4.09)
L8
1	1.03 ${10}^{13}$	397	1.03 ${10}^{11}$	815 (818)	318 (377)	70	20	10.86 (7.44)	9.30 (7.00)
2	1.00 ${10}^{13}$	393	1.03 ${10}^{11}$	956 (1070)	368 (505)	80	31	12.99 (7.44)	12.46 (11.35)
3	6.72 ${10}^{12}$	339	7.22 ${10}^{10}$	507 (593)	200 (246)	48	24	4.43 (3.71)	3.64 (2.82)
4	3.93 ${10}^{12}$	288	6.70 ${10}^{10}$	306 (349)	126 (162)	27	13	7.03 (5.93)	6.32 (4.78)
L11
1	1.03 ${10}^{13}$	397	9.86 ${10}^{10}$	361 (369)	188 (230)	63	20	12.08 (7.83)	9.30 (7.00)
2	1.02 ${10}^{13}$	397	1.15 ${10}^{11}$	442 (486)	224 (286)	73	29	18.27 (12.85)	24.68 (12.26)
3	6.25 ${10}^{12}$	334	6.66 ${10}^{10}$	238 (239)	115 (148)	42	21	5.12 (4.17)	5.38 (3.43)
4	3.93 ${10}^{12}$	288	6.25 ${10}^{10}$	140 (168)	67 (111)	25	13	6.41 (5.19)	7.38 (4.37)

Table 2. Table 2: Best-fit mass function parameters from this and previous works (Schneider et al., 2012 ; Lovell et al., 2014 ) , following the parametrisation from Eq. 2 .

Mass function parameters
model	source	$\gamma$	$\beta$	$M_{\mathrm{hm}}\left[M_{\odot }\right]$
L8-DMO	this work	0.53	-1.3	$1.28\times {10}^8$
L11-DMO	this work	0.27	-1.3	$8.25\times {10}^8$
L8-HYDRO	this work	0.35	-1.3	$1.28\times {10}^8$
L11-HYDRO	this work	0.18	-1.3	$8.25\times {10}^8$
WDM(th. rel.)	Schn+12	1	-1.16	-
WDM(th. rel.)	Lov+14	1	-1.3	-
WDM(th. rel.)	Lov+14(sub)	2.7	-0.99	-

Table 3. Table 3: Average parameters – mean μ 𝜇 \mu and dispersion σ 𝜎 \sigma – of the Gaussian best-fit to the substructure convergence distribution from the 100 and 10 arcsec maps (see Figure 9 ) of the whole sample.

model	$\mu$ (100)	$\sigma$ (100)	$\mu$ (10)	$\sigma$ (10)
CDM	-4.27	0.67	-4.22	0.78
L8	-4.62	0.94	-4.63	0.95
L11	-4.89	1.10	-5.45	1.30

Equations

$$\frac{n_{\mathrm{WDM}}}{n_{\mathrm{CDM}}}={\left(1+\frac{M_{\mathrm{hm}}}{M_{\mathrm{sub}}}\right)}^{\beta },$$

$$\frac{n_{\mathrm{WDM}}}{n_{\mathrm{CDM}}}={\left(1+\gamma \frac{M_{\mathrm{hm}}}{M_{\mathrm{sub}}}\right)}^{\beta },$$

$$\eta (r)=\frac{\mathrm{d}\ln \rho }{\mathrm{d}\ln r}(r)\propto r^{1/\alpha }.$$

$$\rho (r)={\rho }_{-2}\exp \left\{-2\alpha \left[{\left(\frac{r}{r_{-2}}\right)}^{1/\alpha }-1\right]\right\}.$$

$$\Psi (\theta )=\frac{D_{\mathrm{LS}}}{D_{\mathrm{L}}{D}_{\mathrm{S}}}\frac{2}{c^2}\int \Phi (D_{\mathrm{L}},\theta ,z)\mathrm{d}z,$$

$${\Delta }_x\Psi (x)=2\kappa (x).$$

$$\kappa (x)=\frac{\Sigma (x)}{{\Sigma }_{\mathrm{crit}}},\mathrm{with}{\Sigma }_{\mathrm{crit}}=\frac{c^2}{4\pi G}\frac{D_{\mathrm{S}}}{D_{\mathrm{LS}}{D}_{\mathrm{L}}},$$

Full text

††pubyear: 2019††pagerange: The lensing properties of subhaloes in massive elliptical galaxies in sterile neutrino cosmologies–References

1 Introduction

Understanding and unveiling the nature of dark matter is one of the most long-standing challenges in modern astrophysics. According to the standard $\mathrm{\Lambda}$ Cold Dark Matter ($\mathrm{\Lambda}$CDM) model, dark matter constitutes the vast majority of the matter content in the Universe and, together with dark energy, accounts for 95 per cent of the total energy budget (Planck Collaboration et al., 2014). This model has been successful in explaining many aspects of structure formation and evolution, as well as in reproducing the density fluctuations in the early Universe with great accuracy (e.g. Bennett et al., 2013; Planck Collaboration et al., 2014). However, dark matter models are still untested at the small non-linear scales, due to numerical and observational limitations: a number of unsolved discrepancies exist between N-body simulations based on cold dark matter (CDM) and observations, such as the ‘too-big-to-fail’ and the ‘core-cusp’ problems (Klypin et al., 1999; Boylan-Kolchin et al., 2009; Bullock & Boylan-Kolchin, 2017). At the more fundamental level, a dedicated campaign to identify new fundamental particles that fit the requirements for supersymmetric dark matter has not yielded any definitive detections, either indirectly (e.g. the review of Gaskins, 2016), directly (e.g. Aprile et al., 2018) or via collider searches (ATLAS Collaboration, 2018; CMS Collaboration, 2018). Thus the combined tension with both astrophysical and particle physics experimental results necessitates the investigation of alternatives.

Some particle physics models solve the dark matter problem using particles that either evade, or even explain, the problems outlined above. Resonantly produced sterile neutrino dark matter is of particular interest for a series of reasons. It forms part of a well-developed extension of the standard model, called the neutrino minimal standard model ($\nu$MSM), which in addition to providing a plausible dark matter candidate also generates a mechanism to effect neutrino oscillations and baryogenesis (Asaka & Shaposhnikov, 2005; Laine & Shaposhnikov, 2008; Boyarsky et al., 2009). It explains all of these phenomena by introducing just three new, extra particles, which is the minimum required to explain neutrino oscillations and dark matter simultaneously, and could be tested/confirmed by the proposed SHiP experiment (Alekhin et al., 2016).

Sterile neutrinos are expected to decay into X-rays at a rate that is accessible to constraint/detection by X-ray observatories, particularly if the sterile neutrino mass is larger than 2 keV (Abazajian et al., 2001a; Abazajian et al., 2001b). An unexplained X-ray line detected at an energy of 3.55 keV in stacked observations of galaxy clusters (Bulbul et al., 2014), M31 (Boyarsky et al., 2014), the Galactic Centre (Boyarsky et al., 2015; Hofmann & Wegg, 2019) and the Milky Way (MW) halo outskirts (Boyarsky et al., 2018; Cappelluti et al., 2018) counts the decay of a 7.1 keV sterile neutrino among its possible sources. As one of the most promising, although not uncontentious (see Anderson et al., 2015; Jeltema & Profumo, 2016; Ruchayskiy et al., 2016), indirect dark matter detection signals, it constitutes a viable dark matter candidate, and is especially well suited to further study because the particle physics parameters that determine the X-ray decay signal also set the structure formation properties.

A set of complementary probes of the sterile neutrino dark matter cosmology have been derived from the fact that the sterile neutrino particle behaves as warm dark matter (WDM). The term WDM refers to a family of models, including thermal relics and sterile neutrinos, in which there is a cut-off in the linear matter power spectrum at dwarf galaxy scales for viable models; the cut-off scale is at least in part set by the the mass of the WDM particle, and in the specific case of resonantly produced sterile neutrinos can be specified uniquely, notwithstanding systematic uncertainty in the particle physics calculations, by the measured X-ray decay rate and emission energy (see Abazajian & Kusenko, 2019 for a discussion of additional sterile neutrino models). It has been tested with observations of local dwarf galaxies (e.g. Polisensky & Ricotti, 2011; Lovell et al., 2016; Schneider et al., 2017; Cherry & Horiuchi, 2017), Lyman-$\alpha$ forest measurements (Viel et al., 2005, 2013; Iršič et al., 2017; Garzilli et al., 2018), and reionization constraints (Bose et al., 2016) particularly for low-mass sterile neutrinos; these methods have proven effective at ruling out most resonantly produced sterile neutrino dark matter with a particle mass $<6$ keV.

From the observational point of view, techniques such as the detection of (sub)haloes through gravitational lensing (Vegetti & Koopmans, 2009a, b; Vegetti et al., 2010b; Vegetti et al., 2010a; Vegetti et al., 2012, 2014; Nierenberg et al., 2014; Hezaveh et al., 2016; Hsueh et al., 2019) or through gaps in the MW stellar streams (Erkal et al., 2016; Amorisco et al., 2016) are very promising for the understanding of the nature of dark matter. Gravitational lensing is sensitive to the whole mass distribution within the lens galaxies and along the line of sight (Despali et al., 2018; Vegetti et al., 2018): low-mass and non-luminous haloes can be detected via their gravitational effect on the observed lensed images.

The aim of this paper is to study the properties of $z~{}\sim 0.2$ giant elliptical galaxies and their subhalo populations in sterile neutrino models and make a comparison with CDM, with a particular focus on the implications for strong lensing. For this purpose, we ran high-resolution, hydrodynamical zoom-in simulations of analogues of lens galaxies at $z=0.2$ – the mean redshift of the SLACS lens sample (Bolton et al., 2006). This is the first time that hydro simulations have been run of massive ellipticals with sterile neutrinos, while previous works focused on dwarf galaxies and Local Group analogues (Lovell et al., 2017; Bozek et al., 2019).

This paper is structured as follows: in Section 2 we describe the simulations and the halo selection; in Section 3 we present the properties of the subhalo population in CDM and in the two sterile neutrino models, such as the parametrisation of the subhalo mass function, the subhalo density profiles and radial distribution. In Sections 4 and 5 we determine whether these dark matter models can be discriminated using gravitational lensing, by looking at different estimators. In Section 4 we compute the lensing signal of haloes and subhaloes by creating maps of lensing convergence and mock images, with two purposes: comparing the convergence distribution in Section 4.1, and measuring the subhalo power spectrum in Section 4.2. Finally, in Section 5 we use the formalism from Li et al. (2017) and Despali et al. (2018) to calculate the expected number of low-mass ($[10^{7}-10^{10}]$ $\mathrm{M_{\odot}}$) haloes as a function of mass, lens and source redshift, and observational sensitivity in order to determine if we can discriminate between these dark matter models with current or future lens samples. We summarize our findings, discuss the implications for substructure lensing in detail and draw our conclusions in Section 6.

2 Simulations

In this section we summarize the halo selection and the details of our simulations.

Our halo selection is derived from the common membership of two parent samples from the Ref-L100-N1504 EAGLE simulation halo catalogue (Schaye et al., 2015)111The EAGLE catalogues are publicly available at http://icc.dur.ac.uk/Eagle/database.php (McAlpine et al., 2016). Galaxies morphologies computed by Thob et al. (2019) are also available in the database.. The first parent sample is that of Despali & Vegetti (2017), who made a selection of haloes based on a series of properties that were important for matching observed lens galaxies. These halo properties were the mass (defined as $M_{\mathrm{200c}}$ with respect to the critical density), stellar mass, stellar effective radius and velocity dispersion: observed lens galaxies at $z=0.2$ (such as the SLACS sample) follow a linear relation in the stellar effective radius – total stellar mass plane, as derived by Auger et al. (2010). Moreover, the selected galaxies were required to have a disc-to-total mass ratio lower than 0.5 in order to select early-type (ETG) morphologies. The second parent sample was assembled by Oppenheimer et al. (2016) from the same cosmological simulation box. They identified ten group-mass haloes that were suitable for resimulation using zoomed initial conditions, i.e. isolated and of the right mass for their science goals, and subsequently resimulated these at 8-times higher mass resolution than the parent simulation volume. Four of these ten haloes were found to also be a part of the Despali & Vegetti (2017) selection. These were selected for resimulation, creating the first sample of massive elliptical galaxies simulated with sterile neutrino cosmologies and baryonic physics.

The sterile neutrino dark matter cosmology provides for a rich variety of linear matter power spectrum shapes and cut-offs, many of which have already been ruled out by a combination of MW-M31 satellite counts (Polisensky & Ricotti, 2011; Kennedy et al., 2014; Lovell et al., 2014; Lovell et al., 2016) and Lyman-$\alpha$ forest studies (Viel et al., 2005, 2013; Iršič et al., 2017). It is not feasible to run enough simulations to probe in full all of the available parameter space; we therefore take our inspiration from the recently detected 3.55 keV line, one explanation of which is the decay of a sterile neutrino with a mass of 7.1 keV and a neutrino mixing angle $\sin^{2}{2\theta}=[2,8]\times 10^{-11}$. This range corresponds to a production lepton asymmetry, $L_{6}=[9,11.2]$ (Lovell et al., 2016). The $L_{6}=9$ model produces a matter power spectrum significantly warmer than that of CDM, and $L_{6}=11.2$ warmer again than $L_{6}=9$. The coldest matter power spectrum available for a sterile neutrino is attained for $L_{6}=8$, which was found to be consistent with the dark matter decay interpretation of the 3.55 keV line reported by Bulbul et al. (2014) and Boyarsky et al. (2014) but was later ruled out by more stringent limits from Ruchayskiy et al. (2016). In order to take account of any possible systematic uncertainty in the calculation of the matter power spectra for these models, we adopt $L_{6}=8$ rather than $L_{6}=9$ as the colder limit for our study, and adopt $L_{6}=11.2$ for the warmer limit. We refer to these models hereafter as L11 and L8. In both cases the sterile neutrino mass is 7.0 keV rather than 7.1 keV; we do not expect the results for a 7.1 keV particle would be significantly different from those derived here.

The simulations are performed in exactly the same way as their CDM counterparts (Oppenheimer et al., 2016). The galaxy formation model is that of the EAGLE project (Schaye et al., 2015; Crain et al., 2015), which features cooling, star formation, stellar evolution and feedback, associated with the formation of stars and the growth of black holes. We use the version of the model in which the model parameters are recalibrated for the dark matter particle mass of $1.21\times 10^{6}\mathrm{M_{\odot}}$, known as RECAL, as do the equilibrium-chemistry simulations in Oppenheimer et al. (2016). The code is a heavily modified version of P-Gadget3 (Springel et al., 2008a) and uses a pressure-entropy formulation of SPH (Hopkins, 2013). Haloes are identified using the friends-of-friends (FoF) algorithm and are subsequently split into subhaloes using the subfind halo finder (Springel et al., 2001). The cosmological parameters are consistent with the constraints by the Planck Collaboration et al. (2014): $\Omega_{0}=0.307$, $\Omega_{\Lambda}=0.693$, $\Omega_{\mathrm{b}}=0.04825$, $h_{0}=0.6777$, $\sigma_{8}=0.8288$ and $n_{\mathrm{s}}=0.9611$. The only manner in which these simulations differ from their CDM versions is that the initial conditions have been remade with the L11 and L8 matter power spectra, which are shown in fig. 1 of Lovell et al. (2017); like the CDM simulations, our runs were simulated with mass resolution eight times higher than the simulation box from which they were sourced. Moreover, for the sterile neutrino models, we ran both full-hydro and dark matter-only (DMO) versions for each of the four volumes.

Figure 1 shows images of one of our haloes at $z=0.2$, which present qualitatively the differences between the CDM, L8 and L11 models at different scales. In the top row we show projections of the dark matter distribution that are 2 physical Mpc on a side of the CDM (left-hand column), L8 (middle column) and L11 (right-hand column) halo 3 runs. The second row shows instead the central 200 pkpc. In the third row the colour scale represents the lensing convergence – defined as the projected mass density distribution normalized to the critical density, see Section 4 – in a box of 80 pkpc ($\simeq$20 arcsec) on a side, created by using all of the matter within $r_{\mathrm{200c}}$, i.e. baryonic and dark matter combined. Finally, the last row shows mock lensing images, created by placing a source at $z=0.6$: lensed images appear in the very central part, on a scale of the order of $\simeq$2 arcsec/$\simeq$ 10 kpc – or $<3$ per cent of the virial radius.

The lowest mass subhalo – i.e. the bound subfind mass, $M$ – that is resolved in our simulations is $M\simeq 3\times 10^{7}\mathrm{M_{\odot}}\mathrm{h^{-1}}$ (20 particles, which is the default minimum for subfind), but a robust identification is possible only at $M\geq 10^{8}\mathrm{M_{\odot}}\mathrm{h^{-1}}$. At this scale, WDM simulations suffer from spurious fragmentation of filaments, due to the presence of a resolved wavenumber cut-off in the initial power spectrum. These fragments may then be identified as subhaloes by the subhalo finder, artificially increasing the abundance of the low mass objects. This is a purely numerical problem, which is not solved sufficiently by increasing the spatial resolution, but which might be alleviated by future N-body codes that use phase space smoothing techniques (Angulo et al., 2014) or instead apply an adaptive softening (Hobbs et al., 2016). At present, the only solution is to exclude them from the analysis using empirical criteria: in this work we follow the method developed by Lovell et al. (2014), which enables us to eliminate spurious subhaloes based on their peak mass and their Lagrangian shape in the initial conditions, with specific thresholds for each sterile neutrino model (see their work for details on how these are derived). In what follows, we always make use of subhalo catalogues from which the spurious subhaloes have been removed. In all cases, the total mass fraction in spurious subhaloes with respect to $M_{200c}$ is less than 0.0002 and thus does not affect the results.

2.1 Properties of the host haloes

As can be seen from Table 1, the properties of the main haloes and also of the central galaxies are very similar in all three models (CDM, L8 and L11). The values of $M_{200c}$ and the total stellar mass differ by only a few per cent, in agreement with the results from Lovell et al. (2017). The main difference between these dark matter models is indeed the number of subhaloes, as we will discuss in the next section. However, we caution that any change in the galaxy formation model will change both the overall matter distribution (Springel et al., 2018) and thus the disruption rate of satellites (Richings et al., 2018), which has to be borne in mind when interpreting our results. Recently, Ludlow et al. (2019) showed that small galaxies and the inner parts of galaxies in hydrodynamical simulations might be affected by numerical artefacts due to the energy transfer from dark matter to stars when two species of particles with unequal mass are present (equipartition), as is the case in our runs. This has an effect on galaxy sizes and implies that the softening is not a reliable estimate of the minimum reliable size. In principle, this kind of numerical effect could have an effect on the estimate of the lensing signal from haloes and subhaloes, which originates in the central parts of galaxies, as well as on the estimate of density profiles and in general halo properties close to the centre. The softening length of our simulations is $\epsilon=0.350$ physical kpc and thus an effective resolution limit is $2.8\times\epsilon=0.98$ pkpc. We mark these values in all the relevant figures in order to delimit the parameter space where the results have to be interpreted with particular care.

In the main body of the paper we focus on the comparison between CDM and sterile neutrino models, in particular on the relative numbers of subhaloes and thus their relative impact on the number of satellites: any numerical artefacts due to the equipartition effect likely affect all three sets of simulations to a similar degree, and thus we expect that our results on the relative contribution of subhaloes are not strongly affected.

Mukherjee et al. (2019) studied the effect of different feedback mechanisms on the lensing signal of ETGs. They used a suite of $\Lambda$CDM hydrodynamical simulations, run with different variations of the EAGLE galaxy formation model, to measure the density slope at the Einstein radius and the projected mass-density relation, and how these depend on different implementations of the stellar and AGN feedback. They argue that, through a detailed comparison with observed samples, it is possible to use gravitational lensing to constrain feedback models. Here we can complement their findings, identifying that, since the main halo properties are not systematically affected by the change in the dark matter model (see Table 1), the lensing properties of host haloes cannot necessarily be used to distinguish CDM from sterile neutrino models (at fixed baryonic physics model).

Another property that we can measure is the dark matter fraction of the haloes – which has been measured in hydrodynamical CDM simulations (Schaller et al., 2015; Lovell et al., 2018) – in order to see how it differs in sterile neutrino models. Figure 2 shows the dark matter fraction, as a function of distance from the centre, for the four simulated haloes (lines of different colour). For each halo, the dark matter mass fraction is hardly distinguishable in the three models (solid, dashed and dotted lines of the same colour) for $r>0.05r_{\mathrm{200c}}$. A small suppression ($\simeq$ 20 per cent), indicating a smaller dark matter fraction at the very centre of the haloes in both sterile neutrino models, is apparent in the mean distribution; however, we caution that we do not have enough statistics to provide reliable estimates, given the halo-to-halo fluctuation. We can conclude that the overall mass distribution of the four haloes is similar in CDM and sterile neutrino models. This is consistent with the fact that their properties, such as total mass, stellar mass and radius, agree very well (see Table 1). Moreover, it supports the fact that the mass accretion is very similar in the three models, but in CDM a higher percentage of the accreted mass is in the form of clumpy structures. As shown by Schaller et al. (2015), the dark matter fraction increases with halo mass in the central region of the halo, which explains the differences between the four haloes.

3 Subhalo properties

In this section we analyse four dark matter-dependent subhalo population properties that could play a role in setting the gravitational lensing signal of subhaloes: the mass functions, the spatial distribution, the density profiles and the dark matter fraction in subhaloes.

3.1 Mass function

The formation of low-mass ($<10^{10}$ $\mathrm{M}_{\odot}$) dark matter haloes is suppressed due to the cut-off in the initial matter power-spectrum, and the subhalo number density may experience an additional reduction due to the interaction with the host halo, which could feasibly vary between models due to changes in the mass–concentration relation. Nevertheless, previous works have assumed that (to the first order approximation) the degree of suppression is the same or similar for isolated haloes and subhaloes. In particular, Schneider et al. (2012) parametrised the ratio between thermal relic WDM and CDM isolated halo number density as a function of mass as

[TABLE]

where $M_{\mathrm{hm}}$ is the half-mode mass scale as defined by the ratio of the WDM linear matter power spectrum relative to CDM, and $\beta$=-1.16. For the case of MW subhaloes, Lovell et al. (2014) found a better fit with $\beta=-1.3$ or with a slightly different parametrisation:

[TABLE]

where $\beta=-0.99$ and $\gamma=2.7$. It is important to point out that most studies that try to distinguish between different dark matter models using the subhalo counts (e.g. Birrer et al., 2017; Li et al., 2017; Despali et al., 2018; Vegetti et al., 2018; Hsueh et al., 2019) use the parametrisation from Schneider et al. (2012) for any kind of WDM model. This choice is motivated by the fact that recent works focused on thermal relic WDM candidates (e.g. Viel et al., 2005; Lovell et al., 2012; Lovell et al., 2014; Bose et al., 2016), for which this parametrisation was originally developed. As discussed in Section 2, sterile neutrino models are intrinsically more complex, given that they are characterized by the sterile neutrino mass $m_{s}$ and the lepton asymmetry $L_{6}$. The combination of the two determines $M_{\mathrm{hm}}$, while for thermal relic candidates this is determined by the particle mass alone. $M_{\mathrm{hm}}$ has a value of $1.28\times 10^{8}\mathrm{M}_{\odot}$ ($8.25\times 10^{8}\mathrm{M}_{\odot}$) for LA8 (LA11). For the Viel et al. (2005) fit to the thermal relic power spectrum these $M_{\mathrm{hm}}$ correspond to thermal relic particle masses of 4.3 keV and 2.5 keV for LA8 and LA11 respectively.

Note that for the purpose of the fitting formula in Lovell et al. (2014), we should multiply the ‘true’ values of $M_{\mathrm{hm}}$ by the $\gamma$ parameter to derive an ‘effective’ half-mode mass; in this case our LA8 and LA11 models correspond to ‘effective’ thermal relic masses of 6.0 and 4.2 keV respectively. However, it is also the case that the sterile neutrino matter power spectra show a greater variety of cut-off gradients than thermal relic particles, and therefore even knowing $M_{\mathrm{hm}}$ is not sufficient to define the subhalo mass function uniquely. Moreover, we stress that different particle physics codes predict different sterile neutrino momentum distribution functions and consequently different matter power spectra. We performed an approximate re-calculation of our matter power spectra using the public particle physics code published by Venumadhav et al. (2016) together with the class Boltzmann solver (Blas et al., 2011), which returned $M_{\mathrm{hm}}$ for LA8 (LA11) that were $\sim 6$ ($\sim 3$) times higher than those found for our spectra; we defer a detailed comparison to future work.

We fit the ratio between the subhalo mass functions in our simulations with the functional form from Equation (2). Figure 3 shows the subhalo mass function from our DMO (top panel) and hydro (middle panel) runs. The three black curves show the best fit to the subhalo mass function for each dark matter model. In order to obtain a good fit, $\gamma$ has to change for each sterile neutrino model: the best fit values are summarized in Table 2. The CDM mass function with the same normalization and with slope ($\alpha=0.9$ for the dark matter-only run and $\alpha=0.85$ for the hydro run, as derived in Despali & Vegetti, 2017 for the EAGLE hydro simulation) is shown by the dotted line. The lower panel of Figure 3 shows the ratio with respect to the CDM DMO scenario. In the same panel, the grey circles and triangles show the subhalo mass function predicted according to the mass function fit by Schneider et al. (2012), which underestimates the number of subhaloes for our models by up to 30 per cent 222The difference between our sterile neutrino models and the fit from Schneider et al. (2012) might reside in the combination of different aspects: $(i)$ the shallower power-spectrum cut-off in sterile neutrino models with respect to thermal relic WDM; $(ii)$ differences between the halo and the subhalo mass functions; $(iii)$ numerical effects due to different ways of removing the spurious subhaloes (of secondary importance).. The same parametrisation works well for the subhalo population at other redshifts (not shown here), once one takes into account the evolution of the normalisation of the mass function with time (Giocoli et al., 2008).

The resolution of our simulations allows us to probe the subhalo mass function reliably only down to $M\simeq 10^{8}\mathrm{M_{\odot}}\mathrm{h^{-1}}$, where it is suppressed by $\simeq$ 50 (80) per cent in the L8 (L11) model relative to the CDM subhalo mass function but is still mostly flat. However, extending the functional forms from Figure 3 (solid and dashed black lines) to lower masses would result in a sharp decrease at $M<10^{7}\mathrm{M_{\odot}}\mathrm{h^{-1}}$. It remains to be shown definitively whether this drop-off rate describes the sterile neutrino cosmology accurately, and confirmation will require still higher resolution simulations.

Finally, it is important to remember that the number of luminous satellites is another important probe, since any viable dark matter model must be able to reproduce their observed abundance. As shown in Lovell et al. (2016), some sterile neutrino models are already ruled out by satellite counts in MW-like haloes. Given that the total mass of our haloes is larger than for the MW, it is worth investigating the abundance of satellites when they are separated into luminous and dark sub-populations. In Table 1 we list the total number of subhaloes that have stellar masses $M_{*}>10^{6}\mathrm{M_{\odot}}\mathrm{h^{-1}}$ (lowest resolved stellar mass) or $M_{*}>10^{7}\mathrm{M_{\odot}}\mathrm{h^{-1}}$: it is easy to see that the total number are slightly lower in the WDM models, but comparable to CDM333We use these stellar masses only in order to provide a rough estimate of dark vs. luminous satellites; as these masses are at the resolution limit of our simulations, the exact values should not be used for interpreting our results any further.. This again supports the importance of strong lensing as a method to discriminate between different dark matter models, given its ability to detect the total mass and thus the majority of WDM subhaloes.

3.2 Radial distribution

As shown by Bose et al. (2017), the subhalo radial distribution might be different in WDM models with respect to CDM. In particular, even though less numerous, subhaloes in WDM scenarios can be more centrally concentrated than in CDM. In the left-hand panel of Figure 4 we plot the mean radial number density of subhaloes as a function of radius for our three scenarios, rescaled to the mean number density within $r_{200c}$ in each model. The solid, dashed and dotted curves show the best-fit Einasto profiles to the mean number densities. As expected, the subhalo population is more centrally concentrated than in CDM. For all three cases $r_{-2}\simeq 0.54$ $r_{200c}$, while the logarithmic slope is different for each model and decreases for warmer models. If the sterile neutrino models are normalized by the CDM mean values (red dashed and dotted lines in both panels), it is apparent how the subhalo population is suppressed in warmer models.

However, since gravitational lensing is sensitive to projected quantities, we also need to test to what degree this difference in the three dimensional radial distribution is preserved in projection. In the right-hand panel of Figure 4, we show the projected number density of subhaloes as a function of radius, normalized by its value at $r_{200c}$, by averaging over three different projections for each halo. We see that, even though the normalized number density is still higher for the sterile neutrino models, the mean distributions flatten close to the centre (i.e. inner most $\simeq$30 kpc) and the central number densities are of the same order of magnitude for the three models.

3.3 Density profiles

Another recognised effect of WDM is to produce shallower density profiles for the low-mass subhaloes, in the same mass range in which the number of structures is already suppressed (e.g. Lovell et al., 2012; Ludlow et al., 2016). It has been shown (Springel et al., 2008b) that subhaloes in CDM can be well described both by NFW (Navarro et al., 1996) and Einasto (Einasto, 1965) profiles. Here, we use the second since it provides a more flexible parametrisation which might be needed in warmer models.

Thus, we select the subhaloes within $r_{200c}$ and fit their profiles with Einasto profiles, characterized by a power-law logarithmic slope:

[TABLE]

In this model, the integrated density profile is commonly written as

[TABLE]

We perform the fit by optimizing for the three parameters ($\rho_{-2},r_{-2},\alpha$) for each subhalo, where the first two are the density and the radius at which $\rho(r)\propto r^{-2}$ and $\alpha$ defines the steepness of the power-law. Instead of fitting the individual subhalo profiles – which could be more prone to uncertainties – we stack subhaloes of the same mass and we fit the median profile.

$\rho_{-2}$ is 20 (40) per cent lower in the L8 (L11) model, while $r_{-2}$ is roughly 10 (20) per cent larger at fixed $M_{\mathrm{sub}}$. We estimate the subhalo concentration $c=r_{\mathrm{max}}/r_{-2}$, i.e. the ratio between the subhalo size (or the distance of the farthest particle) and $r_{-2}$. $r_{\mathrm{max}}$ is on average the same in the three models, as a function of mass and, at fixed mass, of the same order of magnitude of the virial radius.

In Figure 5, we show the ratio of the concentration–mass relations of sterile neutrino models relative to CDM. The blue and red lines show the ratio for the L8 and L11 scenarios for each mass bin, while the mean suppression obtained by fitting individual profiles is shown by the dotted lines.

We then generate the density profile that would be predicted from the medians of these distributions for four bins in subhalo mass in which we have enough statistics. In Figure 6 we plot the ratio between the predicted density profile in the L8 (upper panel) and L11 (lower panel) models and CDM. The median parameters from the sterile neutrino models yield central densities ($r<1$ kpc) that are systematically lower than the CDM ones, to 30 and 50 per cent at the very centre, consistently for all mass bins.

3.4 Dark matter fraction in subhaloes

Many previous works studied the dark matter fraction in subhaloes (Gao et al., 2004a; Springel et al., 2008a; Xu et al., 2015; Despali & Vegetti, 2017). In WDM models, the halo concentrations are lower (Macciò et al., 2013; Ludlow et al., 2016, – and as in the previous section) and the number of dark subhaloes is suppressed, which can lead to a lower fraction of the total halo mass being located in dark matter subhaloes.

In Section 2 we discussed the total dark matter fraction of the host haloes; we now calculate the dark matter fraction in subhaloes. Its value depends not only on the distance from the centre, but also on halo mass (Gao et al., 2004b; Springel et al., 2008a). For this reason, and since our haloes span almost one order of magnitude in mass, we show the ratio of the sterile neutrino and CDM fractions instead of the absolute value, and present the results in Figure 7. As expected from the suppressed mass function, in sterile neutrino models the fraction is lower than in CDM. The suppression rises from an order of magnitude in the centre to roughly 20 per cent at the virial radius. As shown in previous works (e.g. Springel et al., 2008a; Despali & Vegetti, 2017), the number density of subhaloes increases towards the centre, in a way that is well represented by an Einasto profile (Ludlow et al., 2009), more slowly than the total dark matter density and thus most subhaloes are found in the outer parts of the halo. Given that these are present in similar numbers in all the three runs, the dark matter fraction in subhaloes is similar in the outskirts of the halo. On the other hand, the low-mass subhaloes, which can be found closer to the centre, are more abundant in CDM than in sterile neutrino models, and thus the fraction in subhaloes is suppressed towards the centre, even though the number densities of subhaloes centrally is similar between CDM and WDM.

The two sterile neutrino models give similar predictions and cannot be distinguished with only four simulated systems. However, Lovell et al. (2014) also reported a similar suppression for the slightly different thermal relic WDM model in the MW halo context. Finally, we measure the projected subhalo mass fraction and the projected number of subhaloes as a function of the distance from the centre, thus bringing together the findings from previous subsections. In this case, we use different projections for each halo in order to gain more statistics. Xu et al. (2015) and Despali & Vegetti (2017) have shown that the projected number density of subhaloes as a function of radius is constant at small distances from the centre ($r<0.3r_{\mathrm{200c}}$), while it slowly decreases outside. Moreover, the projected number density in CDM is inversely proportional to the subhalo mass, i.e. low-mass subhaloes are more abundant than high mass haloes. This behaviour is obviously modified in sterile neutrino models. Figure 8 shows the projected number density of subhaloes per kpc*-2* as a function of distance from the centre, focusing on a specific bin in subhalo mass in each panel. The projected density of sterile neutrino models (dashed and dotted lines for L8 and L11, respectively) are of the same order of the CDM ones (solid lines) for the two higher mass bins ($M=10^{9-9.5}\mathrm{M_{\odot}}$), while they are 50 per cent lower for $M=10^{8-8.5}\mathrm{M_{\odot}}$ and 70 per cent lower for the lowest bass bin available in our simulations ($M=10^{7.5}\mathrm{M_{\odot}}$).

4 Lensing properties of subhaloes

In this section, we analyse the lensing properties of the subhalo population in detail. We look at the lensing signal directly by creating mock images and lensing convergence maps, and subsequently measuring the distribution of subhalo convergence and the corresponding power spectrum.

We use the lensing code GLAMER (Metcalf & Petkova, 2014; Petkova et al., 2014) to run ray-tracing simulations and generate lensing maps of the simulated haloes.

For any extended distribution of matter, the lensing effect depends on the projected mass distribution and so is characterized by the effective lensing potential, which is in turn obtained by projecting the three-dimensional potential $\Phi$ on the lens plane:

[TABLE]

where $D_{\text{L}}$, $D_{\text{S}}$ and $D_{\text{LS}}$ are, respectively, the angular diameter distance of the observer to the lens, the observer to the source, and between the lens and the source; $c$ is the speed of light. The gradient of the lensing potential gives the scaled deflection angle, while the Laplacian gives the convergence $\kappa$:

[TABLE]

The lensing convergence is defined as a dimensionless surface density and so effectively corresponds to a scaled projected mass density, characterising the lens system. It can be written as

[TABLE]

where $\Sigma_{\mathrm{crit}}$ is the critical surface density. The value of the lensing convergence in practice determines by how much the background sources appear magnified on the lens plane. For non-spherical mass distributions, the lensed images of the source are also stretched and distorted along privileged directions by the shear, $\gamma$. Finally, the critical curves identify the location on the lens plane in which the magnification is ideally infinite (see the white contours in Figure 1 for an example).

GLAMER calculates the deflection angles, shear and convergence by a tree algorithm, representing each simulation particle with a B-spline in three dimensions as is commonly done in smooth particle hydrodynamics (SPH) simulations. The size of the particles is set to the distance to the $N$th nearest neighbour where $N_{\mathrm{smooth}}$ can be adjusted; in this work $N_{\mathrm{smooth}}=64$. This smoothing scheme provides higher resolution where the particles are dense and higher resolution is therefore justified, and lower resolution where the particles are sparse and shot noise would otherwise be a problem. We set the lens and source redshift respectively to $z_{\mathrm{l}}=0.2$ and $z_{\mathrm{s}}=0.6$. For the sterile neutrino runs, we exclude the particles belonging to spurious subhaloes.

We ran the ray-tracing code through 200 random projections for each halo, in order to obtain a good statistical sample. We look at the lensing convergence generated by subhaloes ($\kappa_{\mathrm{sub}}$). We exclude those subhaloes that have a significant stellar component ($M_{*}>10^{7}\mathrm{M_{\odot}}\mathrm{h^{-1}}$ - see Section 3.1), as we are interested in the dark subhaloes that cause perturbations to the lensing signal but cannot be detected through their light. This is in practice equivalent to removing the same number of massive subhaloes from the maps and has only minor consequences for the distribution of convergence, since it only influences the high-$\kappa_{\mathrm{sub}}$ tail, but not the shape or mean of the distribution.

In particular, we focus on maps of 100 arcsec on a side (thus containing the whole radius of the haloes in projection) and of 10 arcsec on a side, slightly larger than the typical lensing images. The two sets of maps have the same resolution of 0.1 arcsec (similar to that of HST data) and the same depth along the line of sight (i.e. $2\times r_{\mathrm{200c}}$). In this way we want to see how much of the ‘true’ signal from subhaloes can be extracted from individual lensing images, which enclose a region that is much smaller than the virial radius on the plane of the sky (see Figure 1 for an example of the relevant scales). The top panels of Figure 9 show an example of CDM total convergence maps of the two considered fields-of-view. In the left-hand panel, the subhaloes can be identified by eye, while the central region on the right is dominated by the main galaxy and only one extra clump is clearly visible, proving once more the challenge of finding subhaloes with lensing.

We proceed by analysing the distribution of the subhalo lensing convergence and its power spectrum. Previous works have tried to describe the effect of subhaloes on lensed arcs as the presence of residuals (i.e. the difference between the observed lensed images and those predicted by the best lens model) in the form of a Gaussian random field (Bayer et al., 2018; Chatterjee & Koopmans, 2018). Moreover, a number of recent works focused on the subhalo power spectrum, measured from the lensing convergence, using analytical (sub)haloes or numerical simulations comparing cold, warm (thermal relic) models (Cyr-Racine et al., 2018; Brennan et al., 2018) or also self-interacting dark-matter models (Díaz Rivero et al., 2018).

Previous works analysed the dependence of the subhalo power spectrum on both the maximum subhalo mass and the size of the halo sample (and number of projections). While they use analytical density profiles, here we are instead intrinsically limited by the resolution and the number of the simulations. Thus, we can neither characterize fully the variance due to differences between the haloes, nor explore the effect of the substructure population below $M\simeq 10^{8}\mathrm{M_{\odot}}\mathrm{h^{-1}}$ reliably. Moreover, the sterile neutrino models used in this work are colder than the model used in their work and so the difference in the power spectra is unlikely to be as pronounced. Nevertheless, we can explore the variance due to different projections and – since the same haloes have been run with different dark matter models – estimate how the power spectrum changes in sterile neutrino models. Also, previous works have not explored the effect of a smaller field-of-view on the power-spectrum measurement. We discuss the this approach and its limitations in the following paragraphs.

4.1 Lensing convergence of subhaloes

The middle and bottom panels of Figure 9 show the logarithmic distribution of $\kappa_{\mathrm{sub}}$, for the most and least massive haloes in our sample respectively. The points show the median distribution of each model, while the shaded regions enclose the 0.25 and 0.75 quantiles.

In all of the panels, the solid lines show the best fit Gaussian distributions to the data, which reproduce the median well at the high-$\kappa$ end and around the peak, but fail at the low-$\kappa$ tail in all cases. This implies that the overall subhalo convergence distribution is well represented – in particular in the CDM case – by a Log-normal function in real space. By comparing the two panels, it is interesting to note that, while the variation between different projections becomes much larger in the 10 arcsec maps, the position of the peak remains the same for the CDM and L8 models and changes only slightly for L11 (Table 3). The two haloes differ by almost one order of magnitude in mass (see Table 1): the difference in the subhalo convergence distribution is however much smaller, indicating that the measure is quite stable with mass. In Table 3, we list the mean values for the mean and dispersion of the best-fit Gaussian through the whole sample for both sets of maps.

Since subhaloes are a diffuse component, and not as localized as the main galaxy convergence, it could still be possible to measure the true shape of their convergence distribution from small fields-of-view, assuming that the main lens convergence can be measured correctly. However, since the subhalo convergence is not observed directly, but reconstructed as a byproduct of the observational analysis by comparing the real data with the reconstructed lens model, it is also an approximation to assume that any small scale feature is caused by subhaloes (Bayer et al., 2018; Vegetti et al., 2014; Ritondale et al., 2019b).

4.2 Convergence power spectrum

In general, the convergence power spectrum $P(\ell)$ of subhaloes is characterized by a number of features: first of all, the normalization depends on the overall convergence in subhaloes, which is by construction the main difference between the three models that we consider. As already shown in Brennan et al. (2018), the normalization of the subhalo power-spectrum is dominated by the highest mass subhaloes and thus removing them increases the relative difference between different models, as happens when excluding the ‘luminous’ subhaloes from our calculation. As a test, we computed the power spectra using maps that contain all the subhaloes and we also find that in this case the median $P(\ell)$ is essentially identical across the three models. Moreover, since CDM and WDM models differ the most at the lower mass end of the subhalo mass function, it is a natural consequence that the relative difference increases when considering only low-mass subhaloes.

Figure 10 shows the power spectra measured within 100 (left) and 10 arcsec (right), calculated using a module from MOKA (Giocoli et al., 012a). Each thin coloured line shows the measurement for one projection/halo, while the solid, dashed and dotted lines show the median values. The lower panels show the ratio of each median power-spectrum to the CDM equivalent. The discriminating-power of the $P(\ell)$ seems to suffer from the limited field of view more than the convergence distribution does, as the measurements from different 10 arcsec maps scatter much more than those obtained from the 100 arcsec maps.

We find that while the large scale modes ($\ell<0.03~{}\mathrm{arcsec}^{-1}$) are substantially identical for all the projections and for all three models, the power spectrum varies significantly at intermediate and small scales. The mean power spectrum is lower for the sterile neutrino models with respect to CDM, showing a signature of the lack of small scale structures. We note that not only is the normalization different in the three scenarios, but also the slope of the power-spectrum for $k\geq 0.05$. The high-$k$ slope is determined by the density profile and thus the difference between the models is an indication of the different properties of the subhalo profiles.

However, the scatter from different projections is larger than the difference between the medians and values from the different models overlap, as in previous works (Brennan et al., 2018). This suggests that it is virtually impossible to distinguish between the three models by analysing only a few lensing images. The scatter is smaller in the right-hand panel of Figure 10 than in the smaller field of view, but it is relevant in both cases.

Our findings suggest that, using the power spectrum, it might be possible to distinguish between CDM and very warm WDM models, but not the relatively cold sterile neutrino models presented in this work. As the differences increase at progressively lower masses, higher resolution simulations would be required to measure the power spectrum of subhaloes with $M\leq 10^{8}\mathrm{M_{\odot}}$, which might be more significant. At the same time, we anticipate that future observations, that will allow us to detect those low masses – and at the same time identify with more precision the most massive subhaloes, and subsequently to include them in the lens model – will prove more important for power-spectrum studies than current data sets.

5 Detection expectations

As shown in Despali et al. (2018) and Li et al. (2017), subhaloes are only a fraction of the total number of low-mass haloes that contribute to the convergence in the strong lensing images: isolated field haloes, located along the line of sight, constitute at least 50 per cent of the total number of perturbers for any combination of source and lens redshift. Thus, they play a fundamental role in constraining the dark matter model and have to be taken into account. Here we provide approximate detection expectations through gravitational lensing including both field haloes and subhaloes. We assume that the halo mass function suppression in sterile neutrino models is the same as that of subhaloes, as in many previous works (Schneider et al., 2012; Lovell et al., 2014; Despali et al., 2018; Vegetti et al., 2018). We leave an explicit verification of this assumption to a follow-up paper. We use the Sheth & Tormen (1999) CDM halo mass function parametrisation, with the best fit parameters appropriate for the Planck cosmology as calculated by Despali et al. (2016). We then use the best fit parameters $\gamma$ and $\beta$, obtained by fitting Equation 2 to the ensemble of the sterile neutrino runs (see Figure 3 and Table 2) and we use them to estimate the total expected number of perturbers, i.e. low-mass subhaloes and field haloes with masses $M_{\mathrm{min}}\leq M\leq 10^{10}\mathrm{M_{\odot}}$. $M_{\mathrm{min}}$ is defined according to the lowest detectable mass from observational data, as detailed in the following paragraph.

We provide predictions for a representative sample of observed lenses, which reproduce HST data in terms of point-spread-function, resolution and signal-to-noise level appropriate for single-orbit observations - or artificially improved versions of the same data, as discussed below. We consider a sample of $\simeq 30$ lenses, with lens and source redshifts in the range $0.2\leq z_{\mathrm{l}}\leq 0.5$ and $1\leq z_{\mathrm{s}}\leq 2.5$. This resembles the configurations of state-of-the art samples, such as the SLACS (Bolton et al., 2006; Vegetti et al., 2014, 2018) and BELLS GALLERY lenses (Shu et al., 2016; Ritondale et al., 2019a, b). The mean lowest detectable mass in the considered sample is $<M_{\mathrm{low}}\simeq 5\times 10^{9}\mathrm{M_{\odot}}$, while its minimum is $M_{\mathrm{low,min}}\simeq 10^{8}\mathrm{M_{\odot}}$ (possible only in less than one per cent of the relevant area – the distribution lowest detectable mass pixel-by-pixel for realistic data can be see in Ritondale et al., 2019b).444We remind the reader that the data sensitivity in this case corresponds to the lowest detectable mass of (sub)haloes that can influence the surface brightness distribution of the lensed arcs. In particular, a detection is generally considered reliable if the lens model including a secondary clump is preferred to the smooth lens model with a difference in log evidence of $\Delta\log E\geq 50$ (Vegetti & Koopmans, 2009a; Vegetti et al., 2010b, 2012, 2014). The exact value of the lowest detectable mass varies pixel-by-pixel, depending on the position of the secondary pertuber in the image plane.

In order to make some predictions and overcome the current data quality limitations, we look at the distribution of detectable subhaloes and field haloes for each sample, given three different mean values for the lowest detectable mass: the realistic one ($5\times 10^{9}\mathrm{M_{\odot}}$) and two improved values, one and two orders of magnitude lower ($5\times 10^{8}\mathrm{M_{\odot}}$ and $5\times 10^{7}\mathrm{M_{\odot}}$ ). These would in practice correspond to HST data with a higher signal-to-noise or longer exposure times, given that the resolution is kept fixed. Figure 11 shows what are the masses most likely to be detected in each model (different panels) and for these three cases (different lines). In the first panel, the CDM distribution moves almost self-similarly to lower masses when $M_{\mathrm{low}}$ decreases, consistently with the constant slope in the halo and subhalo mass function. On the other hand, for sterile neutrino models (second and third panels), masses below $M\simeq 10^{8}\mathrm{M_{\odot}}$ are strongly suppressed and thus the peak remains above this value in all cases; moreover, the shape of the distribution for lower $M\simeq 10^{8}\mathrm{M_{\odot}}$ changes, with a longer and lower low-mass tail. At the same time, while the three solid histograms are practically identical, the dotted ones present clear differences: at this sensitivity, the three models could be distinguished. This is consistent with the results from Vegetti et al. (2018) and Ritondale et al. (2019b), who found that with currently data it is not possible to put stringent constrains on the nature of dark matter: they can only safely exclude sterile neutrino models with $m_{s}<0.8$ $\mathrm{keV}$, equivalent to $M_{\mathrm{hm}}>10^{12}\mathrm{M_{\odot}}$.

Increasing the number of lenses – at fixed data quality – would result in a linear increase in the total number of detections, but would not alter the distributions in Figure 11. On the other hand an improvement in the lowest detectable mass – i.e. moving to lower values of $M_{\mathrm{low}}$ – would allow us to reach the regime where the three models are different from each other, and thus distinguishable. We plan to extend these results, accurately modelling the effect of resolution and signal-to-noise, in a follow-up paper.

6 Discussion and conclusions

In this work, we used zoom-in simulations with cold or sterile neutrino dark matter to explore the properties of the subhalo population and to characterize the differences between these models in the framework of strong gravitational lensing. We considered two different sterile neutrino models, which could be related to the production of the 3.55 keV unexplained X-ray line observed in clusters of galaxies (Bulbul et al., 2014), M31 (Boyarsky et al., 2014), the Galactic Centre (Boyarsky et al., 2015) and the Milky Way halo outskirts (Boyarsky et al., 2018; Cappelluti et al., 2018). In both models, the particle mass is $m_{s}$=7 keV, but the lepton asymmetry has two different values – $L_{6}=8,11.2$ – which result in a different suppression of low-mass structures (see Figures 1 and 3). We resimulated four haloes hosting massive elliptical galaxies - i.e. typical lens galaxies at redshift $z\simeq 0.2-0.5$ - which have never studied in sterile neutrino hydrodynamical simulations. Here we summarise our findings and discuss the implications for gravitational lensing studies.

Gravitational lensing, and in particular the gravitational imaging technique (Vegetti & Koopmans, 2009a) has the unique ability of detecting dark low-mass haloes, which are not visible as satellites of the main galaxy (or as low-mass galaxies along the line of sight). Through this technique it is possible to detect individual haloes (Vegetti et al., 2010b, 2012; Hezaveh et al., 2016), but also to put constraints on the halo and subhalo mass functions (Vegetti et al., 2014, 2018; Ritondale et al., 2019b). For this reason, the results from simulations and theoretical model are of fundamental importance to interpret observational results.

We started by looking at the properties of the four main haloes in the three different models. We find that the main halo properties, such as total and stellar mass, virial and effective radius, are very similar in the three models. Moreover, the dark matter fraction only differs in the innermost part of the halo ($r\leq 0.05r_{\mathrm{vir}}$), where it is 20 per cent lower in sterile neutrino models. However, with only four haloes, we do not have enough statistics to model this effect precisely. The most important difference remains the number of dark low-mass subhaloes, which are suppressed in sterile neutrino models, while the number of luminous satellites is similar (see Table 1).

Then, we studied the properties of the substructure population, in terms of mass function (Sec. 3.1), spatial distribution (Sec. 3.2 and 3.4) and density profile (Sec. 3.3). Sterile neutrino models are more complex than other WDM candidates, due to the combined effect of $m_{s}$ and $L_{6}$ changing both the matter power spectrum cut-off wavenumber and slope in a non-trivial fashion. For this reason, a two parameter fit is needed to accurately reproduce the subhalo mass function. In particular, the best-fit expression from Schneider et al. (2012) underestimates the number of low-mass subhaloes in our simulations. We provide the best fit parameters for the L8 and L11 subhalo mass functions in Table 2. The fact that the subhalo mass function parametrisation depends on the specific sterile neutrino model (as Equation 2) has an important implication: even measuring the subhalo mass function directly from observational data, as in Vegetti et al. (2018), might not be sufficient to constrain the sterile neutrino model correctly, given the degeneracy between $M_{\mathrm{hm}}$ and $L_{6}$. However, some assumptions can be made to reduce the parameter space: for example, restricting the sterile neutrino mass to 7.1 keV, as we do here, would allow us to put constraints on $L_{6}$ for this particular model, or the other way around. In this perspective, a more precise parametrisation (such as the values in Table 2) will help the comparison with observational results.

The normalization of the subhalo mass function scales with redshift and host halo mass, but the exact number of subhaloes can still vary from system to system. It has been shown that strong lensing studies that try to constrain the projected dark matter fraction in subhaloes $f_{\mathrm{sub}}$ (within a certain radius) and the WDM half-mode mass $M_{\mathrm{hm}}$ simultaneously (e.g. Vegetti et al., 2018; Gilman et al., 2019; Hsueh et al., 2019) find that these two quantities are correlated: similar results can be obtained by increasing $f_{\mathrm{sub}}$ in a warmer model, or decreasing $f_{\mathrm{sub}}$ in a colder dark matter model, as the two contributions balance. Previous observational (Vegetti et al., 2014, 2018) and numerical (Xu et al., 2015; Despali & Vegetti, 2017) works have measured $f_{\mathrm{sub}}$, finding consistent values for CDM of the order of $f_{\mathrm{sub}}\simeq 0.004-0.006$ for subhaloes of mass $10^{6}\leq M_{\mathrm{sub}}\leq 10^{9}\mathrm{M_{\odot}}$. The presence of baryonic physics suppresses the number of subhaloes, by different amounts depending on the details of the implementation, which can help to exclude unrealistic physical models (Despali & Vegetti, 2017). Given that the same baryonic physics model has been used for all the runs, its effect is of the same order (Lovell et al., 2017). In this work we calculated the projected number of subhaloes as a function of the distance from the centre, both for CDM and sterile neutrino models, as has been done previously for CDM alone (Xu et al., 2015; Despali & Vegetti, 2017). These are flat (as a function of radius) in the inner $\simeq$30 kpc and lower in the sterile neutrino models, with a difference from CDM that increases at low subhalo masses, i.e. for $M\leq 10^{8.5}\mathrm{M_{\odot}}$. These results are broadly consistent with previous works based on WDM simulations, both in the form of sterile neutrinos (Lovell et al., 2016) and thermal relics (Schneider et al., 2012; Lovell et al., 2012; Lovell et al., 2014; Ludlow et al., 2016). However, none of these has addressed specifically the host haloes of massive early-type galaxies in sterile neutrino hydrodynamical simulations, as we do in this work. We then measured the dark matter fraction in subhaloes as a function of radius, finding that it is lower towards the centre for the sterile neutrino models than in CDM, similarly to what was found by Lovell et al. (2014) for thermal relic WDM models.

Finally, we measured the subhalo density profiles and calculated the best-fit Einasto profile. We found that the profiles are shallower towards the centre in sterile neutrino models, leading to central densities suppressed by 10 to 50 per cent. For this reason, the subhalo concentrations are lower than in CDM by up to 30 per cent, for subhaloes with masses $M\geq 3\times 10^{8}\mathrm{M_{\odot}}$.

In Section 4, we looked at the differences in the lensing signal of subhaloes between the three models by analysing the distribution of subhalo lensing convergence and its power spectrum. We created convergence maps by ray-tracing through 200 random projections for each halo; for this, we used the particles belonging to non-luminous (and non-spurious) subhaloes. The subhalo population contributes to the total lensing convergence at different scales and in general the subhalo convergence is qualitatively well fitted by a Log-Normal distribution – with departures from gaussianity at the edges – in all models. In Table 2, we list the best-fit parameters of the Gaussian distributions that fit the subhalo convergence in logarithmic space. These distributions describe the subhalo convergence well both when the measurement is done with the subhalo population of the entire halo (i.e. with convergence maps of $\simeq$100 arcsec on a side), and when considering only the inner part of the halo in projection ($\simeq$ 10 arcsec), closer to the location of the lensed images.

Similar information can be obtained from the power spectrum of subhalo convergence: the total power is lower in sterile neutrino models (Figure 10), as well as the slope at $k\geq 0.1$ $\mathrm{arcsec}^{-1}$: this is due to the different properties of the subhalo profile and to the different relative contribution of high and low mass subhaloes, since the former are present in similar numbers. One limitation of the power spectrum approach is that the variation between different projections is larger – or at best of the same order of magnitude – of the intrinsic differences between the power spectrum in the three models. However, the mean and median power spectra show a consistent difference, both when looking at the whole subhalo population and at a 10 arcsec field of view close to the halo centre. In this last scenario, the mean CDM power-spectrum remains consistent with the measurement done on the population within the virial radius, while the sterile neutrino models present a much larger variation and a lower mean – simply due to the lack of structures on small scales in part of the projections. The power-spectrum is dominated by the high-mass subhaloes and thus removing them would increase the difference between the model. In order to measure the power due to the smallest structures alone, in the analysis of observational data one needs to identify the more massive subhaloes and explicitly include them in the model. Moreover, higher resolution observations are needed to probe the low-mass end of the mass function, where the models differ significantly from each other and to overcome the scattering due to different projections and a limited field-of-view. At the same time, the resolution of our simulations does not allow us to measure the power spectrum of subhaloes with mass $M\leq 10^{8}\mathrm{M_{\odot}}$ reliably.

Finally, in Section 5 we follow the approach from Despali et al. (2018) and Li et al. (2017) and provide detection expectations for a realistic sample of observed lenses, which resembles the configuration of the SLACS (Bolton et al., 2006; Vegetti et al., 2014) and BELLS GALLERY (Shu et al., 2016; Ritondale et al., 2019a, b) lenses. We consider both the current data-quality, in terms of resolution and signal-to-noise – with an average detection limit around $M\simeq 5\times 10^{9}\mathrm{M_{\odot}}$ – and artificially improved data which would allow to detect haloes and subhaloes with masses one or two order of magnitude lower. We assume that both the halo and subhalo mass function are suppressed similarly, following Equation 2; we used the best-fit parameters from Table 2 and calculate the number of effective perturbers (both subhaloes and isolated haloes along the line of sight) for the two samples, following the approach from Despali et al. (2018). We calculate the distribution of detectable masses, comparing its peak and shape across the three scenarios. We find that, at the current data sensitivity, it is not possible to put significant constraints on the nature of dark matter or distinguish between CDM and the two sterile neutrino models considered here, consistentl with previous works (Vegetti et al., 2018; Ritondale et al., 2019b). However, if the data sensitivity increased by one or two order of magnitude in mass – i.e. if we could detect (sub)haloes of mass $M\simeq 5\times 10^{8}\mathrm{M_{\odot}}$ or $M\simeq 5\times 10^{7}\mathrm{M_{\odot}}$ – the same samples would allow to distinguish between CDM and sterile neutrino models.

In our simulations, we have used the same baryonic physics model in all the runs. In this way, we can focus solely on the effect of sterile neutrinos on the abundance of small-scale structures, since the effect of a fixed baryonic model on the abundance of subhaloes is very similar in the three cases (as shown also in Lovell et al., 2017). However, it is important to point out that there is a large amount of uncertainty in the effect of different baryonic physics models on the mass function of subhaloes (Despali et al., 2017) and on their survival rates (Kelley et al., 2019). This causes a degeneracy between the effect of baryons and WDM, which would require a full set of hydro simulations that combines variations in both parameter spaces to be fully understood. Future observations with instruments such as E-ELT and VLBI will have sufficient sensitivity to probe the halo mass function below $<10^{7}M_{\odot}$; in this regime the effect of different dark matter models is expected to play the dominant role in shaping the mass function, and thus the uncertainties in the baryonic physics modelling will no longer be an issue.

To conclude, in this work we presented the analysis of hydroynamical simulations with sterile neutrino dark matter. The chosen models – in both cases sterile neutrino with a mass of $m_{s}$ = 7 keV – are among the few alternative dark matter models with a possible observational motivation (i.e. the unexplained 3.5 keV line Bulbul et al., 2014; Boyarsky et al., 2015; Cappelluti et al., 2018). They are, however, colder than most of the warm dark matter models explored so far in simulations (Schneider et al., 2012; Lovell et al., 2012; Lovell et al., 2014; Li et al., 2016). We provide detailed parametrisation of the differences between these sterile neutrino models and CDM. We conclude that the difference in the number of low-mass subhaloes leaves a clear signature on the distribution of lensing convergence and in the subhalo power spectrum. These are, however, not yet detectable with current samples of observed lenses. Future observations, more numerous and with a higher sensitivity (which could be provided by an improved signal-to-noise, longer observational times or a higher resolution) will provide a promising way of discriminating between CDM and sterile neutrino models.

Acknowledgements

We thank the EAGLE collaboration for allowing the use of the EAGLE model in these simulations and Joop Schaye for useful comments. This work was carried out on the Dutch National e-Infrastructure with the support of SURF Cooperative. MRL is supported by a COFUND/Durham Junior Research Fellowship under EU grant 609412, and also acknowledges support by a Grant of Excellence from the Icelandic Research Fund (grant number 173929$-$051). 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). RAC is a Royal Society University Research Fellow. Support for BDO was provided through the NASA ATP grant NNX16AB31G. This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K00042X/1, ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure. We thank Ben Metcalf for giving us access to the lensing code GLAMER and Carlo Giocoli for providing the MOKA libraries. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration, 2013), of the matplotlib (Hunter, 2007) and NumPy packages.