Probing sub-galactic mass structure with the power spectrum of surface-brightness anomalies in high-resolution observations of galaxy-galaxy strong gravitational lenses – I. Power-spectrum measurement and feasibility study | Monthly Notices of the Royal Astronomical Society

ABSTRACT

While the direct detection of the dark-matter particle remains very challenging, the nature of dark matter could be possibly constrained by comparing the observed abundance and properties of small-scale sub-galactic mass structures with predictions from the phenomenological dark-matter models, such as cold, warm, or hot dark matter. Galaxy-galaxy strong gravitational lensing provides a unique opportunity to search for tiny surface-brightness anomalies in the extended lensed images (i.e. Einstein rings or gravitational arcs), induced by possible small-scale mass structures in the foreground lens galaxy. In this paper, the first in a series, we introduce and test a methodology to measure the power spectrum of such surface-brightness anomalies from high-resolution Hubble Space Telescope (HST) imaging. In particular, we focus on the observational aspects of this statistical approach, such as the most suitable observational strategy and sample selection, the choice of modelling techniques, and the noise correction. We test the feasibility of the power-spectrum measurement by applying it to a sample of galaxy-galaxy strong gravitational lens systems from the Sloan Lens ACS Survey, with the most extended, bright, high-signal-to-noise-ratio lensed images, observed in the rest-frame ultraviolet. In the companion paper, we present the methodology to relate the measured power spectrum to the statistical properties of the underlying small-scale mass structures in the lens galaxy and infer the first observational constraints on the sub-galactic matter power spectrum in a massive elliptical (lens) galaxy.

gravitational lensing: strong, methods: statistical, galaxies: structure, dark matter, cosmology: observations

1 INTRODUCTION

Over the last 40 yr, studies of the spatial mass distribution inside galaxies have provided valuable insights into the complex processes of galaxy formation and evolution. Most importantly, the internal mass density profiles of spiral galaxies inferred from their kinematics (Bosma 1978; Rubin, Ford & Thonnard 1978) have led to the hypothesis of a hitherto unknown dominant non-baryonic matter component, referred to as dark matter, which nowadays constitutes a crucial pillar of the concordance dark energy plus cold dark matter (ΛCDM) cosmological model. According to this model and the associated hierarchical structure formation scenario, the early gravitational collapse of dark matter into haloes has created the potential wells necessary for the baryonic gas to cool and condense, finally leading to the formation of the observable galaxies (Blumenthal et al. 1984; White & Frenk 1991; Gao et al. 2007).

Despite this essential role of dark matter in the cosmological structure formation process, its nature and properties remain unknown. The standard CDM paradigm is still challenged by various alternative models, such as warm dark matter (WDM; see e.g. Bode, Ostriker & Turok 2001; Lovell et al. 2014) or self-interacting dark matter (see e.g. Spergel & Steinhardt 2000; Tulin & Yu 2018). These have been proposed in an attempt to explain the striking discrepancy between the number of dwarf satellite galaxies observed in the Local Group and the corresponding predictions from ΛCDM-based simulations (i.e. the Missing Satellites Problem, see e.g. Klypin et al. 1999; Moore et al. 1999; Diemand, Kuhlen & Madau 2007; McConnachie 2012; Drlica-Wagner et al. 2015; Nierenberg et al. 2016; Dooley et al. 2017).

In general, the abundance of such sub-galactic mass structures is determined by the free-streaming length of dark matter in the early Universe. This, in turn, depends on the microscopic properties of the dark-matter particles, such as the particle mass or the strength of the particle–particle interactions. Whereas the standard CDM model predicts galaxy-size haloes to be inhabited by an abundant population of mass structures, the alternative models with less massive, or self-interacting dark-matter particles significantly suppress the formation of sub-galactic mass structures, especially in the low-mass regime below ∼10⁸M_⊙ (see e.g. Bullock & Boylan-Kolchin 2017). Hence, while the direct detection of the dark-matter particle remains challenging (see e.g. Bertone & Tait 2018), its properties could be constrained based on the observed abundance and properties of low-mass structures in a representative sample of galaxies.

However, detecting such low-mass sub-galactic structures beyond the Local Group is a demanding undertaking. They are generally thought to be dark-matter dominated or even purely dark (i.e. devoid of any stars) and, thus, intrinsically invisible. Even if massive enough to form stars, they might be too faint to be observed directly at cosmological distances with the currently available instruments. For this reason, constraints on the properties of sub-galactic mass structures at cosmological distances have been so far inferred mainly from high-resolution observations and modelling of galaxy-scale strong gravitational lenses (e.g. Mao & Schneider 1998; Metcalf & Madau 2001; Dalal & Kochanek 2002; Vegetti, Czoske & Koopmans 2010a; Vegetti et al. 2010b; Nierenberg et al. 2014; Birrer, Amara & Refregier 2017; Gilman et al. 2018; Ritondale et al. 2019; Gilman et al. 2020; Hsueh et al. 2020).

In particular, the phenomenon of galaxy-galaxy strong gravitational lensing makes it possible to detect mass structures in galaxies that, fortuitously, happen to lie along the same line of sight and act as a strong gravitational lens on another galaxy located at a larger distance. Mass structures in the lens galaxy (and possible line-of-sight haloes) induce perturbations to the otherwise smooth lensing potential. These, in turn, perturb the deflection angles of light rays crossing the lens plane in proximity to the mass structures. Thus, even if the structures were purely dark, their gravitational signatures might be observable in the form of the resulting anomalies in the surface-brightness distribution of the extended lensed images (i.e. Einstein ring or gravitational arcs), measured with respect to the best-fitting smooth-lens model (Blandford, Surpi & Kundić 2001; Koopmans 2005; Rau, Vegetti & White 2013). One of the most successful methods utilizing this effect to search for individual mass structures in (massive elliptical) lens galaxies is the gravitational-imaging technique (Koopmans 2005; Vegetti & Koopmans 2009). The application of this technique to deep Hubble Space Telescope (HST) imaging has so far resulted in a detection of two dark-matter subhaloes with the mass of 3.5 × 10⁹M_⊙ and 1.9 × 10⁸M_⊙ at the redshift z = 0.2 and z = 0.881, respectively, with the latter one being the smallest and most distant galactic substructure discovered up to now beyond the local Universe (Vegetti et al. 2010c, 2012).

In order to investigate less massive sub-galactic mass structures, numerously predicted in ΛCDM-based cosmological simulations, Bus (2012), Hezaveh et al. (2016), Diaz Rivero, Cyr-Racine & Dvorkin (2018), and Chatterjee & Koopmans (2018) proposed a complementary statistical approach. Instead of individual massive subhaloes in the lens galaxy, represented by localized potential corrections, the statistical approach models the entire population of small-scale sub-galactic mass structures as Gaussian-random-field (GRF) potential perturbations superposed on the best-fitting smoothly varying lensing potential. In the framework of the theoretical formalism proposed by Chatterjee & Koopmans (2018), both the potential perturbations and the collectively induced surface-brightness anomalies in the lensed images are quantified in terms of their power spectra and related to each other. Successful tests of this approach on mock lensed images, presented by Chatterjee & Koopmans (2018) and Chatterjee (2019), suggest that it might be possible to infer observational constraints on the power spectrum of small-scale mass structures in a (massive elliptical) lens galaxy from the power spectrum of the resulting surface-brightness anomalies in the extended lensed images of the background source galaxy.

This paper is the first in a series of papers aimed at investigating the potential and feasibility of applying this power-spectrum approach to real observational data. The goal of the present paper is to introduce and test the methodology to reliably extract the power spectrum of surface-brightness anomalies from high-resolution HST-imaging of galaxy-galaxy strong gravitational lenses. In particular, we focus on the observational aspects of the power-spectrum measurement, such as the most suitable observational strategy, the sample selection, the choice of modelling techniques, and the noise correction. In the companion paper, Bayer et al. (2023), hereafter referred to as Paper II, we extend the methodology in order to relate the measured power spectrum of surface-brightness anomalies in the lensed images to the statistical properties of the underlying small-scale mass structure in the lens galaxy and infer the first observational constraints on the matter power spectrum in a massive elliptical (lens) galaxy. Future research will apply this approach to a larger sample of lens systems and compare the results with predictions from hydrodynamical simulations, which might eventually allow us to distinguish between the alternative dark-matter models.

The paper is structured as follows. In Section 2, we formalize the concept of surface-brightness anomalies in extended lensed images of a galaxy-galaxy strong gravitational lens system. Section 3 moves on to describe our observational strategy, sample selection, and the imaging data. In Section 4, we present our methodology to measure the power spectrum of surface-brightness anomalies caused by small-scale mass structures in the lens galaxy. Section 5 demonstrates the feasibility of our approach in recovering mock surface-brightness anomalies from simulated lensed images mimicking real observations. Finally, Section 6 provides conclusions and implications for further work.

For a consistent comparison of the inferred smooth lens models with earlier studies by Vegetti et al. (2014), throughout this paper, we assume the following cosmology: H₀ = 73 km s⁻¹Mpc⁻¹, Ω_M = 0.25, and |$\Omega _\Lambda = 0.75$|⁠. Given this cosmology, 1 arcsec corresponds to ∼4 kpc at the redshift of the studied lens galaxies (z_L ∼ 0.2–0.4).

2 SURFACE-BRIGHTNESS ANOMALIES IN EXTENDED LENSED IMAGES

In this section, we first discuss the concept of the hypothetical surface-brightness anomalies that would emerge in the extended lensed images of a background source galaxy as a result of small-scale density fluctuations in the foreground lens galaxy. Subsequently, we elaborate on the observational and modelling challenges that need to be circumvented in order to accurately measure such anomalies.

Let us consider a galaxy-galaxy strong gravitational lens system with extended lensed images described by the surface-brightness distribution |$I({{\boldsymbol x}})$| as a function of the position |${{\boldsymbol x}}$| in the lens plane. The spatial configuration of the lens system is parametrized by the angular diameter distances from the observer to the foreground lens galaxy D_d, from the observer to the background source galaxy D_s and from the lens to the source galaxy D_ds. Following the convention of strong gravitational lensing, we express the surface mass density |$\Sigma ({{\boldsymbol x}})$| of the lens galaxy (including the possible line-of-sight haloes) in units of the critical surface mass density:

$$\begin{eqnarray} \Sigma _{\rm {cr}} = \frac{c^{2}}{4 \pi G} \frac{D_{s}}{D_{d} D_{ds}} \end{eqnarray}$$

(1)

to obtain the commonly used (dimensionless) convergence:

$$\begin{eqnarray} \kappa ({{\boldsymbol x}})=\Sigma ({{\boldsymbol x}})/\Sigma _{\rm {cr}}. \end{eqnarray}$$

(2)

Furthermore, we define the lensing potential, i.e. the gravitational potential of the lens galaxy projected along the line of sight:

$$\begin{eqnarray} \psi ({{\boldsymbol x}})=\frac{1}{\pi }\int _{\mathbb {R}^2}\boldsymbol {{{\boldsymbol dx}}^{\prime }} \kappa ({{\boldsymbol x}}^{\prime })\ \rm {ln}\mid {{\boldsymbol x}}-{{\boldsymbol x}}^{\prime }\mid , \end{eqnarray}$$

(3)

which is related to the convergence by the Poisson’s equation:

$$\begin{eqnarray} \nabla ^2 \psi ({{\boldsymbol x}}) = 2 \kappa ({{\boldsymbol x}}). \end{eqnarray}$$

(4)

The associated (scaled) deflection-angle field:

$$\begin{eqnarray} \boldsymbol {\alpha } ({{\boldsymbol x}})= \nabla \psi ({{\boldsymbol x}}) \end{eqnarray}$$

(5)

determines a mapping between the positions |${{\boldsymbol x}}$| and |${{\boldsymbol y}}$| in the lens- and the source plane, respectively, which is encapsulated in the lens equation:

$$\begin{eqnarray} {{\boldsymbol y}}({{\boldsymbol x}}) = {{\boldsymbol x}}- \boldsymbol {\alpha }({{\boldsymbol x}}). \end{eqnarray}$$

(6)

This mapping together with the principle of surface-brightness conservation in strong gravitational lensing:

$$\begin{eqnarray} I({{\boldsymbol x}}) = S({{\boldsymbol y}}({{\boldsymbol x}})) \end{eqnarray}$$

(7)

builds the foundation for numerical grid-based smooth-lens-modelling codes (e.g. the adaptive grid-based Bayesian lens-modelling code by Vegetti & Koopmans 2009, used in this work) that allow one to simultaneously reconstruct the best-fitting smooth (parametric) lensing potential |$\psi _{M}({{\boldsymbol x}})$| and the unlensed intrinsic surface-brightness distribution of the source galaxy |$S({{\boldsymbol y}})$|⁠.

A discrepancy between the surface-brightness distribution of the observed lensed images |$I({{\boldsymbol x}})$| and the prediction from the best-fitting smooth-lens model |$I_M({{\boldsymbol x}})$| might point towards the presence of mass structure in the lens galaxy. As can be seen from equations (5) and (6), a deviation of the true lensing potential |$\psi ({{\boldsymbol x}})$| from the best-fitting smooth lensing potential |$\psi _{M}({{\boldsymbol x}})$| modifies the mapping between the lens- and the source plane. This, in turn, results in a surface-brightness change |$\delta I({{\boldsymbol x}})$|⁠, such that:

$$\begin{eqnarray} \delta I({{\boldsymbol x}}) = I({{\boldsymbol x}})-I_M({{\boldsymbol x}}) = S\big ({{\boldsymbol x}}-\nabla \psi ({{\boldsymbol x}})\big)-S\big ({{\boldsymbol x}}-\nabla \psi _M({{\boldsymbol x}})\big) \end{eqnarray}$$

(8)

(Blandford et al. 2001; Koopmans 2005). In what follows, we refer to |$\delta I({{\boldsymbol x}})$| as surface-brightness anomalies.

In reality, the extraction of such surface-brightness anomalies from the imaging of real lens systems is complicated by the following issues.

First, observational effects make it impossible to measure the true surface-brightness distribution of the lensed images |$I({{\boldsymbol x}})$|⁠. Instead, the lensed images are blurred by the convolution with the point-spread function (PSF) and the pixellation of the imaging. Moreover, they are affected by the presence of the observational noise and subject to data processing, for example drizzling of the raw data to obtain the final science image (Gonzaga et al. 2012).
Secondly, the reconstructed unlensed surface-brightness distribution of the source galaxy cannot be assumed to perfectly represent the reality. This is due to a degeneracy between the perturbative lensing effect of mass structure in the lens galaxy and the intrinsic surface-brightness fluctuations in the source galaxy. To mitigate this problem, the adaptive grid-based Bayesian smooth-lens-modelling code by Vegetti & Koopmans (2009), used in this work, applies a regularization of the source reconstruction by penalizing solutions with overly strong surface-brightness fluctuations in the source galaxy, as discussed by Warren & Dye (2003) and Koopmans (2005). However, the imposed level of regularization itself is optimized for in the lens-modelling procedure, which might either suppress or enhance the true surface-brightness fluctuations in the source as a result of an over or underregularized source reconstruction, respectively. In other words, the reconstructed smooth-lens model might potentially ‘absorb’ the effect of mass structure into spurious source structure or vice versa – the surface-brightness anomalies due to mass structure might be artificially enhanced if the source reconstruction is overregularized (i.e. too smooth).
Thirdly, the effect of strong gravitational lensing is sensitive to the total mass present in the line-of-sight along which the gravitational lens is observed. Thus, the investigated surface-brightness anomalies might arise not only from mass structure in the lens galaxy but also from possible line-of-sight haloes (Li et al. 2016; Despali et al. 2018). While a detection along the line of sight is valuable for its own sake, it makes the interpretation of the results more complicated.

Hence, one of the main challenges in our approach is to reliably extract the true surface-brightness anomalies |$\delta I({{\boldsymbol x}})$| [as defined in equation (8)] from the smooth-lens-model residuals, taking into consideration all the effects discussed above.

3 OBSERVATIONAL STRATEGY AND DATA

We perform this pilot study based on our HST/WFC3/F390W-observations (Program 12898, Koopmans 2012) of 10 lens systems with highly structured star-forming lensed galaxies selected from the SLACS Survey (Bolton et al. 2008). In this section, we first motivate the choice of the ultraviolet band and discuss our selection criteria for the SLACS sub-sample. Subsequently, we elaborate on the undersampling problem, the dithering strategy, and the data reduction.

3.1 Selection of the observational filter

The level of surface-brightness anomalies caused by the presence of mass structures in the foreground lens galaxy depends not only on the substructure mass but also on the gradient (i.e. level of variations) in the intrinsic surface-brightness distribution of the lensed source galaxy itself. More specifically, if the population of mass structures in the lens galaxy is represented by a potential-perturbation field |$\delta \psi ({{\boldsymbol x}})$| and the intrinsic surface-brightness distribution of the source galaxy is described by |$S({{\boldsymbol y}})$|⁠, then the level of the resulting surface-brightness anomalies |$\delta I({{\boldsymbol x}})$| can be computed as the inner product of the respective gradient fields:

$$\begin{eqnarray} \delta I({{\boldsymbol x}}) = -\nabla S({{\boldsymbol y}}({{\boldsymbol x}})) \cdot \nabla \delta \psi ({{\boldsymbol x}}) \end{eqnarray}$$

(9)

(Blandford et al. 2001; Koopmans 2005).

Hence, the surface-brightness anomalies caused by a given population of mass structures in the lens galaxy can be enhanced by a high level of variations in the intrinsic surface brightness of the source galaxy.

This motivates our choice of the ultraviolet band and the selection of lens systems with highly structured star-forming lensed galaxies. While the optical and infrared bands used in earlier gravitational-imaging studies of SLACS lenses (e.g. Vegetti et al. 2014) capture mostly the smooth old stellar populations, the compact star-forming regions prominent in the ultraviolet band are expected to enhance the surface-brightness gradients in the lensed galaxies, allowing us to improve the sensitivity of our approach to low-mass structures in the lens galaxies (see also Ritondale et al. 2019). We choose to carry out our observations using the Wide Field Camera 3 (WFC3) onboard the HST, which offers the highest resolution and signal-to-noise ratio currently available in the rest-frame ultraviolet (HST/WFC3/F390W, Program 12898; Koopmans 2012).

3.2 Sample selection

For our HST/WFC3/F390W-observations, we select a sub-sample of 10 lens systems from the SLACS Survey (Bolton et al. 2008), which is the so-far largest homogeneous sample of galaxy-galaxy strong gravitational lenses, comprising more than a hundred lens systems. Most of them consist of a massive early-type lens galaxy aligned with a blue star-forming source galaxy, in compliance with our observational strategy. The already existing multicolour HST-imaging data, as well as the extensive lens-modelling and kinematic studies (e.g. Auger et al. 2009; Vegetti et al. 2014) and the availability of the spectroscopic redshifts for both the lens and the source galaxies, make the SLACS lenses an excellent choice to test our methodology.

We intentionally exclude all SLACS systems with late-type lens galaxies to avoid possible degeneracies between the star formation or dust extinction in the lens galaxy and in the lensed images. For all remaining lens systems, we obtain the SDSS spectra and determine the [OII]-flux of the source galaxies, which is assumed to be a good proxy for their star formation rate. Furthermore, we retrieve the already existing one-orbit HST/ACS/F814W-observations and calculate the average surface brightness within the most compact area of the lensed images containing half of the total lensed flux. This can be regarded as the effective surface brightness of the lensed images. A good correlation between the measured [OII]-flux and the surface brightness of the lensed images in F814W indicates that the latter is related to the young stellar populations in the source galaxy (but is more diffuse). Thus, by selecting SLACS lens systems with the brightest lensed images in the F814W-filter, we at the same time select lens systems with highly star-forming source galaxies, in accordance with our observational strategy.

Our final target list comprising 10 SLACS lens systems with the brightest and most extended lensed images is presented in Table 1 together with the basic astrometric and spectroscopic properties: the location on the sky, the spectroscopic redshifts of the lens and the source galaxies, and the stellar velocity dispersion in the lens galaxy (Auger et al. 2009). To improve the quality of the new F390W-imaging data in comparison to the already existing observations in the F814W-filter, we require an average signal-to-noise ratio of 10 in the pixels covering the lensed images and calculate the observing time accordingly. We quote the number of HST-orbits devoted to the observations of each lens system in the last column of Table 1.

Table 1.

Open in new tab

Astrometric and spectroscopic properties of the selected SLACS sub-sample observed with the Hubble Space Telescope in the rest-frame ultraviolet (HST/WFC3/F390W): the target names, location on the sky, spectroscopic redshifts of the lens z_l and the source galaxy z_s, the velocity dispersion of the lens galaxy σ_l (from Auger et al. 2009), and the number of HST-orbits devoted to the observations of each lens system. For the purpose of this study, we perform the lens modelling for only four most suitable lens systems from this sample – SDSS J0252+0039, SDSS J0737+3216, SDSS J1430 + 4105, and SDSS J1627–0053 – characterized by a relatively simple geometry and showing no substantial galaxy-core residuals in the galaxy-subtracted images.

Target lens galaxy	RA	Decl.	z_l	z_s	\|$\sigma _{l} \rm [km/s]$\|	Orbits
SDSS J0252 + 0039	02 52 45.21	+00 39 58.40	0.280	0.982	164 ± 12	2
SDSS J0737 + 3216	07 37 28.45	+32 16 18.60	0.322	0.581	338 ± 16	2
SDSS J0903 + 4116	09 03 15.19	+41 16 09.10	0.430	1.065	223 ± 27	6
SDSS J0912 + 0029	09 12 05.31	+00 29 01.20	0.164	0.324	326 ± 12	5
SDSS J0956 + 5100	09 56 29.78	+51 00 06.60	0.241	0.470	334 ± 15	2
SDSS J0959 + 0410	09 59 44.07	+04 10 17.00	0.126	0.535	197 ± 13	1
SDSS J1430 + 4105	14 30 04.10	+41 05 57.10	0.285	0.575	322 ± 32	2
SDSS J1627–0053	16 27 46.45	–00 53 57.60	0.208	0.524	290 ± 14	3
SDSS J1630 + 4520	16 30 28.16	+45 20 36.30	0.248	0.793	276 ± 16	5
SDSS J2341 + 0000	23 41 11.57	+00 00 18.70	0.186	0.807	207 ± 13	2

Target lens galaxy	RA	Decl.	z_l	z_s	\|$\sigma _{l} \rm [km/s]$\|	Orbits
SDSS J0252 + 0039	02 52 45.21	+00 39 58.40	0.280	0.982	164 ± 12	2
SDSS J0737 + 3216	07 37 28.45	+32 16 18.60	0.322	0.581	338 ± 16	2
SDSS J0903 + 4116	09 03 15.19	+41 16 09.10	0.430	1.065	223 ± 27	6
SDSS J0912 + 0029	09 12 05.31	+00 29 01.20	0.164	0.324	326 ± 12	5
SDSS J0956 + 5100	09 56 29.78	+51 00 06.60	0.241	0.470	334 ± 15	2
SDSS J0959 + 0410	09 59 44.07	+04 10 17.00	0.126	0.535	197 ± 13	1
SDSS J1430 + 4105	14 30 04.10	+41 05 57.10	0.285	0.575	322 ± 32	2
SDSS J1627–0053	16 27 46.45	–00 53 57.60	0.208	0.524	290 ± 14	3
SDSS J1630 + 4520	16 30 28.16	+45 20 36.30	0.248	0.793	276 ± 16	5
SDSS J2341 + 0000	23 41 11.57	+00 00 18.70	0.186	0.807	207 ± 13	2

Table 1.

Open in new tab

Target lens galaxy	RA	Decl.	z_l	z_s	\|$\sigma _{l} \rm [km/s]$\|	Orbits
SDSS J0252 + 0039	02 52 45.21	+00 39 58.40	0.280	0.982	164 ± 12	2
SDSS J0737 + 3216	07 37 28.45	+32 16 18.60	0.322	0.581	338 ± 16	2
SDSS J0903 + 4116	09 03 15.19	+41 16 09.10	0.430	1.065	223 ± 27	6
SDSS J0912 + 0029	09 12 05.31	+00 29 01.20	0.164	0.324	326 ± 12	5
SDSS J0956 + 5100	09 56 29.78	+51 00 06.60	0.241	0.470	334 ± 15	2
SDSS J0959 + 0410	09 59 44.07	+04 10 17.00	0.126	0.535	197 ± 13	1
SDSS J1430 + 4105	14 30 04.10	+41 05 57.10	0.285	0.575	322 ± 32	2
SDSS J1627–0053	16 27 46.45	–00 53 57.60	0.208	0.524	290 ± 14	3
SDSS J1630 + 4520	16 30 28.16	+45 20 36.30	0.248	0.793	276 ± 16	5
SDSS J2341 + 0000	23 41 11.57	+00 00 18.70	0.186	0.807	207 ± 13	2

Target lens galaxy	RA	Decl.	z_l	z_s	\|$\sigma _{l} \rm [km/s]$\|	Orbits
SDSS J0252 + 0039	02 52 45.21	+00 39 58.40	0.280	0.982	164 ± 12	2
SDSS J0737 + 3216	07 37 28.45	+32 16 18.60	0.322	0.581	338 ± 16	2
SDSS J0903 + 4116	09 03 15.19	+41 16 09.10	0.430	1.065	223 ± 27	6
SDSS J0912 + 0029	09 12 05.31	+00 29 01.20	0.164	0.324	326 ± 12	5
SDSS J0956 + 5100	09 56 29.78	+51 00 06.60	0.241	0.470	334 ± 15	2
SDSS J0959 + 0410	09 59 44.07	+04 10 17.00	0.126	0.535	197 ± 13	1
SDSS J1430 + 4105	14 30 04.10	+41 05 57.10	0.285	0.575	322 ± 32	2
SDSS J1627–0053	16 27 46.45	–00 53 57.60	0.208	0.524	290 ± 14	3
SDSS J1630 + 4520	16 30 28.16	+45 20 36.30	0.248	0.793	276 ± 16	5
SDSS J2341 + 0000	23 41 11.57	+00 00 18.70	0.186	0.807	207 ± 13	2

3.3 Observations and data reduction

Our HST/WFC3/F390W-observations were performed between 2013 January 26 and September 16 (Program 12898, Koopmans 2012). Raw HST-images are generally known to be undersampled, which means that the (blurred) point sources are not covered by enough pixels to precisely sample the PSF. Whereas a well-sampled image would have at least two pixels across the full width at half-maximum (FWHM) of the PSF (Nyquist limit), the pixel width of the WFC3 (0.04 arcsec on a side) is relatively large in comparison to the FWHM of the WFC3/UVIS optical performance at 390 nm (0.07 arcsec). In order to alleviate this undersampling problem and more optimally benefit from the superb resolution of the HST-optics, we apply the standard dithering strategy, i.e. we obtain multiple dithered exposures of each target object by slightly shifting the telescope pointing every time the next exposure is taken. Besides the improvement in the PSF-sampling, this dithering technique makes it possible to compensate for cosmic rays, possible dead pixels or columns, and the flat-field effects.

We retrieve the dithered exposures from the Mikulski Archive for Space Telescopes archive¹ in the form of pipeline-pre-processed flat-field calibrated FITS files (flt.fits files) and make use of the Variable-Pixel Linear Reconstruction algorithm (Fruchter & Hook 2002), informally known as Drizzle, to combine them into the final science images of our sample. We perform the drizzling in an automatic way by means of the astrodrizzle task from the drizzlepac package (Gonzaga et al. 2012) in the default configuration. In Section 4.5.2, we investigate the impact of different drizzling settings on our final results in comparison to this reference configuration.

Next, for each lens system, we generate an image cutout centred on the brightest pixel of the lens galaxy, with the side length approximately equal to four Einstein radii, as presented in Fig. 1. Additionally, for the purpose of illustration, Fig. 2 depicts the obtained F390W-imaging for the lens system SDSS J1430 + 4105 in comparison to the archival multiband observations in F606W, F814W, and F160W. The colour-composite image, created using the stiff ² software, combines the F390W-imaging with the archival F814W and F160W-observations.

Figure 1.

HST/WFC3-imaging of the selected SLACS sub-sample in the rest-frame ultraviolet (F390W). Each of these images depicts a massive elliptical galaxy acting as a strong gravitational lens on a star-forming background source galaxy. The source galaxies appear lensed into high-signal-to-noise-ratio gravitational arcs or (in some cases) an almost complete Einstein ring. The image cutouts are centred on the brightest pixel of the lens galaxy and have the side length approximately equal to four Einstein radii of the respective lens system.

Open in new tab Download slide

Figure 2.

A high-resolution multiband view of the Einstein ring in the strong gravitational lens system SDSS J1430 + 4105. From left to right: our new HST/WFC3-observations in F390W, the archival HST-imaging in F606W, F814W, and F160W, and the colour-composite image combining the multiband photometry in F390W, F606W, and F814W obtained using the stiff software. All images are oriented with north pointing up and east pointing left.

Open in new tab Download slide

Finally, we use the tinytim ³ (Krist, Hook & Stoehr 2010) software to model the PSF of the HST/WFC3/F390W-optics in the central pixel of each image cutout. For simplicity, we assume that the PSF in all the other pixels of the considered image cutouts does not deviate significantly from the central pixel. In our PSF-models, we account for the different spectral types of the lens galaxies, which we estimate based on the magnitudes in F555W, F614W, F814W, and F160W, as inferred by Auger et al. (2009).

4 METHODOLOGY

In this section, we present our methodology to measure the power spectrum of surface-brightness anomalies in high-resolution HST-observations of galaxy-galaxy strong gravitational lens systems.

4.1 Analysis synopsis

Our procedure consists of the steps outlined below.

Modelling and subtraction of the lens-galaxy light by means of galfit (Peng et al. 2002) or, alternatively, the b-spline algorithm (Bolton et al. 2006), see Section 4.2;
Smooth lens modelling using the adaptive and grid-based Bayesian lens-modelling technique by Vegetti & Koopmans (2009), see Section 4.3;
Statistical quantification of the residual surface-brightness fluctuations in the lensed images in terms of the azimuthally averaged power spectrum, see Section 4.4;
Estimation of the noise power spectrum based on blank-sky fields (modified to account for the additional Poisson noise present in the imaging of a lens system), see Section 4.5;
Noise-bias correction to reveal the power spectrum of surface-brightness anomalies, see Section 4.6.

In the following sections, we discuss the individual steps in more detail and illustrate them with examples from the analysis of our SLACS sub-sample.

4.2 Lens-galaxy subtraction

Before proceeding with the lens modelling, here we discuss how to reliably estimate and correct for the flux contribution from the lens galaxy in the pixels overlapping with the lensed images. In a typical SLACS lens system, the Einstein radius (∼1 arcsec) is comparable to half the effective radius of the (massive elliptical) lens galaxy (Auger et al. 2009). This means that the lensed images in our sample are projected on the very inner region (∼5–10 kpc from the centre) of the respective lens galaxy and are, thus, contaminated with its light. Since massive elliptical galaxies are empirically known to be characterized by a smoothly varying light distribution and a very regular isophotal structure, a common practice to deal with this overlap is to fit the surface brightness of the lens galaxy with a parametric model and subtract it from the observed image. To this end, we apply and compare the performance of two techniques that have been successfully used in earlier studies of SLACS lenses – the radial b-spline algorithm (Bolton et al. 2006) and the empirical galaxy-fitting code galfit (Peng et al. 2002).

We perform the b-spline modelling using the implementation by Bolton et al. (2006). This technique allows one to find the best-fitting coefficients b_mk and c_mk of the surface-brightness distribution I(R, ϕ) parametrized as follows:

$$\begin{eqnarray} I(\rm {R},\theta) = \sum _{m,k} \big (b_{mk} \rm {cos}(m\phi) + c_{mk} \rm {sin} (m\phi)\big) f_{k}(\rm {R}), \end{eqnarray}$$

(10)

where the radial dependence is modelled with a piecewise (linear, quadratic, or cubic) polynomial function f_k(R) defined on a chosen set of radial intervals k and the angular dependence ϕ is fitted with a chosen number of multipole orders m. We set the radial interval breakpoints every 0.2 arcsec and fit f_k(R) with piecewise-defined cubic polynomials. To model the angular dependence, we begin with the default configuration including the m = 0 (monopole), m = 1 (dipole), and m = 2 (quadrupole) modes. If substantial angular structure is still found in the residual (data-model) image, we iteratively add further multiple orders, i.e. m = 4 (octopole) and higher even terms, until the reduced χ²-statistic is minimized.

As an alternative, we apply the empirical galaxy-fitting technique galfit (Peng et al. 2002) and fit the lens light with a Sérsic profile (Sérsic 1963), which is empirically known to provide a good fit to the surface-brightness distribution of observed massive elliptical galaxies. We iteratively add more Sérsic components if justified by the residual pattern or a high value of the reduced χ²-statistic. In some cases, we improve the model by additionally fitting the disciness/boxiness of the isophotes.

In both approaches, we exclude from the fit all pixels that, for various reasons, should not be taken into account during the galaxy-fitting procedure. More specifically, we mask out all pixels overlapping with the lensed images, nearby satellite galaxies, stellar streams and other astronomical objects not associated with but close in projection to the lens galaxy. In order to generate the mask, we first use the ds9 software⁴ to manually outline these features with a polygon. Subsequently, we determine the set of pixels contained inside this polygon by means of the ds9poly and the fillpoly software (freely available on the galfit webpage⁵). Finally, we make use of the Pyraf-task badpiximage to create a fits-image representing the mask. This initial mask is in some cases adjusted in the course of the galaxy-fitting procedure in order to exclude additional faint features revealed in the residual image.

In each case, the best-fitting model of the surface-brightness distribution in the lens galaxy is finally interpolated over the masked regions and subtracted from the original image. As an example, in Fig. 3, we present the entire procedure of the lens-galaxy subtraction for the lens systems SDSS J0737 + 3216. The figure shows the original HST-image, the applied mask, the best-fitting galfit and b-spline models, and the respective galaxy-subtracted images.

Figure 3.

Lens-galaxy subtraction for SDSS J0737 + 3216 performed using alternatively galfit or the b-spline algorithm. Top row: the observed HST-image in the F390W-band (left-hand panel) and the mask applied to exclude the lensed images from the fit (right-hand panel). Middle 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). Bottom row: the best-fitting b-spline model (left-hand panel) and the corresponding residual image (right-hand panel). The measured flux is expressed in units of electrons per second.

Open in new tab Download slide

As a result, we find that for most of the lens systems in our sample b-spline provide a very good fit to the central region of the lens galaxy; however, the inferred global surface-brightness model has an irregular shape, deviating from our empirical expectations for massive elliptical galaxies. galfit, on the other hand, provides more realistic models, but it very often yields significant galaxy-core residuals, even when fitting multiple Sérsic components. Thus, both techniques need to be applied with caution. The flexibility of b-spline allows to fit features that are not well described by standard parametric functions, but it might also lead to unrealistic models. To the contrary, fitting empirically based parametric functions yields a reasonable solution in most cases but might result in a poor fit if the galaxy deviates from the assumed typical morphology. In Section 4.4.2, we investigate the effect of galfit and b-spline on the resulting power spectrum of the residual image and show that the two lens-galaxy-subtraction techniques lead to almost identical results in this respect. Since our objective is to estimate the surface-brightness contribution of the lens galaxy in the region overlapping with the lensed images and not necessarily to obtain the best estimate of the overall surface-brightness distribution, we choose galfit as the preferred method for the purpose of our analysis.

Fig. 4 shows the final galaxy-subtracted images, obtained using galfit. We note the presence of significant galaxy-core residuals in the majority of the images, which might point towards a cusp or central star formation in the (massive elliptical) lens galaxy (see e.g. Kaviraj et al. 2011). A viable way of mitigating this issue would be to either model the galaxy core separately and only then search for the best-fitting global model or to entirely exclude it from the fit. For the purpose of this study, however, we exclude these problematic lens systems from further analysis and perform the lens modelling only for the remaining four systems – SDSS J0252+0039, SDSS J0737+3216, SDSS J1430 + 4105, and SDSS J1627–0053 – with a relatively simple geometry and no substantial galaxy-core residuals.

Figure 4.

The lens-galaxy-subtracted images for all lens systems in our SLACS sub-sample, obtained using galfit. The significant galaxy-core residuals in some of the images might point towards the presence of a cusp or the central star formation in the (massive elliptical) lens galaxy. For the purpose of this study, we exclude these from further analysis and perform the lens modelling only for the remaining four lens systems having a relatively simple geometry and no substantial galaxy-core residuals.

Open in new tab Download slide

4.3 Smooth lens modelling

We apply the adaptive grid-based Bayesian smooth lens modelling technique by Vegetti & Koopmans (2009) to model each selected lens system under the tentative assumption that the mass in the lens galaxy is distributed smoothly. More specifically, we assume its surface mass density to be well described by the power-law elliptical mass distribution model (PEMD; Barkana 1998), with the convergence parametrized according to the convention used in Vegetti & Koopmans (2009):

$$\begin{eqnarray} \kappa (x,y) = \frac{b \ (2-\frac{\gamma }{2}) \ q^{\gamma -3/2}}{2(x^{2}q^{2}+y^{2})^{(\gamma -1)/2}}. \end{eqnarray}$$

(11)

The model parameters are the lens strength b, the (minor to major) axis ratio q and the (3D) mass-density slope γ (γ = 2 in the isothermal case). Moreover, the mass model is rotated and translated to fit the position angle of the major axis θ (measured with respect to the original telescope rotation) and the centroid location in the lens plane x₀ and y₀. In addition, we model the lensing effect of possible companion objects in the vicinity of the lens galaxy as an external shear field characterized by the shear strength Γ and its position angle Γ_θ.

The lens-modelling code allows us to find the best-fitting parameter values of this PEMD-plus-external-shear macro model and, simultaneously, reconstruct the unlensed pixellated surface-brightness distribution of the source galaxy (on an adaptive grid in the source plane), which, combined together most accurately, reproduce the observed lensed images. However, there are several alternative ways to perform the mapping of pixels and flux values between the lens- and the source plane. First, the resolution of the pixellated source reconstruction can be chosen by setting the value of a parameter referred to as n that determines the linear size of a square in the lens plane out of which only the central 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 reconstruction grid in the source plane, while n = 1 corresponds to casting back every single pixel. Note, however, that all pixels are used in the comparison between the data and the model, independently of the chosen n. Secondly, it is possible to apply different forms of source-grid regularization, such as an adaptive or non-adaptive, variance, gradient, or curvature regularization (see e.g. Suyu et al. 2006). The optimal choice of the source-grid resolution and the form of regularization depends on the level of structure in the source galaxy as well as on the signal-to-noise ratio of the data and is usually made based on the highest value of the marginalized Bayesian evidence.

We obtain smooth-lens models with the highest Bayesian evidence when choosing the highest resolution (n = 1) and an adaptive gradient regularization for all four analysed lens systems (but see Section 4.4.3 for a discussion on the overfitting problem). Table 2 presents the inferred parameter values of the best-fitting PEMD-plus-external-shear macro models, in comparison to the earlier F814W and F555W reconstructions carried out by Vegetti et al. (2014; where available). Figs 5–8 depict the respective modelled data, the inferred best-fitting model, and the reconstructed surface-brightness distribution of the source galaxy on an adaptive grid in the source plane. Moreover, each of these panels shows the resulting residual image representing the deviation of the observed lensed images from the best-fitting model. The key idea of our approach is that, apart from noise, these residual surface-brightness fluctuations might be caused by perturbations in the lensing potential due to small-scale mass structure in the lens galaxy.

Figure 5.

Best-fitting PEMD plus-external-shear macro model of SDSS J0252 + 0039 in the HST/WFC3/F390W-filter, inferred by means of the adaptive and grid-based Bayesian lens-modelling technique of Vegetti & Koopmans (2009). Top row: the lens-galaxy subtracted image overlaid with a mask, used as input for the smooth lens modelling (left-hand panel) and the reconstructed smooth-lens model of the lensed images (right-hand panel). Bottom row: the unlensed surface-brightness distribution of the source galaxy reconstructed on an adaptive grid in the source plane (more specifically a Delaunay triangulation regridded to pixels for plotting purposes) with the source-grid resolution n = 1, i.e. casting back every pixel from the image plane to the source plane (left-hand panel), and the residual (data-model) image showing the remaining surface-brightness fluctuations that are not explained by the smooth-lens model (right-hand panel). The reconstructed parameter values of the best-fitting smooth lensing potential as well as the chosen modelling options can be found in Table 2.

Open in new tab Download slide

Figure 6.

Idem as in Fig. 5 for the lens system SDSS J0737 + 3216.

Open in new tab Download slide

Figure 7.

Idem as in Fig. 5 for the lens system SDSS J1430 + 4105.

Open in new tab Download slide

Figure 8.

Idem as in Fig. 5 for the lens system SDSS J1627–0053.

Open in new tab Download slide

Table 2.

Open in new tab

Parameter values of the best-fitting PEMD plus-external-shear macro model for four lens galaxies from our SLACS sub-sample, selected for further analysis due to a relatively simple geometry and no substantial galaxy-core residuals, in comparison to earlier F814W and F555W reconstructions performed by Vegetti et al. (2014). The models are based on lens-galaxy-subtracted HST/WFC3/F390W-images obtained using alternatively galfit or b-spline. The mass distribution of the lens galaxy is modelled with the following set of free parameters: the lens strength b, the position angle θ (with respect to the original telescope rotation), the axis ratio q, the (3D) mass-density slope γ, the external shear strength Γ, and its position angle Γ_θ. Additionally, we quote the source-grid resolution n and the type of source-plane regularization (Reg.) that lead to models with the highest marginalized Bayesian evidence. G indicates the gradient regularization, C the curvature regularization and adp an adaptive regularization. The typical statistical errors of the parameter values are of the order 10⁻¹ for the angles and 10⁻³ for the remaining parameters.

Lens system	Filter	Lens-galaxy subtraction	\|$b [\rm {arcsec}]$\|	θ (deg.)	q	γ	Γ	Γ_θ (deg.)	n	Reg.
J0252 + 0039	F390W	Galfit	0.996	150.1	0.978	2.066	−0.015	81.4	1	G_adp
	F390W	b-spline	0.996	149.2	0.978	2.066	−0.015	81.4	1	G_adp
	F814W	b-spline	1.022	26.2	0.943	2.047	0.009	101.8	1	G_adp
J0737 + 3216	F390W	Galfit	0.926	66.1	0.862	2.110	0.066	72.8	1	G_adp
	F390W	b-spline	0.933	66.9	0.869	2.102	0.062	72.2	1	G_adp
	F814W	b-spline	0.951	78.3	0.705	2.066	0.050	100.8	1	C
	F555W	b-spline	0.951	77.2	0.709	2.073	0.052	102.5	2	G_adp
J1430 + 4105	F390W	b-spline	1.527	85.4	0.647	2.073	0.029	139.2	1	G_adp
	F814W	b-spline	1.484	61.5	0.710	2.048	0.051	128.6	1	C
J1627–0053	F390W	b-spline	1.217	16.9	0.856	2.006	0.008	101.3	1	G_adp
	F390W	Galfit	1.218	17.0	0.855	2.007	0.008	101.5	1	G_adp
	F814W	b-spline	1.229	14.3	0.912	1.998	0.004	80.0	2	C_adp
	F555W	b-spline	1.212	14.2	0.869	2.058	0.014	87.8	2	C_adp

Lens system	Filter	Lens-galaxy subtraction	\|$b [\rm {arcsec}]$\|	θ (deg.)	q	γ	Γ	Γ_θ (deg.)	n	Reg.
J0252 + 0039	F390W	Galfit	0.996	150.1	0.978	2.066	−0.015	81.4	1	G_adp
	F390W	b-spline	0.996	149.2	0.978	2.066	−0.015	81.4	1	G_adp
	F814W	b-spline	1.022	26.2	0.943	2.047	0.009	101.8	1	G_adp
J0737 + 3216	F390W	Galfit	0.926	66.1	0.862	2.110	0.066	72.8	1	G_adp
	F390W	b-spline	0.933	66.9	0.869	2.102	0.062	72.2	1	G_adp
	F814W	b-spline	0.951	78.3	0.705	2.066	0.050	100.8	1	C
	F555W	b-spline	0.951	77.2	0.709	2.073	0.052	102.5	2	G_adp
J1430 + 4105	F390W	b-spline	1.527	85.4	0.647	2.073	0.029	139.2	1	G_adp
	F814W	b-spline	1.484	61.5	0.710	2.048	0.051	128.6	1	C
J1627–0053	F390W	b-spline	1.217	16.9	0.856	2.006	0.008	101.3	1	G_adp
	F390W	Galfit	1.218	17.0	0.855	2.007	0.008	101.5	1	G_adp
	F814W	b-spline	1.229	14.3	0.912	1.998	0.004	80.0	2	C_adp
	F555W	b-spline	1.212	14.2	0.869	2.058	0.014	87.8	2	C_adp

Table 2.

Open in new tab

Lens system	Filter	Lens-galaxy subtraction	\|$b [\rm {arcsec}]$\|	θ (deg.)	q	γ	Γ	Γ_θ (deg.)	n	Reg.
J0252 + 0039	F390W	Galfit	0.996	150.1	0.978	2.066	−0.015	81.4	1	G_adp
	F390W	b-spline	0.996	149.2	0.978	2.066	−0.015	81.4	1	G_adp
	F814W	b-spline	1.022	26.2	0.943	2.047	0.009	101.8	1	G_adp
J0737 + 3216	F390W	Galfit	0.926	66.1	0.862	2.110	0.066	72.8	1	G_adp
	F390W	b-spline	0.933	66.9	0.869	2.102	0.062	72.2	1	G_adp
	F814W	b-spline	0.951	78.3	0.705	2.066	0.050	100.8	1	C
	F555W	b-spline	0.951	77.2	0.709	2.073	0.052	102.5	2	G_adp
J1430 + 4105	F390W	b-spline	1.527	85.4	0.647	2.073	0.029	139.2	1	G_adp
	F814W	b-spline	1.484	61.5	0.710	2.048	0.051	128.6	1	C
J1627–0053	F390W	b-spline	1.217	16.9	0.856	2.006	0.008	101.3	1	G_adp
	F390W	Galfit	1.218	17.0	0.855	2.007	0.008	101.5	1	G_adp
	F814W	b-spline	1.229	14.3	0.912	1.998	0.004	80.0	2	C_adp
	F555W	b-spline	1.212	14.2	0.869	2.058	0.014	87.8	2	C_adp

Lens system	Filter	Lens-galaxy subtraction	\|$b [\rm {arcsec}]$\|	θ (deg.)	q	γ	Γ	Γ_θ (deg.)	n	Reg.
J0252 + 0039	F390W	Galfit	0.996	150.1	0.978	2.066	−0.015	81.4	1	G_adp
	F390W	b-spline	0.996	149.2	0.978	2.066	−0.015	81.4	1	G_adp
	F814W	b-spline	1.022	26.2	0.943	2.047	0.009	101.8	1	G_adp
J0737 + 3216	F390W	Galfit	0.926	66.1	0.862	2.110	0.066	72.8	1	G_adp
	F390W	b-spline	0.933	66.9	0.869	2.102	0.062	72.2	1	G_adp
	F814W	b-spline	0.951	78.3	0.705	2.066	0.050	100.8	1	C
	F555W	b-spline	0.951	77.2	0.709	2.073	0.052	102.5	2	G_adp
J1430 + 4105	F390W	b-spline	1.527	85.4	0.647	2.073	0.029	139.2	1	G_adp
	F814W	b-spline	1.484	61.5	0.710	2.048	0.051	128.6	1	C
J1627–0053	F390W	b-spline	1.217	16.9	0.856	2.006	0.008	101.3	1	G_adp
	F390W	Galfit	1.218	17.0	0.855	2.007	0.008	101.5	1	G_adp
	F814W	b-spline	1.229	14.3	0.912	1.998	0.004	80.0	2	C_adp
	F555W	b-spline	1.212	14.2	0.869	2.058	0.014	87.8	2	C_adp

4.4 Power-spectrum analysis of the residual surface-brightness fluctuations

The residual images of all investigated lens systems reveal surface-brightness fluctuations that cannot be explained by the assumed smooth PEMD-plus-external-shear macro model of the mass distribution in the lens galaxy, as demonstrated in Figs 5–8. In this section, we estimate the variance of these residuals as a function of their spatial scale, i.e. the power spectrum, following the approach proposed by Bus (2012) and Chatterjee & Koopmans (2018). Moreover, we investigate how the different choices made in the process of the lens-galaxy subtraction and the smooth lens modelling affect the measured residual power spectrum.

4.4.1 Power-spectrum measurement

The goal of this power-spectrum analysis is to decompose the residual surface-brightness fluctuations δI(x) into modes with different length scales λ, expressed in terms of the corresponding wavenumbers k. We note that in this study we follow the convention in which the wavenumber is equal to the reciprocal length scale:

$$\begin{eqnarray} k\equiv \lambda ^{-1} \end{eqnarray}$$

(12)

and is measured in arcsec⁻¹.

We calculate the residual power spectrum for each modelled lens system individually, within the respective mask covering the lensed images. For this, we set the flux values of all pixels located outside the mask to zero and compute the 2D discrete Fourier transform of the masked residual image using the python package numpy.fft. ⁶ The squared magnitudes of the obtained (complex-valued) Fourier coefficients, assigned to the individual pixels of the Fourier-transformed residual image, yield the 2D power spectrum of the residual surface-brightness fluctuations. We further assume the residuals to be isotropic and average this 2D power spectrum along a set of 10 equidistant concentric annuli covering the full Fourier-transformed image. The resulting 1D azimuthally averaged power spectrum P(k) constitutes the final statistic allowing us to perform a statistical comparison of the residual surface-brightness fluctuations resulting from different models.

4.4.2 Effect of the lens-galaxy subtraction

In order to assess whether the choice of the lens-galaxy-subtraction technique, i.e. galfit or b-spline, might be a source of a systematic bias, we calculate the azimuthally averaged power spectrum of the residual images after the lens-galaxy subtraction using both methods for SDSS J0252+0039, SDSS J0737+3216, and SDSS 1627 + 0053. We note that we were not able to obtain a good galfit model for SDSS J1430 + 4105 and, thus, the system is omitted from this comparative power-spectrum analysis. The outcome of this test, presented in Fig. 9, shows that these two galaxy-subtraction techniques lead to almost identical power spectra. As expected, by removing a smooth large-scale surface-brightness component, the lens-galaxy subtraction reduces the measured power spectrum only on the largest spatial scales (i.e. smallest k-values).

Figure 9.

Effect of the lens-galaxy subtraction, performed using alternatively galfit or b-spline, on the power spectrum of surface-brightness fluctuations in the galaxy-subtracted images for SDSS J0252+0039, SDSS J0737+3216, and SDSS 1627 + 0053.

Open in new tab Download slide

We conclude that if the power-spectrum analysis is performed in the ultraviolet band, where the surface-brightness distribution of elliptical galaxies peaks strongly in the centre and decreases very quickly towards the outskirts, the choice of the galaxy-subtraction technique does not appreciably affect the final results. In a different band, however, especially in the infrared where the lens galaxy dominates the surface brightness in the region overlapping with the lensed images, the choice of the lens-light model might significantly alter the analysis outcome. We plan to investigate this issue in our future paper.

4.4.3 Effect of the smooth lens modelling

As stated in Section 4.3, we obtain smooth-lens models with the highest Bayesian evidence when the lens modelling is performed with the highest resolution (n = 1) and an adaptive gradient source-grid regularization. However, it turns out that this choice leads to a power spectrum of the residual surface-brightness fluctuations lying at or below the noise level for all four modelled lens systems (see Section 4.5 for the estimation of the noise power spectrum). This means that all surface-brightness anomalies due to the hypothetical small-scale mass structures in the lens galaxy and even a substantial fraction of the background noise have been modelled as spurious structure in the intrinsic surface-brightness distribution of the source galaxy. This problem is generally known as overfitting (i.e. modelling of the inherent noise present in the data) and arises when the number of free parameters in a model is much larger than the number of the imposed constraints. In our study, it might also arise from an incorrect or incomplete macro model.

In order to investigate this issue, we perform tests with a lower resolution and different forms of the source-grid regularization. As an example, Fig. 10 shows the effect of these different options on the power spectrum of the residual surface-brightness fluctuations measured in the lensed images of the lens system SDSS J0737 + 3216. As is apparent from this figure, the computed residual power spectrum lies at or below the noise level for all models with n = 1 or n = 2, irrespective of the chosen form of regularization. More generally, we find that for almost all investigated lens systems even the choice of n = 2 still leads to the residuals lying below the noise level. From this, we conclude that if the smooth lens modelling in the U-band is carried out with the highest resolution (i.e. n = 1 or n = 2) and a relatively tight mask, as in earlier F814W and F555W reconstructions (see Table 2), the inversion problem to be solved is underconstrained and degenerate.

Figure 10.

Power spectrum of the residual surface-brightness fluctuations remaining in the lensed images of SDSS J0737 + 3216 after the lens-galaxy subtraction and the smooth lens modelling for varying levels of the source-grid resolution and different types of the source-grid regularization. In the legend, G indicates the gradient regularization, C the curvature regularization, ada an adaptive and nada a non-adaptive source-plane regularization. The symbols n1 and n2 stand for the source-grid resolution corresponding to n = 1 and n = 2, respectively. The computed power spectrum lies below the noise level for all models with the highest resolution (n = 1 and n = 2), irrespective of the chosen form of regularization.

Open in new tab Download slide

On the other hand, higher n-values might deteriorate the sampling of the lensed images and, consequently, diminish the accuracy of the source reconstruction. Hence, a balance needs to be found between the possible over and underfitting. In the present study, we mitigate the overfitting problem by lowering the resolution of the source reconstruction even further to n = 3 (i.e. only one pixel out of each contiguous 3 × 3-pixel area is used to create the reconstruction grid in the source plane) while keeping fixed the best-fitting parameter values of the smooth lensing potential inferred with the highest resolution (n = 1). As is shown in Fig. 11 for the lens system SDSS J0252 + 0039, this approach allows us to prevent overfitting and leads to the residual power spectrum lying at or above the noise level for all k-bins. In Section 5, we additionally demonstrate that the choice of n = 3 enables us to successfully recover the known true surface-brightness anomalies in a mock lens system mimicking SDSS J0252 + 0039.

Figure 11.

Power spectrum of the residual surface-brightness fluctuations remaining in the lensed images of SDSS J0252 + 0039 after the smooth lens modelling. The lens modelling is performed alternatively with the highest source-grid resolution (n = 1), using either an adaptive (magenta line) or a non-adaptive (yellow line) source-grid regularization, or with a lower source-grid resolution (n = 3; blue line) chosen for this study to prevent overfitting. As an alternative solution, the red line shows the effect of using a larger mask in the smooth lens modelling combined with n = 1 and a non-adaptive source-grid regularization (for consistency reasons, the power spectrum is calculated within the original tight mask). The sky-background power spectrum (black line) is estimated based on a sample of 20 blank-sky regions (overlaid with the original tight mask) in the proximity to the lens. The estimated total noise power spectrum (green line) accounts for the additional flux-dependent (Poisson) shot noise in the science image.

Open in new tab Download slide

Alternatively, the application of a larger mask, including more noise-dominated pixels in the smooth-lens-modelling procedure (in combination with n = 1 and a high non-adaptive source-grid regularization), might offer another solution to alleviate the overfitting problem in the source reconstruction. As apparent from Fig. 11, both options lead to almost identical power spectra of the residual surface-brightness fluctuations in the lensed images (consistently calculated within the original tight mask), except for a small difference in the lowest analysed k-bin. Our preliminary tests confirm that this does not affect our final results (i.e. exclusion probabilities of the matter-power-spectrum models, see Paper II) significantly. However, since the respective lens models are obtained using different masks and, thus, cannot be considered as inferred from the same data set, a proper comparison in terms of the Bayesian evidence is not possible. We are planning to investigate this alternative in more detail in Paper III of this series (Bayer et al., in preparation).

4.5 Noise power spectrum analysis

Besides possible surface-brightness anomalies due to mass structure in the lens galaxy, the residual surface-brightness fluctuations remaining in the lensed images after subtraction of the best-fitting smooth lens model are partially caused by the observational noise. In this section, we estimate the noise contribution to the measured power spectrum of the residual surface-brightness fluctuations.

4.5.1 The noise-sigma maps

Formally, the observational noise in our HST-imaging can be thought of as a random field. Each pixel is assigned a random variable representing the flux noise, with the expectation value equal to zero (after sky subtraction) and a flux-dependent variance. Two main contributions to this variance are the random fluctuations of the sky background and the Poisson-distributed photon-shot noise from the observed lens system. These are independent random processes, thus the total variance of the observational noise |$\sigma _{n}^2$| in a given pixel can be expressed as the sum of the two variance components:

$$\begin{eqnarray} \sigma _{n}^2 = \sigma _{\rm sky}^2 + \sigma _{P}^2. \end{eqnarray}$$

(13)

We approximate the standard deviation of the sky-background fluctuations σ_sky in the analysed images by the standard deviation of the flux values measured in a sample of close-by blank-sky cutouts σ_sky = 0.002 e⁻ sec⁻¹. The variance of the photon-shot noise |$\sigma _{P}^2$|⁠, on the other hand, is a priori not known but, following the Poisson distribution, equal to the expected flux. The latter can be (due to the large number of counts in our imaging) approximated by the number of electrons per second N measured in a given pixel (after the sky-background subtraction), weighted by the inverse-variance weight W from the weight map computed in the process of drizzling:

$$\begin{eqnarray} \sigma _{P}^2 = N/W. \end{eqnarray}$$

(14)

The estimated standard deviation of the total observational noise σ_n in all individual pixels of our HST-images is finally presented in the form of noise-sigma maps. As an example, Fig. 12 illustrates this procedure for the lens system SDSS J0737 + 3216.

Figure 12.

The process of generating the noise-sigma map for the observed HST/WFC3/F390W-image of the lens system SDSS J0737 + 3216: the drizzled science image with photon counts N (left-hand panel), the inverse-variance weight map W from the drizzling procedure (middle panel), and the resulting noise-sigma map σ_n [see equation (13), right-hand panel]. The standard deviation of the sky background σ_sky is approximated by the value measured in a sample of close-by blank-sky cutouts with the same size as the science image.

Open in new tab Download slide

A noise-sigma map provides a complete description of the noise properties in an image, provided that the random flux fluctuations in the different pixels are statistically independent from each other. However, as we discuss in Section 4.5.2, drizzled images are known to show noise correlations between adjacent pixels, which requires a more thorough noise analysis.

4.5.2 Noise correlations due to drizzling

Despite the fact that the individual pixels in raw HST/WFC3/F390W-images can (ideally) be considered independent, drizzled images are known to show noise correlations. In the process of drizzling, pixels from multiple dithered exposures are aligned and mapped (or informally drizzled) onto a common output grid, based on the relative shift and rotation (i.e. dither) of the respective exposure. The flux of each input pixel is then redistributed over all overlapping output pixels (according to the fractional overlap), which introduces correlations between adjacent pixels in the final drizzled image (for a detailed discussion on noise correlations in drizzled images see Casertano et al. 2000).

In order to investigate the noise-correlation pattern in our data, we create a sample of drizzled blank-sky cutouts located in the proximity to each analysed lens system and quantify their statistical properties in terms of the azimuthally averaged power spectrum. For consistency reasons, we match the size of these blank-sky fields to the size of the respective science image. As an example, Fig. 13 depicts one of these drizzled blank-sky fields, located in proximity to the lens system SDSS J0252 + 0039, in comparison to a realization of uncorrelated Gaussian noise with the same variance of the flux values. Whereas the latter represents the true statistically independent fluctuations of the sky background, the drizzled blank-sky image exhibits a distinct blotchy correlation pattern.

Figure 13.

Noise correlations introduced by drizzling (default configuration): a drizzled blank-sky cutout observed with HST/WFC3/F390W (left-hand panel) versus a realization of an uncorrelated Gaussian noise representing the true statistically independent fluctuations of the sky background (right-hand panel). Both images have the same total variance of flux values, however, the drizzled image shows a distinct patchy correlation pattern.

Open in new tab Download slide

Fig. 14 presents the mean power spectrum measured in a sample of 20 such blank-sky fields located in proximity to SDSS J0252 + 0039, in comparison to a sample of 20 realizations of uncorrelated Gaussian noise with the same total variance. Whereas the power spectrum of the uncorrelated Gaussian-noise realizations is flat (as expected), the power spectrum of the drizzled blank-sky cutouts is scale-dependent and carries a signature of the correlation pattern imposed by the drizzling procedure. The variance in the drizzled images is larger than in the Gaussian-noise realizations on large spatial scales but smaller on small spatial scales (below 0.1 arcsec or 2.5 pixels).

Figure 14.

Azimuthally averaged power spectrum of the sky background in drizzled HST/WFC3/F390W-images: the mean azimuthally averaged power spectrum measured in a sample of 20 drizzled (default configuration) blank-sky regions located in proximity to SDSS J0252 + 0039 (green line) versus the flat power spectrum measured in a sample of 20 mock realizations of uncorrelated Gaussian noise (see Fig. 13 for an example) with the same variance as in the observed blank-sky cutouts (black line).

Open in new tab Download slide

However, a crucial feature of the Drizzle algorithm is the possibility to improve the spatial resolution and reduce noise correlations in the final drizzled image by simultaneously decreasing the pixel scale of the output grid and shrinking the input pixels before mapping them onto the finer output grid (Fruchter & Hook 2002). The pixel size of the output grid is controlled by the final pixscale parameter, whereas the size of the shrunken input pixels, called drops, is varied by means of the final pixfrac parameter. The latter sets the ratio between the linear size of the drop and the original input pixel. The flux of each drop is then redistributed among overlapping output pixels with a weight proportional to the overlap. In comparison to the default configuration, in which both the output pixel and the drop have the same size as the original input pixels, i.e. final pixscale = 0.04 arcsec (for HST/WFC3/UVIS-imaging) and final pixfrac = 1, the flux of the shrunken pixels is redistributed among fewer output pixels, which could help reduce the noise correlations between adjacent pixels.

To test this possibility of reducing the noise correlations in our HST/WFC3/F390W images, we select a blank-sky region in proximity to SDSS J0252 + 0039 and perform a power-spectrum analysis of the random surface-brightness fluctuations for different values of the drizzling parameters final pixscale and final pixfrac. We decrease the final pixscale gradually, from the original value of 0.04 arcsec to 0.033, 0.025 and, finally, 0.02 arcsec. The final pixfrac can be, in principle, varied between 0 (equivalent to sampling with a delta function) and 1 (drop size equal to the original pixel size), but we follow the recommended practice and set the drop size such that it is in each case slightly larger than the output pixels. Fig. 15 presents the resulting azimuthally averaged power spectra, in comparison to the default drizzling configuration.

$Effect of the drizzling parameters final pixscale (linear size of the output pixels) and final pixfrac (size of the drops) on the azimuthally averaged power spectrum of surface-brightness fluctuations in a selected blank-sky region located in proximity to SDSS J0252 + 0039 (see left-hand panel of Fig. 13). The default configuration corresponds to final pixscale = 0.04 arcsec and final pixfrac = 1 (red line).$

Figure 15.

Effect of the drizzling parameters final pixscale (linear size of the output pixels) and final pixfrac (size of the drops) on the azimuthally averaged power spectrum of surface-brightness fluctuations in a selected blank-sky region located in proximity to SDSS J0252 + 0039 (see left-hand panel of Fig. 13). The default configuration corresponds to final pixscale = 0.04 arcsec and final pixfrac = 1 (red line).

Open in new tab Download slide

As can be seen from Fig. 15, lowering the final pixscale (i.e. decreasing the output pixel size) and the final pixfrac (i.e. shrinking the input pixels) when combining the dithered exposures does not allow us to substantially reduce the noise correlations i.e. flatten the noise power spectrum. Moreover, it has a significant effect on the measured power spectrum only in the highest-k bins in which the residual surface-brightness fluctuations revealed in the modelled lens systems have the lowest amplitudes and are very close to the noise level. We stress that any choice of the final pixscale and final pixfrac values allows a valid analysis only if applied to both the science image and the blank-sky cutouts that are used for the estimation of the noise power spectrum. Taking into account that the choice of a lower output pixel scale, while maintaining the same field of view, would substantially increase the number of pixels in the analysed images and, thus, the computational effort of our study (especially the lens modelling), we conclude that the default configuration of the drizzling procedure (i.e. both the input and output pixel size equal to the original pixel size) is a suitable choice for our analysis.

4.5.3 Effect of charge-transfer inefficiency

Due to a gradual degradation process of the HST/WFC3/UVIS-CCDs, the analysed images are additionally affected by the charge-transfer inefficiency (see e.g. Baggett, Gosmeyer & Noeske 2015). This is caused by the radiation damage in space, which leads to defects in the silicon lattice of the CCDs and the formation of spurious trails in the observed images (Massey et al. 2014).

In order to investigate the impact of this issue on the noise properties in our imaging, we perform the drizzling procedure of a selected blank-sky cutout located in vicinity to SDSS J0252 + 0039 using the charge transfer efficiency (CTE) corrected flat-field-calibrated exposures (flc.fits files) and compare the resulting image to the corresponding image based on the default flt.fits files. A careful visual inspection of these two images, presented in Fig. 16, leads to the conclusion that the CTE correction allows us to reduce the level of random surface-brightness fluctuations on the largest spatial scales (smallest k-values). This effect becomes even more apparent in Fourier space. As can be seen from Fig. 17, the CTE correction results in a significant reduction of the noise variance on the largest considered spatial scales (more specifically, the power for k_min = 0.88 arcsec⁻¹ corresponding to the spatial length scale λ = 1.13 arcsec or ∼25 pixels is roughly 40 per cent lower after the CTE correction). However, it does not affect the measured power spectrum on smaller spatial scales (corresponding to higher k-values).

Figure 16.

The effect of the CTE correction: a drizzled HST/WFC3/F390W blank-sky cutout located in proximity to SDSS J0252 + 0039 before (top panel) versus after the CTE correction (bottom panel). The correction allows us to significantly reduce the level of random surface-brightness fluctuations on large spatial scales, which becomes more apparent in Fourier space, see Fig. 17.

Open in new tab Download slide

Figure 17.

The effect of the CTE correction on the azimuthally averaged power spectrum of random surface-brightness fluctuations in a selected HST/WFC3/F390W blank-sky cutout located in proximity to SDSS J0252 + 0039 (corresponding to Fig. 16).

Open in new tab Download slide

While we recommend the use of CTE-corrected flat-field-calibrated exposures (flc.fits files) in future research, due to a substantial variance reduction of the sky-background fluctuations on the largest spatial scales, in the present paper and the accompanying Paper II, we omit the CTE correction and proceed using the standard flt.fits files. We stress that the choice of either of the two options allows a valid analysis as long as both the science image and the blank-sky cutouts used to estimate the noise power spectrum are created in a consistent way (based either on the flt.fits or flc.fits files). Nevertheless, in order to test the impact of the CTE correction on our final results (i.e. exclusion probabilities of the matter-power-spectrum models presented in Paper II), we additionally compute the exclusion probabilities while excluding the largest considered spatial scales and find that this does not significantly affect the derived constraints. We plan to use CTE-corrected images and study this effect in more detail in the analysis of the Jackpot gravitational lens system SDSS J0946 + 1006, which will be presented in Paper III of this series (Bayer et al., in prep).

4.5.4 Estimation of the total noise power spectrum

Here, we estimate the total noise power spectrum in our HST-imaging, which incorporates the combined effect of the random sky-background fluctuations, the noise correlations introduced in the drizzling procedure, and the flux-dependent photon-shot noise from the observed lens system. To this end, we generate a set of scaled sky-background cutouts, located in proximity to each analysed lens system, which can be seen as mock realizations of the total noise in the corresponding science image.

To create the scaled sky-background cutouts, we make use of the set of drizzled blank-sky cutouts generated in Section 4.5.2. These can be treated (after subtracting the mean value) as realizations of a Gaussian random field with the expected value equal to zero and a constant standard deviation over the entire field of view, which we approximate by the standard deviation of the measured flux values. We first divide the blank-sky cutouts by this standard deviation to convert them to the standard normal distribution (i.e. the expected value equal to zero and the standard deviation equal to one) and, subsequently, multiply them by the respective noise-sigma map of the science image (see Section 4.5.1) to incorporate the flux-dependent photon-shot noise from the observed lens system. The average power spectrum measured in the set of these scaled sky-background cutouts, located in proximity to the respective lens system, constitutes our best estimate of the total noise power spectrum in the observed science image.

For a proper comparison with the measured power spectrum of the residual surface-brightness fluctuations in the lensed images (see Section 4.4.3 and Fig. 11), we perform the power-spectrum analysis of the (original and scaled) sky-background cutouts following the same methodology. In particular, before computing the azimuthally averaged power spectrum as specified in Section 4.4.1, we consistently overlay each sky-background cutout with the same mask outlining the lensed images as used in the analysis of the respective science image and set the remaining pixel values to zero.

As an example, Fig. 11 shows the mean power spectrum of the sky-background fluctuations measured in the sample of 20 masked blank-sky regions in proximity to the lens system SDSS J0252 + 0039 and the total noise power spectrum estimated based on the corresponding scaled blank-sky cutouts, as discussed above. We stress that the deviation between the noise power spectra shown in Figs 11 and 14 is due to the difference in the applied window function; that is, in Fig. 11, the power spectrum is computed based on masked cutouts for a consistent comparison with the residual power spectrum measured within the mask, while in Fig. 14 the cutouts are unmasked to investigate the pure effect of drizzling.

4.6 Power spectrum of surface-brightness anomalies due to small-scale mass structures in the lens galaxy

As a final step of our methodology, we perform a noise correction of the residual surface-brightness fluctuations and infer a (conservative) upper-limit constraint on the power spectrum of surface-brightness anomalies induced in the lensed images by the hypothetical small-scale mass structures in the lens galaxy.

As thoroughly discussed in Section 4.4.3, we obtain our best estimate for the power spectrum of the residual surface-brightness fluctuations in the lensed images after the subtraction of the best-fitting smooth lens model inferred with a lower resolution (corresponding to n = 3) and using the original tight mask (blue line in Fig. 11). This choice allows us to prevent the overfitting problem as well as mitigate the degeneracy between the surface-brightness anomalies due to mass structures in the lens galaxy and the intrinsic surface-brightness fluctuations in the source galaxy itself. As can be seen from Fig. 11 for the lens system SDSS J0252 + 0039, the measured residual power spectrum exceeds in this case the noise power spectrum in the five lowest-k bins (largest considered spatial scales) ranging from k_min = 0.88 to k_max = 7.95 arcsec⁻¹, while it is at the noise level for all higher-k modes.

Assuming that the perturbations due to mass structure and the observational noise are statistically independent, we consider the corresponding power spectra to be additive. Under this assumption, we can simply subtract the estimated total noise power spectrum from the residual power spectrum. The procedure is illustrated in Fig. 18 for the lens system SDSS J0252 + 0039, where the difference of these two power spectra corresponds to the red line (shown only in the range of k-modes for which the residual exceeds the noise level). This noise-corrected power spectrum of the residual surface-brightness fluctuations in the lensed images constitutes our upper-limit constraint on the power spectrum of surface-brightness anomalies due to small-scale mass structures in the lens galaxy and is the final outcome of the methodology introduced in this paper. In the companion Paper II, we intend to extend our methodology and relate this measurement to the statistical properties of the underlying small-scale mass structures (more specifically, the sub-galactic matter power spectrum) in the massive elliptical lens galaxy.

Figure 18.

Upper-limit constraints on the power spectrum of surface-brightness anomalies due to small-scale mass structures in the lens galaxy SDSS J0252 + 0039: the power spectrum of the residual surface-brightness fluctuations after the smooth lens modelling with n = 3 (blue line), the mean power spectrum of the sky background in a sample of 20 blank-sky regions located in proximity to the lens (black line), the estimated total noise power spectrum including the photon-shot noise (green line), and the noise-corrected power spectrum of the residual surface-brightness fluctuations, constituting our upper limit on the power spectrum of surface-brightness anomalies due to small-scale mass structures in the lens galaxy (red line).

Open in new tab Download slide

5 PERFORMANCE TEST

Finally, in this section, we test the performance of the introduced methodology in recovering the true power spectrum of mock surface-brightness anomalies from a simulated image mimicking real HST/WFC3/F390W-observations of the lens system SDSS J0252 + 0039, in which the lens galaxy is perturbed by small-scale sub-galactic mass structures with known statistical properties.

We model these hypothetical small-scale mass structures in the lens galaxy as a realization of GRF potential perturbations |$\delta \psi _{\rm {GRF}}({{\boldsymbol x}})$| superposed on a PEMD-plus-external-shear smooth lensing potential. Following Chatterjee & Koopmans (2018), we assume |$\delta \psi _{\rm {GRF}}({{\boldsymbol x}})$| to be fully characterized by a power-law power spectrum:

$$\begin{eqnarray} P_{\delta \psi }(k) = A \times k^{-\beta } \end{eqnarray}$$

(15)

with two free parameters, i.e. the variance of the potential perturbations |$\sigma ^2_{\delta \psi }$| and the power-law slope β. The power spectrum obeys the following normalization condition:

$$\begin{eqnarray} A \Big (\sigma ^2_{\delta \psi }, \beta , L \Big) = \frac{L^2 \sigma ^2_{\delta \psi }}{\sum _{k_{x}}\sum _{k_{y}} \Big (\sqrt{k_{x}^2+k_{y}^2} \Big) ^{-\beta }}, \end{eqnarray}$$

(16)

where the wavenumbers k_x, k_y correspond to the reciprocal wavelength of the associated harmonic waves e^{−2πik · x} propagating in the x and y direction in the Fourier representation of the GRF and L indicates the side length of the analysed image measured in arcsec, see Paper II for a more thorough discussion of the applied formalism.

To simulate the effect of such potential perturbations on the lensed images of SDSS J0252 + 0039, we add a realization of |$\delta \psi _{\rm {GRF}}({{\boldsymbol x}})$| to the best-fitting PEMD-plus-external-shear smooth lensing potential, inferred for this system in Section 4.3, and repeat the lensing operation of the reconstructed source galaxy. We set |$\sigma ^2_{\delta \psi } = 4 \times 10^{-4}$| and β = 4, such that the power spectrum of the induced mock surface-brightness anomalies resembles the residual power spectrum revealed in the real system. To account for the observational effects, we convolve this simulated image with the Tiny-Tim PSF of the HST/WFC3/F390W-optics and add a realistic noise realization. For simplicity, this noise realization is generated based solely on the noise-sigma map (see Section 4.5.1) and, thus, does not reflect the noise correlations found in the drizzled images (see Section 4.5.2). This simplifying assumption is justified by a low level of noise compared to the surface-brightness anomalies induced by the small-scale mass structures.

We perform smooth lens modelling of this simulated image using the same methodology that was applied to the real observed data, within the same mask outlining the lensed images. The modelling is carried out without re-optimizing for the parameter values of the best-fitting smooth lensing potential. By doing so, we assume that the parametric lensing potential can be reconstructed accurately and focus instead on investigating the degeneracy between the anomalies caused by the small-scale mass structures in the lens galaxy and the intrinsic surface-brightness fluctuations in the source galaxy itself. We will test the validity of this assumption in Paper III of this series (Bayer et al., in preparation). In the current modelling procedure, we apply an adaptive source-grid regularization and varying levels of the source-grid resolution, i.e. the number of pixels cast back from the lens plane to the source plane, corresponding to n = 1, 2, 3, and 4 (see Section 4.3). In each case, we determine the resulting power spectrum of the residual surface-brightness fluctuations remaining in the lensed images after the subtraction of the best-fitting smooth lens model and compare it with the known true power spectrum of the imposed surface-brightness anomalies.

The results of this performance test are presented in Fig. 19. From this figure, it can be seen that the power spectrum of the residual surface-brightness fluctuations lies significantly below the noise level when the modelling is performed with the highest source-grid resolution (n = 1; i.e. each pixel is cast back from the lens plane to the source plane). As in the analysis of the real system, this overfitting can be explained by the absorption of the induced surface-brightness anomalies, and partially even the observational noise, in the source structure. However, this degeneracy can be alleviated by lowering the resolution (i.e. choosing higher n-values) of the adaptive source grid, which leads to a better agreement between the reconstructed and the true residual power spectrum. A comparison of the power spectra corresponding to n = 3 and n = 4 suggests that convergence is reached for n = 3 and lowering the source-grid resolution even further does not allow us to thoroughly suppress this degeneracy (but would lead to a considerably less accurate source reconstruction). The absorption of the potential perturbations into the source structure persists on the smallest considered k-scales.

$Performance test of the introduced methodology in recovering the power spectrum of mock surface-brightness anomalies induced in the lensed images of SDSS J0252 + 0039 by a realization of Gaussian random field potential fluctuations $\delta \psi _{\rm {GRF}}({{\boldsymbol x}})$ with known statistical properties. The plot depicts the true azimuthally averaged power spectrum of the induced surface-brightness anomalies (blue line), the power spectrum of the overlaid white-noise realization (black line), the power spectrum of these two effects combined (i.e. the power spectrum to be recovered, red line), and the power spectra of the residual surface-brightness fluctuations actually recovered from the perturbed lensed images after the subtraction of the best-fitting smooth lens model obtained with varying source-grid resolution (n = 1, 2, 3, and 4; green lines from bottom to top).$

Figure 19.

Performance test of the introduced methodology in recovering the power spectrum of mock surface-brightness anomalies induced in the lensed images of SDSS J0252 + 0039 by a realization of Gaussian random field potential fluctuations |$\delta \psi _{\rm {GRF}}({{\boldsymbol x}})$| with known statistical properties. The plot depicts the true azimuthally averaged power spectrum of the induced surface-brightness anomalies (blue line), the power spectrum of the overlaid white-noise realization (black line), the power spectrum of these two effects combined (i.e. the power spectrum to be recovered, red line), and the power spectra of the residual surface-brightness fluctuations actually recovered from the perturbed lensed images after the subtraction of the best-fitting smooth lens model obtained with varying source-grid resolution (n = 1, 2, 3, and 4; green lines from bottom to top).

Open in new tab Download slide

All in all, based on the results of this performance test with a realistic mock lens, we conclude that our approach allows us to recover the true power spectrum of mock surface-brightness anomalies when the smooth lens modelling is performed with n = 3 as the most suitable source-gird resolution for the given data quality and the choice of a relatively narrow mask.

6 CONCLUSIONS AND OUTLOOK

In this paper, the first in a series, we have introduced and tested a novel methodology to reliably measure the power spectrum of surface-brightness anomalies in extended lensed images of galaxy-galaxy strong gravitational lens systems. To illustrate our approach, we have applied it to a SLACS sub-sample observed with HST/WFC3 in the ultraviolet and discussed the modelling challenges. Finally, as a proof of concept, we have demonstrated the feasibility of the introduced methodology by recovering the true power spectrum of mock surface-brightness anomalies from simulated lensed images mimicking real HST/WFC3/F390W-observations of the lens system SDSS J0252 + 0039.

One of the main challenges in the power-spectrum measurement turned out to be the degeneracy between the surface-brightness anomalies due to the presence of mass structures in the lens galaxy and the intrinsic surface-brightness fluctuations in the source galaxy itself. While this degeneracy is less problematic in the case of the direct detection of individual subhaloes with masses above the detection limit, as in Vegetti et al. (2014), this issue requires a more careful consideration in the power-spectrum approach. Our test on simulated lensed images has shown that the degeneracy can be alleviated by performing the smooth lens modelling with a lower source-gird resolution to prevent overfitting.

In the companion Paper II (Bayer et al. 2023), our main objective is to extend the introduced methodology such that the estimated power spectrum of surface-brightness anomalies in the lensed images can be traced back to the statistical properties of the underlying small-scale mass structures in the lens galaxy. With this goal in mind, we carry out a systematic study of mock surface-brightness anomalies induced by GRF potential perturbations with varying statistical properties and compare the results to the real measurement performed in the present paper. For a pilot application of the extended methodology, we choose one of the lens systems from the investigated SLACS sub-sample, SDSS J0252 + 0039, due to its simple geometry and a high signal-to-noise ratio of the lensed images. As a final result of the combined analysis, we infer the first observational constraints on the matter power spectrum in a massive elliptical (lens) galaxy.

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

http://archive.stsci.edu/hst/search.php

https://www.astromatic.net/software/stiff

http://www.stsci.edu/hst/observatory/focus/TinyTim

http://ds9.si.edu/site/Home.html

http://users.obs.carnegiescience.edu/peng/work/galfit/galfit.html

https://docs.scipy.org/doc/numpy/reference/routines.fft.html

References

Auger

M. W.

Treu

Bolton

A. S.

Gavazzi

Koopmans

L. V. E.

Marshall

P. J.

Bundy

Moustakas

L. A.

2009

ApJ

705

1099

Month:	Total Views:
May 2023	6
June 2023	26
July 2023	10
August 2023	8
September 2023	4
October 2023	16
November 2023	16
December 2023	22
January 2024	14
February 2024	18
March 2024	15
April 2024	7
May 2024	15
June 2024	15
July 2024	30
August 2024	11
September 2024	16
October 2024	21
November 2024	17
December 2024	6
January 2025	11
February 2025	13
March 2025	10
April 2025	9

Article Contents

Probing sub-galactic mass structure with the power spectrum of surface-brightness anomalies in high-resolution observations of galaxy-galaxy strong gravitational lenses – I. Power-spectrum measurement and feasibility study

ABSTRACT

1 INTRODUCTION

2 SURFACE-BRIGHTNESS ANOMALIES IN EXTENDED LENSED IMAGES

3 OBSERVATIONAL STRATEGY AND DATA

3.1 Selection of the observational filter

3.2 Sample selection

3.3 Observations and data reduction

4 METHODOLOGY

4.1 Analysis synopsis

4.2 Lens-galaxy subtraction

4.3 Smooth lens modelling

4.4 Power-spectrum analysis of the residual surface-brightness fluctuations

4.4.1 Power-spectrum measurement

4.4.2 Effect of the lens-galaxy subtraction

4.4.3 Effect of the smooth lens modelling

4.5 Noise power spectrum analysis

4.5.1 The noise-sigma maps

4.5.2 Noise correlations due to drizzling

4.5.3 Effect of charge-transfer inefficiency

4.5.4 Estimation of the total noise power spectrum

4.6 Power spectrum of surface-brightness anomalies due to small-scale mass structures in the lens galaxy

5 PERFORMANCE TEST

6 CONCLUSIONS AND OUTLOOK

ACKNOWLEDGEMENTS

DATA AVAILABILITY

Footnotes

References

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only