© ESO 2022

1. Introduction

The Hubble constant H₀ is a key parameter for understanding the past and future evolution of our Universe. Achieving a subpercent accuracy on its value is essential for exploring cosmological models beyond the standard Λ cold dark matter (CDM). Several teams have achieved high-precision measurements using different methods (e.g. Abbott et al. 2018; Riess et al. 2019; Freedman et al. 2020; Pesce et al. 2020; Philcox et al. 2020; Planck Collaboration VI 2020; Schombert et al. 2020; Wong et al. 2020; Blakeslee et al. 2021). In recent years, the H₀ inference has raised tension between early-time probes, such as measurements from the cosmic microwave background and baryonic oscillations, and late-time probes using distance-ladder techniques (Verde et al. 2019). The time-delay technique, which consists of inferring the ‘time-delay distance’ from the modelling of strongly lensed quasars for which time-delays have been measured (Refsdal 1964) is among the late-time probes, and confers the advantage of providing a distance estimate not anchored on a local distance ladder. However, the accuracy of this technique is currently based on an a priori knowledge of early-type galaxies and of their mass density profiles (e.g. Birrer et al. 2016; Ding et al. 2021a).

The main lensing galaxy mass distribution is often assumed to follow an elliptical power-law density profile or to be the sum of a baryonic component and dark matter elliptical density component. The role of the assumption of a radial density profile in the time-delay distance inference has been extensively discussed in the literature (Schneider & Sluse 2013, 2014; Xu et al. 2016; Kochanek 2020, 2021; Birrer et al. 2020). Less attention has been given to quantifying the impact of an angular change of the density profile on the distance inference. Suyu et al. (2010) quantified the impact of angular structures on H₀ for the time-delay lens system B1608+656 by allowing for a pixelized perturbation of an elliptical density profile. These authors showed that, for that system, perturbations of < 2% of the lens potential were required, yielding a change in H₀ of less than 1%. This result motivated the use of mostly elliptical mass distribution in subsequent time-delay studies by the H0LICOW and STRIDES Collaborations. Nevertheless, a more generic investigation of the role of angular structure on H₀ is still missing. In the present work, we focus on specific angular structures: Fourier-like/multipole perturbations to the elliptical shape of the density profile. The study of early-type galaxies shows that their isophotes are well described by ellipses with Fourier deviations, also called multipoles, of orders three and four (e.g. Rest et al. 2001; Hao et al. 2006; Pasquali et al. 2006; Krajnović et al. 2013; Mitsuda et al. 2017). While the third-order components are generally of small amplitude, those of order four can be substantial, giving rise to clear discy or boxy shapes of the light profile. To our knowledge, the imprint left by multipoles present in the lensing galaxy on extended lensed images has not been discussed in detail in the literature, nor has its impact on H₀.

Due to lensing cross-section, massive early-type galaxies constitute the vast majority of the lens population discovered so far. They are commonly well modeled with density profiles that are perfectly elliptical. Apart from nearby galaxies explicitly included in the model, the only source of angular perturbation generically considered is tidal perturbation –also known as external shear– caused for example by other galaxies along the line of sight towards the lens. Those perturbations of the gravitational potential ψ are quadrupolar, of the form ψ ∝ cos(2χ), where χ is the azimuthal angle (Trotter et al. 2000; Kochanek 2006; Chu et al. 2013). Higher order multipoles, corresponding to ψ ∝ cos(mχ), are often ignored as they are expected to displace lensed image positions by Δθ < 0.003″ (Sect. 4 in Kochanek 2006). This forecast however relies on a rather small multipole amplitude (a₄/θ_E ∼ 0.01, where θ_E is the Einstein radius of the considered lens), and on the assumption that the light provides a good proxy of the perturbation existing in the total mass.

Only a few studies have included an investigation of the constraints offered by lensing observables on those multipoles. Trotter et al. (2000) used very long baseline interferometry (VLBI) observations of multiple structures in MGJ0414+0534 to constrain high-order multipoles (m = 3, 4) in that system. Claeskens et al. (2006) also used the compact star formation regions observed in J1131−1231 to constrain m = 3, 4 terms of the multipole expansion. In both cases, the physical meaning of those multipoles was undetermined as either their amplitude or direction could hardly fit any reasonable physical model. The role of multipoles in introducing scatter around the fundamental surface of quad lenses (Gomer & Williams 2018) was recently investigated, revealing that the multipoles may not be sufficient to explain the observed deviation of lensed quasar positions from this surface. The presence of high-order multipoles can also modify flux ratios between quasar images (e.g. Möller et al. 2003; Winn et al. 2003; Keeton et al. 2003, 2005). Those perturbations to the flux ratios introduced by macro-structures in the density profile need to be accounted for when using flux ratio anomalies to constrain the halo mass function (e.g. Mao & Schneider 1998; Xu et al. 2015; Gilman et al. 2019, 2021; Hsueh et al. 2020). Observational constraints on the amplitude of the high-order multipoles associated to the total mass distribution would therefore be highly valuable, as they could be used to quantify the impact of low-amplitude boxyness and discyness on flux ratio anomalies.

With the advent of deep, high-angular-resolution images of lensed systems, as achievable with the Hubble Space Telescope (HST) or with adaptive optics systems, it is now possible to detect extended images of distant lensed galaxies, even those hosting a quasar (e.g. Bolton et al. 2006; Lagattuta et al. 2012; Suyu et al. 2017). The Einstein rings are known to provide numerous constraints on the mass model and the extended source shape (Kochanek et al. 2001). The only systems for which the shape of the Einstein ring has been used to constrain the amplitude of multipoles are SDSS J0924+0219, HE 0435−1223, B1938+666, and PG 1115+080 (Yoo et al. 2005, 2006). For each galaxy, the deviation of the mass density from ellipsoidal shape because of fourth-order multipoles was found to be consistent with zero, and only upper limits on a₄ could be derived. Extended lensed images are now commonly reproduced by lensed models down to the noise, even when the density profile of the lensing galaxy is modelled as an elliptical power-law density profile (Birrer et al. 2016; Shajib et al. 2019).

The primary aim of this paper is to assess the detectability of multipolar components, and our second goal is then to quantify the extent to which they modify the H₀ inference. By quantifying the detectability of multipoles in high-resolution lensed images, we intent to find out if multipole-component-induced deformations are effectively negligible in those systems or have been absorbed by the model. In the latter case, we quantify the deviation of the time-delay distance induced by the fitted elliptical model. While this is the main motivation of this work, the presence of multipoles in the total mass profile is also important for our general understanding of galaxies and of the ability of interactions between baryons and dark matter to shape the galaxy density profile. It is also relevant for dark matter studies based on flux ratio anomalies (Hsueh et al. 2017, 2020; Gilman et al. 2020a,b, 2021; Nierenberg et al. 2020): the ability to constrain the presence of (dark) multipole-like structures in the lens from the analysis of the lensed images would enable further tightening of the constraints on the dark matter properties inferred from flux ratio anomalies.

In Sect. 2, we introduce our general methodology, the mathematical definitions of the main lens mass profiles used in our analysis, and explain mock-image creation and fitting procedure. In Sect. 3, we systematically test the detectability of multipoles depending on different parameters and discuss their detectability in real cases. In Sect. 4, we then quantify the effect of the presence of multipoles on H₀ inference at a single-galaxy level and at a population level and assess their potential influence in cosmological lensing analyses. We summarise and conclude in Sect. 5.

For all calculations presented in this paper involving non-angular quantities, we assume a flat ΛCDM cosmological model with Ω_m = 0.3, Ω_Λ = 0.7, and H₀ = 70 km s⁻¹ Mpc⁻¹.

2. Methodology

To quantify the detectability of multipoles from strongly lensed images, we designed a controlled set of simulated lensed systems containing multipoles and modelled them using state-of-the-art techniques. Specifically, we emulate mock HST-like data of a quasar and its host lensed by a singular isothermal ellipsoid (SIE) with fourth-order multipoles and shear. Those multipoles are the dominant non-zero multipole contributions that control boxyness or discyness of galaxies (see Fig. 1). We then fit the mock images with an SIE and shear (i.e. without fourth-order azimuthal structures).

Fig. 1. Isodensity contours for an SIE mass profile (grey dashed line) and an SIE + fourth-order multipole mass profile (black plain line). The multipole component has a strength a4 = 0.01 (top), a4 = 0.03 (middle), or a4 = 0.05 (bottom) and is oriented along the SIE main axis to create a discy profile (left), or at 45° to the SIE main axis to create a boxy profile (right).

With this setup, we then analyse the effect of different factors on the quality of the retrieved fit. A good fit means that the multipoles remain unnoticed in imaging data of a lens system. This is the situation where multipoles may introduce unnoted systematic errors into the analysis, especially into the H₀ inference. Conversely, a poor fit would be an indication of modelling inadequacy. We assume that any structure detected in the residuals¹ would be identified as a need to include multipoles in the lens mass model. However, we do not explore the ability of macro-models to effectively recover the input multipoles. Instead, we seek the limit at which visible residual patterns become present over the noise. This depends on the mass and source model assumptions and on the data quality. Specifically, we explore the role of the shear, the lens axis-ratio, the additional freedom granted by relaxing the slope constraint in the fit –and the influence of the radial mass profile in that case–, the signal-to-noise ratio (S/N), the multipole strength, and the freedom allowed in the quasar host galaxy morphology during the fit.

We finally look at the impact of unnoticed multipolar components on time-delay cosmography. Focusing on our previous mock systems in cases where multipoles remain undetected, we first explore the H₀ recovery for those specific systems. We then extend the study by constructing a sample of lenses that emulates a realistic population of lenses and look at the combined H₀ inference for different S/N levels.

2.1. Lens mass profile

To construct mock images, one needs to assess a lensing mass profile. As we are interested in the effect of azimuthal structures such as boxyness and discyness in the lensing galaxy, we first select a radially simple mass profile. The density profile of lensing galaxies is generally well approximated by a mass distribution that follows a power-law density profile and some external shear (Suyu et al. 2009, 2010; Treu & Marshall 2016; Yıldırım et al. 2020; Wong et al. 2020). The mean slope of the profile has been found to be close to isothermal (Courbin & Minniti 2002; Koopmans et al. 2003; Suyu et al. 2010), such that the singular isothermal ellipsoid (SIE) provides a simple but realistic choice of density profile for the lens². In addition, we add a multipole perturbation to the ellipticity and include a shear term to account for the tidal perturbation of the potential, such as the one due to line-of-sight mass distribution. Below, we describe the lensing potential associated to each of these components. The deflection angle and Hessian matrix can be found analytically by differentiation of the potential.

Following Keeton & Kochanek (1998), the potential of our SIE is given by the potential of a non-singular isothermal ellipsoid (NIE) with a virtually null core size. This translates into a potential –for an SIE aligned with the x-axis– at position θ = (x, y):

$$\begin{array}{c}\hfill {\psi }_{\mathrm{SIE}}(\mathit{\theta })=x{\alpha }_x+y{\alpha }_y,\end{array}$$

with

$$\begin{array}{cc}& {\alpha }_x(\mathit{\theta })=\frac{b}{\sqrt{1-q^2}}arctan\left(\frac{\sqrt{1-q^2}x}{\sqrt{q^2{x}^2+y^2}}\right),\hfill \\ \hfill & {\alpha }_y(\mathit{\theta })=\frac{b}{\sqrt{1-q^2}}\mathrm{arctanh}\left(\frac{\sqrt{1-q^2}y}{\sqrt{q^2{x}^2+y^2}}\right),\hfill \end{array}$$

where q is the axis ratio, and b is the scale radius and is equal to ${\theta }_{\mathrm{E}}\sqrt{q}$, where θ_E is the Einstein radius. The difference in potential between two SIEs with different ellipticities is a quadrupole moment.

The multipole has several definitions in the literature. We chose to use the same convention as Xu et al. (2015) and Keeton et al. (2003). The lensing potential at a position θ = (x, y) = (θ cos ϕ, θ sin ϕ) is given by:

$$\begin{array}{c}\hfill {\psi }_{\mathrm{multipole}}(\mathit{\theta })=\theta \frac{1}{1-m^2}{a}_{\mathrm{m}}cos(m(\varphi -{\varphi }_{\mathrm{m}})),\end{array}$$

where m is the order of the multipole, a_m is the multipole strength, and ϕ_m is the main axis orientation angle. The fourth-order multipole is an octupolar moment of the potential. It introduces a ‘discy’ shape when ϕ₄ = 0 and a ‘boxy’ shape when ϕ₄ = 45° (see Fig. 1). In the sample studied by Hao et al. (2006), 70% of elliptical galaxies display a multipole a₄ < 0.02 (see Fig. 2). We note that other authors (e.g. Bender et al. 1988, 1989; Hao et al. 2006; Penoyre et al. 2017; Frigo et al. 2019) prefer the normalised a₄ and b₄ convention, using (a₄/a)_conv and (b₄/a)_conv where a is the length of the semi-major axis of the reference ellipse, (a₄/a)_conv represents the cosine deformations of that ellipse, and (b₄/a)_conv represents the sine deformations. This normalised convention is more appropriate when dealing with specific isophotes at different radii, while our a₄ convention is particularly useful for the creation of a general lens density profile that can be defined independently of the main lens shape. Conversion between the two conventions is explained in Appendix B.2 of Xu et al. (2015) and is summarised here:

$$\begin{array}{cc}& a_4=\sqrt{{(a_4/a)}_{\mathrm{conv}}^2+{(b_4/a)}_{\mathrm{conv}}^2}\times \frac{{\theta }_{\mathrm{E}}{(}^{″})}{\sqrt{q}}\hfill \end{array}$$(1)

$$\begin{array}{cc}& {\varphi }_4=\frac{1}{4}\times arctan\left(\frac{{(b_4/a)}_{\mathrm{conv}}}{{(a_4/a)}_{\mathrm{conv}}}\right)·\hfill \end{array}$$(2)

Fig. 2. Distribution of fourth-order multipoles in the light of local elliptical galaxies from Hao et al. (2006). The two panels display the distribution for the normalised (a4/a)conv and (b4/a)conv convention (left) and our a4 convention used in this paper (right), for θE = 2″ and q = 0.8 (see Eqs. (1) and (2)).

As we only consider pure boxyness and discyness in the following analysis (i.e. we always have (b₄/a)_conv = 0), the backward conversion is simply

$$\begin{array}{c}\hfill {(a_4/a)}_{\mathrm{conv}}=\frac{\pm a_4\times \sqrt{q}}{{\theta }_{\mathrm{E}}{(}^{″})},\end{array}$$(3)

which is positive in the discy case and negative in the boxy case.

The shear is defined by the lensing potential at position θ = (x, y):

$$\begin{array}{c}\hfill {\psi }_{\mathrm{shear}}(\mathit{\theta })=\frac{1}{2}({\gamma }_1{x}^2+2{\gamma }_{2}xy-{\gamma }_1{y}^2),\end{array}$$

where γ₁ and γ₂ are the components of the complex shear. The shear strength ${\gamma }_{\mathrm{ext}}=\sqrt{{\gamma }_1^2+{\gamma }_2^2}$ and its orientation ${\varphi }_{\mathrm{ext}}=\frac{1}{2}arctan({\gamma }_2/{\gamma }_1)$. The shear is a quadrupole moment of the potential.

For our baseline model, we consider θ_E = 2″ and q = 0.8 for the SIE. We orient the SIE with an angle of 22° with respect to the x-axis. The orientation of the multipole ϕ₄ is equal to either 0° or 45° in a reference frame aligned with the main lens, and this will create either discy or boxy mass profiles. The strength of the multipole is a₄ = 0.01 by default. We fix the shear amplitude to 0.05 and its orientation to 30° with respect to the SIE main axis. The shear is therefore neither oriented on a multipole axis nor on the SIE main axis. No spurious artificial shear due to truncation of pixelated mass maps is introduced (Van de Vyvere et al. 2020) because we have an analytical description of each component of the lens mass profile.

While one can argue that using an Einstein radius of 2″ is on the upper edge of the distribution of Einstein radii in observed lensing systems, it is particularly relevant to use such an θ_E in our study to highlight the effect of multipoles. As explained in Sect. 3.8, for a fixed exposure time and source magnitude and shape, a smaller Einstein radius would yield a smaller signal-to-noise ratio in the ring and therefore reduce our ability to explore the lens properties that most impact the multipole detectability.

2.2. Mock images creation

The modelling of extended lensed images requires high-spatial-resolution data, which in general are provided by space-based or by ground-based AO systems (e.g. Lagattuta et al. 2010; Chen et al. 2016; Suyu et al. 2017). In this work, we chose to emulate mock³ HST images observed with the WFC3 camera through the filter F160W. This choice is guided by the data used in many time-delay cosmography studies (Suyu et al. 2017), and in other studies using high-quality mock lens images (e.g. Ding et al. 2021a; Park et al. 2021; Wagner-Carena et al. 2021). The near-infrared (NIR) band confers the advantage that, in general, it provides the highest host brightness (see e.g. Ross et al. 2009), and is a priori the most favourable band with which to analyse extended lensed features. In addition, we also consider a transparent lens, which is equivalent to assuming that the lens light is perfectly subtracted and provides a negligible contribution to the photon noise in the ring. Due to the lack of lens light, this setup slightly increases the S/N in the ring compared to a real NIR HST image. On the other hand, state-of-the-art lens modelling is based on multi-band HST data. This increases the effective signal in the ring. Moreover, the bluer bands have higher angular resolution and better sampling of the PSF due to smaller pixels, even if they are not able to fully balance the lack of host flux at those wavelengths. The relative advantage of the bluer bands is also that the lens light drops off more quickly, and therefore blends less with the ring. This last benefit is counterbalanced in our setup by the use of a transparent lens. Therefore, our single-band setup should, at least partially, compensate for the lack of bluer bands and give a reasonable indication of the detectability of multipole-induced deformation of the ring.

We used the same PSF as the one used in the Time Delay Lens Modeling Challenge (TDLMC, Ding et al. 2021a) for the Rung 2 and 3, created from the drizzling of eight PSFs extracted from real HST images. We use a pixel size of 0.08″, typical of drizzled NIR HST images. We include photon and read-out noise (RON). The zero-point in the AB system is 25.96 mag. Our baseline setup corresponds to an exposure of 5400 s. The ‘sky’ brightness, that is zodiacal light and Earth shine, is fixed at 22.3 mag arcsec⁻² and the RON is set to 21 e⁻. Those values are based on the WFC3 documentation Dressel (2012) and concur with other works creating typical HST WFC3/F160W images (Ding et al. 2021a; Park et al. 2021; Wagner-Carena et al. 2021).

We consider a simple source light distribution constituted from an AGN on top of a circular host following a Sersic profile. The choice of a circular source is based on the fact that using a more complicated source with a principal axis may produce results that will be highly dependent on its orientation, while the goal of this analysis is to create mock systems where the azimuthal structure of the lens is the only source of potential bias. We set the contrast between the QSO and its host based on the work of Dunlop et al. (2003) and Jahnke & Wisotzki (2003). This corresponds to a host that is brighter than the QSO by 0.5 mag. For our baseline model, the unlensed quasar apparent brightness is set to 21 mag and the unlensed host is set to 20.5 mag. This yields a mean magnitude of the lensed QSO of 18.5 mag. The host galaxy morphology is a circular Sersic with a half-light radius R_sersic = 0.1″ and a Sersic index n_sersic = 3. The source is placed at a redshift z_s = 2.

The lensing galaxy density profiles are described in the previous section. The lens redshift is z_l = 0.271. For each mock lens density profile, we only consider quad configurations. We focus on four-lensed-images configurations for two reasons: first, cosmographic studies prioritise such systems as they provide multiple delays. Second, dark matter investigations based on flux ratio anomalies currently only rely on samples of quads. More specifically, we consider four quad configurations of lensed images by positioning the source at various locations with respect to the inner caustic. The four sources and associated image positions are illustrated in Fig. 3. The source is either close to the long axis cusp (cusp_l), close to the short axis cusp (cusp_s), close to a fold caustic, or centred (cross image configuration). Specifically, we locate the source at 80% of the distance between the centre and the cusp and fold caustics. The source for the cross-configuration is at a 5% distance from the centre in the direction of the long axis cusp. The lenstronomy software⁴ is used to compute the mock images (Birrer & Amara 2018).

Fig. 3. Illustration of the four lensing configurations used in this work. The inner caustic (bold) and associated critical curve are represented in black. The four different source positions are shown with star symbols and the corresponding image positions are drawn with triangles. Each configuration is displayed in a single colour (see legend).

The time-delays are calculated within lenstronomy. We fix the time-delays to the ‘ground truth’ from the macro model but we consider a 2% uncertainty on the delay with a lower boundary of 1 day when performing the lens modelling.

2.3. Modelling setup

Once the mock image is created, we model it with a standard modelling technique constrained by the pixels inside a mask which omits all pixels that are below twice the background noise level. We use lenstronomy and its particle swarm optimisation (PSO, Kennedy & Eberhart 1995; Shi & Eberhart 1998) to find the model parameters that maximise the likelihood. The likelihood contains terms associated to the image residuals, to the time-delays, and to the position of the four ray-traced back point-source positions in the source plane (see Birrer & Amara 2018). This last likelihood term allows efficient sampling by relaxing the constraint that all point sources should ray-trace back to the exact same position whilst still demanding a sufficiently accurate matching of the lens equation. The model considered is identical to the input mass profile except that no multipoles are included: the lens model is an SIE + shear, the source model is a circular Sersic, and the four point-sources are modelled in the image plane with their positions in the source plane matching the source centre. After the PSO has converged, a Markov chain Monte Carlo (MCMC) using emcee (Goodman & Weare 2010; Foreman-Mackey et al. 2013) is performed to sample the posterior probability distribution.

3. Multipole detectability

We first study the base case described in the previous section. We set the multipole to a₄ = 0.01 and ϕ₄ = 0° with respect to the SIE main axis for discy cases and ϕ₄ = 45° with respect to the SIE main axis for boxy cases; see top panel of Fig. 1. This multipole level of boxyness and discyness is typical of what is observed in the light profile of elliptical galaxies (Fig. 2). The resulting fits are shown in Fig. 4. The residuals and χ² displayed are for the best model, that is, the model with the highest likelihood. The H₀ results are the median value of the MCMC sampling with the error corresponding to the 0.16 and 0.84 quantile. Regardless of the image configuration and multipole angle, we see an alternation of positive and negative residuals beyond the 3σ level when running azimuthally along the ring. Those correlated structures are sufficiently large and numerous to yield a reduced χ² much larger than 1. The impact on H₀ inference is discussed in Sect. 4.

Fig. 4. Results of the base case mock images created and fitted as described in Sect. 2. Those systems are also displayed in red in Table 1. Top 2 rows: mock images made from one boxy galaxy for four configurations and residuals associated to the fit for each image with χ2 and H0 inference. Bottom 2 rows: same as in the top rows but for the discy lensing galaxy.

Our baseline setup shows that the imprint of multipoles on lensed images would be ubiquitous. We now systematically investigate how deviation from the baseline assumption impacts those results. We quantify the influence of the following parameters:

– Shear: we investigate the extreme case where the shear amplitude of the mock system is equal to zero.

– Lensing galaxy ellipticity: we consider three values representative of existing lensing galaxy samples, q = (0.7, 0.8, 0.9).

– Assumption on the fitted lens mass profile: we allow the power-law slope to deviate from the isothermal case in the modelling while keeping the input slope isothermal.

– Assumption on the fitted lens mass profile when the mock lensing galaxy is not a SIE: we extend the previous point by using a mock lensing galaxy with varying profile slope with radius, following a baryon + dark matter prescription, and allow the fitted power-law slope to deviate from the isothermal case.

– S/N: we test three different S/Ns by varying the exposure time and/or the source and quasar brightness.

– Multipole strength: we vary a₄ between 0.01 and 0.05. We note that comparison of those values to multipole amplitudes found in the literature may require the use of (3).

– Assumption of the source morphology: we include shapelets in the source reconstruction during the modelling.

The following sections discuss the results obtained for each setup change. Table 1 summarises the main results of those tests, providing for each case the value of the reduced χ² for the best model and inferred value of H₀ from the MCMC sampling (median, 0.16 quantile, and 0.84 quantile). Two subsamples of mock systems appear multiple times in the table: our baseline case (red) and a lower S/N (3000 s exposure time with a bright source) case with a₄ = 0.01 (green). Our base case is the ideal case for which multipoles would leave ubiquitous imprints in the lensed images (see Fig. 4); it is used as the starting point to compare the influence of different parameters such as the shear, the galaxy ellipticity, and the S/N. Because imprints from multipoles are ubiquitous with the fiducial case, we have built a lower S/N system to study other factors (see Fig. 5). This lower S/N setting is afterwards used for studying the influence of the multipole strength and the addition of source structure (i.e. shapelets) in the modelling. A specific pair of tests illustrating how the source can absorb perturbations introduced by multipoles is shown in Fig. 6. The corresponding cases are outlined in purple in Table 1.

Fig. 5. Same as Fig. 4, but with an exposure time of 3000 s instead of 5400 s and the unlensed quasar and extended source fainter by 1.75 mag. Those results are also displayed in green in Table 1. We note that the mask, which is defined by a threshold of twice the background noise, may sometimes have holes inside the Einstein ring due to the lack of lensed signal at those locations.

3.1. Influence of the shear

We compare the modelling of our base case, which is characterised by γ_ext = 0.05 and an orientation of 30° with respect to the SIE major axis, with systems generated without shear. We find that the absence of shear does not significantly change the goodness of the fit, apart from in two specific cases.

As shown in Table 1, the general trend is the inability of the shear to absorb deformations introduced by multipoles. This is expected as a fourth-order multipole (i.e. m = 4) is an octupolar moment while the shear is a quadrupolar one. However, two mock systems stand out: the discy ones with cusp configurations. For those two cases, we find that good residuals can be achieved. This may be explained by the specific symmetry: the multipole, the SIE, and the source position are aligned, thus creating images in a relatively narrow region with conserved axial symmetry (for the potential, magnification, and deflection). It is thus easier to mimic the effect of input multipole with an SIE+shear in this small symmetric region. In addition, for those two systems, we observe that the fit of an SIE on an SIE+multipole is not equivalent to fitting an SIE+shear on SIE+shear+multipole. The former is not presented in the table but is a good fit, while the latter is a poor fit. This result may appear surprising at first sight as the same additional component (the shear) is added to the mock system and to the model. However, while the lensing potential, deflection, and convergence of the different components sum up linearly, the magnification is not linearly changed, thus explaining the difference between sheared and shearless cases.

3.2. Influence of the lens ellipticity

Early-type lensing galaxies generally have axis ratios in the range q ∈ [0.65 − 0.95] (e.g. Sluse et al. 2012; Rusu et al. 2016; Shajib et al. 2019). We therefore tried three different axis ratios for the lensing galaxy: a more elliptical galaxy with q = 0.7, the base case with q = 0.8, and a rounder galaxy with q = 0.9. We find that the rounder the galaxy, the poorer the fit (see Table 1). This means that the detectability of multipoles is increased when the galaxy is rounder.

The main cause of this effect may be the following: as the ellipticity increases (i.e. q decreases), the size of the inner caustic increases, such that the magnification of a source located at the same relative position with respect to the caustic gets smaller. For a fixed exposure time, this yields a lower effective S/N in the ring, reducing our ability to detect multipole imprints in the lensed host images.

Apart from this general trend, a careful look at Table 1 shows that for a small number of cases the fit slightly improved with the increase in q. This may be due to the specific realisation of the noise, and the slight dependence of the results on the position of the source with respect to the caustics. We note that to test the impact of q, we compared systems holding the value of a₄ to 0.01. According to Eq. (3), this implies a different value of a₄/a for the different q. Keeping a₄/a constant instead of a₄ constant does not change the increase in χ² with the axis ratio q.

3.3. Influence of a relaxed slope constraint

Above, we model the lens with an isothermal mass model, which is our fiducial lensing galaxy profile. However, a power law with a free slope is often used in the modelling of real lensing systems. We therefore followed the same fitting procedure using a power-law elliptical mass density profile (PEMD; Barkana 1998) with an additional step in the fitting that allows the slope to be free to vary. A PEMD is characterised by a 3D density profile $\rho (r)\propto r^{-{\gamma }^{\prime }}$. This model is equivalent to an SIE for a power-law index γ′ = 2. See Appendix B.1 for more details.

As can be seen in Table 1, a slightly lower χ² is obtained when relaxing the slope. However, changing the slope of the modelled profile is not sufficient to absorb the residual patterns due to the multipole. There is therefore no evidence for a degeneracy between the modelled mass density slope of the lens and the input fourth-order multipoles.

The impact of input angular structures on a fitted monopole was recently discussed by Kochanek (2021). He explains how angular structures impact the Einstein ring and, if fit with too simple a model, can either produce large χ² or drag the radially fitted model to an incorrect slope. Here we have a true lens mass profile that displays boxyness or discyness and fit it with an elliptical power law and shear. Because the multipoles only weakly change the mass within the Einstein radius, and octupole perturbations of the potential cannot be mimicked by a combination of quadrupolar components (i.e. ellipticity and shear), the freedom allowed to the fitted slope of the power law is not sufficient to absorb the multipole effects.

3.4. Influence of a relaxed slope constraint when the mock lens follows a composite mass profile

While we observe in the previous section that relaxing the constraint on the fitted slope does not significantly improve the fit, one may wonder whether or not this is due to the fact that the same family of models, i.e. power law, is used for both the mock image creation and the fit, thus not allowing any possible degeneracy linked to the slope. To test this hypothesis, we created mock lensing galaxies with a composite mass profile comprised of a baryonic Chameleon (Dutton et al. 2011) profile and a Navarro-Frenk-White (NFW; Navarro et al. 1996) dark matter component. The Einstein radius of these new lensing galaxies is chosen to be 2″, as used for our fiducial mock systems. We parameterized our mass profiles such that the lensing galaxy is similar to those observed in massive ellipical lensing galaxies, with a dark matter fraction of 0.32 within the Einstein radius, with an axis ratio q = 0.8 constant with radius, and a slope close to isothermal at the Einstein radius. The complete set of characteristics of our mock profile is given in Appendices B.2 and B.3.

This experiment being mainly an expansion of Sect. 3.3, we display the results of this test separately in Table 2. The first row of Table 2 contains the results when no multipole is present in the composite lens mass; this is used as a fiducial case to ensure the reference χ² is close to one⁵. As already observed, when the input mass profile follows a SIE, the freedom allowed in the fitted slope only slightly improves the χ². On the other hand, we also observe that the typical χ² obtained with the power-law model is similar when the truth is a composite (Table 2, rows 2–3) or a SIE (Table 1, rows 6–7).

The ability of the power-law model to fit the composite mock system when a₄ = 0 demonstrates that the two models are degenerate. The addition of multipoles to the mock lenses leaves the radial profile and Einstein radius unchanged, such that the same region of the radial profile is probed by the host galaxy images. The similarity of the χ² obtained with the composite and power-law mock mass profiles is consistent with that interpretation. On the other hand, the similarity between the χ² for both the composite and SIE mock systems (with and without multipoles) informs us that the role of the radial mass profile is less important for χ² than effects such as modelled source shape freedom (Sect. 3.7) or other input characteristics such as the S/N (Sect. 3.5).

We note that for more complex azimuthal structures, the choice of the radial profile may not be ignored. For example, if we consider a varying ellipticity with radius but the multipole kept constant, we do not exclude the possibility of greater interplay between the radial profile and the azimuthal structures: the imprint of the multipoles on the lensed ring gets stronger for a rounder galaxy, as seen in Sect. 3.2. As a consequence, a mock galaxy displaying varying slope and ellipticity with radius could potentially behave differently from what we have seen here under multipolar perturbations. However, this level of complexity is beyond the scope of this paper.

3.5. Influence of signal-to-noise ratio (S/N)

Summarising the S/N of images of a lensed system with a single number is a difficult task, as the point source and the extended images need to be accounted for. We therefore opted for a more qualitative criterion based on the exposure length and the brightness of the source. Our base case has an exposure time of 5400 s and a very bright source which has an intrinsic brightness of 20.5 mag for the host and 21 mag for the quasar. This translates into integrated lensed magnitudes around 17 (17) mag for the four images of the host (quasar). Considering each image separately, the typical mean lensed image magnitude is 18.5 mag. We then reduce the exposure time to 3000 s to match the exposure obtained typically with one orbit of HST. For that exposure time, we consider two intrinsically fainter sources: a first source that has a mean observed magnitude of typically 20.25 mag and another one of 21 mag. Those systems with lower S/N correspond to unlensed quasar brightnesses of respectively 22.75 mag and 23.5 mag, and are referred to as ‘bright’ and ‘faint’ source cases below. While this latter case corresponds to a particularly faint quasar image, the corresponding host brightness is sometimes observed in lensed systems (Ding et al. 2021b).

From Figs. 4 and 5 and Table 1, we see that there is a clear decrease in the prevalence of the residuals as the data quality decreases. Indeed, the highest S/N fits (our base case, in red in the table) display ubiquitous patterns in the residuals and the χ² is always above 1.2. For the 3000 s exposure time and bright source case (green in the table), we predominantly see cases for which the multipole imprint is hidden in the noise, but there are a few systems with χ² close to 1.2 with light imprints. This is no longer the case for the faint source cases, for which the χ² is consistently lower than 1.1.

The exposure time and source brightness are not the only factors that control the S/N, but they provide a reasonable observation-based proxy of the data quality. Other factors, such as the Einstein radius, the contrast between quasar and host, and so on, also influence the S/N, as discussed in Sect. 3.8.

3.6. Influence of the multipole strength

To characterise the dependence of the fit on the input multipole strength, we created mock mass profiles for four additional values of a₄ ∈ [0.01, 0.05]. Figure 1 provides a visual representation of the induced deviation from a pure elliptical profile. As the residuals for a₄ = 0.01 for our baseline setup have already been detected, we chose to vary the strength of the multipole for our second setup corresponding to an exposure time of 3000s and a bright source (see Sect. 3.5). The latter case with a₄ = 0.01 is shown above in Fig. 5 and the green line in Table 1. If we convert the multipole convention of Hao et al. (2006) (Eq. (3)) to our notation, we see that 70% of the galaxies of their sample have a₄ < 0.02 and only 4% of their sample have a₄ > 0.05 (see Fig. 2). The low multipole cases (a₄ ≤ 0.02) are therefore the most likely to be representative of a real sample of galaxies, while the case a₄ = 0.05 should be seen as a rare extreme situation that could plausibly happen in a large sample (typically 100) of time-delay lenses. We may note that boxy and discy galaxies are not evenly distributed, with 64% of the sample by Hao et al. (2006) being discy. However, the relative rate of boxyness or discyness is also found to depend on the galaxy mass or brightness (Hao et al. 2006; Kormendy et al. 2009; Mitsuda et al. 2017).

As one would expect, it is less easy to achieve a good fit as the multipole strength increases (see Table 1). More precisely, if we assume a threshold of χ² = 1.2 to distinguish between poor and good residuals, multipoles get detected above a₄ = 0.03 for discy galaxies while a₄ = 0.02 is the detection limit for boxy galaxies.

3.7. Influence of complexity in the modelled source

When modelling real systems, the extended source is hardly ever modelled as one simple light profile (e.g. Birrer et al. 2016, 2019). Most of the time, the source is either pixelised with regularisation applied (Brewer & Lewis 2008; Suyu & Halkola 2010; Nightingale & Dye 2015; Galan et al. 2021) or analytical with a main simple profile plus a set of basis functions, such as shapelets, that are designed to capture morphological host complexity beyond axisymmetry (Birrer et al. 2015, 2016, 2019; Birrer & Amara 2018; Shajib et al. 2019). Above, we model the source with a circular Sersic profile because we know that this is its true input profile. However, more flexibility in the source should be allowed, as in state-of-the-art modelling of strongly lensed quasars. Due to our choice of modelling software, we chose to add shapelets to the source (as e.g. Birrer et al. 2016, 2019; Shajib et al. 2019). We performed an extra step in the fitting procedure of all the systems presented in the previous sections, adding shapelets to the source with maximal polynomial order n_max = 5.

As shown in Table 1, the improvement when shapelets are added in the source is ubiquitous. When the fit is already good without shapelets, that is, typically cases with a₄ = 0.01 for the bright source, there is only a slight improvement of the χ². When the fit without shapelets displays large residual features, that is, typically when χ² > 1.2, the improvement introduced by shapelets is more significant, as illustrated in Fig. 6. The retrieved source is generally more structured than the input circular one but is not unphysical (e.g. the source surface brightness remains positive).

Fig. 6. Illustration of the influence of complexity in the modelled source. Top: mock image created for 3000 s exposure time and a bright source. The lens has the same characteristics as for the base case (see Sect. 2.1), but is discy with a4 = 0.04 in a fold configuration. Bottom: fit of the above image with (right) and without (left) adding shapelet reconstruction for the source during the fitting process. These results are also displayed in purple in Table 1.

The approximate degeneracy between the source morphology and the main lens potential is likely an expression of the source position transformation (SPT; Schneider & Sluse 2014; Unruh et al. 2017; Wagner 2018; Wertz et al. 2018). Additional complexity in the source can compensate for complexity in the lensing potential, which has further-reaching effects in general than the specific exploration of multipoles in this work. In practice, the complexity of the source is unknown, but it could be informed using a sophisticated prior that emulates the known morphologies of host galaxies. As an example, the dark spots in the source reconstructions for our fits may be mathematically acceptable (in the sense that they have non-negative brightness values), but the human eye can recognise them as being unusual and may suspect that they are instead absorbing the complexity of the potential. With no straightforward way to quantify and implement this human intuition (or certainty as to whether it can always be trusted), identification of whether or not the source reconstruction corresponds to the intended quantity is a difficult task.

With our tests, we find that shapelet reconstruction can absorb some of the multipole complexity. However, one could also wonder whether or not allowing even more freedom in the source could absorb any additional imprint of lens complexity. While answering this question is not the main purpose of this work, we concisely explored two different ways of allowing more complexity in the source: allowing ellipticity in the Sersic source, and allowing shapelets with a maximal polynomial order of eight instead of five (both maximal orders were investigated in Birrer et al. 2019). After testing a few ‘extreme cases’ with χ² > 1.3, we found that the ellipticity in the source cannot compensate the lensing potential complexity. The improvement of the fit thanks to this new degree of freedom is negligible. On the other hand, shapelets up to order n_max = 8 yield a substantial improvement of the χ² compared to the fits obtained with n_max = 5. Indeed, for the last row of Table 1, five of the eight mock images can be fitted with χ² < 1.2 with n_max = 8, while only two of the eight good fits can be obtained with n_max = 5. This improvement in the image residuals comes at a cost in the source plane, namely the appearance of clumpy structures in the source that may hopefully be identified as unrealistic (see Appendix A). In summary, the freedom in the source reconstruction helps to absorb some of the multipole complexity, but it seems that the associated distortion of the source may not systematically mimic structures that exist in real galaxies. Future work is needed to more thoroughly explore the constraining power and limitations of an arbitrary number of shapelets.

3.8. Other influences

Apart from the factors explored in the previous sections, there are yet other ingredients that could influence the detectability of multipoles, such as the Einstein radius, the source intrinsic shape, and the lens light contamination. We did not systematically quantify their impact, but we explain their influence qualitatively.

As the Einstein radius gets smaller, the different images of the same background source are brought closer together and can easily overlap, especially in cusp cases. The extended source images consist of fewer pixels and the quasar images, with the PSF spreading their light and covering the ring, dominate most of the flux. As the multipolar component is visible in residuals mainly thanks to the host flux, having less host pixels reduces the ability to detect imprints from multipole perturbations. This effect can even be stronger than a change of S/N as performed in Sect. 3.5. However, precisely quantifying the role of the Einstein radius remains difficult as the detectability of multipolar structures in the ring becomes much more dependent on the lensed-images configuration, source shape, PSF extension, and contrast for smaller Einstein radii.

Similarly, if the contrast increases between the extended host flux and the quasar flux (i.e. the host is fainter for a given quasar magnitude), the multipolar components are less detectable at fixed exposure time. In addition, the source radial profile also plays a role: as the host displays a shallower profile (decrease of n_sersic), the total brightness is the same but the host is more rapidly lost in the noise. This therefore decreases the total S/N in the host flux and the multipole is less detectable. Following the same lines, if a source is larger (i.e. increase of R_sersic) and keeps the same magnitude in the source plane, the host image spreads over more pixels with less flux in each, also decreasing the S/N in the ring. We did not explore the influence of the source azimuthal shape (e.g. an elliptical source). We expect the effect to be less straightforward than a change in n_sersic or R_sersic. We anticipate a complex dependence of the multipole detectability on the orientation of the source relative to the caustics.

In our analysis, we considered a transparent lens (i.e. no lens light). However, the presence of the lens light may have two different effects. First, the lens light and the source images can be blended, making it difficult to resolve the two components: part of the flux attributed to the lens might in fact come from the host images perturbed by multipoles in the lens mass. To avoid blending issues between source and lens light, it is common to use images of the same lens system in different spectral bands (Shajib et al. 2019). We did not simulate any lens light in our analysis and therefore did not test the efficiency of multiband imaging to mitigate the potential issue arising from blending either. Second, the lens light in itself could potentially help in detecting multipole presence as multipolar components in the mass are linked to those in the light. The detectability of multipole in the lens light profile is discussed in more detail in Sect. 3.9. The above list may not be exhaustive, but encompasses the main sources of effects that could affect the results presented in this work.

3.9. Detectability in real lens systems

We see in the previous sections that several factors influence the goodness of the fit of our mock images. For the highest quality data, the imprint of the multipole on the lensed images is ubiquitous. However, there are circumstances where the multipole imprint can remain hidden, such as if shear is absent and the main axes of both the ellipticity and multipole are aligned with the source position. Multipole perturbations are also more difficult to detect in galaxies with a higher ellipticity (i.e. small axis ratio). Leaving the mass profile slope to vary improves the residuals but only slightly, even when the true input mass profile slope varies with radius. The S/N on the other hand is crucial for detecting multipoles: with low S/N, the multipole contribution can be totally concealed in the noise. Of course, it also depends on the strength of the multipole contribution and on the freedom allowed in modelling the source, for example using a shapelet basis on top of an analytical model. Other factors such as the contrast, the source radial profile, or the Einstein radius also add more complexity to the problem.

Overall, the different factors that can increase or decrease the detectability of multipoles combine together and create a unique lensing system each time. It is not possible to characterise each system with a single number that would reflect its robustness to multipole detectability because it depends on so many parameters. However, it is possible to compare a system to our three different S/N levels. For our three S/N levels as defined by exposure time and source brightness, we essentially find three tiers of multipole detectability: 1. multipole will never go unnoticed in the residuals, 2. multipole can go unnoticed if it is at a low level or absorbed through freedom in the source structure, and 3. multipole will nearly always go unnoticed. However, exposure time alone is not a complete descriptor of these thresholds, because a lens with for example a lower exposure time but a brighter extended source relative to the quasar (or larger Einstein radius, as another example) will mimic a higher level of multipole detectability. Consequently, in order to accurately quantify the detectability of multipoles in a given system, the system must be precisely emulated so as to classify it in one of the described tiers.

A complementary path that can be followed to assess the role of multipoles in the model is to search for evidence of Fourier modes in the luminosity profile of the lens. Mitsuda et al. (2017) and Pasquali et al. (2006) demonstrated that multipolar perturbation of the light profiles can be measured in elliptical galaxies at redshifts similar to those of lensing galaxies. It is a priori less clear whether or not such a measurement is feasible in lensing galaxies for which the lens light is blended with the light from the lensed images. To assess the feasibility of such a measurement in lensed systems, we used F814W HST and K-band AO images of the lensed quasar QSO J1131–1231 (Sluse et al. 2003; Suyu et al. 2013; Chen et al. 2016). While the large Einstein ring (θ_E ∼ 1.8″) and bright lensing galaxy may ease the detection of multipoles in the lens, the lensed images are particularly bright with features extending down to the very inner region of the lens (see e.g. Fig. 4 of Claeskens et al. 2006). This makes this system a good test case. We performed multipole measurements on the images subtracted from the best ring model published by Suyu et al. (2013) and Chen et al. (2016) using two different techniques. First, we fitted ellipses with fourth-order Fourier components using the Python package photutils. The results obtained on HST and AO data are compatible with each other for radii 0.4 < θ < 0.9″. They suggest a mean value of a₄ < 0.02. At radius θ < 0.4″, there is an insufficient number of pixels to enable any sensible measurement. At radius θ > 0.9″, the residuals from the PSF and/or lower S/N in the lens prevented the algorithm from converging. For AO data, measurements were possible up to θ = 2″, confirming a₄ ≤ 0.02. We used a second method consisting in fitting the full lens intensity profile with a Sersic model using galfit (Peng et al. 2010), but including a fourth-order Fourier mode to describe the azimuthal profile of the lens. While the fit of a single Sersic profile provides a poor description of the profile in the inner regions of the lens (see e.g. Claeskens et al. 2006), the amplitude a₄ = 0.01 agrees for the two data sets. However, the same is not true for ϕ₄ which is offset in the two data sets by 30 deg. We may note that these values of a₄ do not appear unrealistic from visual inspection of the images, as a₄ > 0.04 are easily detectable by eye in good S/N data, but no deviation from a perfect ellipse is apparent. In summary, this test demonstrates that a measurement of a₄ in the light of lensing galaxies may be feasible in high-resolution images, but more work is required to assess the amplitude of systematic errors in such a measurement.

Finally, we note that we did not test the inclusion of multipoles in the lens mass model. This is an alternative path that could also be explored, but is beyond the scope of this paper.

4. Impact on H₀ inference

We show in the previous section that azimuthal perturbations of the potential, if present, could remain unnoticed in lensed images. Importantly, those perturbations can be absorbed in the lens and/or source parameters, which then differ from their true values. In this section, we quantify their impact on the value of H₀ inferred from those (biased) models. While the goodness of the fit was the main topic of the previous section, we discuss here those results in terms of the H₀ inference.

4.1. At a single-galaxy level

During our systematic tests shown in Table 1, the inferred H₀ values range typically between 60 and 80 km s⁻¹ Mpc⁻¹ while the fiducial value is 70 km s⁻¹ Mpc⁻¹. We disregard the results obtained for models yielding a poor χ², as we consider that those models would not be included in a subsequent analysis. We also disregard the tests from Table 2 because the results are redundant with those of Sect. 3.3 already present in Table 1. Figure 7 summarises the inferred H₀ values obtained for χ² < 1.2. All the cases listed in Table 1 fulfilling this criterion are included but we exclude duplicates. Specifically, we consider only once the coloured rows that appear multiple times in the table and only the shapelet reconstruction for the cases with 3000 s exposure time and a bright source. In our sample, the discy and boxy galaxies do not behave in the same way: the discy galaxy models display values of H₀ of around 70 km s⁻¹ Mpc⁻¹ and above; while the boxy ones show H₀ values with an extended tail below the fiducial H₀. If we consider only fits with χ² < 1.1, the same trends are observed. As our sample is neither homogeneous nor highly representative of a galaxy population, Fig. 7 is not the distribution of H₀ we would have derived if the value of a₄ of the lensing galaxies were randomly drawn from the distribution observed in a population of elliptical galaxies. However, each H₀ inference is representative of a possible elliptical galaxy.

Fig. 7. Median H0 inference for mock systems from Table 1 with χ2 < 1.2 (see text for more information). The histogram for boxy galxies is displayed in blue, that for discy ones is in green, and the sum of both is in red. The histograms are representative of the retrieved H0 values obtained in our systematic tests that encompass multipole amplitudes of plausible lenses. However, this does not represent an inferred H0 distribution from a realistic population of lenses (see Fig. 9).

We may now wish to quantify the impact of multipoles on H₀ when considering a population of lens systems. For this purpose, we first need to quantify the probability of having lensing galaxies deviating from a pure elliptical profile. For this purpose, we need first to identify if there is any evidence for a dependence on a₄ on the galaxy luminosity, mass, environment, and/or whether a₄ varies as function of the galacto-centric radius and of the redshift. We look at those various aspects in the following (Sects. 4.2 and 4.3).

4.2. The a₄ distribution for elliptical galaxies

Hao et al. (2006) studied the light of 847 galaxies at redshift z < 0.05. They found that boxy galaxies tend to be bigger and brighter than discy ones and that discy galaxies are more abundant, with abundance ratios depending on the magnitude cutoff. The boxy galaxies are also more likely to be found in denser environments compared to discy ones. This corroborates the findings of Bender et al. (1988, 1989) in several ways, except that the latter found equivalent abundances of boxy and discy galaxies. A search for a redshift evolution of a₄ was performed by Pasquali et al. (2006) and Mitsuda et al. (2017). Pasquali et al. (2006) studied a small sample of 18 galaxies with 0.5 < z < 1.0 and found more discy galaxies but suggested that this might be due to their limited sample. They corroborate previous findings that discy galaxies have higher ellipticities while boxy galaxies are more luminous, wider, and less elliptical in general. Pasquali et al. (2006) also found that the galaxies in their sample are similar to the ones at z = 0 and only a possible trend in a₃ with redshift was observed. Mitsuda et al. (2017) extended this work by studying 130 galaxies at z ∼ 1 and 355 at z ∼ 0. These authors found, independently of the redshift, that galaxies showing a₄ contributions are very common and that discy shapes are favoured for galaxies with log(M_*/M_⊙) < 11.5, while the boxy galaxies start to dominate at higher masses. Mitsuda et al. (2017) also looked at the a₄ radial profile of the disc-shaped NGC 4697 in more detail: in the innermost region, the a₄ parameter is more uncertain and in some filters switches to a boxy profile. For the other parts of the galaxy, the a₄ is more or less constant with radius. The change of a₄ with galacto-centric radius has been studied more extensively by several authors (e.g. Rest et al. 2001; Kormendy et al. 2009; Krajnović et al. 2013). The first two authors studied the largest sample of elliptical galaxies. They find a great diversity in a₄ variation with radius. A profile can be boxy in the centre and discy in the outskirts, or be relatively constant or even display a monotonic increase or decrease. Discyness or boxyness can therefore exist at the regions encompassed by typical Einstein rings, that is mostly the region within 1−2 effective radii. In summary, at redshifts, masses, and ellipticities of galaxies considered for lensing, the presence of galaxies with nonzero boxyness or discyness is common.

As the mass is not directly observable, we need to rely on numerical simulations to quantify the presence of such structures in the mass distribution. No authors have yet reported the presence of a₄ components in the total mass distribution, that is both stars and dark matter. However, the stellar mass dominates at radii ⪅R_Ein, and the studies focusing on multipolar components in the sellar mass are thus relevant for our purpose. Penoyre et al. (2017) studied stellar mass distribution in the Illustris simulations (Vogelsberger et al. 2014) and found a few galaxies with clear boxyness and discyness but could not conclude on precise a₄ measurements due to too low spatial resolution. Frigo et al. (2019) presented ‘zoom’ simulations and studied the boxyness and discyness of their sample more thoroughly. They did not find clear boxyness in their sample. They found systematic discyness if the galaxies are viewed edge-on. When seen at random orientations, the galaxies are still systematically discy but with lower a₄ values. The AGN feedback in their simulations tends to create less azimuthal structures at z = 0 while the AGN influence is not predominant at z = 1 where both AGN and no-AGN cases presented discyness, especially for more elliptical galaxies.

In summary, we do not find any evidence for multipoles to be less prominent in galaxies characteristic of lensing galaxies, namely galaxies with high mass and z > 0.2, than in close-by galaxies. In conjunction, we see that for a single galaxy, realistic amplitudes of the multipoles can introduce a bias on H₀ (Fig. 7). In the following section, we quantify the impact of multipoles on H₀ derived from a population of lensing galaxies.

4.3. At the population level

In this section, we use our mock systems and a plausible distribution of a₄ to quantify their influence on the inferred distribution of H₀ from an ensemble of lens systems. First, we repeat the H₀ inference for values of a₄ in the range [0.00, 0.05]⁶, a fixed axis ratio q = 0.8, and the three different levels of S/N obtained by modifying the exposure time and source brightness (see Sect. 3.5). This way, we expect to cover a representative range of S/N achievable in observations. As explained in Sect. 3.9, those three S/N levels are more broadly representative of three tiers of multipole detectability (almost always detectable; sometimes detectable; almost never detectable).

We then proceed to the fit of those mock systems with an SIE+shear model, as outlined in Sect. 2.3. When the fit is not good enough, that is when χ² > 1.2, we add shapelets in the source light. Results are displayed in Fig. 8 (considering only cases with χ² < 1.2). As already seen in the previous section, poor fits are generally obtained for the highest S/N considered, such that only a few cases appear on the figure. The inferred value of H₀ tends to be underestimated when the galaxy is boxy and overestimated when discy. Moreover, the lower the S/N, the larger the spread in the H₀ inference. We note that the typical uncertainties reported in Fig. 8 being the mean of the 68% credible intervals, not only strongly depend on the S/N but also on the configurations: cross configurations yield uncertainties roughly twice as large compared to the other configurations. Thus, having more (less) good fits for crosses will lengthen (shorten) the final typical error bars.

Fig. 8. H0 inference from mock systems created with different a4 and different S/N (yellow circles: 3000 s faint source; green squares: 3000 s bright source; red triangles: 5400 s very bright source). The boxy and discy cases are in the left and right panels, respectively. The typical uncertainties representing the 0.16 and 0.84 quantiles are displayed in the legend. The fiducial H0 is 70 km s−1 Mpc−1 and is represented with the blue dashed line.

To forecast the distribution of values of H₀ for a population of discy and boxy galaxies similar to known early-type galaxies, we combine the results from Fig. 8, weighting each data point based on the distribution of a₄ derived by Hao et al. (2006) (Fig. 2). The distribution of a₄ is discretized into six bins: the bin with a₄ ≤ 0.005 contains 10% of the sample, the 0.005 < a₄ ≤ 0.015 bin contains 46%, the 0.015 < a₄ ≤ 0.025 bin takes 21%, the one with 0.025 < a₄ ≤ 0.035 contains 10%, the 0.035 < a₄ ≤ 0.045 bin contains 6%, and the 0.045 < a₄ ≤ 0.055 bin contains 3%. We neglect the rare systems with a₄ > 0.055. For each bin, the probability density function (PDF) is the average of the different PDFs associated to each point on the graph showing a₄ versus H₀ (Fig. 8) for a given type (boxy or discy) and a given a₄ value. Those PDFs are then weighted with the probability of encountering the a₄ value and by the fraction of discy and boxy galaxies assumed in the sample. The final PDF is thus

$$\begin{array}{c}\hfill \mathrm{PDF}(H_0)={\displaystyle \sum _{i=1}^6\mathrm{PDF}(i)\times w(i),}\end{array}$$(4)

where i is the index corresponding to the a₄ bin, PDF(i) is the associated PDF for each a₄ bin and type (boxy or discy), and w(i) is the weight associated to each a₄ bin multiplied by the relative fraction of boxy and discy galaxies. In Hao et al. (2006), discy galaxies represent 64% of the sample while 36% of it are more boxy. However, Mitsuda et al. (2017) show that the fraction of boxy/discy galaxies depends on the stellar mass of those galaxies. When log(M_*/M_⊙) > 11.5, boxy and discy galaxies are equally observed. The results for the two relative distributions of boxy and discy galaxies are reported in solid lines in Fig. 9 and the median and 68% credible intervals are written in Table 3 in the H₀ columns. Those distributions correspond to the underlying H₀ distributions from which samples are drawn. In other words, from a sample of n lensing galaxies, one is expected to infer n independent measurements of H₀ that would follow the observed distribution.

Fig. 9. Underlying distribution of H0 for our base case (5400 s, very bright source) (red), and for 3000 s bright (green) and faint (yellow) sources. Each panel corresponds to two different fractions of boxy and discy galaxies (left: 36% boxy, 64% discy, as reported in Hao et al. 2006; right: 50% boxy, 50% discy). The corresponding median and 68% credible interval are displayed below in solid lines. The dotted lines are the 68% credible interval for the multipole-free distribution. It can be seen in this figure, as well as in Table 3, that the presence of multipoles becomes a significant source of uncertainty for lower S/N.

Because boxy and discy galaxies have opposite impacts on the inferred value of H₀, the final distribution is overall unbiased at the subpercent level. This result is independent of the value of the S/N considered. However, by increasing the fraction of boxy galaxies, we see that the distribution is slightly driven to lower values of H₀ by less than 1%. The larger scatter on H₀ when the S/N decreases is mainly due to the larger uncertainty for a fixed value of a₄ which is stronger for lower S/N, especially for high multipolar deformations.

Beyond bias, a noise term due to multipole presence is expected. If we consider that a multipole introduces an uncertainty on the measurement of H₀ that is independent from the other sources of uncertainty, that is when no multipoles are present in the lensing galaxy, we can calculate the multipole contribution to the total uncertainty budget. Indeed, we have

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{multipole}}^2={\sigma }_{\mathrm{tot}}^2-{\sigma }_{\mathrm{no}\mathrm{multipole}}^2,\end{array}$$(5)

where σ_tot is the total uncertainty quoted in Table 3, σ_multipole is the contribution from the multipoles, and σ_{no multipole} is the uncertainty when modelling images with a multipole-free lensing galaxy.

The σ_{no multipole} is obtained by summing the PDFs of the a₄ = 0.00 case with the different configurations weighted by the occurrence of each lens configuration in the calculation of the total PDFs (in Table 3 or Fig. 9). As a cross would have an uncertainty larger than a cusp or a fold, it is important to weight the configurations to calculate a multipole-free PDF that mimics the one with multipoles because, for example, for the highest S/N, not all configurations yield a sufficiently good fit. We find, for equal occurrence of boxyness and discyness (right panel of Fig. 9 and Table 3), that the multipoles introduce around 2% random uncertainty on the H₀ distribution (i.e. σ_multipole/H₀ = 0.02) for the highest S/N considered, and above 5% for the lower S/N. If instead we consider a distribution of multipole values similar to Hao et al. (2006) (left panel of Fig. 9 and Table 3), the relative uncertainty on H₀ due to multipoles drops below 2% in the highest S/N case, and is between 4.5% and 5% for the two other S/N cases. In the lower S/N cases, multipoles constitute the dominant source of uncertainty on H₀ (i.e. σ_multipole > σ_{no multipole}).

We note that we did not weight our sample by the occurrence of the different configurations in real lensing systems. By systematically trying to fit the four configurations, that is the two cusps, the fold, and the cross, we aim at encompassing all possibilities. However, in real cases, some configurations might be more likely than others, and those occurrences may also change from one sample to another. We did not take into account this effect but we are confident that the effect is negligible in terms of systematic bias at a population level. However, it may modify the exact value of the uncertainty on H₀.

4.4. Impact for the TDCOSMSO analysis

The TDCOSMO/H0LiCOW Collaboration combined the results of seven lenses to make their H₀ inference (Wong et al. 2020; Birrer et al. 2020). Those seven lenses were imaged in one or more spectral bands, usually including the F160W/WFC3 HST band with exposure times varying between 2000 s and 26000 s. Those NIR images are a priori the most favourable for analysing extended lensed features (see discussion in Sect. 2.2), which are the most prone to deformation from boxyness or discyness. Their Einstein radii range between 0.8″ and 2″ (Birrer et al. 2020). The host galaxies are sometimes significantly brighter than the quasar (e.g. the contrast in RXJ1131, whose host is brighter than the quasar by 0.9 mag, which is stronger than what we considered in this study, i.e. a 0.5 mag difference) and sometimes the quasar dominates the source flux (e.g. WFI2033 with its quasar being 1.5 mag brighter than the host; Ding et al. 2021b). Overall, the TDCOSMO systems are generally very bright, with very high-quality data. For most of the TDCOSMO lenses, we expect large-amplitude multipoles to be detectable if present, and thus to have an unbiased H₀ inference already at a single-galaxy level. However, a detailed analysis for each lens is needed to assess the exact detectability of multipoles and its potential impact on the H₀ inference for a given lens. Doing so is possible by creating a mock of each system with increasing values of a₄, as done in Sect. 2, to see the multipole threshold at which patterns in the residuals are detectable. By combining this with an analysis of the light of the lensing galaxy, it is possible to assess a range of plausible values of a₄, and to quantify how they may impact H₀ inference.

As a population, the seven lenses should yield an unbiased combined inference of H₀: their sample is drawn from distributions centred on the true H₀ at a subpercent level (see Fig. 9 and Table 3). Indeed, considering the TDCOSMO data quality, we conceive that for most lenses, the multipoles are detectable (typically, red or sometimes green curve in Fig. 9). Hypothetically, assuming the worst-case scenario –that the multipole were not ever detectable (yellow curve in Fig. 9)– as long as the lenses are typical of the elliptical population (i.e. with balanced quantity of boxy and discy galaxies and with a typical distribution of a₄ as in Hao et al. 2006), the bias on H₀ would still be below the percent level (albeit with increased scatter).

We note that in the TDCOSMO sample, the different single-galaxy inferences are combined by multiplying their PDFs to find the underlying true H₀ and the error on this measurement. The distributions we present in this work (Fig. 8) are the underlying distributions of H₀ from which a sample will be drawn. Thus, the uncertainties reported in Table 3 are the underlying distribution width, not the error one would get by combining x independent H₀ inferences from x lenses (the combined inference being increasingly precise as more lenses are included). The uncertainties on the combined H₀ inference quoted in the TDCOSMO analysis (Wong et al. 2020; Birrer et al. 2020) already take into account any broadening due to the possible presence of multipoles.

5. Conclusion

State-of-the-art lens modelling generally assumes that the projected mass density of lensing galaxies follows a perfect ellipse or a combination of elliptical baryonic and dark components. However, boxy or discy isophotal profiles, also referred to as fourth-order multipoles, are observed in massive elliptical galaxies. While the amplitude of those multipoles is unknown for the total mass distribution, lensing galaxies are baryon-dominated in their inner regions. It therefore seems reasonable to assume that the results of studies investigating boxyness and discyness in the light of galaxies are also applicable to the mass of lensing galaxies.

In this study, we investigated how deviations in galaxy morphology from a pure ellipse modify extended lensed images and become detectable when considering full lens image information. For this purpose, we built an ensemble of mock lensed images of a quasar and its host galaxy quadruply imaged due to strong lensing by a galaxy that is almost a singular isothermal ellipsoid (SIE): to mimic discy or boxy mass density distributions, we perturbed this density profile by adding fourth-order multipoles with amplitudes matching those detected in the light. Image quality typical of the one achieved with the WFC3 camera on board the Hubble Space Telescope was considered. We then fitted the mock images without any multipolar components.

As our base case, we create a mock image of a very bright QSO+host source with a long exposure time similar to what has been achieved for time-delay cosmography lensed systems. The source we consider is a circular Sersic lensed by an SIE with an Einstein radius of 2″, with moderate ellipticity (q = 0.8) and small multipole perturbation a₄ = 0.01. This multipole amplitude is barely detectable by eye, and is typical of what is generally measured in the light profile of elliptical galaxies (a₄ < 0.02 is expected for 70% of the systems according to Hao et al. 2006).

The fitting of this base case showed substantial residuals, suggesting that multipoles, when present, leave ubiquitous imprints in extended strongly lensed images. However, we show that various parameters may influence that detection:

• When the SIE, the multipole, the source position, and the shear are aligned, the imprint of the multipoles on the ring can be absorbed in the modelling by a combination of ellipticity and shear.

• Rounder lensing galaxies will create images that are more sensitive to the presence of multipoles.

• The freedom allowed in the lens mass model with a varying power-law slope only slightly decreases the patterns present in the residuals, independently of the family of models used for the input mass profile.

• Allowing for a structured source luminosity by means of shapelets can absorb substantial residuals coming from the multipoles.

• Naturally, the S/N plays a role: the lower the S/N, the more easily patterns in the residuals are hidden.

• Other parameters, such as the Einstein radius and the host–quasar brightness contrast, may decrease(increase) the detectability of multipoles depending on the host–quasar image blending.

We looked at the impact of the multipole on the H₀ retrieval when the model residual is compatible with the noise. We see that boxy and discy galaxies do not have the same impact on H₀: the boxy galaxies tend to bias H₀ low, while the discy galaxies mostly bias H₀ high. In both cases, H₀ can be biased by several percent for a given system. However, if we consider a sample of lenses with multipolar contributions drawn in a way that is representative of what is observed in the light of galaxies, high multipolar contributions are less probable and the mix of both boxy and discy galaxies tends to reduce the bias to a level of below 1%. However, we note that a specific selection of lenses may bias the final H₀ inference if for example only boxy galaxies are used. The S/N of the images will mainly broaden the final H₀ distribution.

On a lens-by-lens basis, an upper limit on the bias on H₀ cannot be given without performing a simulation that thoroughly emulates the specific characteristics of a given lens. For a given modelled lensing system meant to be used for cosmological inference on its own, this can be performed by following the methodology presented here, creating a mock of the system and adding multipoles to it in order to assess the threshold at which multipoles are detectable in the residuals. In addition, one can also look at the lensing galaxy light profile. Multipoles in the light can in principle be used as a proxy to the presence of multipoles in the total mass, with the caveat that the presence of multipoles in the dark matter is unsettled.

However, instead of looking at the bias for each system individually, one can consider a population approach. If using a wide sample of lensing galaxies displaying multipolar distribution representative of what is observed in the light of elliptical galaxies, we show that the combined H₀ should not be biased by the presence of multipole structure. This is also true for the TDCOSMO data. Our results suggest that thanks to the high quality of those data, and the ability of the models to reproduce them down to the noise, the bias on the combined H₀ inference from the seven lenses of TDCOSMO is expected to be below the 1% level. Furthermore, the uncertainties quoted in published analyses already take into account any slight broadening caused by to multipolar components.

Beyond the cosmology inference, this study also shows that, with high-S/N data, the multipoles have an imprint on the image of the extended source. Those observational constraints may be instrumental to better quantifying the contribution of multipoles to flux-ratio anomalies in lensed quasars (Mao & Schneider 1998; Kochanek & Dalal 2004; Gilman et al. 2019). This work might therefore help in extending existing studies of multipolar-like azimuthal perturbation as the source of flux ratio anomalies (e.g. Keeton et al. 2003; Xu et al. 2015) by calculating the effective contribution of the multipoles to the flux-ratio anomalies based on the perturbations seen in the ring. However, the accuracy with which multipoles can be recovered from the lens models has not yet been tested, and is left for future work.

Acknowledgments

The authors thank Daniel Gilman, Martin Millon, Veronica Motta, Sampath Mukherjee, Anna Pasquali, Sherry Suyu, Tommaso Treu, Jenny Wagner, Liliya L.R. Williams, and Kenneth Wong for their useful remarks on this work. This work makes use of lenstronomy v1.5.1 (Birrer & Amara 2018) and of the following Python packages: Python v3.6.5 (Oliphant 2007; Millman & Aivazis 2011), Astropy (Astropy Collaboration 2013, 2018), Numpy (van der Walt et al. 2011), Scipy (Virtanen et al. 2020), Matplotlib (Hunter 2007), Pandas (McKinney 2010; The pandas development team 2020), Photutils (Bradley et al. 2020). This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (COSMICLENS: grant agreement No. 787886). This programme is supported by the Swiss NationalScience Foundation (SNSF). D.D.X. acknowledges Tsinghua University Initiative Scientific Research Program ID 2019Z07L02017 and the Chinese Academy of Sciences (CAS) project no. 114A11KYSB20170054. G.V. has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodovska-Curie grant agreement No. 897124.

References

Appendix A: Source reconstruction with shapelets

In Sect. 3.7, we model mock images with additional freedom in the source shape model. We present in Fig. A.1 the source reconstruction obtained when modelling mock systems created with a₄ = 0.05 (last row of Table 1) with a source more complex than a smooth circular Sersic. For this purpose, we added shapelets to the source light distribution. We present the results obtained with two different maximal polynomial orders n_max of the shapelets. We see that more complex clumpy structures appear in the source with the increase of n_max. Those complex structures create lensed images that mimic those produced by multipoles, sometimes even achieving a reasonably good fit. The change in the source morphology may however hardly mimic structures seen in real galaxies. Models that predict such a source morphology would probably be considered as a symptom of a shortcoming in the macro-model of the lensing galaxy. We suggest that in the future, a machine-learning-based strategy could help to identify such features as artefacts caused by an incorrect macro model in an objective manner.

Fig. A.1. Shapelet reconstruction of the source for four mock systems with a4 = 0.05, 3000 s exposure time and a bright source. The displayed configurations are either fold configuration (left) or cusp_s configuration (right). The input lensing galaxy is either boxy (top) or discy (bottom). For each mock system, we consider two fitting setups: including shapelets in the source model with maximal polynomial order nmax = 5 or nmax = 8. The associated residuals and source shape reconstruction are displayed. For good models (i.e. χ2 < 1.2), the morphology of the source gets less realistic as nmax increases.

Appendix B: Other mass profiles

In Sect. 3.3–3.4, mass models deviating from the base setup are used. A power-law model is used in the fitting procedure in Sect. 3.3–3.4. In addition, a composite model made of a baryonic component represented by a Chameleon profile and a dark matter component emulated with a NFW profile are used to create the mock lensing galaxy in Sect. 3.4. Hereafter, we specify the equations associated to those mass profiles.

B.1. Power-law model

The convergence of an elliptical power-law profile is given by:

$$\begin{array}{c}\hfill \kappa ({\theta }_1,{\theta }_2)=\frac{3-{\gamma }^{\prime }}{2}{\left(\frac{{\theta }_{\mathrm{E}}}{\sqrt{q{\theta }_1^2+{\theta }_2^2/q}}\right)}^{{\gamma }^{\prime }-1},\end{array}$$(B.1)

where (θ₁, θ₂) are the coordinates in a reference frame oriented along the main axis of the power-law ellipsoid, γ′ is the slope of the profile (γ′ = 2 corresponds to the isothermal case), θ_E is the Einstein radius, and q is the axis ratio.

B.2. Chameleon model

The Chameleon mass model (Dutton et al. 2011; Suyu et al. 2014, Gomer et al. in prep.) is a profile with finite mass equal to the difference between two non-singular isothermal ellipsoids (NIEs) with the same normalisation but different core radii. It also mimics a Sersic profile for Sersic indexes roughly between 1 and 4. This profile is mostly used to represent baryonic matter in mass as well as in light. One formulation of convergence of the Chameleon profile is the following:

$$\begin{array}{cc}\hfill \kappa ({\theta }_1,{\theta }_2)=& \frac{L_0}{1+q}(\frac{1}{\sqrt{{\theta }_1^2+{\theta }_2^2/q^2+4w_c^2/{(1+q)}^2}}\hfill \\ \hfill & -\frac{1}{\sqrt{{\theta }_1^2+{\theta }_2^2/q^2+4w_t^2/{(1+q)}^2}}),\hfill \end{array}$$(B.2)

where (θ₁, θ₂) are the coordinates in a reference frame oriented along the main axis of the Chameleon ellipsoid, L₀ is a normalisation factor, q is the axis ratio, w_c is a proxy for the core radius of the first NIE, and w_t is a proxy for that of the second NIE.

In Sect. 3.4, we use the following values: q = 0.8, w_c = 0.024, w_t = 2.8, and L₀ such that the deflection angle at 1″ (if the Chameleon were circular) is 1.6″. Such a Chameleon profile mimics a Sersic profile of index n = 4 and effective radius R_eff = 2.3″. The radial profile of this Chameleon is presented with the dotted orange curve in Fig. B.1.

Fig. B.1. Radial convergence profile of the composite model used in Sect. 3.4. The baryonic component emulated by a Chameleon profile (dotted orange), the dark matter component rendered with a NFW profile (dashed green), and the sum of both components (dash-dotted blue) are represented. The vertical solid black line indicates the Einstein radius of the total composite mass profile.

B.3. NFW model

The Navarro-Frenk-White (Navarro et al. 1996) mass profile is the most commonly used profile to represent dark matter halos. The three-dimensional density of such a circular profile is given by:

$$\begin{array}{c}\hfill \rho (r)=\frac{{\rho }_0}{(r/r_s){(1+r/r_s)}^2},\end{array}$$(B.3)

where r is the radius, ρ₀ the central density, and r_s the scale radius. We refer to Golse & Kneib (2002) for a more complicated formulation of such a profile with ellipticity introduced in the lensing potential, as used in our simulations with lenstronomy.

In Sect. 3.4, we use the following values: ρ₀ = 2.3 × 10¹⁵ M_⊙/Mpc³, R_s = 0.071 Mpc, and the ellipticity is chosen such that the axis ratio in equal to 0.8 in convergence. Physically, this profile has a virial radius R₂₀₀ = 0.42 Mpc and a virial mass M₂₀₀ = 1.1 × 10¹³ M_⊙. The radial profile of this elliptical NFW is drawn in dashed green at Fig. B.1.

All Tables

All Figures

	Fig. 1. Isodensity contours for an SIE mass profile (grey dashed line) and an SIE + fourth-order multipole mass profile (black plain line). The multipole component has a strength a4 = 0.01 (top), a4 = 0.03 (middle), or a4 = 0.05 (bottom) and is oriented along the SIE main axis to create a discy profile (left), or at 45° to the SIE main axis to create a boxy profile (right).
In the text

	Fig. 2. Distribution of fourth-order multipoles in the light of local elliptical galaxies from Hao et al. (2006). The two panels display the distribution for the normalised (a4/a)conv and (b4/a)conv convention (left) and our a4 convention used in this paper (right), for θE = 2″ and q = 0.8 (see Eqs. (1) and (2)).
In the text

	Fig. 3. Illustration of the four lensing configurations used in this work. The inner caustic (bold) and associated critical curve are represented in black. The four different source positions are shown with star symbols and the corresponding image positions are drawn with triangles. Each configuration is displayed in a single colour (see legend).
In the text

	Fig. 4. Results of the base case mock images created and fitted as described in Sect. 2. Those systems are also displayed in red in Table 1. Top 2 rows: mock images made from one boxy galaxy for four configurations and residuals associated to the fit for each image with χ2 and H0 inference. Bottom 2 rows: same as in the top rows but for the discy lensing galaxy.
In the text

	Fig. 5. Same as Fig. 4, but with an exposure time of 3000 s instead of 5400 s and the unlensed quasar and extended source fainter by 1.75 mag. Those results are also displayed in green in Table 1. We note that the mask, which is defined by a threshold of twice the background noise, may sometimes have holes inside the Einstein ring due to the lack of lensed signal at those locations.
In the text

Fig. 6. Illustration of the influence of complexity in the modelled source. Top: mock image created for 3000 s exposure time and a bright source. The lens has the same characteristics as for the base case (see Sect. 2.1), but is discy with a4 = 0.04 in a fold configuration. Bottom: fit of the above image with (right) and without (left) adding shapelet reconstruction for the source during the fitting process. These results are also displayed in purple in Table 1.

In the text

Fig. 7. Median H0 inference for mock systems from Table 1 with χ2 < 1.2 (see text for more information). The histogram for boxy galxies is displayed in blue, that for discy ones is in green, and the sum of both is in red. The histograms are representative of the retrieved H0 values obtained in our systematic tests that encompass multipole amplitudes of plausible lenses. However, this does not represent an inferred H0 distribution from a realistic population of lenses (see Fig. 9).

In the text

Fig. 8. H0 inference from mock systems created with different a4 and different S/N (yellow circles: 3000 s faint source; green squares: 3000 s bright source; red triangles: 5400 s very bright source). The boxy and discy cases are in the left and right panels, respectively. The typical uncertainties representing the 0.16 and 0.84 quantiles are displayed in the legend. The fiducial H0 is 70 km s−1 Mpc−1 and is represented with the blue dashed line.

In the text

Fig. 9. Underlying distribution of H0 for our base case (5400 s, very bright source) (red), and for 3000 s bright (green) and faint (yellow) sources. Each panel corresponds to two different fractions of boxy and discy galaxies (left: 36% boxy, 64% discy, as reported in Hao et al. 2006; right: 50% boxy, 50% discy). The corresponding median and 68% credible interval are displayed below in solid lines. The dotted lines are the 68% credible interval for the multipole-free distribution. It can be seen in this figure, as well as in Table 3, that the presence of multipoles becomes a significant source of uncertainty for lower S/N.

In the text

Fig. A.1. Shapelet reconstruction of the source for four mock systems with a4 = 0.05, 3000 s exposure time and a bright source. The displayed configurations are either fold configuration (left) or cusp_s configuration (right). The input lensing galaxy is either boxy (top) or discy (bottom). For each mock system, we consider two fitting setups: including shapelets in the source model with maximal polynomial order nmax = 5 or nmax = 8. The associated residuals and source shape reconstruction are displayed. For good models (i.e. χ2 < 1.2), the morphology of the source gets less realistic as nmax increases.

In the text

	Fig. B.1. Radial convergence profile of the composite model used in Sect. 3.4. The baryonic component emulated by a Chameleon profile (dotted orange), the dark matter component rendered with a NFW profile (dashed green), and the sum of both components (dash-dotted blue) are represented. The vertical solid black line indicates the Einstein radius of the total composite mass profile.
In the text