© The Authors 2022

1. Introduction

One of the most groundbreaking findings during the last century was the discovery of a mass excess in nearby galaxies that could not be explained by the amount of ordinary (luminous) matter that was found to exist in those galactic systems. The observation of flat rotation curves in spiral galaxies (e.g., Rubin et al. 1980; Bosma 1981; Corbelli & Salucci 2000), the discrepancy between observed velocity dispersion measurements and those predicted by the virial theorem in elliptical galaxies and globular clusters (e.g., Zwicky 1933; Faber & Jackson 1976), and the presence of collapsed structures at high redshift (e.g., Gunn & Gott 1972) were incontrovertible evidence for the existence of a new extraordinary form of matter in the Universe, called dark matter (DM), which neither emits nor absorbs radiation. The DM hypothesis has led to several correct predictions and has explained many observational discrepancies that had emerged in the past from the comparison of observational data with the first cosmological scenarios that were dust-radiation-only models. Despite the great success of the DM model, its nature remains unknown, making it one of the most fundamental unsolved questions in physics.

The most widely accepted scenario for the origin of DM is the so-called cold dark matter (CDM), a part of the standard ΛCDM cosmological model, that has been remarkably successful in explaining the properties of a wide range of large-scale observations, including the accelerating expansion of the Universe (Perlmutter et al. 1999), the power spectrum of the cosmic microwave background (CMB), (Page et al. 2003), and the observed abundances of different types of light nuclei (Cyburt et al. 2016). However, the ΛCDM paradigm still presents some discrepancies with observations, mostly at small scales. Such small-scale challenges include, among others, the “cusp-core” problem, the missing satellites problem, the too-big-to-fail problem, and the angular momentum catastrophe (for a review, see Bullock & Boylan-Kolchin 2017; see also Perivolaropoulos & Skara 2022).

An appealing solution to those problems is to modify the intrinsic properties of DM particles. During the past few years, numerous DM alternatives and ΛCDM extensions have been proposed by several authors, with the purpose to address some of the ΛCDM challenges. One of the most promising DM alternatives is the warm dark matter (WDM) model (e.g., Viel et al. 2005; Lovell et al. 2012), where particles have a rest mass on the order of a few keV, such as sterile neutrinos or thermal relics, that had non-negligible velocities at early times. Another very popular DM scenario is the self-interacting dark matter (SIDM) model where particles interact with each other (e.g., Spergel & Steinhardt 2000) having non-negligible cross sections, on the order of ∼1 cm² g⁻¹ (e.g., Zavala et al. 2013). Other more exotic DM alternatives include the following: ultra-light axion dark matter (Schwabe et al. 2016), dark atoms (for a review, see Cline 2022), and fuzzy dark matter (e.g., Kulkarni & Ostriker 2022). Although these latter models are not examined here, the toolkit we have developed can be straightforwardly adapted to any DM model for which the redshift-dependent mass function and the density profile of haloes and subhaloes can be calculated.

The properties of the DM particle affect the formation of DM structures on all scales, their stability, as well as their evolution in time. In addition, the fundamental attributes of DM particles modify the primordial power spectrum describing the initial overdensity seeds of cosmological structures. So, differences in the intrinsic DM particle properties between different models are expected to lead to measurable deviations in the resulting mass function of collapsed objects. For instance, models that include light particles, such as WDM, feature a sharp cutoff in the differential halo mass function below a critical mass scale, which depends on DM particle mass: the mildly relativistic velocities of WDM particles in the early Universe led to small-scale density fluctuations being washed out (free streaming; e.g., Melott & Schramm 1985; Viel et al. 2005). The density profile of DM haloes also turns out to be noticeably different from model to model. For example, virialized haloes made of WDM particles typically have lower central densities with respect to CDM haloes of the same mass, by virtue of their generally later formation epochs (e.g., Lovell et al. 2012).

The study of DM haloes below subgalactic scales turns out to be particularly crucial for the exploration of the nature of DM. Nevertheless, it is extremely challenging to detect such haloes directly, in order to either measure their number density in the Universe or examine their internal structure, since they might not even form galaxies due to their small size. So, the only possible way to explore and study them is through gravitational effects.

One of the most promising methods of detecting subgalactic DM haloes is strong gravitational lensing, where light that passes near a massive object (the lens) is being deflected, traveling a longer path than it would in the absence of the gravitational potential of the lens (e.g., Weinberg 1972). As a result, when a compact background source (for instance, a radio loud quasar) emits radiation with the lens being in between the source and the observer and close enough to the line of sight, then the path of the light is affected strongly, resulting in the emergence of multiple images of the background source on sky with different magnifications (e.g., Vegetti et al. 2012), provided its projected surface density exceeds a threshold. This effect is commonly known as strong lensing (see for example, Wright & Brainerd 2000). In the special case where the source displays intrinsic variability, observable time delays between the different images (pulses) may occur (Zackrisson & Riehm 2010).

Gravitational lensing can be used to detect compact objects (COs) that could not be detected otherwise, such as primordial black holes (PBHs) or dense DM haloes. Press & Gunn (1973) introduced the idea of assessing the cosmological abundance of COs through their strong gravitational lensing effect on distant background sources. They demonstrated that the cosmological mass density of COs can be constrained by deriving the fraction of lensed radio sources. Later on, Wilkinson et al. (2001) carried out a search for milli-lenses (gravitational-lensing images with milli-arcsec separations) in Very Long Baseline Interferometry (VLBI) observations of a sample of 300 compact radio sources, but no lensed systems in the mass range ∼10⁶ M_⊙ to ∼10⁸ M_⊙ were found. Their negative result allowed them to place an upper limit Ω_CO ≲ 0.01 (95% confidence) on the cosmological density of COs in this mass range, concluding that the contribution of a primordial supermassive BHs population to the dark matter content of the Universe is negligible. The currently ongoing Search for MIlli-LEnses (SMILE) project (Casadio et al. 2021) expands the search for milli-lenses in the range ∼10⁶ M_⊙ to ∼10⁹ M_⊙, to a complete sample of ∼5000 radio-loud sources using VLBI data.

Motivated by the potential of the SMILE project, in this work we develop a novel method to exploit its upcoming results with the purpose to derive constraints on the nature of DM and discriminate between currently viable DM scenarios. Our approach is based on the concept of the lensing optical depth, representing the probability for an observed source to be gravitationally lensed by a foreground mass distribution. The prescription for the implementation of this method can be found in Zackrisson & Riehm (2007). Recently, several authors have followed similar approaches to place limits on the abundance of PBHs, using Fast Radio Bursts (FRBs), (e.g., Leung et al. 2022; Zhou et al. 2022a,b), Gamma-ray Bursts (GRBs), (e.g., Kalantari et al. 2021), afterglows of GRBs (Gao et al. 2022), and compact radio sources (Zhou et al. 2022c), while others explore the diffractive lensing of gravitational waves (GWs) emanating from binary black hole mergers by small DM haloes to probe the nature of DM (see Guo & Lu 2022). Here, we pursue the possibility of subgalactic DM haloes acting as gravitational milli-lenses.

We derive the expected number of milli-lenses in the source sample of the SMILE project for various DM models by calculating the milli-lensing optical depth as a function of the source’s redshift. This in turn depends on the halo mass function, as well as on the projected surface mass density. Both of these physical quantities have noticeable differences between various scenarios, and hence the milli-lensing optical depth exhibits differences between DM models.

The layout of this paper is as follows. In Sect. 2 we describe our calculation of the milli-lensing optical depth. In Sect. 3 we discuss the analytic descriptions we use for the structure of DM haloes for various cosmological DM scenarios, and their corresponding mass functions. In Sect. 4 we present the results of our calculations, which we discuss in Sect. 5

2. Lensing probabilities

The principal result of any survey for lensing systems in the observable Universe is the number of confirmed lensed images in a complete sample of sources (e.g., Myers et al. 2003; Browne et al. 2003). To maximize the constraining power of this product, we have to connect it to theoretical models that predict the expectation value of lensing events taking into account the differences in the abundance and density profile of DM haloes between various DM models. The most straightforward way to achieve this is to compute the lensing optical depth for any given DM scenario.

The lensing optical depth depends strongly on the mass function of gravitational lenses and on the surface density profile of each halo, which essentially is related to the density profile. It also depends on the cosmology. In this paper, we fix the cosmological parameters to be H₀ = 100 h km s⁻¹ Mpc⁻¹, h = 0.7, Ω_m = 0.3, Ω_Λ = 0.7, n = 0.97, δ_c(0)=1.674, and σ₈ = 0.8. The overall results, however, are not sensitive to small variations in these parameters.

2.1. Milli-lensing optical depth

In order for our results to be applicable to the SMILE project, we are interested in lenses that produce multiple images with angular separation on the order of milli-arcseconds (milli-lenses). Thus, we focus on lenses of masses (10⁶ − 10⁹) M_⊙, which result in lensed images of angular separation that lie on the range ∼(3 − 100) mas, considering both the lens and the source to be at cosmological distances. For the calculation of the milli-lensing optical depth, we adopt the prescription of Zackrisson & Riehm (2007). We treat the lens as a massive object of mass M_l with an angular Einstein radius

$$\begin{array}{c}\hfill {\beta }_E=\sqrt{4{\displaystyle \frac{G{M}_l}{c^2}}{\displaystyle \frac{D_{\mathit{ls}}}{D_{\mathit{ol}}{D}_{\mathit{os}}}}},\end{array}$$(1)

where D_os,  D_ls, and D_ol are the angular-diameter distances from the observer to the source, from the lens to the source, and from the observer to lens, respectively, with the lens being at redshift z while the source is located at redshift z_s. D_AB can be written as

$$\begin{array}{c}\hfill D_{\mathit{AB}}(z_A,z_B)={\displaystyle \frac{c}{1+z_B}}{\displaystyle {\int }_{z_A}^{z_B}}{\displaystyle \frac{\mathrm{d}z}{H(z)}},\end{array}$$(2)

where H(z) is the Hubble parameter,

$$\begin{array}{c}\hfill {\left({\displaystyle \frac{H(z)}{H_0}}\right)}^2={\mathrm{\Omega }}_m{(1+z)}^3+{\mathrm{\Omega }}_{\mathrm{\Lambda }},\end{array}$$(3)

with Ω_m,  Ω_Λ referring to the present values of the density parameters for matter and dark energy, respectively.

The milli-lensing optical depth for a source at redshift z_s is given by

$$\begin{array}{c}\hfill \tau (z_s)={\displaystyle {\int }_0^{z_s}}\left|{\displaystyle \frac{c\mathrm{d}t}{\mathrm{d}z}}\right|\mathrm{d}z{\displaystyle {\int }_{{10}^6{M}_{\odot }}^{{10}^9{M}_{\odot }}}\sigma (M_l,z,z_s){\displaystyle \frac{\mathrm{d}n(M_l,z,z_s)}{\mathrm{d}{M}_l}}\mathrm{d}{M}_l,\end{array}$$(4)

where σ(M_l, z, z_s) is the lensing (effective) cross section

$$\begin{array}{c}\hfill \sigma (M_l,z,z_s)\equiv \pi {\beta }_E^2{D}_{\mathit{ol}}^2={\displaystyle \frac{4\pi G{M}_l}{c^2}}{\displaystyle \frac{D_{\mathit{ol}}{D}_{\mathit{ls}}}{D_{\mathit{os}}}},\end{array}$$(5)

and dn/dM_l is the differential lens mass function

$$\begin{array}{c}\hfill {\displaystyle \frac{\mathrm{d}n(M_l,z,z_s)}{\mathrm{d}{M}_l}}={\displaystyle \frac{\mathrm{d}M(M_l,z,z_s)}{\mathrm{d}{M}_l}}{\displaystyle \frac{\mathrm{d}n(M,z)}{\mathrm{d}M}}.\end{array}$$(6)

In Eq. (6), dn/dM is the differential halo mass function (see Sect. 3.2). We caution the reader of the two different masses entering Eq. (6): the lens mass M_l, and the halo mass M. These two are not in general the same because only the part of the halo in which the projected surface mass density exceeds the critical strong lensing threshold can act as a gravitational lens. The critical surface density value for a source at redshift z_s undergoing strong gravitational lensing by a foreground DM halo (lens) at redshift z is

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_{\mathit{cr}}(z,z_s)\equiv {\displaystyle \frac{M_l}{\sigma (M_l,z,z_s)}}={\displaystyle \frac{c^2}{4\pi G}}{\displaystyle \frac{D_{\mathit{os}}}{D_{\mathit{ol}}{D}_{\mathit{ls}}}}\gtrsim {10}^9{M}_{\odot }\mathrm{k}p{c}^{-2}.\end{array}$$(7)

Therefore, in order to calculate the halo mass M for given M_l,  z, and z_s, we demand a solution of the equation

$$\begin{array}{c}\hfill \mathrm{\Sigma }(M;M_l,z)={\mathrm{\Sigma }}_{\mathit{cr}}(z,z_s),\end{array}$$(8)

where Σ is the projected halo surface density described extensively in Sect. 3. Solving this equation numerically, we obtain M(M_l, z, z_s). We use the central finite difference approximation to estimate the derivative dM/dM_l, which appears in Eq. (6). Given that the halo surface mass density is obtained after an integration of the density profile, it is clear that the results will differ significantly from model to model, since each DM scenario predicts a different density profile.

2.2. Expectation value of lensing events

Once we obtain the milli-lensing optical depth, we evaluate the expectation number of lensing events in the SMILE source sample using

$$\begin{array}{c}\hfill N_l={\displaystyle \sum _{i=1}^{N_{\mathrm{sources}}}}1-exp(-\tau (z_{s,i})).\end{array}$$(9)

For τ ≪ 1, we can approximate Eq. (9) by

$$\begin{array}{c}\hfill N_l\approx {\displaystyle \sum _{i=1}^{N_{\mathrm{sources}}}}\tau (z_{s,i}).\end{array}$$(10)

The source sample, as well as their corresponding redshifts are described next.

2.3. SMILE sample

The source sample considered in this study is the one of SMILE¹: a complete sample built starting from the complete sample used in the Cosmic Lens All-Sky Survey (CLASS; Myers et al. 2003; Browne et al. 2003), the most successful search to date for gravitational lens systems at galactic scales using radio frequencies. The complete sample of 11 685 sources presented in CLASS is drawn from two other catalogs: the 5 GHz GB6 catalog (Gregory et al. 1996), and the 1.4 GHz NVSS catalog (Condon et al. 1998). The CLASS catalog contains sources from declination 0° to 75°, with a minimum flux density of 30 mJy at 5 GHz, flat spectral index (< 0.5) between 1.4 and 5 GHz, and Galactic latitude (|b|≥10°). The complete sample of 11685 sources has been initially followed up in CLASS with low resolution Very Long Array (VLA) observations at 8 GHz. The SMILE started from the complete sample in CLASS and selected sources with total flux density at 8 GHz ≥ 50 mJy. The 4968 sources that satisfy such a requirement make a complete sample of flat spectrum sources at declination [0°, +75°].

In order to obtain redshift measurements for sources in the SMILE sample, we used the Optical Characteristics of Astrometric Radio Sources (OCARS) catalog (Malkin 2018), containing redshift measurements of a large number of radio sources observed in different VLBI astrometry programs. Of the 4968 sources in SMILE, 2781 have an optical counterpart within 3 arcsec with redshift measurements in OCARS. For the remaining sources we searched for an optical counterpart within 3 arc seconds, with known redshift, in NED². In total, we collected redshifts for ∼2/3 of sources in SMILE. For the remaining ∼1/3, we followed a rather conservative approach generating redshift measurements by randomly selecting values from the known redshift sample. Their distribution is shown in Fig. 1.

Fig. 1. Redshift distribution for sources used in this study. Cyan solid line represents the redshift distribution of sources with known redshift, whereas golden dashed line stands for the distribution of the randomly selected redshift measurements from the known redshift sample.

3. DM haloes & mass functions

3.1. Halo size & structure

The internal structure of dark matter haloes affects the lensing optical depth, since the threshold for strong gravitational lensing is associated with the projected surface mass density, which in turn is related to the shape of the density profile. For a detailed review of the various mass densities see Zavala & Frenk (2019), while a thorough comparison among various density profiles can be found in Merritt et al. (2006) work.

An important finding of the past decades is that spherically averaged DM density profiles in N-body cosmological simulations have a universal form (Navarro et al. 1997). Such density profiles are described by a simple functional form characterized by only two free parameters. The first one is the concentration parameter, denoted by c_Δ, which quantifies how concentrated the mass is toward the center of the halo. The other one is the characteristic radius, r_s, which determines the distance from the center above which the density profile becomes steeper, that is to say quantifies roughly the size of the core. These two parameters are related to each other through

$$\begin{array}{c}\hfill c_{\mathrm{\Delta }}={\displaystyle \frac{r_{\mathrm{\Delta }}}{r_s}},\end{array}$$(11)

where r_Δ is the virial radius.

Since the distribution of mass is continuous, the boundary of a halo cannot be defined precisely. So, another major challenge is to come up with a robust method of determining the size of a halo uniquely. Thus far, numerous papers that deal with this problem have been published by several authors (Cole & Lacey 1996; White 2001; Cuesta et al. 2008; Zavala & Frenk 2019). In general, the radius of a halo can be defined through the overdensity parameter, Δ(z), which in principle depends on the cosmology (Bryan & Norman 1998; Tinker et al. 2008; Naderi et al. 2015; Seppi et al. 2021). In particular, it represents the radius where the mean interior density is Δ(z) times the critical density of the Universe ρ_cr(z), namely

$$\begin{array}{c}\hfill {\displaystyle \frac{3}{4\pi r_{\mathrm{\Delta }}^3}}{\displaystyle {\int }_0^{r_{\mathrm{\Delta }}}}\rho (\mathit{r})d^3\mathit{r}=\mathrm{\Delta }{\rho }_{\mathit{cr}},\end{array}$$(12)

where the critical density is given by

$$\begin{array}{c}\hfill {\rho }_{\mathit{cr}}(z)={\displaystyle \frac{3H^2(z)}{8\pi G}}={\rho }_{cr,0}{\left({\displaystyle \frac{H(z)}{H_0}}\right)}^2,\end{array}$$(13)

with G being the Newtonian gravitational constant, while ρ_cr, 0 accounts for the critical density of the Universe at redshift z = 0, and H(z) is given in Eq. (3).

The halo mass M_Δ, which is the mass contained within a sphere of radius r_Δ, is given by

$$\begin{array}{c}\hfill M_{\mathrm{\Delta }}=\mathrm{\Delta }{\displaystyle \frac{4\pi }{3}}{r}_{\mathrm{\Delta }}^3{\rho }_{\mathit{cr}},\end{array}$$(14)

and as a result the halo radius can also be written as

$$\begin{array}{c}\hfill r_{\mathrm{\Delta }}(M_{\mathrm{\Delta }},z)={\left({\displaystyle \frac{3M_{\mathrm{\Delta }}}{4\pi \mathrm{\Delta }{\rho }_{\mathit{cr}}(z)}}\right)}^{1/3}.\end{array}$$(15)

However, the most commonly used way to determine the halo’s size is to consider that the overdensity parameter Δ is fixed and equal to 200, since it turns out to be a rather convenient way to define the boundary of the halo and simplifies the calculations (e.g., Cole & Lacey 1996). Taking this fact into account, we fix the overdensity to be Δ = 200, throughout this paper, and therefore the halo mass is M₂₀₀ (hereafter, M), the halo radius is r₂₀₀, and the concentration is c₂₀₀ (hereafter, c).

Other quantities used extensively below are the enclosed mass, M_enc, the projected surface mass density, Σ, and the lens mass, M_l. The enclosed mass is defined as the mass that is contained within a sphere of radius r, namely

$$\begin{array}{c}\hfill M_{\mathrm{enc}}(r)=4\pi {\displaystyle {\int }_0^r}{r}^{\prime 2}\rho (r^{\prime })\mathrm{d}{r}^{\prime },\end{array}$$(16)

where we assume that we deal with spherically symmetric objects. The projected surface density is derived simply by the integration of the mass density along the line of sight (Wright & Brainerd 2000; Mo et al. 2010; Dhar & Williams 2010; Retana-Montenegro et al. 2012; Lapi et al. 2012)

$$\begin{array}{c}\hfill \mathrm{\Sigma }(s)={\displaystyle {\int }_{-\infty }^{\infty }}\rho (\mathit{r})\mathrm{d}z,\end{array}$$(17)

where $s=\sqrt{r^2+z^2}$ is the projected radius (orthogonal to the line of sight) relative to the center.

Having defined the projected surface density we can infer the gravitational lens mass M_l, by carrying out an integration of the Σ(s) over the disk which has radius s (see for example Eq. (41) in Retana-Montenegro et al. 2012)

$$\begin{array}{c}\hfill M_l(s)=2\pi {\displaystyle {\int }_0^s}x\mathrm{\Sigma }(x)\mathrm{d}x.\end{array}$$(18)

This quantity is the mass contained within an infinite cylinder of radius s in which the mass distribution is characterized by the mass density profile ρ(r). Equation (8) cannot be solved independently, but must simultaneously satisfy Eq. (18), owing to the fact that one of the input parameters in Eq. (8) is the lens mass. Therefore, the relation between lens mass and halo mass that is needed in the computation of the lensing optical depth is derived only after solving this coupled nonlinear system.

3.2. Mass function

The computation of the lensing optical depth depends on the number density of lensing objects, which in our case are DM haloes. So, we need to obtain a formula that determines the distribution of virialized DM haloes (in the field) per volume element and per halo mass at a given redshift z, in order to estimate the lensing optical depth. This problem has been addressed by several authors in the past, either analytically (Press & Schechter 1974; Bond et al. 1991; Pavlidou & Fields 2005) or numerically (Jenkins et al. 2001; Tinker et al. 2008). All these works are based on the spherical collapse scenario (e.g., Gunn & Gott 1972; Naderi et al. 2015). Improvements using ellipsoidal collapse do exist (e.g., Sheth et al. 2001), but we do not use them in this paper.

3.2.1. CDM mass function

The differential halo mass function of CDM haloes (the number of haloes with mass between the range M and M + dM per proper volume at a given redshift) reads (Press & Schechter 1974, see also Appendix A)

$$\begin{array}{c}\hfill n(M,z)\equiv {\displaystyle \frac{\mathrm{d}N}{\mathrm{d}M}}=\sqrt{{\displaystyle \frac{2}{\pi }}}{\displaystyle \frac{{\rho }_m(z)}{M}}{\displaystyle \frac{{\delta }_c(z)}{{\sigma }_M^2}}\left|{\displaystyle \frac{\mathrm{d}{\sigma }_M}{\mathrm{d}M}}\right|exp[-{\displaystyle \frac{{\delta }_c^2(z)}{2{\sigma }_M^2}}],\end{array}$$(19)

where M refers to the halo mass, z is the redshift, ρ_m(z) is the mean matter density of the Universe at redshift z, and δ_c(z) denotes the overdensity of a structure collapsing at redshift z linearly extrapolated to the present. Moreover, σ_M is the rms of the density field smoothed on scale M (see Appendix A).

3.2.2. WDM mass function

Although the halo mass function for WDM haloes shows a similar behavior to the CDM one on galaxy clusters scales, it exhibits a cutoff below the dwarf galaxy scale, owing to the free streaming of WDM particles in the early Universe. The most commonly applied method to derive the WDM halo mass function is the development of N-body simulations that evolve the primordial density field perturbations in time, leading to collapsed DM haloes (e.g., Schneider et al. 2012; Bose et al. 2016; Lovell 2020a,b). In this study, we choose to use the numerical fit offered by Lovell (2020b) to take into account the cutoff in the mass function of WDM haloes with respect to the CDM one. The halo mass function in the case of WDM is given by

$$\begin{array}{c}\hfill n_{\mathrm{WDM}}(M,z)=n_{\mathrm{CDM}}(M,z){(1+{\left({\displaystyle \frac{\alpha M_{\mathrm{hm}}}{M}}\right)}^{\beta })}^{\gamma },\end{array}$$(20)

where M_hm is a characteristic mass scale (the half-mode mass), while α,  β,  γ are parameters of the fit that have been found to be 2.3,  0.8,   − 1, respectively. The half-mode mass is associated with the free-streaming length which in turn is related to the rest mass of the WDM particle.

Here, we are interested in exploring the case where the WDM is made of collision-less particles (thermal relics) of mass m_WDM = 3.3 keV. In this scenario, the theoretical value for the half-mode mass is M_hm ≃ 2 × 10⁸ M_⊙ (e.g., Bose et al. 2016). This value for the half-mode mass coincides with the one for the well-motivated sterile-neutrino model in which particles are assumed to have a rest mass equal to 7 keV and lepton asymmetry number L₆ = 8.66 (see Bose et al. 2016). Sterile neutrinos are part of the neutrino Minimal Standard Model (νMSM; Boyarsky et al. 2009), which is a simple extension to the Standard Model of particle physics. It has been introduced to explain the unidentified 3.53 keV X-ray line observed recently in galaxies (see for example, Boyarsky et al. 2014) by considering this line to be the decay signal of those 7 keV sterile neutrinos.

Given that the cutoff in the halo mass function of these two WDM models is determined by the same half-mode mass and that the internal structure of haloes is identical, the inferences of this work concerning the 3.3 keV thermal relic WDM particle will also be valid for the 7 keV sterile neutrinos model. In Fig. 2, we display the differential halo mass function for the CDM model and for the WDM one investigated here as a function of the halo mass for various redshifts. Solid lines correspond to the proper density of haloes of mass M at different epochs (redshifts) divided by the present critical density of the Universe, whereas dashed lines represent the same quantities, but for WDM. A major difference between the CDM mass function and the mass function of WDM is that the latter one exhibits a cutoff at the dwarf-galaxy scales (∼10⁹ M_⊙). This distinctive feature of the halo mass function in the WDM scenario has a considerable implication to the ability of DM haloes to act efficiently as gravitational milli-lenses on background sources, because the mass function is directly involved in the calculation of the milli-lensing optical depth (Eq. (4)).

Fig. 2. Comparison of the CDM differential halo mass function with the one of the WDM model for various redshifts. Solid lines refer to CDM, while dashed lines correspond to WDM. The vertical axis corresponds to the halo mass while the y-axis shows the proper density of haloes of mass M normalized to the present value of the critical density of the Universe.

3.3. CDM halo density profile

In this study, we employ the Navarro-Frenk-White (NFW) profile for the description of the mass distribution within CDM haloes (Navarro et al. 1995, 1996, 1997)

$$\begin{array}{c}\hfill \rho (r,r_s,c,z)={\rho }_{\mathit{cr}}(z){\displaystyle \frac{{\varphi }_c}{(r/r_s){(1+r/r_s)}^2}},\end{array}$$(21)

where r_s is the characteristic radius, while ϕ_c is calculated by

$$\begin{array}{c}\hfill {\varphi }_c={\displaystyle \frac{200}{3}}{\displaystyle \frac{c^3}{ln(1+c)-c/(1+c)}},\end{array}$$(22)

with c being the concentration parameter, which is not mass independent, but correlates strongly with the halo mass, as well as with the redshift, following a simple scaling law (see for example, Bullock et al. 2001; Neto et al. 2007; Prada et al. 2012; Dutton & Macciò 2014; Klypin et al. 2016; Shan et al. 2017; Ragagnin et al. 2019, 2021). The NFW profile predicts a cuspy halo’s center, since the mass density goes as ∼r⁻¹ near the center of the halo.

Using Eq. (21) along with Eq. (16), we obtain the enclosed mass

$$\begin{array}{c}\hfill M_{\mathrm{enc}}(M,z,r_{200}(M,z))={\displaystyle \frac{M}{ln(1+c)-c/(1+c)}}F(r,c,r_{200}),\end{array}$$(23)

where

$$\begin{array}{c}\hfill F(r,c,r_{200})=ln(1+\mathit{cr}/r_{200})-{\displaystyle \frac{\mathit{cr}/r_{200}}{1+\mathit{cr}/r_{200}}}.\end{array}$$(24)

From Eqs. (17) and (21), we obtain the surface density

$$\begin{array}{c}\hfill \mathrm{\Sigma }(s,M,z,r_{200}(M,z))=2{\displaystyle \frac{{\varphi }_c{\rho }_{\mathit{cr}}(z)r_{200}{c}^{-1}}{c^2{(s/r_{200})}^2-1}}\mathcal{S}(s,c,r_{200}(M,z)),\end{array}$$(25)

where we have defined for convenience

$$\begin{array}{c}\hfill \mathcal{S}(s,c,r_{200}(M,z))=\{\begin{array}{cc}1-{\displaystyle \frac{{cos}^{-1}(r_{200}/cs)}{\sqrt{c^2{(s/r_{200})}^2-1}}}\hfill & \mathrm{if}s>r_s\hfill \\ 1-{\displaystyle \frac{{cosh}^{-1}(r_{200}/cs)}{\sqrt{1-c^2{(s/r_{200})}^2}}}\hfill & \mathrm{if}s<r_s\hfill \end{array}.\end{array}$$(26)

From Eq. (18), the lens mass then is

$$\begin{array}{c}\hfill M_{\mathrm{lens}}(s,M,z,r_{200}(M,z))={\displaystyle \frac{M}{ln(1+c)-c/(1+c)}}\mathcal{H}(s,c,r_{200}),\end{array}$$(27)

where

$$\begin{array}{cc}& \mathcal{H}(s,c,r_{200})=ln(cs/2r_{200})\hfill \\ \hfill & +\{\begin{array}{cc}{\displaystyle \frac{2}{\sqrt{{(cs/r_{200})}^2-1}}}arctan\sqrt{{\displaystyle \frac{cs/r_{200}-1}{cs/r_{200}+1}}}\hfill & \mathrm{if}cs>r_{200}\hfill \\ 1\hfill & \text{if}cs=r_{200}\hfill \\ {\displaystyle \frac{2}{\sqrt{1-{(cs/r_{200})}^2}}}\mathrm{arctanh}\sqrt{{\displaystyle \frac{1-cs/r_{200}}{cs/r_{200}+1}}}\hfill & \mathrm{if}cs<r_{200}\hfill \end{array}.\hfill \end{array}$$(28)

3.4. WDM halo density profile

Warm dark matter is made of particles that had non-negligible thermal velocities at early times. This major difference is expected to have an impact on the concentration of mass near the center, but not in the shape of the distribution of mass within a halo. Indeed, the density profile in WDM models can be well described by a NFW profile (see e.g., Lovell et al. 2014; Bose et al. 2016). However, DM haloes consisting of WDM are typically formed at smaller redshifts with respect to the formation of CDM haloes. This difference affects the concentration of the halo, which generally reflects the density of the Universe at the epoch of halo formation (for a detailed discussion see, Schneider et al. 2012; Bose et al. 2016; Zavala & Frenk 2019). Therefore, we assume that the mass distribution in WDM haloes is consistent with the NFW profile, but more fuzzy, that is to say less concentrated around the center. In order to calculate the concentration parameter for WDM haloes as a function of the mass and the redshift we use the findings of Bose et al. (2016). They offer the following simple functional form for the concentration parameter

$$\begin{array}{c}\hfill {\displaystyle \frac{c_{\mathrm{WDM}}}{c_{\mathrm{CDM}}}}={(1+{\gamma }_1{\displaystyle \frac{M_{\mathrm{hm}}}{M}})}^{-{\gamma }_2}{(1+z)}^{\beta (z)},\end{array}$$(29)

where γ₁ = 60, γ₂ = 0.17, and β(z)=0.026z − 0.04. M_hm is the half-mode mass and in this work is set to be M_hm = 2 × 10⁸ M_⊙ corresponding either to the model of thermal relics WDM particles of rest mass m_WDM = 3.3 keV or to the 7 keV sterile neutrinos model (an extension to the Standard model). Due to their smaller concentrations, WDM haloes will be less likely to exceed the lensing surface-density threshold, resulting in a lower milli-lensing optical depth.

3.5. SIDM halo density profile

3.5.1. SIDM core-like halo center

The SIDM model was originally introduced by Spergel & Steinhardt (2000) to explain observations of central densities in galaxies within the Local Group. Since then, numerous authors have argued that the self-interaction of particles leads to a core-like profile rather than a cusp-like profile. Such a feature could alleviate the “cusp-core” problem arising for CDM, and hence is considered a well motivated DM alternative.

There are two kinds of SIDM theories. In the first case the scattering rate per particle, Γ(r), is velocity independent, which implies that the ratio of the effective cross section, σ, to the dark matter particle’s mass, m, is constant (e.g., Rocha et al. 2013; Elbert et al. 2015). In the second scenario, the scattering rate is velocity dependent and falls rapidly as the velocity increases (e.g., Zavala et al. 2013). For a recent discussion on SIDM models, as well as on the observational constraints on the self-scattering cross section, see Tulin & Yu (2018).

In this work, we assume that the scattering rate per particle is velocity independent and has the following form

$$\begin{array}{c}\hfill \mathrm{\Gamma }(r)\propto \rho (r)(\sigma /m)v_{\mathrm{rms}}(r),\end{array}$$(30)

where ρ(r) is the DM mass density at radius r, while v_rms is the rms speed of dark matter particles. We consider a typical value for the ratio σ/m ∼ 1 cm² g⁻¹, since SIDM models with smaller values, on the order of 0.1 cm² g⁻¹, are very similar to the CDM models even on scales smaller than dwarf galaxies and cannot produce detectable deviations from CDM predictions (Zavala et al. 2013). On the other hand, higher values of the cross section per mass, ∼10 cm² g⁻¹, have already been ruled out by cluster observations (see e.g., Dawson et al. 2012).

For the structure of SIDM haloes, we again assume spherical symmetry, but now we use a core-like profile. In fact, the mass density is well approximated by the Burkert profile (see Burkert 1995), which also has two free parameters and is given by the following formula

$$\begin{array}{c}\hfill {\rho }_B(r,r_b,{\rho }_b)={\displaystyle \frac{{\rho }_b}{(1+r/r_b)(1+{(r/r_b)}^2)}},\end{array}$$(31)

where r_b is the scale (core) radius, while ρ_b is the central density. As in the NFW profile, the free parameters of the Burkert profile scale with the halo mass. In the special case where σ/m ∼ 1 cm² g⁻¹, Rocha et al. (2013) have provided a couple of simple scaling laws that connect both the r_b and ρ_b with the halo mass, using data from N-body simulations. These relations (Eqs. (17) and (20) in Rocha et al. 2013) are given below

$$\begin{array}{c}\hfill {\displaystyle \frac{r_b}{1\mathrm{kpc}}}=2.21{\left({\displaystyle \frac{M_{\mathrm{vir}}}{{10}^{10}{M}_{\odot }}}\right)}^{0.43},\end{array}$$(32)

$$\begin{array}{c}\hfill {\displaystyle \frac{{\rho }_b}{M_{\odot }/p{c}^3}}=0.029{\left({\displaystyle \frac{M_{\mathrm{vir}}}{{10}^{10}{M}_{\odot }}}\right)}^{-0.19}.\end{array}$$(33)

Since these relations have been derived using the virial mass, M_vir, instead of M₂₀₀ which we have employed throughout this work, we have to rescale the density profile to be consistent with Eq. (12). Equation (31) can be recast as

$$\begin{array}{c}\hfill {\rho }_B(r,M)=\mathcal{A}(M){\displaystyle \frac{{\rho }_b}{(1+r/r_b)(1+{(r/r_b)}^2)}},\end{array}$$(34)

where

$$\begin{array}{c}\hfill \mathcal{A}(M)={\displaystyle \frac{M{r}_b^{-3}}{2\pi {\rho }_b}}{[ln(1+{\displaystyle \frac{r_{200}}{r_b}})+{\displaystyle \frac{1}{2}}ln(1+{\displaystyle \frac{r_{200}^2}{r_b^2}})-{tan}^{-1}\left({\displaystyle \frac{r_{200}}{r_b}}\right)]}^{-1},\end{array}$$(35)

and now we can use the halo mass M ≡ M₂₀₀ instead of the virial mass in Eqs. (32) and (33). The term $\mathcal{A}(M)$ has been derived by requiring the mean density inside a sphere of radius r₂₀₀ to be 200ρ_cr (see Eq. (12)).

Using Eqs. (16) and (34), we obtain for the enclosed mass

$$\begin{array}{c}\hfill M_{\mathrm{enc}}(r)=\mathcal{A}(M)\pi {\rho }_b{r}_b^3[ln(1+{\displaystyle \frac{r^2}{r_b^2}})+2ln(1+{\displaystyle \frac{r}{r_b}})-2{tan}^{-1}\left({\displaystyle \frac{r}{r_b}}\right)].\end{array}$$(36)

The surface density cannot in general be derived analytically and therefore we have to perform the integration numerically. In the special case where s = 0, that is for the column density through the line of sight, the integration that returns the surface density Σ_B(0) yields the closed-form expression

$$\begin{array}{c}\hfill {\mathrm{\Sigma }}_B(s=0,M)=\mathcal{A}(M){\displaystyle \frac{\pi }{2}}{\rho }_b{r}_b,\end{array}$$(37)

where the index B indicates that this surface density arises from the Burkert profile. The surface mass density is maximized when s = 0, since ρ(r) is a monotonically decreasing function of r, and as a result for a SIDM halo of given mass, the maximum value of the surface density is determined by Eq. (37). This feature is of great importance in strong gravitational lensing where the surface density must exceed a critical threshold in order to significantly bend a light ray.

In order to have a qualitative picture of the differences between the three mass density profiles mentioned above, in Fig. 3, we apply them to a DM halo of mass M = 10⁸ M_⊙ at redshift z = 0. For SIDM particles the density profile near the center is flat, while for CDM particles the profile near the center goes as r⁻¹, since we have used the NFW profile. For WDM particles the profile is NFW-like but with larger characteristic radius than in the CDM reflecting the fact that the concentration in the WDM scenario is smaller than the one in CDM.

Fig. 3. Comparison of mass density profiles for a given halo of mass M = 108 M⊙ at redshift z = 0 for various dark matter scenarios. The concentration for the CDM case is given by Eq. (38).

3.5.2. SIDM core collapse

Even though most of SIDM models are in favor of a less dense core-like halo center, the strong self-interaction developed between particles in the innermost region of the halo might have important implications in the dynamical evolution of the halo. In one scenario, strong self-interactions between particles induce a negative heat capacity, eventually leading to the formation of a dense central core in the inner part of the halo (see e.g., Yang & Yu 2021, 2022). Yang & Yu (2021) demonstrated that such a scenario can be successful in explaining the observational excess of small-scale gravitational lenses in galaxy clusters reported in Meneghetti et al. (2020). They exploited the fact that at late stages of the gravothermal evolution of a halo composed of SIDM, the core might undergo gravothermal collapse, resulting in a highly dense halo center, thereby increasing its lensing effect on background sources compared to CDM haloes.

In the most extreme case, the collapsed core can further contract, eventually leading to the formation of a supermassive black hole (SMBH) at the halo center. This scenario was firstly proposed and studied extensively by Feng et al. (2021; see also Feng et al. 2022) as a possible mechanism to explain the existence and origin of SMBHs at high redshifts (z ∼ 6 − 7). Essentially, SIDM offers a natural mechanism for triggering dynamical instability, a necessary condition to form a black hole. This scenario can be tested and well-constrained through milli-lensing, since the central SMBH can effectively act as a strong gravitational lens and produce multiple images of a compact background source.

Given that studies dealing with the core collapse scenario do not provide an exact formula for the final mass distribution of DM inside the collapsed halo, we shall restrict ourselves in investigating here only the latter, most extreme, scenario of core collapse where the formation of a SMBH takes place from the gravothermal collapse of the core. The exploration of this model yields an upper limit on the expectation value of lensing events in the SMILE source sample in the case of the SIDM scenario.

4. Results

4.1. CDM

4.1.1. CDM: model A

We start by investigating the CDM scenario using a concentration-mass relation derived from N-body simulations. We employ the relation given in Ragagnin et al. (2019) to determine the dependence of the concentration parameter c on the redshift z, as well as on the halo mass M:

$$\begin{array}{c}\hfill c(M,z)=6.02{\left({\displaystyle \frac{M}{{10}^{13}{M}_{\odot }}}\right)}^{-0.12}{\left({\displaystyle \frac{1.47}{1+z}}\right)}^{0.16}.\end{array}$$(38)

In Fig. 4 we plot with a blue solid line the milli-lensing optical depth obtained for this c(M, z). The value of the milli-lensing optical depth is well below ∼10⁻⁴ implying that even a sample of ten thousands distant (z ∼ 5) compact sources is highly unlikely to produce at least one lensing event. Indeed, performing the summation in Eq. (10) over all sources involved in the SMILE sample, we end up with the value ⟨N_exp⟩≃1.5 × 10⁻³, which makes detection of a milli-lens improbable.

Fig. 4. Lensing optical depth as a function of the source redshift for different dark matter scenarios.

4.1.2. CDM: model B

As a limiting case of the possible effect of the concentration-mass relation on our results, we also test a power-law extrapolation to lower masses of the empirical (fitted from observations rather than simulations) c − M relation shown in Fig. 13 of Prada et al. (2012)

$$\begin{array}{c}\hfill logc(M,z)=4.23-0.25log(M/M_{\odot })-0.16log\left({\displaystyle \frac{1+z}{1.47}}\right).\end{array}$$(39)

Regarding the dependence on redshift, we consider that it is identical to Eq. (38), but we stress that most studies suggest a weak dependence of the concentration on redshift, so even if we slightly modify the last term in Eq. (39) associated with the redshift dependence, the overall results do not change noticeably. In practice, the concentration parameter is set by the halo mass. We note that Eq. (39) predicts higher values of the concentration parameter with respect to the ones inferred from N-body simulations. Although this c − M relation has been derived from galaxy cluster observations and might overestimate the c parameter of haloes on subgalactic scales, recently Şengül & Dvorkin (2022) investigated the strong lens system JVAS B1938+666, concluding that subgalactic DM haloes can be highly concentrated (c ≈ 60), in line with Eq. (39).

Using this relation in Eq. (4), we obtain the green dash-dotted line in Fig. 4, showing the milli-lensing optical depth as a function of the source redshift. Subsequently, using Eq. (10) to compute the expectation value of lensing events in the source sample of SMILE, we obtain ⟨N_exp⟩≃1.2. This value deviates remarkably from the one corresponding to model A (see Sect. 4.1.1), demonstrating that the concentration-mass relation plays a crucial role in the process of strong gravitational milli-lensing and can thus be strongly constrained with milli-lensing observations. This value also places an upper limit on the expectation number of detected milli-lenses in the SMILE’s source sample, in the case where the properties of DM particles are in line with the framework of the CDM model.

A comparison between the two concentration-mass relations related to the CDM scenario for redshift z = 0 can be found in Fig. 5. The concentration-mass relation given by Eq. (38), (Model A) is displayed with a blue solid line, while the green dash-dotted line stands for the c(M) considered in model B (i.e., Eq. (39)).

Fig. 5. Concentration-mass relation at z = 0 for different dark matter scenarios.

4.2. SIDM

In the SIDM model, there are two possibilities, which lead to quite different internal structure of haloes. The first corresponds to haloes described by a core-like density profile, while the latter refers to haloes of collapsed cores.

4.2.1. SIDM core-like halo center

The standard scenario is the one where the inner part of haloes is characterized by a core-like profile yielding the projected surface mass density of Eq. (37). Since the surface mass density is maximized at the center of the halo (as long as the density profile is a decreasing function of r), if the central region does not exceed the critical threshold for strong lensing, then the halo will not act as a strong lens. Using Eqs. (37) and (7), we conclude that the halo mass of a SIDM halo must be ≳10¹⁴ M_⊙, for the surface mass density at the center to exceed the strong lensing threshold. However, this mass scale corresponds to galaxy clusters and therefore no trustful inferences can be done without including the effect of strong lensing due to the presence of baryons. The main finding is that SIDM-only subgalactic haloes cannot produce milli-lensing images since they are not dense enough to satisfy the strong lensing criterion.

4.2.2. SIDM core collapse: model C

The second scenario related to SIDM haloes is based on the gravothermal core collapse process that might take place in the inner parts of SIDM haloes (see Sect. 3.5.2). Assuming that the halo in the beginning was described by a NFW profile with the concentration-mass relation to given by Eq. (38), we can calculate the mass enclosed inside a projected disk with radius equal to the scale radius r_s. Then, we consider the most extreme case where the entire core collapses into a very small but extremely dense core that eventually results in the formation of a compact object. This collapsed core is the part of the halo that can produce strong gravitational lensing of light emitted by background sources.

In Fig. 4, we show with the black dotted line the milli-lensing optical depth in the case of SIDM core collapse. Having obtained the optical depth, we carry out the sum shown in Eq. (10) over the redshifts of the SMILE project sources and find the value ⟨N_exp⟩≃13.

4.3. WDM: model D

In order to derive the milli-lensing optical depth for the scenario of WDM, where particles are supposed to have a rest mass m_WDM = 3.3 keV (thermal relic) or be sterile neutrinos with a rest mass equal to 7 keV, we take into account that the halo mass function is different from that in CDM, and so we adopt the fit offered by Lovell (2020b; Eq. (20)). Regarding the density profile we again employ the NFW one, but with the concentration parameter to be given by Eq. (29). This concentration-mass relation is however a fit that relates the concentration of WDM haloes with the one corresponding to CDM haloes, so we simply consider that the concentration-mass relation of CDM haloes is the one shown in Eq. (38) and in such a way we obtain a formula for the concentration of WDM haloes as a function of the halo mass and redshift. In Fig. 5, we display the concentration-mass relation which corresponds to WDM haloes at z = 0 with a red dashed line.

Performing the integration of Eq. (4), we find the milli-lensing optical depth in the case of WDM, shown in Fig. 4 with the red dashed line. Combining this result with Eq. (10), we compute the expectation number of detected WDM milli-lenses, obtaining ⟨N_exp⟩≃1.1 × 10⁻³. It is therefore extremely unlikely to detect any milli-lenses with SMILE if DM is in the form of WDM.

5. Discussion and conclusions

In this work we have explored the ability of subgalactic DM haloes to act as milli-lenses on background sources resulting in multiple images of the same source with angular separation on the order of milli-arcseconds, considering different DM models. We have developed a semi-analytical method to estimate the expectation value of detected milli-lenses in several DM scenarios, computing the lensing optical depth. We have modeled the number density and internal structure of haloes using either (semi)analytical calculations or fits to N-body simulation results, depending on the DM model. We have restricted ourselves in applying the point-mass lens approximation to infer the lens mass, imposing the effective surface threshold criterion for strong lensing to connect the lens mass to the halo mass. Finally, we used the milli-lensing optical depth in each scenario to calculate the expectation number of detected milli-lenses in the source sample of the SMILE project.

We found that the probability of strong milli-lensing by DM haloes strongly depends on the model, being regulated by the properties of DM particles which dictate the inner structure of haloes, as well as their number density. We have shown that even within the CDM model, the lensing optical depth is quite sensitive to variations in the concentration-mass relation, in agreement with Amorisco et al. (2022), leading to very different expectation values of detected milli-lenses in the SMILE source sample. Milli-lensing observations might therefore enable us to constrain the concentration-mass relation down to subgalactic mass scales.

In addition, we have demonstrated that DM scenarios which are in favor of core-like density profiles, such as the SIDM one investigated here, are unlikely to produce milli-lenses because they predict haloes with low-density centers. However, our method allows to also probe scenarios like core collapse which enhance considerably the probability of milli-lensing.

Finally, we have shown that haloes consisting of WDM lead to an extremely small milli-lensing optical depth due to their combination of low concentration and mass-function cutoff. Even if a steeper concentration-mass relation (such as Eq. (39)) is used, the cutoff in the number density below subgalactic scales still prevent the milli-lensing optical depth from increasing significantly. Therefore, the detection of milli-lenses would provide definitive evidence against the WDM model and more generally models that exhibit a cutoff in their halo mass function affecting the 10⁶ − 10⁹ M_⊙ mass scales.

In Fig. 6, we summarize our results, plotting the expectation value of detected milli-lenses for all models investigated in this study. The blue point corresponds to the CDM model A (Sect. 4.1.1), while the green point refers to the CDM model B (Sect. 4.1.2). The black point corresponds to the core collapse SIDM scenario C (Sect. 4.2.2) and the red point to the WDM model D (Sect. 4.3). Even among the limited number of DM models studied here, milli-lensing observations of source samples comparable to that of SMILE hold significant discriminating power.

Fig. 6. Expectation number of detected milli-lenses in the SMILE project. The error-bars have been derived assuming Poisson-like error in the calculation of the expectation values.

Even though subgalactic haloes are expected to be almost empty of baryons, one source of uncertainty in our work might be the fact that we ignore the overall effect of baryons in the internal structure of haloes, which in principle can alter the ability of those haloes to act as milli-lenses. Although such systems are DM dominated, the baryons might play a crucial role in the strong lensing and therefore it is left for a future investigation, since the purpose of this paper is mostly to highlight the point that milli-lensing observations can be used to constrain the nature of dark matter and further discriminate between currently viable models.

Another possible source of uncertainty arises from the fact that some DM haloes might host a SMBH at their center, which would contribute significantly in the strong lensing signal from those haloes, thereby modifying our results. However, neither the fraction of haloes that host such objects at their centers nor the accurate relation between SMBH and halo mass are known. In addition, there is still much discussion on explaining the existence of unexpectedly large SMBHs in luminous quasars observed at high redshifts, which might challenge current theoretical models. In general milli-lensing experiments, and in particular the upcoming results of the SMILE project (see Sect. 2), could be exploited to test theoretical models of SMBH formation in the early Universe through their imprints on milli-lensing signals, and thus milli-lensing surveys might have strong discriminatory power in this context. Nevertheless, this issue requires very careful and detailed treatment in modeling the connection of DM haloes to SMBHs that are embedded inside them, and therefore will be addressed in a future paper. Hence we do not take this possibility into account and ignore the effect of SMBHs in this study.

Here, we assumed a certain redshift distribution for the ∼1/3 of sources for which no redshift measurements were currently available. Given the lack of knowledge about the distances of these sources, we adapted a rather plausible redshift distribution, which was chosen to be similar to the one of the sources with known redshift. This conservative choice partially affects the results, although not our qualitative conclusions. For instance, if we considered that a considerable fraction of sources with missing redshift measurements are among the weakest and most distant sources, then the expectation number of detected milli-lenses would increase. In this case, the constraints computed in this study would be lower limits.

It should be noted, however, that the qualitative findings in this work are insensitive in variations of the redshift distribution corresponding to the fraction of sources with unknown distances. As a matter of fact, even if we considered that all sources of unknown redshift (∼1600) had a redshift at about 5, then as can be inferred from Fig. 4, the lensing optical depth for the models (A) & (D), which predict zero milli-lenses, would still be well-below ∼10⁻⁴, therefore could not yield measurable deviations in the expectation values, which would again vanish for these two models. Nevertheless, better estimates of the expectation values of detected milli-lenses can be only obtained by increasing the number of sources with measured redshifts in the future.

In this study, we have concluded that: 1) the source sample included in the SMILE project is sufficiently large to enable inferences about the nature of DM; 2) WDM haloes are highly unlikely to produce even a single strong milli-lensing event in the source sample of SMILE; 3) SIDM haloes can only act as strong milli-lenses in the case where self-interactions trigger the core collapse mechanism, leading to highly dense cores; and 4) the ability of CDM subgalactic haloes to act as milli-lenses strongly depends on the mass-concentration relation. Finally, we have shown that if CDM is indeed the relevant model for describing the properties of DM particles, then milli-lensing observations will enable us to further constrain the relationship between concentration and halo mass down to subgalactic mass scales.

Acknowledgments

We thank the anonymous referee for several insightful comments and suggestions that helped us improve this manuscript. N.L. would like to thank Hai-Bo Yu for fruitful discussions and comments related to the core collapse scenario in the context of SIDM. N.L. and K.T. acknowledge support by the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme under grant agreement No. 771282. V.P. acknowledges support by the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the “First Call for H.F.R.I. Research Projects to support Faculty members and Researchers and the procurement of high-cost research equipment grant” (Project 1552 CIRCE), and by the Foundation of Research and Technology – Hellas Synergy Grants Program through project MagMASim, jointly implemented by the Institute of Astrophysics and the Institute of Applied and Computational Mathematics. CC acknowledges support by the European Research Council (ERC) under the HORIZON ERC Grants 2021 programme under grant agreement No. 101040021. K.T. acknowledges support from the Foundation of Research and Technology – Hellas Synergy Grants Program through project POLAR, jointly implemented by the Institute of Astrophysics and the Institute of Computer Science.

References

Appendix A: Press-Schechter mass function

In this appendix, we offer a short prescription for calculating the halo mass function. For more details on the subject, see Mo et al. (2010). Using the Press-Schechter formalism (Press & Schechter 1974), the differential halo mass function is expressed as

$$\begin{array}{c}\hfill n(M,z)\equiv {\displaystyle \frac{\mathit{dN}}{\mathit{dM}}}=\sqrt{{\displaystyle \frac{2}{\pi }}}{\displaystyle \frac{{\rho }_m(z)}{M}}{\displaystyle \frac{{\delta }_c(z)}{{\sigma }_M^2}}\left|{\displaystyle \frac{d{\sigma }_M}{\mathit{dM}}}\right|exp[-{\displaystyle \frac{{\delta }_c^2(z)}{2{\sigma }_M^2}}],\end{array}$$(A.1)

where ρ_m(z)=Ω_m(1 + z)³, δ_c(z) is the overdensity of a structure collapsing at redshift z linearly extrapolated to the present epoch, and σ_M accounts for the linear rms fluctuation (variance) of the density field on scale M. Under the spherical collapse assumption and concordance cosmology, the critical overdensity is given by (e.g., Eq. C30 in Pavlidou & Fields 2005)

$$\begin{array}{c}\hfill {\delta }_c(z)={\displaystyle \frac{1}{D(z)}}{\delta }_c(0),\end{array}$$(A.2)

where D(z) is the normalized linear growth factor, that is D(z = 0)=1, calculated by

$$\begin{array}{c}\hfill D(z)=G[{(2\omega )}^{1/3}/(1+z)]/G\left[{(2\omega )}^{1/3}\right],\end{array}$$(A.3)

while ω = Ω_Λ/Ω_m and

$$\begin{array}{c}\hfill G(u)={\displaystyle \frac{{(2+u^3)}^{1/2}}{u^{3/2}}}{\displaystyle {\int }_0^u}{\displaystyle \frac{y^{3/2}}{{(2+y^3)}^{3/2}}}dy.\end{array}$$(A.4)

The variance σ_M, normalized to be equal to σ₈ when R = 8h⁻¹ Mpc, is given by (see Eq. 10 in Pavlidou & Fields 2005)

$$\begin{array}{c}\hfill {\sigma }_M^2={\sigma }_8^2{\displaystyle \frac{{\displaystyle {\int }_0^{\infty }}dkP(k)W^2(kR(M))k^2}{{\displaystyle {\int }_0^{\infty }}dkP(k)W^2(k8h^{-1}\mathrm{Mpc})k^2}},\end{array}$$(A.5)

where the σ₈ parameter is a direct observable quantity, while W(kR) refers to the window function (filter) and R is the characteristic radius of the filter related to the mass M through

$$\begin{array}{c}\hfill R(M)={\left({\displaystyle \frac{M}{{\gamma }_f{\rho }_m}}\right)}^{1/3},\end{array}$$(A.6)

where γ_f is a parameter depending on the shape of the filter. Here, we employ a sharp in k-space filter given by

$$\begin{array}{c}\hfill W(\overrightarrow{k};R)=\{\begin{array}{cc}1\hfill & \text{if}|\overrightarrow{k}|\le R^{-1}\hfill \\ 0\hfill & \text{if}|\overrightarrow{k}|>R^{-1}\hfill \end{array}.\end{array}$$(A.7)

Given this choice, γ_f becomes equal to 6π² (Mo et al. 2010).

In Eq. (A.5), P(k) stands for the linear matter power spectrum. Invoking linear theory, we can write the matter power spectrum to be

$$\begin{array}{c}\hfill P(k)\propto T^2(k)P_{\mathit{init}}(k),\end{array}$$(A.8)

where P_init is the initial power spectrum proportional to kⁿ with n = 0.97, while T(k) corresponds to the transfer function. Bond & Efstathiou (1984) offer the following simple numerical formula for the transfer function T(k), which is consistent with the ΛCDM model (see also Jenkins et al. 2001)

$$\begin{array}{c}\hfill T(k)={\displaystyle \frac{1}{{[1+{[aq+{(bq)}^{3/2}+{(cq)}^2]}^{\nu }]}^{1/\nu }}},\end{array}$$(A.9)

where q = k/Γ, Γ = Ω_m, 0h, ν = 1.13, a = 6.4 h⁻¹ Mpc, b = 3 h⁻¹ Mpc, and c = 1.7 h⁻¹ Mpc.

All Figures

	Fig. 1. Redshift distribution for sources used in this study. Cyan solid line represents the redshift distribution of sources with known redshift, whereas golden dashed line stands for the distribution of the randomly selected redshift measurements from the known redshift sample.
In the text

	Fig. 2. Comparison of the CDM differential halo mass function with the one of the WDM model for various redshifts. Solid lines refer to CDM, while dashed lines correspond to WDM. The vertical axis corresponds to the halo mass while the y-axis shows the proper density of haloes of mass M normalized to the present value of the critical density of the Universe.
In the text

	Fig. 3. Comparison of mass density profiles for a given halo of mass M = 108 M⊙ at redshift z = 0 for various dark matter scenarios. The concentration for the CDM case is given by Eq. (38).
In the text

	Fig. 4. Lensing optical depth as a function of the source redshift for different dark matter scenarios.
In the text

	Fig. 5. Concentration-mass relation at z = 0 for different dark matter scenarios.
In the text

	Fig. 6. Expectation number of detected milli-lenses in the SMILE project. The error-bars have been derived assuming Poisson-like error in the calculation of the expectation values.
In the text