-
PDF
- Split View
-
Views
-
Cite
Cite
D Bayer, S Chatterjee, L V E Koopmans, S Vegetti, J P McKean, T Treu, C D Fassnacht, K Glazebrook, Probing sub-galactic mass structure with the power spectrum of surface-brightness anomalies in high-resolution observations of galaxy–galaxy strong gravitational lenses. II. Observational constraints on the subgalactic matter power spectrum, Monthly Notices of the Royal Astronomical Society, Volume 523, Issue 1, July 2023, Pages 1310–1325, https://doi.org/10.1093/mnras/stad1402
- Share Icon Share
ABSTRACT
Stringent observational constraints on the subgalactic matter power spectrum would allow one to distinguish between the concordance ΛCDM and the various alternative dark-matter models that predict significantly different properties of mass structure in galactic haloes. Galaxy–galaxy strong gravitational lensing provides a unique opportunity to probe the subgalactic mass structure in lens galaxies beyond the Local Group. Here, we demonstrate the first application of a novel methodology to observationally constrain the subgalactic matter power spectrum in the inner regions of massive elliptical lens galaxies on 1–10 kpc scales from the power spectrum of surface-brightness anomalies in highly magnified galaxy-scale Einstein rings and gravitational arcs. The pilot application of our approach to Hubble Space Telescope (HST/WFC3/F390W) observations of the SLACS lens system SDSS J0252+0039 allows us to place the following observational constraints (at the 99 per cent confidence level) on the dimensionless convergence power spectrum |$\Delta ^{2}_{\delta \kappa }$| and the standard deviation in the aperture mass σAM: |$\Delta ^{2}_{\delta \kappa }\lt 1$| (σAM < 0.8 × 108 M⊙) on 0.5-kpc scale, |$\Delta ^{2}_{\delta \kappa }\lt 0.1$| (σAM < 1 × 108 M⊙) on 1-kpc scale and |$\Delta ^{2}_{\delta \kappa }\lt 0.01$| (σAM < 3 × 108 M⊙) on 3-kpc scale. These first upper-limit constraints still considerably exceed the estimated effect of CDM subhaloes. However, future analysis of a larger sample of galaxy–galaxy strong lens systems can substantially narrow down these limits and possibly rule out dark-matter models that predict a significantly higher level of density fluctuations on the critical subgalactic scales.
1 INTRODUCTION
The dark-energy-plus-cold-dark-matter (ΛCDM) concordance cosmological model is known to successfully reproduce the observed large-scale (larger than ∼1 Mpc) distribution of matter in the Universe (e.g. Vogelsberger et al. 2014; Schaye et al. 2015; Guo et al. 2016; Planck Collaboration VI 2020). However, on smaller galactic and subgalactic scales, theory and observations appear to diverge (see Bullock & Boylan-Kolchin 2017, for a recent review on the small-scale challenges to the ΛCDM paradigm). One of the main discrepancies, known as the Missing Satellites Problem (MSP), lies in the fact that the number of dwarf satellite galaxies observed in the Local Group (∼100; e.g. McConnachie 2012; Drlica-Wagner et al. 2015) is much lower than the numerous abundance of subhaloes populating galactic haloes in ΛCDM-based numerical simulations of the cosmological structure-formation process (e.g. Klypin et al. 1999; Moore et al. 1999; Diemand, Kuhlen & Madau 2007; Nierenberg et al. 2016; Bullock & Boylan-Kolchin 2017; Dooley et al. 2017; Tanaka et al. 2018; Fielder et al. 2019; Navarro 2019).
The currently favoured interpretation of the MSP is that the missing subhaloes do exist but are extremely inefficient at forming stars due to a variety of baryonic processes, such as feedback from massive stars and active galactic nuclei, tidal stripping or photoionization squelching (e.g. Thoul & Weinberg 1996; Bullock, Kravtsov & Weinberg 2000; Somerville 2002; Sawala et al. 2014; Despali et al. 2018; Kim, Peter & Hargis 2018), and thus remain undetectable for conventional imaging surveys. Alternatively, the MSP might point towards dark-matter models with higher thermal velocities in the early Universe (referred to as warm dark matter, Bode, Ostriker & Turok 2001), which predict a suppression of structure formation on the subgalactic scales (see e.g. Menci, Fiore & Lamastra 2012; Nierenberg et al. 2013; Viel et al. 2013; Lovell et al. 2014; Vegetti et al. 2018; Zavala & Frenk 2019; Hsueh et al. 2020, for some recent studies). Lastly, the Local Group might just be a biased environment with less abundant substructures and, thus, not representative of the entire Universe (e.g. Muller et al. 2018). In order to fully clarify these ambiguities and test the aforementioned solutions to the small-scale tensions of the ΛCDM paradigm, it is crucial to constrain the properties of substructures in galaxies beyond the Local Universe on mass scales below that of luminous satellites.
The key techniques to search for the faint, or even truly dark substructures in galaxies at cosmological distances are based on the study of their gravitational imprints on the lensed images of galaxy-scale strong gravitational lens systems. Substructures, if present in the foreground lens galaxy (or along its line of sight), introduce perturbations to the otherwise smooth lensing potential, thus leading to a deviation between the observed surface-brightness distribution of the lensed images and the prediction from the best-fitting smooth-lens model, e.g. in the form of flux-ratio anomalies in multiply imaged gravitationally lensed quasars (e.g. Mao & Schneider 1998; Metcalf & Madau 2001; Dalal & Kochanek 2002; Nierenberg et al. 2014; Birrer, Amara & Refregier 2017; Gilman et al. 2018, 2020; Hsueh et al. 2020) or surface-brightness anomalies in the lensed emission of an extended background galaxy (e.g. Blandford, Surpi & Kundić 2001; Koopmans 2005; Vegetti & Koopmans 2009a; Vegetti, Czoske & Koopmans 2010a; Vegetti et al. 2010b; Rau, Vegetti & White 2013; Vegetti et al. 2014, 2018; Ritondale et al. 2019; Despali et al. 2022; He et al. 2022, 2023; O’Riordan et al. 2023; Wagner-Carena et al. 2023). Such anomalies can be measured, modelled, and traced back to the underlying substructures in the lens galaxy.
In particular, the gravitational-imaging technique developed by Koopmans (2005) and Vegetti & Koopmans (2009a) aims at the detection of individual dark-matter subhaloes in (massive elliptical) lens galaxies at cosmological distances through sophisticated computational modelling of surface-brightness anomalies measured in deep high-resolution imaging of galaxy–galaxy strong gravitational lens systems. The successful detection of two dark-matter subhaloes with masses 3.5 × 109 and 1.9 × 108 M⊙ in lens galaxies at the redshift z = 0.2 and z = 0.9, respectively, reported in Vegetti et al. (2010b, 2012), demonstrates that galactic subhaloes can be identified as localized corrections to a smooth gravitational potential, with the detection threshold depending on the angular resolution and the signal-to-noise ratio of the available data. Recently, Despali et al. (2018) showed that this technique can also be used to detect individual line-of-sight haloes and Vegetti et al. (2018) derived the resulting constraints on the properties of dark matter. However, despite these encouraging results, the current mass-detection threshold of the gravitational-imaging technique (∼108 M⊙ for HST-observations of SLACS lenses; Vegetti et al. 2014, 2018) still lies above the mass regime in which the alternative dark-matter models could be easily distinguished.
In an attempt to detect gravitational signatures of dark-matter subhaloes with masses below this detection threshold, an alternative statistical approach has emerged in the literature over the recent years (Bus 2012; Hezaveh et al. 2016; Chatterjee & Koopmans 2018; Diaz Rivero, Cyr-Racine & Dvorkin 2018a; Díaz Rivero et al. 2018b; Cyr-Racine, Keeton & Moustakas 2019). Within this framework, the low-mass subhaloes in the lens galaxy are viewed as a statistical population and modelled in terms of the subgalactic matter power spectrum. The latter can potentially be constrained from a statistical analysis of the collectively induced surface-brightness anomalies in the extended lensed images of galaxy–galaxy strong gravitational lens systems. Despite the considerable number of recent theoretical studies (e.g. Brennan et al. 2019; Çaǧan Şengül et al. 2020; Galan et al. 2022; Dhanasingham et al. 2023), there has not been a single attempt to apply this approach to real observational data. This paper is part of a series aiming to close this gap.
In the companion paper, Bayer et al. (2023), hereafter referred to as Paper I, we have introduced the methodology to measure the power spectrum of the hypothetical surface-brightness anomalies in high-resolution Hubble Space Telescope (HST) observations of galaxy–galaxy strong gravitational lens systems and applied it to a subsample from the Sloan Lens ACS Survey (SLACS; Bolton et al. 2008). In this study, our aim is to extend the methodology and relate this power-spectrum measurement to the statistical properties of the underlying low-mass structures in the respective lens galaxy. The main challenge when applying this statistical approach to observational data, however, is that in reality the detected surface-brightness anomalies may arise not only from the presence of low-mass dark-matter subhaloes (or field haloes along its line of sight, see e.g. Despali et al. 2018; Çaǧan Şengül et al. 2020; Minor et al. 2021), but from many different kinds of inhomogeneities in the total (dark and baryonic) mass distribution of the lens galaxy, such as e.g. small-scale dark-matter density fluctuations, globular clusters, tidal streams, dynamical distortions in the baryonic mass distribution or edge-on discs (Vegetti et al. 2014; Hsueh et al. 2016; Gilman et al. 2017; Hsueh et al. 2017, 2018; Galan et al. 2022), which are not explicitly included in the smooth-lens model.
Thus, instead of considering the subhaloes only, we set out to constrain the statistical properties of the overall density fluctuations (i.e. deviations from the best-fitting smooth-lens model) in the total mass distribution of (massive elliptical) lens galaxies on the subgalactic 1–10 kpc scales, which we model in terms of Gaussian random field (GRF) potential perturbations (Bus 2012; Hezaveh et al. 2016; Chatterjee & Koopmans 2018). As a proof of concept, we apply the complete methodology to high-resolution HST/WFC3/F390W observations of the SLACS lens system SDSS J0252+0039 (Koopmans 2012). This analysis leads to the first-ever observational constraints on the power spectrum of GRF potential perturbations in a massive elliptical (lens) galaxy and on the corresponding subgalactic matter power spectrum. Future work will compare these constraints with predictions from hydrodynamical simulations (such as Illustris or EAGLE; Vogelsberger et al. 2014; Schaye et al. 2015), which incorporate both dark and baryonic matter.
This paper is structured as follows. Section 2 provides a concise description of the GRF-formalism applied to model the small-scale density fluctuations in the lens galaxy and presents our procedure adopted to uncover, quantify and interpret the associated surface-brightness anomalies in the lensed images. In Section 3, we present the analysed HST observations of our pilot lens system SDSS J0252+0039. Sections 4, 5, and 6 describe the lens-galaxy subtraction, smooth lens modelling and the power-spectrum analysis of the residual surface-brightness fluctuations, respectively. In Section 7, we generate a catalogue of mock perturbed lensed images to be compared with the real observations. The inferred observational upper-limit constraints on the (projected) subgalactic matter power spectrum are presented in Section 8 and compared to predictions from the ΛCDM model in Section 9. The final Section 10 provides conclusions and implications of this work for further research.
For a consistent comparison of the inferred lens model with earlier studies by Vegetti et al. (2014), throughout this work, we assume the following cosmology: H0 = 73 km s−1 Mpc−1, ΩM = 0.25 and |$\Omega _\Lambda = 0.75$|. Given this cosmology, 1 arcsec corresponds to 4.11 kpc at the redshift of the analysed lens galaxy (zL = 0.280) and 7.88 kpc at the redshift of the source galaxy (zS = 0.982).
2 METHODOLOGY
Low-mass structures in lens galaxies are commonly modelled in terms of corrections to the best-fitting smoothly varying parametric model of the lensing potential (e.g. Koopmans 2005; Vegetti & Koopmans 2009a). While previous studies have mainly focused on individual massive dark matter subhaloes (e.g. Vegetti et al. 2012, 2014) and line-of-sight haloes (Despali et al. 2018), here we consider the overall departures from a smooth-lens model, i.e. small-scale density fluctuations in the lens galaxy arising not only from subhaloes or haloes along the line of sight, but from the complex total (dark and baryonic) mass distribution of real (massive elliptical lens) galaxies in general.
In the companion Paper I (Bayer et al. 2023), we have introduced a methodology allowing us to measure the power spectrum of the associated collectively induced surface-brightness anomalies in high resolution HST-observations of galaxy–galaxy strong gravitational lens systems. Here, we extend this methodology and relate the measured power spectrum to the statistical properties of the underlying small-scale density fluctuations in the total mass distribution of the lens galaxy. We model these as Gaussian random field (GRF) potential perturbations superposed on the best-fitting smoothly varying lensing potential, as first suggested by Bus (2012), Hezaveh et al. (2016) and Chatterjee & Koopmans (2018).
In Section 2.1, we first summarize the concept of the GRF potential perturbations in a lens galaxy and the resulting surface-brightness anomalies in the extended lensed images of the background source galaxy. Section 2.2 then outlines the complete procedure to measure the power spectrum of the surface-brightness anomalies, relate it to the statistical properties of the underlying GRF potential perturbations in the lens galaxy and, finally, infer constraints on the subgalactic matter power spectrum.
2.1 GRF potential perturbations in the lens galaxy
We treat the hypothetical small-scale density fluctuations in the lens galaxy as a statistical ensemble and model the associated potential corrections as a homogeneous and isotropic Gaussian random field superposed on a smoothly varying lensing potential (Bus 2012; Hezaveh et al. 2016; Chatterjee & Koopmans 2018; Vernardos & Koopmans 2022). More formally, we assume that the lensing potential of the lens galaxy |$\psi ({{\boldsymbol x}})$| as a function of the position |${{\boldsymbol x}}$| in the lens plane can to the first order be approximated by the sum of a smoothly varying parametric component |$\psi _{0}({{\boldsymbol x}})$| and a Gaussian potential perturbation field |$\delta \psi _{\rm {GRF}}({{\boldsymbol x}})$|:
with |$\langle \delta \psi _{\rm {GRF}}({{\boldsymbol x}})\rangle =0$| and no covariance between |$\delta \psi _{\rm {GRF}}({{\boldsymbol x}})$| and |$\psi _{0}({{\boldsymbol x}})$|.
In this approximation, the resulting deflection angle |$\alpha ({{\boldsymbol x}})$| can be separated into the deflection caused by the smooth component of the lensing potential |$\boldsymbol \alpha _{0}({{\boldsymbol x}})$| and a perturbation |$\boldsymbol {\delta \alpha }_{\rm {GRF}}({{\boldsymbol x}})$| caused by the differential lensing effect of the GRF potential perturbations:
The differential deflection-angle field |$\boldsymbol {\delta \alpha }_{\rm {GRF}}({{\boldsymbol x}})$| can be in this case directly linked to the underlying potential perturbations |$\delta \psi _{\rm {GRF}}({{\boldsymbol x}})$|:
independently of the smooth lensing component |$\psi _{0}({{\boldsymbol x}})$|. Similarly, the associated convergence-perturbation field |$\delta \kappa _{\rm {GRF}}({{\boldsymbol x}})$| (i.e. surface-mass density perturbations in units of the critical surface-mass density for lensing) is directly related to the potential-perturbation field |$\delta \psi _{\rm {GRF}}({{\boldsymbol x}})$| via the Poisson equation:
The perturbative effect of such Gaussian potential fluctuations |$\delta \psi _{\rm {GRF}}({{\boldsymbol x}})$| in the lens galaxy leads to observable surface-brightness anomalies |$\delta I_{\rm {GRF}}({{\boldsymbol x}})$| in the lensed images of the background source galaxy, measured with respect to the smooth-lens model |$I_0({{\boldsymbol x}})$|:
where |$S({{\boldsymbol y(x)}})$| stands for the unlensed intrinsic surface brightness of the source galaxy at the corresponding location in the source plane and the mapping between the source and the lens plane |${{\boldsymbol y(x)}}$| is given by the lens equation:
A crucial feature of a Gaussian random field is that its properties are entirely characterized by the second-order statistics. Hence, the statistical behaviour of the considered Gaussian potential perturbations |$\delta \psi _{\rm {GRF}}({{\boldsymbol x}})$| is fully described by the 2-point correlation function or, alternatively, its Fourier transform – the power spectrum. Following Chatterjee & Koopmans (2018), we assume the power spectrum of |$\delta \psi _{\rm {GRF}}({{\boldsymbol x}})$| to be isotropic and to follow a power law:
The amplitude A is related to the total variance of the potential perturbations |$\sigma ^2_{\delta \psi }$| (integrated between the spatial scales equal to the inverse of the image length and the inverse of the pixel scale, respectively), while the power-law slope β describes the distribution of this variance over the different length scales (i.e. k-modes) and, thus, determines the scale dependence of the hypothetical small-scale structure in the lens galaxy. For a science image with a side length L (measured in arcsec), we determine A by specifying the overall variance of the potential perturbations |$\sigma ^2_{\delta \psi }\equiv \langle (\delta \psi -\langle \delta \psi \rangle)^2 \rangle = \langle \delta \psi ^2 \rangle$| in the considered field of view:
where the integration is performed over the corresponding Fourier grid and the wavenumbers kx and ky, measured in arcsec−1, are calculated according to the convention in which the wavenumber k is equal to the reciprocal wavelength λ:
of the associated harmonic wave e−2πik · x in the Fourier representation of the GRF field. Substituting Pδψ(k) in equation (8) with equation (7) and replacing the integrals by a summation over discrete pixels with the size dkx = dky = L−1 in k-space finally leads to the following normalization condition:
As can be seen from equation (10), this normalization depends not only on the specific combination of β and |$\sigma ^2_{\delta \psi }$|, but also on the field-of-view and the pixel scale of the analysed image.
Using the methodology developed in this work, we seek to derive observational constraints on the variance |$\sigma ^2_{\delta \psi }$| and the slope β of the power-law power spectrum Pδψ(k), as defined in equations (7–10), from the power spectrum of the associated surface-brightness anomalies |$\delta I_{\rm {GRF}}({{\boldsymbol x}})$| measured in high-resolution HST-observations of galaxy–galaxy strong gravitational lens systems. From now on, we refer to a specific combination of the parameters |$\sigma ^2_{\delta \psi }$| and β as a matter-power-spectrum model.
2.2 Analysis overview
Our approach consists of three main components:
First, we extract the surface-brightness anomalies from the extended lensed images of an observed lens system and quantify their statistical properties in terms of the power spectrum PδI(k), see Paper I for a thorough discussion of the methodology and the modelling degeneracies. To this end, we perform the following steps:
HST-observations and data reduction by means of the drizzlepac package (Gonzaga et al. 2012), see Section 3;
modelling and subtraction of the surface-brightness contribution from the foreground lens galaxy using galfit (Peng et al. 2002), see Section 4;
simultaneous reconstruction of the smooth lensing potential |$\psi _{0}({{\boldsymbol x}})$| and the intrinsic surface-brightness distribution of the background source galaxy |$S({{\boldsymbol y}})$|, using the adaptive and grid-based Bayesian lens-modelling technique by Vegetti & Koopmans (2009a), see Section 5;
statistical quantification of the residual surface-brightness fluctuations in the lensed images in terms of the azimuthally averaged power spectrum PδI(k), see Section 6.1;
estimation and correction for the noise power spectrum, including both the sky background and the flux-dependent photon shot noise, to obtain an upper limit on the power spectrum of the surface-brightness anomalies in the observed lensed images; see Section 6.2.
To be conservative in coping with the modelling degeneracies discussed in Paper I, we treat any residual surface-brightness fluctuations remaining in the lensed images after the lens-galaxy subtraction, the smooth-lens modelling and the noise correction as an upper limit on the effect caused by the density fluctuations in the lens galaxy.
Second, we generate a catalogue of mock lensed images perturbed by GRF potential fluctuations |$\delta \psi _{GRF}({{\boldsymbol x}})$| with known statistical properties (i.e. the total variance |$\sigma ^2_{\delta \psi }$| and the power-law slope β of the GRF power spectrum Pδψ(k), as defined in equations 7–10) for a systematic comparison with the observed data. Subsequently, we subtract the best-fitting smooth lens model obtained for data from the mocks and compute the power spectrum of the residual surface-brightness fluctuations; see Section 7.
Third, we compare the power spectra of surface-brightness anomalies in the mock perturbed lensed images to the upper limit inferred from the observational data and estimate the exclusion probability for each considered matter-power-spectrum model, i.e. combination of |$\sigma ^2_{\delta \psi }$| and β, given the observational measurement. Based on these results, we additionally derive upper-limit constraints on the corresponding power spectrum of perturbations in the deflection angle Pδα(k) and the convergence (i.e. projected mass density) Pδκ(k); see Section 8.
To comply with the Gaussianity assumption for the hypothetical overall small-scale inhomogeneities in the total mass distribution of the massive elliptical (lens) galaxies, we intend to apply our method solely to lens systems for which the effect of individual dominant subhaloes and satellite galaxies has already been modelled or excluded in previous studies. Moreover, in the future, we plan to additionally use hydrodynamical simulations to thoroughly test the validity of this assumption.
3 OBSERVATIONS AND DATA REDUCTION
The mass of the smallest structures that can be recovered in a galaxy-scale halo using the method of gravitational imaging (Koopmans 2005; Vegetti & Koopmans 2009a) is limited by the sensitivity of the observational set-up to tiny surface-brightness perturbations in the lensed images. As argued by Blandford et al. (2001), Koopmans (2005) and explicitly demonstrated by Rau et al. (2013), the surface-brightness anomalies induced in the lensed images by a given density fluctuation in the lens galaxy are enhanced when the lensed background source is either compact or highly structured. Taking into account that blue star-forming regions are considerably more structured than the redder old stellar populations, the choice of a lens system with a star-forming source galaxy combined with the selection of an ultraviolet observational filter might substantially improve the overall sensitivity of our approach.
Following this idea, we apply our methodology to HST/WFC3/F390W-imaging data (with the central wavelength of 390 nm, referred to as U-band in the remainder of the paper) of the galaxy–galaxy strong gravitational lens system SDSS J0252+0039 from the SLACS Survey. This lens system consists of a massive elliptical lens galaxy at redshift zl = 0.280 and a blue star-forming source galaxy at redshift zs = 0.982 (Bolton et al. 2008). The choice of SDSS J0252+0039 is motivated by the large mass of its lens galaxy (total mass of 1011.25 M⊙ inside the Einstein radius), the high surface-brightness gradient of the star-forming background galaxy and the relatively simple lens geometry which reduces the uncertainty in the lens modelling (Auger et al. 2009).
As one of the SLACS gravitational lens systems, SDSS J0252+0039 has already been modelled, based on HST-observations at near-infrared (in F160W: Program 11202, PI Koopmans) and optical wavelengths (in F814W: Program 10866, PI Bolton and in F606W: Program 11202, PI Koopmans), see Bolton et al. (2008) and Auger et al. (2009, 2010). In particular, Vegetti et al. (2014) applied the technique of gravitational imaging (Koopmans 2005; Vegetti & Koopmans 2009a) to the HST/ACS data in the V and I bands to search for signatures of individual massive subhaloes in the lens galaxy. No gravitational signatures of localized substructure have been identified above the mass-detection threshold of ∼108 M⊙. Here, we take advantage of the enhanced source gradient in the U band observations and search for the collectively induced gravitational imprints of much smaller density fluctuations.
The HST/WFC3/F390W-observations of SDSS J0252+0039 were carried out on 2013 August 25, (Program 12898; PI Koopmans, Koopmans 2012). We retrieve the acquired eight dithered flat-field calibrated images from the MAST archive.1 and apply the drizzle method (Gonzaga et al. 2012) to combine them into the final science image. For this, we use the astrodrizzle task from the drizzlepac package in its default configuration, with the original HST/WFC3 rotation, the output-pixel size of 0.0396 arcsec and the drop size equal to the original pixel size (final pixfrac = 1). Our analysis is based on an image cutout of 121 by 121 pixels (corresponding to an area of 4.48 arcsec on a side) centred on the brightest pixel of the lens galaxy, which we present in Fig. 1.

HST/WFC3/F390W-imaging of the strong gravitational lens system SDSS J0252+0039: the drizzled science image with the side length of 4.48 arcsec.
A known effect of the drizzling procedure are the correlations between adjacent pixels in the output science image (see Casertano et al. 2000). These correlations are additionally enhanced by the effect of the charge-transfer inefficiency (CTI) in the HST/WFC3/UVIS-CCDs (see e.g. Baggett, Gosmeyer & Noeske 2015). To investigate and possibly alleviate these instrumental effects, we repeat the drizzling procedure with different settings (i.e. different values of the final pixfrac parameter and the output pixel scale) and compare the results to those obtained for the default configuration. Moreover, we perform a comparison between the power spectra obtained by drizzling either the original or the CTI-corrected flat-field calibrated exposures for a sample of blank-sky cutouts (located in vicinity to the lens system). Based on this analysis, discussed in detail in Paper I, we decide to proceed by drizzling the original dithered exposures in the default configuration of the astrodrizzle task.
To account for further observational effects present in the science image, we obtain the point-spread function (PSF) of the HST/WFC3/UVIS optics using the PSF-modelling software tinytim 2 (Krist, Hook & Stoehr 2010), with the approximation of the G8V spectral type for the lens galaxy (based on its known magnitudes in the V , I, and H bands from Auger et al. 2009). Even though the TinyTim-PSF might not be a perfect representation of the real telescope optics, any minor deviations from the true PSF and the assumed spectral-energy distribution of the lens system affect the measured power spectra of the residual surface-brightness fluctuations only on scales below the full width at half maximum (FWHM) of the PSF, corresponding to wavenumbers k ⪆ FWHM−1 = (0.07 arcsec)−1 ≈ 14 arcsec−1, which is beyond the regime considered in this work (see Section 6.2).
Finally, for visual purposes and in order to assess the possible presence of dust in the studied lens system, we use the stiff 3 package to create a colour-composite image of SDSS J0252+0039. The resulting Fig. 2 shows a smooth early-type lens galaxy and uniformly blue lensed images of the background galaxy. A visual inspection of this image does not identify any indications of dust extinction either in the lens or the source galaxy. This conclusion is confirmed by a more quantitative dust analysis, discussed in Appendix A.

Colour-composite image of the strong gravitational lens system SDSS J0252+0039 combining the new HST/WFC3/F390W-photometry (blue) with the already existing HST-observations in the visual (F814W; green) and infrared (F160W; red) bands, obtained using the stiff software.
4 LENS-GALAXY SUBTRACTION
The Einstein radius of SDSS J0252+0039 (∼1 arcsec) is comparable to the effective radius of its massive elliptical lens galaxy. Consequently, the lensed images overlap with the inner region of the lens galaxy, as can be seen in Fig. 1. To correct for the light contribution from the lens galaxy, we model its surface-brightness distribution with a combination of several Sérsic components (Sérsic 1963) – a parametric class of smooth light profiles widely used to model early-type galaxies.
The applied lens-galaxy-subtraction procedure is illustrated in Fig. 3. We perform the modelling by means of the two-dimensional fitting algorithm galfit (Peng et al. 2002) after masking out all pixels covering the lensed images, as shown in the top row. The best-fitting galfit model, shown in the bottom left-hand panel of Fig. 3, consists of two Sérsic components with Sérsic indices n = 4.44 and n = 0.09 and effective radii 1.2 and 0.96 arcsec (corresponding to ∼4.8 and ∼3.8 kpc at the redshift of the lens), respectively. The small value of the second Sérsic index suggests that the light distribution of the lens galaxy can be well described by the standard de Vaucouleurs profile (i.e. Sérsic profile with n = 4) and a diffuse stellar halo. As an alternative, we additionally carry out the lens-galaxy subtraction using the b-spline algorithm (Bolton et al. 2006), but choose to continue our analysis based on the galfit model due to the slightly higher Bayesian evidence of the resulting best-fitting smooth-lens model, see Section 5. The final lens-galaxy-subtracted image is presented in the bottom right-hand panel of Fig. 3. This residual image constitutes our best estimate of the lensed emission from the background galaxy only and is as such used for the smooth-lens modelling.

Lens-galaxy subtraction for SDSS J0252+0039 using galfit. Top row: the drizzled science image in F390W (left-hand panel) and the applied mask, covering the lensed source emission (right-hand panel). Bottom row: the best-fitting galfit model of the surface-brightness distribution in the lens galaxy (left-hand panel) and the lens-galaxy-subtracted residual image (right-hand panel).
5 SMOOTH LENS MODELLING
In search of possible surface-brightness anomalies in the lensed images, we first model the studied lens system assuming that the global mass distribution in the lens galaxy is smooth and to first order well described by the power-law elliptical mass-distribution model (PEMD, Barkana 1998) in an external shear field. This smooth lens modelling is performed by means of the adaptive and grid-based Bayesian lens-modelling technique of Vegetti & Koopmans (2009a) that allows us to find the best-fitting parameter values of the PEMD model for the lens galaxy and, simultaneously, reconstruct the unlensed surface-brightness distribution of the background galaxy on an adaptive grid in the source plane. We follow the parametrization used by Vegetti & Koopmans (2009a) and model the projected mass density of the lens galaxy in terms of the convergence κ:
This is a function of the position in the lens plane (x, y) and has the following parameters – the lens strength b, the position angle of the major axis θ (with respect to the original telescope rotation), the (minor to major) axis ratio q, the (three-dimensional) mass-density slope γ (with γ = 2 for the isothermal case), the centre coordinates of the mass distribution in the lens plane x0 and y0, the external shear strength Γ and its position angle Γθ. This parametrization tries to reduce the degeneracy between the mass enclosed by the lensed images, the axis ratio of the lens galaxy and its mass-density slope.
Several options for the inversion of the lensing operation are available to choose from when using the smooth-lens-modelling code by Vegetti & Koopmans (2009a). To prevent overfitting, the reconstruction of the background galaxy can be regularized by applying an adaptive or non-adaptive, variance, gradient, or curvature source-plane regularization. Furthermore, one can choose the source-grid resolution, which is characterized by the fraction and spacing of pixels that are cast back from the image plane to the source plane when generating a grid for the source reconstruction. More specifically, the source-grid resolution is controlled by the parameter n that sets the linear size of a square in the image plane out of which only one pixel is cast back to the source plane. For example, if n = 3 only one pixel out of each contiguous 3 × 3-pixel area is used to create the Delauney-tesselation grid in the source plane. We note that, even if n > 1, all pixels are still used to calculate the Bayesian evidence and compare the model to the data.
As pointed out by Suyu et al. (2006) and Vegetti et al. (2014), the optimal choice for the regularization and the source-grid resolution depends crucially on the smoothness of the modelled lensed images and may be different for each specific lens system. Thus, the common practice is to perform the smooth lens modelling for different combinations of the above options, and find the optimal settings based on the highest Bayesian evidence. Our tests in this respect show that the chosen source-grid properties hardly affect the best-fitting parameter values of the lensing potential. We find, however, that they have a significant effect on the reconstructed unlensed surface-brightness distribution of the background source galaxy and, consequently, on the level of the residual surface-brightness fluctuations remaining in the lensed images after the best-fitting smooth-lens model has been subtracted.
We achieve the highest value of the marginalized Bayesian evidence when the smooth lens modelling is performed by casting back each pixel from the lens plane to the source plane (referred to as n = 1) and using the adaptive gradient source-grid regularization, see Fig. 4. The best-fitting parameter values of the lensing potential are presented in Table 1. These are in a good agreement with the parameter values inferred by Vegetti et al. (2014) in an earlier analysis of the HST/ACS/F814W-data (I band), except for the apparent discrepancy in the position angle of the major axis θ. This discrepancy can be explained by the rotational invariance owing to the nearly spherical symmetry of the modelled lens galaxy (axis ratio q very close to 1) and the negligible external shear Γ, hardly altering the lensing potential. All in all, the best-fitting parameter values inferred in both bands indicate that SDSS J0252+0039 has a nearly spherical isothermal mass-density profile (γ ≈ 2) with the Einstein radius of ΘE ≈ b ≈ 1 arcsec.

Smooth lens modelling of SDSS J0252+0039 in the U-band (HST/WFC3/F390W) by means of the adaptive and grid-based Bayesian lens-modelling technique of Vegetti & Koopmans (2009a), performed with the source-grid resolution n = 3, i.e. casting back only one pixel out of each contiguous 3 × 3-pixel area from the image plane to the source plane. Top row: the lens-galaxy-subtracted image overlaid with a mask, used as input for the smooth lens modelling (left-hand panel) and the best-fitting smooth-lens model of the lensed images (right-hand panel). Bottom row: the reconstructed intrinsic surface-brightness distribution of the background galaxy (left-hand panel) and the residual image showing the remaining surface-brightness fluctuations in the lensed images, possibly caused by small-scale mass structure in the lens galaxy (right-hand panel).
Parameter values of the best-fitting PEMD model for the lens galaxy SDSS J0252+0039 inferred based on the U-band (F390W) imaging, in comparison to the I-band (F814W) reconstruction by Vegetti et al. (2014). The optimized parameters are the lens strength b, the position angle θ (with respect to the original telescope rotation), the axis ratio q, the (three-dimensional) mass-density slope γ, the external shear strength Γ and its position angle Γθ. In both bands, the reconstruction was carried out using the adaptive gradient source regularization and casting back each pixel (n = 1). The typical statistical errors of our U-band reconstruction are of the order 10−1 for the angles and 10−3 for the remaining parameters.
Filter . | b [arcsec] . | θ [deg] . | q . | γ . | Γ . | Γθ [deg] . |
---|---|---|---|---|---|---|
F390W | 0.996 | 150.1 | 0.978 | 2.066 | −0.015 | 81.4 |
F814W | 1.022 | 26.2 | 0.943 | 2.047 | 0.009 | 101.8 |
Filter . | b [arcsec] . | θ [deg] . | q . | γ . | Γ . | Γθ [deg] . |
---|---|---|---|---|---|---|
F390W | 0.996 | 150.1 | 0.978 | 2.066 | −0.015 | 81.4 |
F814W | 1.022 | 26.2 | 0.943 | 2.047 | 0.009 | 101.8 |
Parameter values of the best-fitting PEMD model for the lens galaxy SDSS J0252+0039 inferred based on the U-band (F390W) imaging, in comparison to the I-band (F814W) reconstruction by Vegetti et al. (2014). The optimized parameters are the lens strength b, the position angle θ (with respect to the original telescope rotation), the axis ratio q, the (three-dimensional) mass-density slope γ, the external shear strength Γ and its position angle Γθ. In both bands, the reconstruction was carried out using the adaptive gradient source regularization and casting back each pixel (n = 1). The typical statistical errors of our U-band reconstruction are of the order 10−1 for the angles and 10−3 for the remaining parameters.
Filter . | b [arcsec] . | θ [deg] . | q . | γ . | Γ . | Γθ [deg] . |
---|---|---|---|---|---|---|
F390W | 0.996 | 150.1 | 0.978 | 2.066 | −0.015 | 81.4 |
F814W | 1.022 | 26.2 | 0.943 | 2.047 | 0.009 | 101.8 |
Filter . | b [arcsec] . | θ [deg] . | q . | γ . | Γ . | Γθ [deg] . |
---|---|---|---|---|---|---|
F390W | 0.996 | 150.1 | 0.978 | 2.066 | −0.015 | 81.4 |
F814W | 1.022 | 26.2 | 0.943 | 2.047 | 0.009 | 101.8 |
However, despite a remarkably good agreement between the model and data, the source reconstruction obtained using n = 1 turns out to be under-regularized, i.e. most of the surface-brightness fluctuations in the lensed images, and partially even the observational noise, are ‘absorbed’ in the source structure. The power spectrum of the residual surface-brightness fluctuations (Section 6) remaining in the lensed images after the subtraction of the best-fitting smooth-lens model lies below the noise power spectrum for both adaptive and non-adaptive source-grid regularization (see Paper I). While a thorough analysis of the degeneracy between the fluctuations in the lensing potential and the intrinsic source structure is the subject of a future work (but see also Chatterjee 2019; Vernardos & Koopmans 2022), we continue our study based on the (more conservative) source reconstruction obtained using a lower grid resolution (n = 3), while keeping fixed the parameter values inferred from the highest-resolution model (n = 1), see Fig. 4 and Table 1. The particular choice of n = 3 is motivated by the results of our tests carried out for a mock lens system mimicking SDSS J0252+0039 (see Paper I). These tests suggest that the degeneracy between the density fluctuations in the lens galaxy and the intrinsic source structure is significantly reduced when the modelling is performed with n = 3 (or higher).
The obtained model of the lensed images (top right-hand panel of Fig. 4), which would be observed if the lensing potential was indeed smooth, is subsequently subtracted from the observed data to uncover residual surface-brightness fluctuations. These indicate a possible deviation of the true lensing potential from the assumed smooth (PEMD) model. As can be seen from the residual image in the bottom right-hand panel of Fig. 4 and the corresponding signal-to-noise ratio image presented in Fig. 5, the revealed residual surface-brightness fluctuations significantly exceed the noise level.

The signal-to-noise ratio of the residual surface-brightness fluctuations remaining in the lensed images of SDSS J0252+0039 after subtraction of the smooth-lens model obtained with the source-grid resolution of n = 3.
6 POWER-SPECTRUM ANALYSIS OF THE SURFACE-BRIGHTNESS ANOMALIES IN THE LENSED IMAGES
In this section, we estimate the power spectrum of surface-brightness anomalies in SDSS J0252+0039, defined as residual surface-brightness fluctuations in the lensed images caused by the presence of density fluctuations in the lens galaxy (i.e. deviations of the real mass distribution from the best-fitting PEMD model). To start with, in Section 6.1, we take into account all the residual surface-brightness fluctuations revealed in the lensed images after the subtraction of the lens-galaxy light and the best-fitting smooth-lens model. We note, however, that such residual fluctuations might originate not only from the possible mass structure in the lens galaxy, but also from other phenomena, for example the observational noise, uncertainties in the PSF-model or intrinsic structure in the lensed source that has not been recovered in the lens-modelling procedure due to the choice of a relatively low source-grid resolution (n = 3). Hence, in order to constrain the surface-brightness anomalies resulting solely from the density fluctuations in the mass distribution of the lens galaxy, it is crucial to quantify the other effects, in particular the observational noise (see Section 6.2).
6.1 Power spectrum of residual surface-brightness fluctuations
To quantify the overall residual surface-brightness fluctuations, we follow the procedure introduced in Paper I and compute the azimuthally averaged power spectrum of the residual image within the mask covering the lensed images (shown in Fig. 3). To achieve this, we set the flux values of the pixels located outside the mask to zero and calculate the two-dimensional discrete Fourier transform (DFT) of the masked residual image, using the Python package numpy.fft 4 The squared magnitude of the Fourier coefficient assigned to each pixel yields the two-dimensional power spectrum. Assuming isotropy of the modelled potential perturbations |$\delta \psi ({{\boldsymbol x}})$|, we average the power-spectrum values along a set of ten equidistant concentric annuli spanning from kmin = 0.88 to kmax = 16.79 arcsec−1 (corresponding to the spatial scales between λmin = 0.22 and λmax = 4.65 kpc at the redshift of the lens galaxy zl = 0.280).
6.2 Noise correction
In order to characterize the noise properties in the analysed HST-image, we create a sample of 20 selected blank-sky regions with the same size (121 by 121 pixels), located in proximity to the lens system. The first rough estimate of the noise level is given by the standard deviation of the flux values in this blank-sky sample: σsky = 0.002 e− sec−1. However, this estimate does not take into account the photon-shot (Poisson-distributed) noise, which depends on the number of detected photons and, consequently, varies from pixel to pixel.
A more precise description of noise properties, including both the sky-background and the photon-shot noise, is provided by the noise-sigma map, which quantifies the standard deviation of noise for each pixel separately. Considering that the Poisson variance of the photon counts is similar to the measured number of photons and the raw HST images of the lens are drizzle-combined using an inverse-variance map weighting, we construct the noise-sigma map for our drizzled HST science image according to the following formula:
where N is the number of photo-electrons per second detected in a particular pixel (after the sky-background subtraction) and W is the weight of this pixel taken from the weight map of our image provided by the drizzling pipeline. This noise-sigma map is also used in the Bayesian smooth-lens modelling procedure, presented in Section 5. Since the Poisson noise approaches Gaussian noise for large number counts, as is the case for the studied image, in the remaining part of our analysis, we approximate the photon-shot noise by an additive Gaussian noise |$\mathcal {N}(0, \sigma _{n}^2)$| with a variance |$\sigma _{n}^2$| adapted to the flux value in a particular pixel.
However, due to noise correlations in the drizzled images and the charge-transfer inefficiency (CTI), discussed in more detail in Paper I, the noise correction of the residual surface-brightness fluctuations requires a more extended approach. To mimic the true noise properties in the science image, we use the selected blank-sky regions to generate a sample of 20 scaled sky-background realizations, which account for both the realistic noise-correlation pattern and the spatially varying flux-dependent photon-shot noise. This is implemented by first dividing the flux values of the blank-sky regions by their standard deviation and subsequently multiplying them by the noise-sigma map of the science image (equation 12). For consistency reasons, these scaled sky-background realizations are finally overlaid with the same mask as the one used in the analysis of the real data (Section 6). The average power spectrum measured in this sample is used as our best estimate for the total noise power spectrum of the science image (see Paper I for a more thorough discussion).
Finally, we use this estimated total noise power spectrum to perform the noise correction of the residual surface-brightness fluctuations. We assume that the observational noise and the potential fluctuations δψ, perturbing the smooth lensing potential, are independent stochastic processes and consider the corresponding power spectra to be additive. This allows us to subtract the estimated noise power spectrum from the power spectrum of the total residual surface-brightness fluctuations, as illustrated in Fig. 6. The difference of these two power spectra is treated in our further analysis as an upper limit to the power spectrum of the surface-brightness anomalies PδI (k) caused by the hypothetical small-scale mass structures in the lens galaxy.

Power spectrum of residual surface-brightness fluctuations remaining in the lensed images of SDSS J0252+0039 after the lens-galaxy subtraction, the smooth lens modelling and the noise correction. The difference between the power spectrum of residual surface-brightness fluctuations after the smooth lens modelling with n = 3 (blue line) and the estimated total noise power spectrum (green line) constitutes our upper limit on the power spectrum of surface-brightness anomalies due to mass structure in the lens galaxy SDSS J0252+0039 (red line).
However, a comparison between the power spectrum of the revealed residual surface-brightness fluctuations with the estimated noise power spectrum shows that for the highest considered k-values (corresponding to scales below three pixels), the residual fluctuations reach the noise level. This indicates that no surface-brightness anomalies are detected on these scales. For this reason, in our further analysis, we consider only the perturbation wave numbers ranging from kmin = 0.88 to kmax = 7.95 arcsec−1, which corresponds to the spatial scales between λmin = 0.52 and λmax = 4.65 kpc at the redshift of the lens galaxy, see Fig. 6. Performing the analysis on scales above three pixels together with our choice of n = 3 in the smooth-lens-modelling procedure allows us to neglect the effects of small errors in the PSF (with FWHM = 0.07 arcsec), the correlations between adjacent pixels due to drizzling and the possible residual errors in the source-light modelling.
7 CATALOGUE OF PERTURBED LENSED IMAGES
In order to interpret the power spectrum of the surface-brightness anomalies measured in the observed image of SDSS J0252+0039, in this section, we generate a catalogue of simulated lensed images which mimic these observations but are perturbed by Gaussian potential perturbations |$\delta \psi _{\rm {GRF}}({{\boldsymbol x}})$| from the power-law power spectrum |$P_{\delta \psi }(k;\sigma ^2_{\delta \psi }, \beta)$| (see Section 2.1) on a grid containing 100 × 100 different combinations of β and |$\sigma ^2_{\delta \psi }$|. The considered values of β are equidistant within the interval [3, 8], whereas |$\sigma ^2_{\delta \psi }$| (within the studied field of view) is varied by obtaining 100 values evenly spaced in the logarithmic range [10−6, 10−1].
For each considered combination of β and |$\sigma ^2_{\delta \psi }$|, we generate a pixelated realization of |$\delta \psi _{\rm {GRF}}({{\boldsymbol x}})$| from the corresponding power spectrum |$P_{\delta \psi }(k;\sigma ^2_{\delta \psi }, \beta)$| (equations 7 and 10) and superimpose it on the best-fitting smooth lensing potential |$\psi _0({{\boldsymbol x}})$| of SDSS J0252+0039, obtained for the observational data in Section 5. We then apply this perturbed lensing potential |$\psi _0({{\boldsymbol x}}) + \delta \psi _{\rm {GRF}}({{\boldsymbol x}})$| to repeat the lensing operation of the previously reconstructed pixellated surface-brightness distribution of the background source galaxy and finally obtain the resulting perturbed lensed image. Fig. 7 shows three examples of these mock images for different values of |$\sigma ^2_{\delta \psi }$| and β, together with the underlying realizations of |$\delta \psi _{\rm {GRF}}({{\boldsymbol x}})$| and |$\delta \kappa _{\rm {GRF}}({{\boldsymbol x}})$|.

Mock surface-brightness anomalies induced in the gravitational arcs of SDSS J0252+0039 by Gaussian potential perturbations |$\delta \psi _{\rm {GRF}}({{\boldsymbol x}})$| in the lens galaxy, fully characterized by the power spectrum |$P_{\delta \psi }(k) \propto \sigma ^2_{\delta \psi } \times k^{-\beta }$|. Left column: unperturbed lensed images, i.e. the reconstructed source galaxy lensed through the best-fitting smooth-lens model (upper panel) and three examples of perturbed lensed images (for |$\sigma ^2_{\delta \psi } = 2.154 \times 10^{-5}$| and β = 5.5, |$\sigma ^2_{\delta \psi } = 2.783 \times 10^{-4}$| and β = 4.25, |$\sigma ^2_{\delta \psi } = 10^{-3}$| and β = 3; from top to bottom) with a realistic noise realization overlaid for visualization purposes. Middle column: the underlying realizations of |$\delta \psi _{\rm {GRF}}({{\boldsymbol x}})$|. Right column: the corresponding convergence (i.e. surface mass density) perturbations |$\delta \kappa _{\rm {GRF}}({{\boldsymbol x}})=\frac{1}{2}\nabla ^2 \delta \psi _{\rm {GRF}}({{\boldsymbol x}})$|.
To obtain the power spectra of surface-brightness anomalies for these perturbed mock lensed images, we follow the methodology applied to the observational data, i.e. we subtract the best-fitting smooth-lens model from each of the mocks and compute the power spectrum of the residual surface-brightness fluctuations within the same mask and using the same set of ten bins for the azimuthal averaging as for the observational data. We stress that in this paper, for simplicity, we do not model the perturbed lensed images individually, but assume that the same best-fitting smooth lens model that was inferred for the data holds for all the perturbed mocks. We plan to investigate the validity of this assumption in our future paper (Bayer et al. in prep). Moreover, to reduce the sample variance, we consider ten different realizations of |$\delta \psi _{\rm {GRF}}({{\boldsymbol x}})$| for each combination of |$\sigma ^2_{\delta \psi }$| and β and average the resulting power spectra over this sample. Fig. 8 presents a few examples of the averaged mock power spectra. A comparison of these to the real measurement, illustrated in Fig. 8, allows us to infer exclusion probabilities for the considered portion of the parameter space spanned by |$\sigma ^2_{\delta \psi }$| and β, which is discussed in the next section.

Upper limit on the power spectrum of surface-brightness anomalies due to mass structure in the lens galaxy SDSS J0252+0039 (red line), corresponding to the red line in Fig. 6, in comparison to a mock catalogue containing power spectra of surface-brightness anomalies induced by Gaussian potential perturbations with known values of the integrated variance |$\sigma ^2_{\delta \psi }$| and the power-law power-spectrum slope β according to equations (7) and (10) (dashed lines). For clarity of presentation, we show only 9 out of 104 mock power spectra in our catalogue.
8 CONSTRAINTS ON THE SUB-GALACTIC MATTER POWER SPECTRUM
In this section, we first discuss our approach to assess whether a particular matter-power-spectrum model (i.e. combination of |$\sigma ^2_{\delta \psi }$| and β) leads to surface-brightness anomalies that are significantly different from the measured residual surface-brightness fluctuations. Subsequently, we present and discuss the resulting constraints on the statistical properties of the potential perturbations in the lens galaxy SDSS J0252+0039 (Section 8.1). Finally, we express these constraints in terms of the variance of the resulting deflection-angle perturbations (Section 8.2), the dimensionless convergence power spectrum |$\Delta ^{2}_{\delta \kappa }(k)$| and the aperture mass on three specific spatial scales: 0.5, 1, and 3 kpc (Section 8.3).
8.1 Exclusion probability of matter-power-spectrum models
A comparison between the power spectra of the mock surface-brightness anomalies and the power spectrum of the actually observed residual surface-brightness fluctuations enables us to determine exclusion probabilities for different combinations of |$\sigma ^2_{\delta \psi }$| and β. As discussed in Section 2.2, we conservatively treat the estimated residual surface-brightness fluctuations as an upper limit to the anomalies caused solely by potential perturbations in the lens galaxy. Consequently, we rule out a particular matter power-spectrum model only if the power spectrum of the resulting surface-brightness anomalies exceeds the measured power spectrum on all considered scales, i.e. in all five analysed k-bins between kmin = 0.88 and kmax = 7.95 arcsec−1.
In order to determine the exclusion probability, we treat the power-spectrum value measured in the data for each individual k-bin as the expected value of a Gaussian random variable |$P^{D}_{\delta I}(k)$| with the standard deviation (corresponding to the error bars in Fig. 8) estimated based on the variance in our sample of mock noise realizations introduced in Section 6.2. Under this assumption, we first compute the probability that the power spectrum of the mock surface-brightness anomalies for a given matter-power-spectrum model exceeds the actually measured power spectrum for each k-bin separately in the following way:
with |$P^{M}_{\delta I}(k)$| denoting the mock power spectrum for a particular choice of |$\sigma ^2_{\delta \psi }$| and β. Assuming that our exclusion-probability estimates for different perturbation wavenumbers k are statistically independent, the final exclusion probability for a particular matter-power-spectrum model can be calculated as the product of the exclusion probabilities in all considered k-bins:
The inferred exclusion probabilities for the considered range of matter-power-spectrum models (i.e. combinations of |$\sigma ^2_{\delta \psi }$| and β) are presented in Fig. 9. The depicted grid contains 50 different values for |$\sigma ^2_{\delta \psi }$|, evenly spaced in the logarithmic range [10−5.0, 10−2.5], and 100 values for β, equidistant in the interval [3, 8]. Colour coded is the inferred exclusion probability of the corresponding matter-power-spectrum models |$P_{\rm excl} \Big (\sigma ^2_{\delta \psi },\beta \Big)$|, defined as the probability that the resulting surface-brightness anomalies exceed the observed ones on all considered spatial scales. The superposed (black) isoprobability line indicates models with exclusion probability larger than 99 per cent.

Exclusion probabilities |$P_{\rm excl} \Big (\sigma ^2_{\delta \psi },\beta \Big)$| for a subset of considered subgalactic matter-power-spectrum models, inferred from our analysis of the lens system SDSS J0252+0039. Each point of this exclusion plot corresponds to a particular combination of the integrated variance |$\sigma ^2_{\delta \psi }$| and the power-law slope β, assumed to fully characterize the power spectrum of the hypothetical GRF potential perturbations |$P_{\delta \psi }(k) \propto \sigma ^2_{\delta \psi } \times k^{-\beta }$| in the lens galaxy. The superposed black isoprobability line indicates models with exclusion probability larger than 99 per cent. The white contour lines connect matter-power-spectrum models with the same variance of the associated deflection-angle perturbations |$\sigma ^2_{\delta \alpha }$| in the analysed field of view.
Based on the exclusion plot presented in Fig. 9, we rule out matter-power-spectrum models with |$\sigma ^2_{\delta \psi }$| exceeding ∼10−2.5 (on the spatial scales between L−1 and pixel scale−1 and within the range of the power-law power-spectrum slope 3 ≤ β ≤ 8) at the 99 per cent confidence level. This corresponds to the upper bound on the standard deviation of such potential fluctuations σδψ ≤ 0.06. The excluded models would lead to power spectra of surface-brightness anomalies that significantly exceed the observed upper limits on all considered scales. Moreover, we find that for a specific value of the integrated variance |$\sigma ^2_{\delta \psi }$|, the exclusion probability depends on the power-spectrum slope β and, thus, on the exact distribution of the integrated variance between the different spatial scales. Shallower slopes, assigning a larger fraction of |$\sigma ^2_{\delta \psi }$| to small scales, i.e. high k-modes, are found more likely to be excluded. This leads to the conclusion that potential fluctuations on smaller spatial scales have a stronger perturbative effect on the lensed images than potential fluctuations on larger scales.
8.2 Upper limits on the deflection-angle perturbations
To explain the asymptotic shape of the exclusion plot, presented in Fig. 9, we convert our original constraints on the power spectrum of the potential perturbations |$P_{\delta \psi }\Big (k, \sigma ^2_{\delta \psi },\beta \Big)$| into constraints on the power spectrum of the corresponding perturbations in the deflection angle |$P_{\delta \alpha } \Big (k, \sigma ^2_{\delta \psi },\beta \Big)$|, making use of the following relation:
According to equation (15), the slope of |$P_{\delta \alpha } \Big (k, \sigma ^2_{\delta \psi },\beta \Big)$| decreases by 2.0 (i.e. becomes shallower) with respect to the slope of |$P_{\delta \psi }\Big (k, \sigma ^2_{\delta \psi },\beta \Big)$|. The total variance in the differential deflection angle |$\sigma ^2_{\delta \alpha }$| over the analysed field of view is obtained by integrating Pδα(k) over all pixels of the two-dimensional Fourier grid:
The result is presented in Fig. 9, where the overlaid (white) isocontours correspond to the same values of |$\sigma ^2_{\delta \alpha }$|. These isocontours almost perfectly follow the overall shape of the exclusion limits, which indicates a strong correlation between the exclusion probability of a matter-power-spectrum model and the total variance of the associated deflection-angle perturbations. This suggests that |$\sigma ^2_{\delta \alpha }$| is the fundamental quantity that determines the level of the resulting surface-brightness anomalies in the lensed images. Consequently, the exclusion probability is almost insensitive to the slope β of the Pδα(k), i.e. the distribution of the total variance in the deflection angle over the different length scales. Based on Fig. 9, we exclude (with 99 per cent probability) a matter-power-spectrum model if the corresponding total variance in the differential deflection field is larger than 6 × 10−3, independently of the slope β. This insight is valuable for a future analysis of additional lens systems, which instead of considering potential perturbations |$\delta \psi _{\rm {GRF}}({{\boldsymbol x}})$| might more efficiently focus on the corresponding deflection-angle perturbations |$\boldsymbol {\delta \alpha }_{\rm {GRF}}({{\boldsymbol x}})$|. The inferred threshold value of |$\sigma ^2_{\delta \alpha }$|, however, might vary between different lens systems and depend on the chosen field of view, the PSF and the signal-to-noise ratio of the analysed image.
For completeness, we derive the corresponding constraints on the dimensionless deflection-power spectrum defined as follows:
for the spatial scales of 0.5, 1, and 3 kpc, see Fig. 10. These exclude |$\Delta ^{2}_{\delta \alpha }(k)$| larger than 0.001 on all considered spatial scales at the 99 per cent confidence level.

Upper-limit constraints on the dimensionless convergence power spectrum |$\Delta ^{2}_{\delta \kappa }$| and the dimensionless differential-deflection power spectrum |$\Delta ^{2}_{\delta \alpha }$| in the lens galaxy SDSS J0252+0039 for three different subgalactic scales. Each point of the exclusion plot corresponds to a particular subgalactic matter-power-spectrum model |$P_{\delta \psi }(k) \propto \sigma ^2_{\delta \psi } \times k^{-\beta }$|, as specified in Fig. 9. The white contour lines connect models with the same value of |$\Delta ^{2}_{\delta \kappa }$| (left-hand panel) and |$\Delta ^{2}_{\delta \alpha }$| (right-hand panel) on a particular scale. Overlaid in black is the 0.99-contour of the exclusion probability. Note that, β represents the power-spectrum slope of the originally investigated potential perturbations |$\delta \psi _{\rm {GRF}}({{\boldsymbol x}})$|. The slopes of the corresponding power spectra of the convergence perturbations |$\delta \kappa _{\rm {GRF}}({{\boldsymbol x}})$| and the deflection perturbations |$\delta \alpha _{\rm {GRF}}({{\boldsymbol x}})$| decrease by 4 and 2, respectively (in terms of the absolute value).
8.3 Upper limits on the convergence power spectrum
Finally, the inferred constraints on the power spectrum of the potential perturbations Pδψ(k) can be translated into constraints on the associated power spectrum of the convergence (i.e. surface mass density) perturbations Pδκ(k). By expressing equation (4) in Fourier space, Pδκ(k) can be related to Pδψ(k) as follows:
For a future comparison with the ΛCDM predictions, we express these constraints in terms of the dimensionless convergence power spectrum:
which quantifies the contribution of a particular length scale λ = k−1 to the total variance of the convergence perturbations.
We determine |$\Delta ^{2}_{\delta \kappa }(k)$| corresponding to each combination of |$\sigma ^2_{\delta \psi }$| and β for three different spatial scales: 0.5, 1 and 3 kpc. The smallest scale of 0.5 kpc corresponds to ∼3 pixels, which is the smallest spatial scale considered in our analysis. The largest considered scale, on the other hand, is limited by the size of the lensed images (gravitational arcs) in SDSS J0252+0039. Fig. 10 shows the resulting contour lines connecting matter-power-spectrum models with the same value of |$\Delta ^{2}_{\delta \kappa }$| inferred on a particular scale, overlaid on the original exclusion plot. Based on Fig. 10, we rule out matter-power-spectrum models with |$\Delta ^{2}_{\delta \kappa }$| larger than 1 on 0.5-kpc scale, larger than 0.1 on 1-kpc scale and larger than 0.01 on 3-kpc scale, at the 99 per cent confidence level.
These results can be additionally interpreted in terms of the standard deviation of the total convergence perturbation |$\sigma _{\delta \kappa }(\lambda) \equiv \sqrt{\Delta ^{2}_{\delta \kappa }(\lambda)}$| within the aperture diameter equal to the considered scale λ in an infinitely large sample of circular regions, randomly chosen in proximity to the Einstein radius. For the lensing-mass distribution of SDSS J0252+0039, we infer the following upper limits on this standard deviation (on the spatial scales between L−1 and pixel scale−1, within the range of the power-law power-spectrum slope 3 ≤ β ≤ 8): σδκ(0.5 kpc) < 1 on 0.5-kpc scale, σδκ(1 kpc) < 0.3 on 1-kpc scale and σδκ(3 kpc) < 0.1 on 3-kpc scale. With the critical surface-mass density for SDSS J0252+0039 Σcr ≈ 4 × 108 M⊙ kpc−2, these constraints can be translated into upper limits on the integrated standard deviation σAM(λ) in the aperture mass (within an aperture of diameter λ in the lens plane) in the inner region of the lens galaxy SDSS J0252+0039: σAM(0.5 kpc) < 0.8 × 108 M⊙, σAM(1 kpc) < 1 × 108 M⊙, and σAM(3 kpc) < 3 × 108 M⊙, at the 99 per cent confidence level.
9 DISCUSSION
In order to compare the derived upper-limit constraints with the predictions from the ΛCDM model, we now provide a simple heuristic estimation of the expected dimensionless convergence power spectrum |$\Delta ^{2}_{\delta \kappa }(\lambda)$| due to CDM subhaloes in the dark-matter halo of SDSS J0252+0039. For the sake of simplicity, in our estimation, we only consider the contribution from CDM subhaloes, neglecting both possible line-of-sight haloes and any other fluctuations in the baryonic or dark-matter distribution of the lens galaxy. Furthermore, we assume the projected substructure mass fraction of 0.005 near the Einstein radius of the lens halo in the subhalo mass range between 4 × 106 and 4 × 109 M⊙ (Vegetti & Koopmans 2009b). Finally, we conservatively treat the subhaloes as point masses and apply the Poisson statistics to analytically predict their abundance from the CDM substructure mass function of the form dN/dM ∝ M−1.9 (Springel et al. 2008), thus neglecting the possible suppression of the substructure population due to baryonic processes (Despali et al. 2018). Taking into consideration that the studied lens galaxy is well described by the Singular Isothermal Sphere model and, thus, its convergence κ ∼ 0.5 in proximity to the Einstein radius, we estimate the following upper limits on |$\Delta ^{2}_{\delta \kappa }(\lambda)$| due to CDM subhaloes: |$\Delta ^{2}_{\rm {CDM}}(0.5 \ \mathrm{kpc})\lt 10^{-3}$|, |$\Delta ^{2}_{\rm {CDM}}(1 \ \mathrm{kpc})\lt 4 \times 10^{-4}$|, and |$\Delta ^{2}_{\rm {CDM}}(3 \ \mathrm{kpc})\lt 10^{-4}$|. While this point-mass approach is solely an approximation, it provides a conservative upper limit on the convergence contribution from CDM subhaloes. Any other choice of the subhalo density profile, for example the more realistic Navarro–Frenk–White (NFW) density profile, would only lower the predicted |$\Delta ^{2}_{\rm {CDM}}(\lambda)$|.
A comparison of this estimation with the inferred observational constraints on |$\Delta ^{2}_{\delta \kappa }$| for a flat convergence power spectrum (corresponding to the power-spectrum slope β = 4 for the potential perturbations; see Fig. 10) leads to the conclusion that the estimated contribution from CDM subhaloes lies significantly below the observational upper limits on all considered scales. This preliminary conclusion does not change even if we take into account that in reality the total number of haloes perturbing the lensed images might be a few times higher than the estimated number of subhaloes, due to the possible presence of unbound haloes along the line of sight (Despali et al. 2018).
We attribute the above discrepancy mainly to the fact that, unlike our heuristic ΛCDM-based predictions, the inferred observational upper-limit constraints refer to the total (dark and baryonic) mass distribution projected along the line of sight, including the complex baryonic and dark matter structure of the lens galaxy (e.g. Gilman et al. 2017; Hsueh et al. 2018) and the line-of-sight haloes. Further research is required to adequately compare these observational limits with hydrodynamical simulations, which model not only the formation of dark-matter haloes and subhaloes, but also include baryonic processes and their effect on the overall mass distribution in galaxies.
10 SUMMARY AND CONCLUSIONS
The alternative dark-matter models and galaxy-formation scenarios predict significantly different levels of mass structure on the subgalactic scales (e.g. Lovell et al. 2014). In this work, we have introduced a novel methodology to observationally constrain the statistical properties of such small-scale density fluctuations in the total projected mass distribution of massive elliptical lens galaxies by means of the power-spectrum analysis of surface-brightness anomalies measured in highly magnified galaxy-scale Einstein rings and gravitational arcs. Our approach is based on the theoretical framework introduced by Chatterjee & Koopmans (2018).
The pilot application of the presented methodology to the lens system SDSS J0252+0039 from the SLACS Survey leads to the following conclusions:
The enhanced intrinsic source-galaxy structure in the analysed U-band data requires a higher source-grid resolution and leads to more severe degeneracies between the source and lens models than it was the case for the I-band data previously modelled by Vegetti et al. (2014). Whereas this degeneracy is less problematic when trying to identify individual subhaloes with masses above the detection limit (see Vegetti et al. 2014), alleviating it is crucial when performing a power-spectrum analysis. In this paper, we have addressed this degeneracy by lowering the resolution of the source reconstruction in the smooth lens modelling to suppress the absorption of the potential perturbations and the observational noise into the source structure. This strategy has been shown to be effective in the performance test of our methodology discussed in Paper I.
Our analysis of SDSS J0252+0039 rules out the presence of Gaussian potential perturbations |$\delta \psi _{\rm {GRF}}({{\boldsymbol x}})$| with the variance |$\sigma ^2_{\delta \psi }$| exceeding ∼10−2.5 at 99 per cent C.L. (on the spatial scales between L−1 = 0.88 arcsec−1 and pixel scale−1 = 16.79 arcsec−1 and within the range of the GRF power-spectrum slope 3 ≤ β ≤ 8).
In order to account for the effect of the chosen field-of-view, we infer the corresponding constraints on the dimensionless convergence power spectrum |$\Delta ^{2}_{\delta \kappa }(\lambda)$| on three different subgalactic scales and rule out matter-power-spectrum models with |$\Delta ^{2}_{\delta \kappa }(0.5 \ \mathrm{kpc})\gt 1$| on 0.5-kpc scale, |$\Delta ^{2}_{\delta \kappa }(1 \ \mathrm{kpc})\gt 0.1$| on 1-kpc scale, and |$\Delta ^{2}_{\delta \kappa }(3 \ \mathrm{kpc})\gt 0.01$| on 3-kpc scale (at the 99 per cent C.L.).
The inferred constraints can be translated into to the following limits on the standard deviation σAM(λ) of the aperture mass (integrated within a cylinder with the respective diameter λ in the lens plane) in proximity to the Einstein radius of the lens galaxy: σAM (0.5 kpc) < 0.8 × 108 M⊙, σAM (1 kpc) < 1 × 108 M⊙ and σAM (3 kpc) < 3 × 108 M⊙ (at the 99 per cent C.L.).
We find that the fundamental quantity that determines the level of surface-brightness anomalies in the lensed images and, thus, the probability of the matter-power-spectrum model exclusion, is the total variance in the differential deflection angle |$\sigma ^2_{\delta \alpha }$| (on the spatial scales between L−1 and pixel scale−1) resulting from the underlying potential perturbations |$\delta \psi _{\rm {GRF}}({{\boldsymbol x}})$|. Consequently, the exclusion probability is nearly insensitive to the slope of the deflection-angle power spectrum Pδα(k), i.e. the distribution of |$\sigma ^2_{\delta \alpha }$| over the different length scales. Based on our analysis, |$\sigma ^2_{\delta \alpha }\lt 6\times 10^{-3}$| within the entire considered range 3 ≤ β ≤ 8. This insight is valuable for our future analysis of further galaxy-scale lens systems, which might be carried out by perturbing the deflection angle |$\alpha ({{\boldsymbol x}})$| instead of the lensing potential |$\psi ({{\boldsymbol x}})$|. The threshold value itself, however, might vary for different lens systems and depend on the chosen field-of-view, the PSF, and the signal-to-noise ratio of the analysed image.
In future work, we intend to investigate the modelling degeneracies in more detail, analyse a larger sample of lens systems, and infer more stringent constraints on the dark-matter and galaxy-formation models by comparing these results to hydrodynamical simulations.
ACKNOWLEDGEMENTS
The authors would like to thank the anonymous reviewer for his/her constructive and valuable comments on this work. We also thank Georgios Vernardos for his suggestions, many of which were very helpful. This study is based on observations made with the NASA/ESA Hubble Space Telescope, obtained from the data archive at the Space Telescope Science Institute. Support for this work was provided by a VICI grant (project number 614.001.206) from the Netherlands Organization for Scientific Research (NWO) and by a NASA grant (HST-GO-12898) from the Space Telescope Science Institute. STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5–26555. DB acknowledges support by the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. TT acknowledges support by the Packard Foundation through a Packard Research Fellowship. CDF acknowledges support from the NSF under grant AST-1715611.
DATA AVAILABILITY
The images and mock data analysed in this work are available from the corresponding author upon reasonable request. The raw HST images are publicly available in the Mikulski Archive for Space Telescopes (MAST).
Footnotes
References
APPENDIX A: DUST ANALYSIS
Among all HST images available for SDSS J0252+0039 (U, V, I, and H bands), the U-band photometry is most sensitive to the presence of dust. Thus, clumpy dust in the lens galaxy could potentially cause the small-scale variations in the surface brightness of the lensed images that have been so far interpreted as arising from density fluctuations in the lensing-mass distribution. In order to investigate this possible degeneracy, we compare coaligned images of SDSS J0252+0039 in the U and I band. For this comparison, the I-band image is drizzled to the reference frame and pixel size of the U-band image. The U-band image, on the other hand, is smoothed with a Gaussian profile to a larger blur (according to the Gaussian cascade smoothing) of the I-band image: the PSF of the F390W-filter is well described by a Gaussian profile with the FWHM 0.07 arcsec, whereas the PSF of the F814W filter has the FWHM 0.1 arcsec. Subsequently, both images are divided by the respective standard deviation of the photo-electron counts in the empty sky (|$0.0012 \, \mathrm{e^-} \, \mathrm{sec^{-1}}$| for the smoothed U band and |$0.0045 \, \mathrm{e^-} \, \mathrm{sec^{-1}}$| for the I-band image drizzled to the pixel scale of the U band) to obtain the signal-to-noise ratio for each pixel. Finally, to assess the effect of dust, we generate a ratio image between the I and the U band and conclude that the variations present across the lensed images can be attributed to differences in the corresponding point-spread functions of the compared filters. An additional visual inspection of the colour image in Fig. 2, based on observations in the UV (F390W), the visual (F814W), and the infrared (F160W) bands, shows no indication for dust extinction either in the lens or the source galaxy.