Galaxy mass profiles from strong lensing I: The circular power-law model

TL;DR

This paper derives a simple, analytic method to determine the slope of galaxy mass profiles in strong lensing systems using only flux and position ratios, with high precision.

Contribution

It introduces an analytic expression linking observable ratios to the mass profile slope, independent of source structure or Einstein angle, and validates it with mock data.

Findings

The slope $\gamma$ is uniquely determined by flux and position ratios.

The uncertainty in $\gamma$ depends only on image position ratio and S/N.

Mock observations agree with theoretical predictions.

Abstract

In this series of papers we develop a formalism for constraining mass profiles in strong gravitational lenses with extended images, using fluxes in addition to positional information. We start in this paper with a circular power-law profile and show that the slope $\gamma$ is uniquely determined by only two observables: the flux ratio $f_1/f_2$ and the image position ratio $\theta_1/\theta_2$ of the two images. We derive an analytic expression relating these two observables to the slope, a result which does not depend on the Einstein angle or the structure or brightness of the source. We then find an expression for the uncertainty on the slope $\sigma_\gamma$ that depends only on the position ratio $\theta_1/\theta_2$ and the total S/N in the images. For example, in a system with position ratio $\theta_1/\theta_2=0.5$, S/N $=100$ and $\gamma=2$ we find that $\gamma$ is constrained to a precision of $\pm0.03$. We then test these results against a series of mock observations. We invert the images and fit an 11 parameter model, including ellipticity and position angle for both lens and source and measure the uncertainty on $\gamma$. We find agreement with the theoretical estimate for all mock observations. In future papers we will examine the radial range of the galaxy over which the constraint on the slope applies, and extend the analysis to elliptical lenses.

Tables

Table 1. Table 1: The parameters, and their priors, for the synthetic observations. The priors on the ellipticity parameters are equivalent to restricting the lens or source position angle to the range 0 ≤ ϕ < π 0 italic-ϕ 𝜋 0\leq\phi<\pi .

	Parameter	Model	Prior limits
	$b$ arcsec	$1$
Lens	$\gamma$	2
	${\epsilon }_x$	0
	${\epsilon }_y$	0
	$\{{\beta }_{S_x},{\beta }_{S_y}\}$ arcsec	$\{0.0-0.8,0.0\}$
Source	$I_0$	$1$
	$n_s$	$2$
	$r_{\mathrm{eff}}$	0.1, 0.2
	${\epsilon }_{S_x}$	0
	${\epsilon }_{S_y}$	0
	S	$100$

Equations

$$\rho (r)\propto r^{-\gamma },$$

$$\Sigma (R)\propto R^{1-\gamma },$$

$$\kappa (\theta )=\frac{\Sigma (D_{\mathrm{d}}\theta )}{{\Sigma }_{\mathrm{crit}}},$$

$${\Sigma }_{\mathrm{crit}}=\frac{c^2}{4\pi G}\frac{D_{\mathrm{s}}}{D_{\mathrm{ds}}{D}_{\mathrm{d}}},$$

$$M(\theta )={\Sigma }_{\mathrm{crit}}{D}_{\mathrm{d}}^2{\int }_0^{\theta }\kappa ({\theta }^{\prime })2\pi {\theta }^{\prime }\text{differential}{\theta }^{\prime }.$$

$$M({\theta }_{\mathrm{E}})={\Sigma }_{\mathrm{crit}}\pi {\theta }_{\mathrm{E}}^2{D}_{\mathrm{d}}^2,$$

$$\kappa (\theta )=\frac{3-\gamma }{2}{\left(\frac{{\theta }_{\mathrm{E}}}{\theta }\right)}^{\gamma -1}.$$

$$\left|\theta -\alpha (\theta )\right|-\beta =0,$$

$$\alpha (\theta )={\theta }_{\mathrm{E}}{\left(\frac{{\theta }_{\mathrm{E}}}{\theta }\right)}^{\gamma -2}.$$

$${\theta }_{\mathrm{E}}^{\gamma -1}=\frac{{\theta }_1+{\theta }_2}{{\theta }_1^{2-\gamma }+{\theta }_2^{2-\gamma }}.$$

$$\beta =\frac{{\theta }_2^{\gamma -1}-{\theta }_1^{\gamma -1}}{{\theta }_1^{\gamma -2}+{\theta }_2^{\gamma -2}}.$$

$${\beta }^{\prime }=\frac{1-{{\theta }_{\mathrm{r}}}^{\gamma -1}}{1+{{\theta }_{\mathrm{r}}}^{\gamma -2}}.$$

$$f_{\mathrm{r}}=\left|\frac{\mu ({\theta }_1)}{\mu ({\theta }_2)}\right|,$$

$$\mu (\theta {)}^{-1}=\left[1-{\left(\frac{{\theta }_{\mathrm{E}}}{\theta }\right)}^{\gamma -1}\right]\left[1+(\gamma -2){\left(\frac{{\theta }_{\mathrm{E}}}{\theta }\right)}^{\gamma -1}\right],$$

$$\theta ={\theta }_{\mathrm{E}}(2-\gamma {)}^{1/(\gamma -1)}.$$

$$f_{\mathrm{r}}=\left|\left[\frac{1-{\theta }_{\mathrm{E}}^{\gamma -1}{\theta }_2^{1-\gamma }}{1-{\theta }_{\mathrm{E}}^{\gamma -1}{\theta }_1^{1-\gamma }}\right]\left[\frac{1+(\gamma -2){\theta }_{\mathrm{E}}^{\gamma -1}{\theta }_2^{1-\gamma }}{1+(\gamma -2){\theta }_{\mathrm{E}}^{\gamma -1}{\theta }_1^{1-\gamma }}\right]\right|.$$

$$f_{\mathrm{r}}={\theta }_{\mathrm{r}}\left|\frac{{{\theta }_{\mathrm{r}}}^{2-\gamma }+(\gamma -2){\theta }_{\mathrm{r}}+\gamma -1}{(\gamma -1){{\theta }_{\mathrm{r}}}^{2-\gamma }+(\gamma -2){{\theta }_{\mathrm{r}}}^{1-\gamma }+1}\right|.$$

$${\mu }_{\mathrm{tot}}=\frac{2(1+{\theta }_{\mathrm{r}})}{1-{\theta }_{\mathrm{r}}},$$

$$S=\frac{f_1+f_2}{\sqrt{a{f}_1+a{f}_2}}.$$

$${\sigma }_y^2=\sum _i^n{\sigma }_{x_i}^2{\left|\text{partialderivative}f{x}_i\right|}^2,$$

$${\sigma }_{f_{\mathrm{r}}}=f_{\mathrm{r}}{\left[{\left(\frac{\sqrt{a{f}_1}}{f_1}\right)}^2+{\left(\frac{\sqrt{a{f}_2}}{f_2}\right)}^2\right]}^{\frac{1}{2}}.$$

$${\sigma }_{f_{\mathrm{r}}}=\frac{1}{S}\left(1+f_{\mathrm{r}}\right)\sqrt{f_{\mathrm{r}}}.$$

$${\sigma }_{{\theta }_k}^2=\frac{r_{\mathrm{eff}}^2}{S_k^2},$$

$${\sigma }_{{\theta }_{\mathrm{r}}}={\theta }_{\mathrm{r}}{\left[{\left(\frac{r_{\mathrm{eff}}}{S_1{\theta }_1}\right)}^2+{\left(\frac{r_{\mathrm{eff}}}{S_2{\theta }_2}\right)}^2\right]}^{\frac{1}{2}}.$$

$${\sigma }_{{\theta }_{\mathrm{r}}}=\frac{1}{S}\frac{r_{\mathrm{eff}}}{{\theta }_{\mathrm{E}}}{\left[\frac{({\theta }_{\mathrm{r}}+1{)}^3({{\theta }_{\mathrm{r}}}^3+1)}{4{\theta }_{\mathrm{r}}}\right]}^{\frac{1}{2}}.$$

$${\sigma }_{\gamma }^2={\left|\text{partialderivative}\gamma f_{\mathrm{r}}\right|}^2{\sigma }_{f_{\mathrm{r}}}^2+{\left|\text{partialderivative}\gamma {\theta }_{\mathrm{r}}\right|}^2{\sigma }_{{\theta }_{\mathrm{r}}}^2.$$

$$F({\theta }_{\mathrm{r}},f_{\mathrm{r}},\gamma )=0.$$

$$\text{partialderivative}\gamma {\theta }_{\mathrm{r}}=-\frac{\text{partialderivative}\ast F{\theta }_{\mathrm{r}}}{\text{partialderivative}\ast F\gamma },\text{partialderivative}\gamma f_{\mathrm{r}}=-\frac{\text{partialderivative}\ast F{f}_{\mathrm{r}}}{\text{partialderivative}\ast F\gamma }.$$

$${\sigma }_{\mathrm{iso}}=\frac{2}{1-{f_{\mathrm{r}}}^2}\sqrt{{\sigma }_{f_{\mathrm{r}}}^2+{\sigma }_{{\theta }_{\mathrm{r}}}^2}.$$

$$\frac{{\sigma }_{{\theta }_{\mathrm{r}}}}{{\sigma }_{f_{\mathrm{r}}}}=\frac{r_{\mathrm{eff}}}{{\theta }_{\mathrm{E}}}{\left[\frac{({{\theta }_{\mathrm{r}}}^3+1)({\theta }_{\mathrm{r}}+1)}{4{{\theta }_{\mathrm{r}}}^2}\right]}^{\frac{1}{2}},$$

Full text

††pubyear: 2019††pagerange: Galaxy mass profiles from strong lensing I: The circular power-law model–A

1 Introduction

Strong gravitational lensing, that is, when a source is multiply imaged by a lensing galaxy, can provide measurements of the masses of galaxies that are much more accurate than measurements from dynamics. Strong lensing gives, principally, the total projected mass interior to the Einstein angle $\theta_{\mathrm{E}}$ (the radius inside which the average density is high enough for multiple imaging to occur). Unfortunately strong lensing, on its own, yields only limited information on the radial profile of the mass – or at least that is the prevailing view. Kochanek (1991) made the first detailed study of the use of strong lensing to measure mass profiles, drawing pessimistic conclusions. Dynamical measurements complement the information from lensing, and, as summarised below, this combination has provided a wealth of results on the mass profiles of galaxies, and their evolution. Two thorough reviews (Kochanek, 2006; Treu, 2010) cover the use of strong lensing for measuring galaxy mass profiles.

The accurate measurement of the mass profile of the lensing galaxy is central to the determination of the Hubble constant from measured time delays (Birrer et al., 2018; Oguri & Kawano, 2003; Rusu et al., 2017; Sluse et al., 2017; Sonnenfeld, 2018; Suyu et al., 2010, 2017; Wong et al., 2017, 2018; Xu et al., 2016). Another application is in the detection of mass substructure in galaxy haloes (Hezaveh et al., 2016; Koopmans, 2005; Nierenberg et al., 2014), where it is important to determine simultaneously the best-fit smooth galaxy mass profile, against which substructure is identified, to avoid biasing the substructure results. Strong lensing mass determinations have also proved useful in studies of the stellar initial mass function (Auger et al., 2010; Brewer et al., 2014; Spiniello et al., 2011; Treu et al., 2010). The measurement of the exponent $\gamma$ of power-law fits to galaxy density profiles, $\rho(r)\propto r^{-\gamma}$, its mean value and scatter, and their evolution, has provided an exacting test for theories of the formation of galaxies (Barnabè et al., 2009, 2011; Bolton et al., 2008; Dye & Warren, 2005; Dye et al., 2008; Koopmans & Treu, 2003; Koopmans et al., 2006; Koopmans et al., 2009; Rusin et al., 2002; Spiniello et al., 2011; Treu & Koopmans, 2002; Treu et al., 2010; Courteau et al., 2014).

Most current methods to measure galaxy mass profiles of distant galaxies combine the lensing observables with dynamical information, ostensibly to break degeneracies which exist in each method alone. Pure dynamics methods suffer from a degeneracy between the measured mass profile and the radial profile of the anisotropy of the stellar orbits (e.g. Courteau et al., 2014; Koopmans & Treu, 2003). Pure lensing methods also suffer from degeneracies in determining the mass profile. The best known of these is the mass-sheet transformation (MST, Falco et al., 1985), but Schneider & Sluse (e.g. 2014) have identified a wider set of degeneracies, known collectively as the source-position transformation, of which the mass-sheet transformation is a particular example. In the original and most widely used lensing+dynamics analysis method (Treu & Koopmans, 2002; Koopmans & Treu, 2003; Koopmans et al., 2006), the only lensing information used is the total mass inside the Einstein angle, which breaks the dynamical mass-anisotropy degeneracy. This is sometimes framed the other way round, that the dynamics breaks the mass-sheet degeneracy. Barnabè & Koopmans (2007) developed a more advanced lensing+dynamics method which combines 2D velocity dispersion information with the full lensing flux information. These lensing+dynamics studies have shown that total mass profiles in early-type galaxies are close to isothermal, i.e. $\gamma\sim 2$, and such a profile slope in a spherical mass distribution gives a flat rotation curve.

Returning to pure strong-lensing analyses, lensing degeneracies may often be avoided by the use of a parametric form for the mass profile, e.g. a power law (Wagner, 2018). The results of such a parametric analysis will be useful to the extent that the model is correct. The power-law profile is a popular model for parameterising the mass profile in the central regions of galaxies, and this has been the preferred model in lensing+dynamics analyses, as well as pure dynamics studies (e.g. Chae et al., 2014) and in simulations (e.g. Wang et al., 2018; Mukherjee et al., 2018; Schaller et al., 2015; Xu et al., 2016).

Although there have been a number of pure lensing analyses that use power-law fits, and that provide good constraints on $\gamma$, rather little attention has been paid to the results. One reason we have identified that may explain this is a general distrust within the lensing community of the use of flux information to provide constraints on the mass profile. This is for historical reasons, because the majority of strong lensing theory was developed for analysing quasar images. Fits to multiply-imaged quasars only ever use the positions of the images, possibly supplemented by time delays, because the flux information is unreliable due to variability, microlensing or substructure. As a consequence the theory of the use of fluxes remains comparatively undeveloped. Without the fluxes there is too little information in images of lensed quasars to constrain the parameters of the mass profile, as demonstrated in the original analysis by Kochanek (1991). Improved constraints require positional information from multiple source systems, with $>4$ images (Cohn et al., 2001; Trotter et al., 2000), but such systems are rare.

For extended sources, where the source is a galaxy rather than a quasar, variability and microlensing may be neglected, so the flux information is useful, provided extinction in the lens may be neglected i.e. the lens may be treated as transparent. Substructure remains an issue – which is why it can be measured (see references above). The images can be highly stretched, providing multiple resolution elements, each yielding flux and position information. Analyses of mass profiles using extended sources (e.g. Rusin et al., 2002; Dye & Warren, 2005; Dye et al., 2008; Hezaveh et al., 2016; Bellagamba et al., 2017; Spingola et al., 2018; Wucknitz et al., 2004) have provided the most precise measurements of $\gamma,$ with uncertainties as small as $0.02$. This is substantially more precise than the original lensing+dynamics method and at least as precise as the updated method mentioned above. A direct comparison of the lensing+dynamics analysis of the Einstein ring 0047-2803 by Koopmans & Treu (2003), with a pure lensing analysis of the same data, was made by Dye & Warren (2005). They found that the pure lensing analysis provided much stronger constraints, and in particular (their fig. 3) it was possible to make an accurate decomposition using a two component stars+dark matter model, and to measure the stellar mass-to-light ratio, removing the strong degeneracies between the two components found in the corresponding lensing+dynamics analysis.

Given that strong gravitational lensing on its own can measure $\gamma$ as precisely as $0.02$, why then is it not the preferred method for measuring galaxy mass profiles? The reason seems to be the belief that the constraints on the slope only apply over the radial range of the images, and this is often only a small fraction of $\theta_{\mathrm{E}}$. This belief is made explicit in Kochanek (2006) which states ‘it is important to remember that the actual constraints on the density structure really only apply over the range of radii spanned by the lensed images’, while Treu (2010) states ‘It should be noted that lensing is mostly sensitive to the projected mass-density slope at the location of the images, rather than the average inside the images’. Numerous other papers contain similar statements, or make clear that the measured value of $\gamma$ quoted only applies at the Einstein angle (e.g. Chae et al., 2014; Hezaveh et al., 2016; Koopmans et al., 2006; Spingola et al., 2018; Suyu et al., 2017; van de Ven et al., 2009). The earliest example of this view that we can find is Kochanek (1995) who states ‘As with all lens systems the mass distribution is strongly constrained only over the multiply imaged region’. We have not found similar earlier quotes, so this may be the origin of the idea. Nevertheless we have not found a calculation which explicitly proves this statement. We will examine this belief closely in later papers in this series.

To summarise, flux information has been underused in the analysis of strongly lensed images of extended sources, and the theory of the constraints on the galaxy mass profile provided by such images is underdeveloped. Our goal in these papers is to produce a detailed understanding of how images of extended sources constrain the mass profile, and to develop a method to determine the lens surface mass density, together with the uncertainties, as a function of radius – in other words to determine where in the profile the constraints come from.

In this first paper we analyse the simple case of a transparent circularly symmetric power-law lens, and develop a general theory of how accurately the power-law slope $\gamma$ may be measured for any image configuration. We treat the case of an isolated galaxy. In general, and especially for high precision work, for example measuring the Hubble constant, contributions from other galaxies nearby or along the line of sight must also be accounted for. Fortunately good solutions to this problem exist: for example McCully et al. (2017) present a comprehensive treatment of this issue. A discussion of the circular power-law lens is provided by Kochanek (2006), but it is limited to consideration of positional information. The circular power-law lens is also addressed by Suyu (2012), although not in a systematic way. Treating the case of an extended source, and inverting a synthetic lensed image, she concludes that ‘the relative thickness of the arcs accurately constrains the lens mass distribution’. In their inversion of the lens SDSS J1148+1930, Bellagamba et al. (2017) draw a similar conclusion, identifying the radial magnification ratio between pairs of images as the relevant quantity (see also Sonnenfeld, 2018). Here we consider the general two-image case for the circular power-law lens. We show how the positional information on its own, or the flux (i.e. magnification) information on its own, measures only a degenerate combination of $\gamma$ and the source position $\beta$, but that the combination of positional and flux information breaks the degeneracies and provides a measurement of both quantities. In later papers we consider more complicated mass models, and extend the analysis to include elliptical lenses.

In Section 2 we present a theoretical analysis of the problem, and derive an expression that relates $\gamma$ to two observables; the ratio of the radial angle of the two images (the position ratio), and the ratio of the fluxes. We go on to derive the uncertainty on the measured value of $\gamma$ as a function of the total signal to noise ratio of the two images and the position ratio. These results are independent of the structure of the source. In Section 3 and Section 4 we create realistic artificial observations for isolated, transparent, isothermal, circular lenses over a range of source positions, and constrain the parameters of a power-law model by inverting the images. The method for creating these observations is detailed in Section 3, and the results are compared to the theoretical estimates in Section 4. Appendix A lists expressions for the correlations between the parameters in the power-law model.

2 The circular power-law lens

The mass density of a singular spherical power-law lens is

[TABLE]

where $r$ is the physical three-dimensional radial coordinate and $\gamma$ is the mass profile slope. We adopt the thin lens approximation (Schneider et al., 1992) which allows us to integrate the above into a two-dimensional surface-mass density

[TABLE]

where $R$ is now a physical two-dimensional radius perpendicular to the line of sight. This is the singular circular power-law lens.We can then scale this distribution by the distances involved and work strictly in angular coordinates as follows. We define the dimensionless surface mass density $\kappa$, also called the convergence, as

[TABLE]

where $\theta$ is the angular radial coordinate. The critical density of the system, $\Sigma_{\mathrm{crit}}$ is

[TABLE]

where $D_{\mathrm{ds}}$, $D_{\mathrm{d}}$ and $D_{\mathrm{s}}$ are the angular diameter distances from the lens (also called the deflector) to the source, from the observer to the lens and from the observer to the source respectively (e.g. Schneider et al., 1992).

The total mass $M$ enclosed by a radius $\theta$ is

[TABLE]

The Einstein angle, $\theta_{\mathrm{E}}$, is the angular radius such that

[TABLE]

or in other words, the radius inside which the average density is $\Sigma_{\mathrm{crit}}$. The appropriately normalised form of $\kappa$ is then

[TABLE]

This power-law profile includes three special cases: $\gamma=1$, $M(\theta)\propto\theta^{2}$ and the distribution is a sheet of uniform density; $\gamma=2$, $M(\theta)\propto\theta$, corresponding to an isothermal sphere; $\gamma=3$, $M(\theta)\propto 1$, equivalent to the lens being a point mass.

The lens system is illustrated in Fig. 1 and is characterised by four parameters: the mass profile slope $\gamma$, the Einstein angle $\theta_{\mathrm{E}}$, the source position $\beta$ and the source flux $f_{\mathrm{S}}$. We make the simplifying assumption that the source size is small relative to $\theta_{\mathrm{E}}$. Then we can derive image positions and magnifications without reference to the structure of the source. We also assume that the centroid of the mass distribution is known.

For $\gamma>2$ there are always two images. For $\gamma<2$ the number of images is either one or three. Third central (demagnified) images are a rarity in lensing observations (e.g. Rusin & Ma, 2001) and we do not take them into account here. In any case, the presence of a third image would only enhance the constraints on $\gamma$ as it provides more information in the image plane while also restricting the model to profiles with $\gamma<2$. Our analysis is therefore limited to the use of two images. The image plane then provides four observable quantities: the two image positions $\theta_{1}$ (inner) and $\theta_{2}$ (outer) and the flux of each image $f_{1}$ and $f_{2}$. With four observables, the four system parameters may be constrained directly.

In the following we show that a measurement of the ratios ${f_{\mathrm{r}}}=f_{1}/f_{2}$ and ${\theta_{\mathrm{r}}}=\theta_{1}/\theta_{2}$ is sufficient to determine the mass profile slope $\gamma$, independent of the other system parameters. We then proceed to derive the uncertainty on this measurement as a function of ${\theta_{\mathrm{r}}}$ and the signal to noise ratio of the observation.

2.1 Image positions

In the following the angles $\beta$ (source position), $\theta$ (image position), and $\alpha$ (reduced deflection angle) are always positive quantities. For a source at some position $\beta$, the lens forms images at $\theta_{1}$ and $\theta_{2}$ where $\theta_{1}<\theta_{\mathrm{E}}<\theta_{2}$. The positions themselves are given by the roots of the lens equation

[TABLE]

and $\alpha(\theta)$ for the power-law profile is given by

[TABLE]

Using the lens equation, and the fact that both images share the same source position, as shown by Kochanek (2006) we can write the Einstein angle as a function only of $\gamma$ and the image positions, yielding

[TABLE]

We could also have eliminated $\theta_{\mathrm{E}}$ rather than $\beta$ when combining Eqs. 8 and 9. In this case one obtains

[TABLE]

We define the reduced source position, $\beta^{\prime}=\beta/\theta_{2}$. Combining this with Eq. 11 we get an expression for $\beta^{\prime}$ in terms of $\gamma$ and ${\theta_{\mathrm{r}}}$:

[TABLE]

We choose $\beta^{\prime}=\beta/\theta_{2}$ as the definition of the reduced source position rather than the natural $\beta/\theta_{\mathrm{E}}$, because the resulting expression Eq. 12 is much simpler.

2.2 Image fluxes

The flux at a given image position is the product of the source flux and the magnification at that position. The flux ratio is then

[TABLE]

where $\mu(\theta)$ is the scalar magnification, provided below. The flux ratio is independent of the source flux. We take the absolute value of the ratio, because the magnification, defined by

[TABLE]

may be negative, depending on the image parity (which is not used). The magnification is plotted for different values of $\gamma$ in Fig. 1. The magnification diverges as $\theta\rightarrow\theta_{\mathrm{E}}$, or similarly as $\beta\rightarrow 0$, where an infinitely magnified circular image is formed at $\theta_{\mathrm{E}}$. For profiles with $\gamma<2$ the magnification also diverges at

[TABLE]

The third, central image is located inside this angular radius.

By combining Eqs. 13 and 14 we have

[TABLE]

Using Eq. 10 we can eliminate $\theta_{\mathrm{E}}$ and by substituting ${\theta_{\mathrm{r}}}$ for the position ratio we finally obtain

[TABLE]

Eq. 17 is the key equation in this paper. The equation may be numerically inverted to provide the relation between $\gamma$ and the two observables ${f_{\mathrm{r}}}$ and ${\theta_{\mathrm{r}}}$, demonstrating that a measurement of the flux ratio and position ratio for the two images provides a measurement of $\gamma$. This result is plotted in Fig. 2. We note that ${\theta_{\mathrm{r}}}={f_{\mathrm{r}}}$ yields $\gamma=2$. This is evident by inserting $\gamma=2$ into Eq. 17, and is as expected for the isothermal sphere. Since $\beta^{\prime}$ depends only on $\gamma$ and ${\theta_{\mathrm{r}}}$ (see Eq. 12) it is also uniquely determined by ${f_{\mathrm{r}}}$ and ${\theta_{\mathrm{r}}}$. Contours of $\beta^{\prime}$ are plotted in Fig. 2 as dashed lines. The use of the coordinates ${f_{\mathrm{r}}}$ and ${\theta_{\mathrm{r}}}$ has eliminated the two system parameters $\theta_{\mathrm{E}}$ and $f_{\mathrm{S}}$, which are irrelevant to the determination of $\gamma$.

It is significant that for $\gamma<2$ the contours of $\gamma$ become vertical at small values of ${\theta_{\mathrm{r}}}$. This corresponds to the point where the inner image disappears, given by Eq. 15, which is the same as setting the denominator of Eq. 17 to zero. As an example, for $\gamma=1.5$ we find this condition holds when ${\theta_{\mathrm{r}}}=3-2\sqrt{2}$.

Fig. 2 makes clear the importance of flux information for measuring the mass profile. Measurement of the image positions alone provides only a degenerate combination of $\gamma$ and $\beta^{\prime}$, equivalent to a vertical slice in Fig. 2. The addition of fluxes yields a single combination of $\gamma$ and $\beta^{\prime}$, equivalent to a point in the figure. As is also evident from the figure, the uncertainties on $\gamma$ and $\beta^{\prime}$ depend on the uncertainties on ${f_{\mathrm{r}}}$ and ${\theta_{\mathrm{r}}}$, as well as on the actual values of these two parameters, i.e. the location in the figure, since the spacing of the contours varies across the plane. Systems approaching a complete ring, where ${\theta_{\mathrm{r}}},{f_{\mathrm{r}}}\rightarrow 1$, can have any value of $\gamma$ and reproduce the same observables, leaving $\gamma$ unconstrained when $\beta^{\prime}\rightarrow 0$.

2.3 Observational uncertainties

We use a simple model for the noise in the two images, that is tractable analytically, and that is a good approximation to the actual situation in images of interest, for example Hubble Space Telescope (HST) images of the BELLS GALLERY sample of Shu et al. (2016). The main result in the current paper is the calculation of the uncertainty on the measured value of $\gamma$ as a function of three quantities: the combined S/N in the two images $S$, the position ratio ${\theta_{\mathrm{r}}}$, and the value of $\gamma$ itself. Since the results are scaled to $S$ the only requirement on the noise model is to correctly represent the ratio of the S/N in the two images.

There will be four main contributions to the uncertainty in a particular image: i) detector noise (read noise and dark current); ii) sky background counts; iii) image counts from the source; and iv) counts at the location of the image from the lensing galaxy, that will have been subtracted. For the first two terms the contribution to the variance is proportional to the size of the image, which is proportional to the flux, because lensing preserves surface brightness. The third term is also proportional to the flux. For the fourth term, although the variance is again proportional to the size of the image, the constant of proportionality will be different for the two images, because the surface brightness of the lensing galaxy will be different at the two image locations. We assume simply that the variance in each image is proportional to the flux in the image. This will be a good approximation when the image counts from the source are greater than the counts from the lensing galaxy. This is true for example for the BELLS GALLERY sample. The model will also be generally valid when the contributions from the lensing galaxy at the locations of the two images are not too different. This will be true for larger values of the position ratio, ${\theta_{\mathrm{r}}}\sim 1$. When the position ratio is small, and the lens galaxy is bright relative to the lensed source, the noise model will be less good. Configurations of this sort are less interesting in practise, because the magnifications are small, so cases with large S/N that are the most useful will be rare. The total magnification, summed for the two images, for the case $\gamma=2$ is given by

[TABLE]

and so for ${\theta_{\mathrm{r}}}\leq 0.4$, $\mu_{\mathrm{tot}}\leq 14/3$.

2.4 Uncertainty on the flux ratio

The noise model is that the variance in an individual image is proportional to the flux, i.e. $\sigma_{f}^{2}=af$ where $a$ is a constant. The summed signal to noise ratio $S$ in our two-image observation is then

[TABLE]

For a quantity $y$ which is a function of more than one parameter $y=f(x_{1},...,x_{n})$ the uncertainty on $y$, $\sigma_{y}$, is given by

[TABLE]

for small values of $\sigma_{x_{i}}$. This gives an expression for the uncertainty on the flux ratio

[TABLE]

Using Eq. 19 to eliminate $a$ gives

[TABLE]

2.5 Uncertainty on the image position ratio

Consider an image at position $\theta_{k}$. The radial size of the image is the source size $r_{\mathrm{eff}}$ multiplied by the radial magnification. For $\gamma=2$ the radial magnification is unity and at other values of $\gamma$ the correction to this is small enough to ignore here. If the image signal is normally distributed with variance $r_{\mathrm{eff}}^{2}$ around its true position and the position is estimated as the centroid of the image then

[TABLE]

where $S_{k}$ is the signal to noise ratio for that image, given by $S_{k}=f_{k}/\sqrt{af_{k}}$ (see Section 2.4). Again using Eq. 20 we add $\sigma_{\theta_{k}}$ for images 1 and 2 in quadrature

[TABLE]

Considering the case where $\gamma=2$ we can use the fact that ${\theta_{\mathrm{r}}}={f_{\mathrm{r}}}$ and $\theta_{2}=2\theta_{\mathrm{E}}/(1+{\theta_{\mathrm{r}}})$ to simplify the above and obtain

[TABLE]

The term in brackets is $\sim 2$ for all but the smallest values of ${\theta_{\mathrm{r}}}$.

2.6 Uncertainty on the slope

According to Eq. 17 and Fig. 2, the slope can be determined via a measurement of both the position ratio and the flux ratio. Using Eq. 20, the uncertainty on that measurement, $\sigma_{\gamma}$ can be estimated by

[TABLE]

Equation 17 can be rewritten as a function $F$ of ${\theta_{\mathrm{r}}}$, ${f_{\mathrm{r}}}$ and $\gamma$, such that

[TABLE]

We can then find the partial derivatives in Eq. 26 by the implicit function theorem,

[TABLE]

The exact forms of the above are, in general, complicated; however, for an isothermal lens ($\gamma=2$) they become quite simple. We use the fact that ${\theta_{\mathrm{r}}}={f_{\mathrm{r}}}$ and by evaluating the derivatives Eq. 26 becomes

[TABLE]

Now consider the relative contributions to the uncertainty from the measurements of flux ratio and position ratio. Using Eqs. 22 and 25 the ratio of the uncertainties is

[TABLE]

which is of the order $r_{\mathrm{eff}}/\theta_{\mathrm{E}}$ at all except the smallest values of ${\theta_{\mathrm{r}}}$. This shows that if the source is substantially smaller than the Einstein radius the uncertainty on the flux ratio dominates that of the position ratio and we can safely ignore the $\sigma_{\theta_{\mathrm{r}}}^{2}$ term in Eq. 29. Finally, then, for the uncertainty on the mass profile slope in a singular isothermal sphere, observed at a S/N of $S,$ we have

[TABLE]

Since galaxies are approximately isothermal, the accuracy with which $\gamma$ can be measured in a circular galaxy using strong lensing is encapsulated by this simple formula.

We have evaluated $\sigma_{\gamma}$ using Eq. 26, and ignoring the $\sigma_{\theta_{\mathrm{r}}}^{2}$ term, for a range of values of $\gamma$. From inspection of Eq. 22 and Eq. 17, it can be seen that generally $\sigma_{\gamma}\propto 1/S$, with a complicated dependence on ${\theta_{\mathrm{r}}}$ and $\gamma$, which simplifies to Eq. 31 for $\gamma=2$. The results for different $\gamma$ and $S$ are plotted in Figure 3. In the upper plot, we fix $S=100$, and plot curves of $\sigma_{\gamma}$ for different $\gamma$. In the lower plot, we fix $\gamma=2$ and vary $S$. In each panel, the dashed line corresponds to Eq. 31, for $S=100$. For an isothermal circular system with $S=100$ and $0.4<{\theta_{\mathrm{r}}}<0.8$, Eq. 31 yields $0.02<\sigma_{\gamma}<0.09$. The lower plot also includes some systems from the BELLS GALLERY (labelled in the figure) as illustrative examples. We selected four near-circular ($\varepsilon<0.15$), two-image systems and measured the S/N according to the method in Section 3.3.

The parameter $S$ separates out the dependence of $\sigma_{\gamma}$ on the quality of the observations, and reveals the dependence of $\sigma_{\gamma}$ on the lensing configuration. For a lens of particular $\gamma$, the only lens variable on which $\sigma_{\gamma}$ depends is ${\theta_{\mathrm{r}}}$, and $\sigma_{\gamma}$ decreases as ${\theta_{\mathrm{r}}}$ decreases, or, equivalently, as the image separation increases. Figure 2 illustrates the behaviour of the constraint in more detail. As ${\theta_{\mathrm{r}}}$ decreases from unity, moving R to L in the figure, two effects combine to improve the constraint: $\gamma$ becomes less sensitive to changes in the observables, illustrated by the increased spacing in the contours; and the flux ratio ${f_{\mathrm{r}}}$, which dominates the error budget, is more precisely measured, illustrated by the shrinking error bars. The result is an improving $\sigma_{\gamma}$ as the image separation increases for a constant signal to noise ratio, an effect which we will continue to see in the more complicated systems considered in later papers. To state this result another way, we can take the reciprocal of Eq. 31 and use the image positions rather than the ratio to obtain

[TABLE]

This shows that the precision to which $\gamma$ is measured is proportional to radial range $\theta_{2}-\theta_{1}$ divided by the geometric mean of the image positions (recall that $\theta_{1}$ and $\theta_{2}$ are both positive quantities), which is very similar to the Einstein angle for small image separations (and is precisely the Einstein angle when $\gamma=3$).

In a population of real lenses $S$ is naturally a function of $\beta$. Moving the source away from the lens axis produces images of lower magnification and therefore smaller $S$ for the same level of noise. It is crucial to account for this in creating mock observations: fixing $S$ ensures that changes in $\sigma_{\gamma}$ across different image configurations are a function only of the physical ability of those configurations to constrain the mass distribution, rather than the quality to which they have been observed.

In summary, measuring the position ratio and flux ratio completely determines the slope in a circular power-law lens, irrespective of the Einstein angle or the structure in the source. Using the scalar magnification, the flux ratio can be written as a function of the image positions, the Einstein angle and the slope. The lens equation can then be used to eliminate the Einstein angle, resulting in an expression for the flux ratio which depends only on the position ratio and the mass profile slope, illustrated in Fig. 2. Estimating the uncertainty on the observables then gives an analytic expression for the uncertainty on the slope (Eq. 26) which has a simple form for an isothermal lens (Eq. 31).

3 SIMULATED OBSERVATIONS

In this section we detail the parameterisation of the lens and the source and describe the procedure for creating and inverting a simulated strong lensing observation. The purpose of the simulated observations is both to verify the previous analytic results, plotted in Fig. 3, and to confirm the simulation and inversion methodology. In this way we will then be able to extend this work to more complicated problems using simulated observations alone. All the simulations use $\gamma=2$, and the noise is scaled such that the total S/N in the two images is $S=100$. If theory and simulation agree, we expect the measured uncertainties to lie along the dashed line plotted in the two panels in Fig. 3.

We define a square image plane covering $6\,\mathrm{arcsec}\,\times\,6\,\mathrm{arcsec}$ with a pixel width of $0.04$ arcsec to mimic HST WFC3 observations. The image plane position vector at the centre of the $i$th pixel is $\bm{\theta}_{i}$ and the surface brightness in this pixel is $s_{i}(\bm{x})$ where $\bm{x}$ is the vector of lens and source parameters. We assume a transparent lens. We have not included the complication of convolution with a point spread function (PSF) in the modelling. In real images the effect of the PSF must be accounted for, but there is no difficulty in principle to correct for the PSF, and the lensing results will not be impacted provided the PSF HWHM is significantly smaller than the unlensed source effective radius, and the images are well sampled.

We define a function $I(\bm{\beta},\bm{x}_{\mathrm{S}})$ that gives the source plane surface brightness for a given source plane position vector $\bm{\beta}$ and a set of source parameters $\bm{x}_{\mathrm{S}}$. The lens equation then gives the surface brightness in the image plane

[TABLE]

where $\bm{\alpha}$ is the vector deflection angle, itself a function of the lens parameters $\bm{x}_{\mathrm{L}}$. The set of all image plane surface brightnesses comprises our model,

[TABLE]

The evaluation of $\mathcal{M}$ at a given $\bm{x}$ has two parts: first find the corresponding source plane coordinates for each image pixel, from the deflection angle; second compute the source surface brightness at those coordinates.

We now detail the parameterisation of the lens and the source. In fitting, as far as possible, it is important not to impose restrictions on the models, other than the assumption of a power law profile. Therefore although both the lens and the source as circular, in fitting we parameterise each as elliptical. For real data it is best to pixelise the source to avoid biasing the results (Nightingale et al., 2018).

3.1 Mass modelling

We use a singular circular power-law galaxy for the lens, but fit a model of a singular power-law ellipsoid (SPLE, Tessore & Metcalf, 2015) with the following parameters; a lensing strength $b$, a mass profile slope $\gamma$, an ellipticity $\varepsilon_{\mathrm{L}}$ and a position angle $\phi_{\mathrm{L}}$. The ellipticity is defined as $\varepsilon=1-q$ where $q$ is the axis ratio (minor/major) of the mass distribution, and the position angle is the anti-clockwise angle from the $x$-axis to the semi-major axis of the mass distribution. The centroid of the mass distribution is assumed known. In a real observation the centroid of the mass would normally be taken as the centroid of the light of the lensing galaxy, which can usually be measured very precisely. Alternatively, the centroid could be parameterised.

The SPLE has homoeoidal, elliptical isodensity contours and its convergence is given by

[TABLE]

where $\theta_{\varepsilon}$ is an elliptical radius defined by

[TABLE]

where $i$ and $j$ are the major and minor axes respectively. In other words, $\theta_{\varepsilon}$ is the semi-minor axis of the ellipse passing through $\bm{\theta}=(\theta_{i},\theta_{j})$.

The lensing strength, $b$, is the elliptical analogue of the Einstein angle. Specifically it is the semi-minor axis of the ellipse where the average interior surface density is $\Sigma_{\mathrm{crit}}$. For a lens with lensing strength $b$ the equivalent circular lens has an Einstein angle $\theta_{\mathrm{E}}=b/\sqrt{q}$. We use the solutions in Tessore & Metcalf (2015) to efficiently calculate $\bm{\alpha}(\bm{\theta})$ for the SPLE.

In summary the mass model has four parameters, $b$, $\gamma$, $\varepsilon_{\mathrm{L}}$, and $\phi_{\mathrm{L}}$.

3.2 Source modelling

We use a circular source with a Sérsic profile to create the simulated observations, but fit an elliptical Sérsic profile for the source surface brightness in the inversion. The elliptical Sérsic profile is given by

[TABLE]

where $I_{0}$ is the surface brightness at the centre of the source, given by $\bm{\beta}_{\mathrm{S}}=\left(\beta_{\mathrm{S}_{x}},\beta_{\mathrm{S}_{y}}\right)$ and $r_{\mathrm{eff}}$ is the effective radius. $\beta_{\varepsilon}(\bm{\beta})$ is an elliptical radius in the source plane given by

[TABLE]

where $q_{\mathrm{S}}$ is the source axis ratio, defined in the same way as the lens axis ratio above, and $\beta_{\varepsilon}$ is appropriately transformed depending on the source position angle $\phi_{\mathrm{S}}$. The source ellipticity is given by $\varepsilon_{\mathrm{S}}=1-q_{\mathrm{S}}$. An expression for the constant $b_{n}$, chosen such that $r_{\mathrm{eff}}$ is the half-light radius of the galaxy, is provided by Ciotti & Bertin (1999).

To correctly sample the source function near its centre we adaptively sub-pixelise the image plane as follows. Equation 33 is evaluated a single time for each pixel. Each pixel is then assigned a level of sub-pixelisation proportional to its surface brightness in the first pass. The mean of the sub-pixel brightnesses then gives the final brightness in each pixel. Up to $100$ sub-pixels are used in the brightest pixels.

In summary the source model has seven parameters, $\beta_{\mathrm{S}_{x}}$, $\beta_{\mathrm{S}_{y}}$, $I_{0}$, $n_{s}$, $r_{\mathrm{eff}}$, $\varepsilon_{\mathrm{S}}$ and $\phi_{\mathrm{S}}$.

3.3 Addition of noise

The process outlined so far produces a model image plane $\mathcal{M}$ for a given set of input parameters $\bm{x}$. We now use this model evaluation to produce a mock observation by adding noise. The noise model, described in Section 2.3, requires that the variance in each of the images is proportional to the flux in the image. This can be achieved most simply by applying uniform noise across the image plane, since image size is proportional to flux. We would like to produce a set of simulated observations of the same $S$. Therefore we need a process for fixing $S$ in a consistent way across different simulated observations.

To do this we add noise $n_{i}$ to each pixel from a normal distribution with zero mean and variance $\sigma_{\mathrm{d}}^{2}$. The new set of surface brightness values with added noise comprise our data $\mathcal{D}$:

[TABLE]

If the $n_{i}$ are distributed around zero with variance $\sigma_{\mathrm{d}}^{2}$ then the total signal to noise ratio in the simulated observation is

[TABLE]

where the sum is restricted to only include pixels that belong to the two lensed images. To find these pixels we convolve the image plane with a Gaussian kernel. The kernel standard deviation scales with source size, and we use $\sigma=0.08$ arcsec (2 pix.) for a source with $r_{\mathrm{eff}}=0.1$ arcsec. After convolving, image pixels are chosen as those with brightness above $2\sigma_{\mathrm{d}}$. We assign a masking variable $m_{i}$ where $m_{i}=1$ for an image pixel, [math] otherwise. The level of noise $\sigma_{\mathrm{d}}$ required to reach a target signal to noise ratio $S_{\mathrm{T}}$ is then given by

[TABLE]

The values of $m_{i}$ depend on $\sigma_{\mathrm{d}}$ so the $\sigma_{\mathrm{d}}$ which sets the correct $S_{\mathrm{T}}$ must be found iteratively.

To decide on a suitable value for the total S/N in our mock observations we measured $S$ in the $2\sigma_{\mathrm{d}}$ mask for $16$ images in the BELLS GALLERY sample (Shu et al., 2016), which were found to have $50<S<500$ with a mean $S=241$. Based on this we chose $S=100$ as representative of the typical quality of images available for fitting lens models.

3.4 Parameter estimation

Recall that our data $\mathcal{D}$ is a set of pixel values $\{d_{i}\}$ where each $d_{i}=s_{i}+n_{i}$ and each $n_{i}$ is drawn from $\mathcal{N}(0,\sigma_{\mathrm{d}})$. With this assumption, the probability that we observe $d_{i}$ in a given pixel is

[TABLE]

where $s_{i}$ can be found by Eq. 33 for a given $\bm{x}$. The product of these probabilities from all $N$ pixels

[TABLE]

gives the probability of observing $\mathcal{D}$ if the true model is $\mathcal{M}$. The model is only a function of the parameters $\bm{x}$ so we can write the above as

[TABLE]

where $C$ is some constant dependent only on $\sigma_{\mathrm{d}}$ and $\chi^{2}$ has the standard definition

[TABLE]

Using Bayes’s theorem,

[TABLE]

we obtain the posterior probability of the model parameters as

[TABLE]

where $\mathrm{Pr}\!\left(\bm{x}|\mathcal{I}\right)$ is the probability of the model having parameters $\bm{x}$ before any observation, also called the ‘prior’ probability of $\bm{x}$. This is conditioned on prior information $\mathcal{I}$ which in this case is the physical range of the parameters. We adopt a prior that is uniform and proper in each parameter within the finite bounds specified by Table 1.

We then marginalise over the parameters $\beta_{\mathrm{S}_{x}}$, $\beta_{\mathrm{S}_{y}}$, $I_{0}$, $n_{s}$, $r_{\mathrm{eff}}$, $\varepsilon_{\mathrm{S}}$, $\phi_{\mathrm{S}}$, $b$, $\varepsilon_{\mathrm{L}}$, and $\phi_{\mathrm{L}}$ to find $\mathrm{Pr}\!\left(\gamma|\mathcal{D},\mathcal{I}\right)$. We define the posterior uncertainty in $\gamma$, $\sigma_{\gamma}$, as the mean distance from the median of $\mathrm{Pr}\!\left(\gamma|\mathcal{D},\mathcal{I}\right)$ to its $16$th and $84$th percentiles. If $\mathrm{Pr}\!\left(\gamma|\mathcal{D},\mathcal{I}\right)$ is normally distributed then $\sigma_{\gamma}$ is equivalent to the $1\sigma$ uncertainty. Using $\sigma_{\gamma}$ to describe the constraint on $\gamma$ in this way is a convenient but incomplete summary of the posterior which will also feature correlations between the parameters. In later parts where we refer to ‘the constraint on the profile slope’ or similar, $\sigma_{\gamma}$ as defined here is the specific quantity we are referring to.

3.5 Implementation and posterior sampling

For a given mock observation we constrain the parameters using a Markov chain Monte Carlo (MCMC) sampler. We use the affine-invariant ensemble sampler introduced by Goodman & Weare (2010) and implemented in the emcee Python package by Foreman-Mackey et al. (2013). The affine-invariance of the sampler makes it very useful for the highly correlated posterior distributions present in lensing problems. It also scales well in parallel, further improving convergence time. The sampler uses the log-posterior

[TABLE]

When sampling, we swap the parameters $\varepsilon$ and $\phi_{\mathrm{L}}$ for $\varepsilon_{x}$ and $\varepsilon_{y}$, defined such that

[TABLE]

which accounts for the unconstrained position angle when the lens is circular (Hezaveh et al., 2016). We take the same approach with the source ellipticity and position angle.

To ensure as short a burn-in time as possible we initialise the walkers in a Gaussian ball around the true parameter values which we already know for our simulated observations. We could opt to initialise the walkers in a uniform distribution across the space defined by Table 1 and achieve the same results, albeit with a longer convergence time. After an initial burn-in phase with a small number of walkers ($\sim 50$) the number of walkers is increased ($\sim 300$), using the final positions of the burn-in walkers as seed locations. Again, this is done in the interest of speed. At the end of the sampling run, walkers that did not converge to the main ensemble are removed, found automatically using a clustering algorithm. Typically $\lesssim 1$ per cent of the walkers are removed in this way.

4 Results

4.1 Summary of simulated observations

The simulated observations all have circular, isothermal lenses with $\theta_{\mathrm{E}}=1.0$ arcsec. The sources are also circular, with $n_{s}=2$, and central surface brightness (which is arbitrary) set to unity. Two sets of simulated observations were created, one with the source effective radius $r_{\mathrm{eff}}=0.1$ arcsec, and one with $r_{\mathrm{eff}}=0.2$ arcsec. For each set the source position was moved along the $x$ axis over the range $0\leq\beta_{\mathrm{S}_{x}}\leq 0.8$ arcsec in steps of 0.01 arcsec. The relation between ${\theta_{\mathrm{r}}}$, ${f_{\mathrm{r}}}$ and $\beta$ for an isothermal lens is

[TABLE]

so the simulated observations cover the range $0.11\leq{\theta_{\mathrm{r}}}\leq 1.0$. All the parameters, and their priors, are summarised in Table 1. Pictures of the models, for a range of ${\theta_{\mathrm{r}}}$ are provided in the top panel of Fig. 4, for the case of $r_{\mathrm{eff}}=0.1$ arcsec.

4.2 Constraints on the mass profile slope

In Fig. 4 we plot the measured value of $\sigma_{\gamma}$ for all the simulated observations, together with the theoretical curve according to Eq. 31. The constraints on the mass profile slope measured from the simulations match the predicted curve very closely, for both values of $r_{\mathrm{eff}}$. This result confirms that the theoretical analysis captures all that is relevant in measuring $\gamma$. The earlier analysis made no assumptions about the source structure, yet we have correctly predicted $\sigma_{\gamma}$ for extended sources. This shows that the theoretical analysis should apply equally to any source brightness distribution, and the unknown structure of the source – if properly treated (see below) – has no impact on the precision of the measurement of $\gamma$. The result also verifies the simulation methodology, including the method for measuring $S$ in the images.

The size of the source has no significant effect on the constraint as long as $r_{\mathrm{eff}}\ll\theta_{\mathrm{E}}$, as is true for both cases above. By assuming that $\sigma_{f_{\mathrm{r}}}$ dominates over $\sigma_{\theta_{\mathrm{r}}}$ in deriving Eq. 31 we removed any dependence of $\sigma_{\gamma}$ on $r_{\mathrm{eff}}$ and the data shows this assumption to be valid. In our mock observations the true (isothermal) slope is found $64.9\%$ of the time to within $1\sigma$ and $94.1\%$ of the time to within $2\sigma$, confirming that the results are unbiased. For real data a Sérsic profile will only be an approximation to the true source profile, and as noted earlier to avoid biased results a pixelised source should be used (Nightingale et al., 2018).

At ${\theta_{\mathrm{r}}}\rightarrow 1$ we should in theory see $\sigma_{\gamma}\rightarrow\infty$ as there is only one infinitely magnified image at $\theta=\theta_{\mathrm{E}}$ (see upper right of Fig. 2). However there is some evidence that the constraint for $r_{\mathrm{eff}}=0.2\,\mathrm{arcsec}$ is slightly better than the theoretical prediction for large values of ${\theta_{\mathrm{r}}}$. We attribute this to the fact that in the simulated observations the Einstein ring has a finite width and the small difference in ${\theta_{\mathrm{r}}}$ and ${f_{\mathrm{r}}}$ between pixels across the extension of the ring is enough to provide a weak constraint on $\gamma$. At the other extreme, as ${\theta_{\mathrm{r}}}\rightarrow 0$, both sets of data start to deviate from the prediction. At ${\theta_{\mathrm{r}}}<0.2$ the term in brackets in Eq. 25 begins to increase and our assumption that the position ratio uncertainty does not contribute to $\sigma_{\gamma}$ starts to break down.

4.3 Parameter constraints and correlations

Analysis of the correlations between the parameters of the fit provides further insight into the lensing properties of the singular power-law galaxy. The form of the posterior probability density functions (PDFs) is similar across the range of systems, containing a complicated set of correlations. Figure 5 shows the PDF for four parameters that display interrelated correlations. Additionally the parameters $n_{s}$ and $I_{0}$ are correlated. This is a generic feature in fitting Sérsic profiles to galaxies, and is unrelated to determination of the lensing parameters.

The correlations plotted may all be understood through the earlier theoretical analysis. We have computed the theoretical slopes of the correlations for the six panels labelled A$-$F in the figure, plotted as dashed lines, and their derivation is explained below. All the computed correlations match the results of the simulations well. The correlation between $\beta_{\mathrm{S}_{x}}$ and $\gamma$ is the degeneracy between these parameters for fixed ${\theta_{\mathrm{r}}}$, corresponding to a vertical line in Fig. 2. The degeneracy is broken by the flux information, and the uncertainty along the line is set by the uncertainty on the flux ratio. Although the correlation is extremely narrow, the axis ratio of the correlation relates to the ratio of the uncertainty on ${\theta_{\mathrm{r}}}$ to the uncertainty on ${f_{\mathrm{r}}}$. For a larger source the correlation is broader.

To derive the slope of the correlation it is sufficient to assume that the positions are measured perfectly. Then we simply differentiate Eq. 11 with respect to $\gamma$, assuming $\theta_{1}$ and $\theta_{2}$ are constant. Since the model is isothermal, we insert $\gamma=2$ into the resulting expression, and substitute $b=\theta_{\mathrm{E}}=(\theta_{1}+\theta_{2})/2$ from Eq. 10 to yield the simple relation

[TABLE]

which is plotted in panel D in Fig. 5. In a similar way we differentiate Eq. 10 with respect to $\gamma$ to derive the correlation between $\gamma$ and $b$

[TABLE]

which is plotted in panel F. Application of the chain rule leads to the correlation between $b$ and $\beta$, plotted in panel B.

In frames A, C and E in Fig. 5 the inferred radius of the source is correlated with the system parameters $\beta$, $b$, and $\gamma$. This arises because of the dependence of magnification on $\gamma$ for fixed position $\theta$, Eq. 14 . Since the inferred source radius is inversely proportional to the square root of the absolute value of the magnification $r_{\mathrm{eff}}\propto\mu^{-1/2}$, then

[TABLE]

From Eq. 14 we can compute $\differential{\mu}/\differential{\gamma}$. Combining these two relations with the derivatives of Eqs. 10 and 11 we can compute all the correlations in frames A, C and E. Where $\mu$ appears in an expression for a correlation, as in these three panels, we compute the gradient at both images and take the mean weighted by image size, proportional to image position for an isothermal lens. For the explicit route to each of the gradients used in the plot see Appendix A.

5 Conclusions

We have presented an analysis of how observations of an extended background source lensed by a singular circular power-law lens can be used to measure the slope $\gamma$ of the density distribution $\rho\propto r^{-\gamma}$ in the lens. Using just the positions of the two images provides only a degenerate combination of $\gamma$ and the source position $\beta$. The flux information breaks this degeneracy. We showed that if the lens is isothermal, with $\gamma=2$, the precision with which $\gamma$ may be measured is given by a simple analytic expression. This precision improves with decreasing image position ratio, or increasing image separation, and increasing signal to noise ratio.

The theoretical analysis was based simply on the measurement of two quantities, the position ratio and the flux ratio, and made no assumption about the source surface-brightness profile. Nevertheless the results were confirmed to apply in the case of a full inversion of lensed images for an extended source. We created synthetic observations for a circular isothermal power-law lens, and a circular source, but we fit with elliptical models for both the lens and source. We recovered unbiased values of the 11 input parameters of the lens and the source, and the results for $\sigma_{\gamma}$ as a function of ${\theta_{\mathrm{r}}}$ match the theoretical curve accurately. The fact that the results of the full inversion, fitting the surface-brightness distribution, agree with the theoretical analysis, which made no assumptions about the source surface-brightness profile, shows that fitting the surface brightness profile in a real case, using e.g. the method of Nightingale et al. (2018), utilises the full positional and flux information in determining the mass slope. In effect a complicated galaxy light profile may be considered as a set of compact sources, each contributing to the measurement of $\gamma$.

These results, especially the strong agreement between the theoretical analysis and the synthetic observations, validate our inversion method and our method for setting the total signal to noise ratio. This is crucial to the work in the rest of this series where purely analytic results will not be possible and we will have to rely on the synthetic observations alone. Having established a theory for the constraints on the profile slope in the simplest case, the rest of the series will address primarily two issues: the radial range of the galaxy over which the constraint on the slope applies, and the strength of constraint available in elliptical lenses.

Acknowledgements

We are grateful to Simon Dye for discussions on aspects of this work, to Yiping Shu for sending us images from the BELLS GALLERY sample, and to the Imperial College Research Computing Service for HPC resources and support. SJW thanks Tom Broadhurst for a conversation long ago, before CMO’R went to school, that sparked interest in this topic. CMO’R is supported by an STFC Studentship.

Appendix A Correlations in Power Law Profiles

We have six correlations between four parameters; $\beta$, $r_{\mathrm{eff}}$, $b$ and $\gamma$. In Section 2 we derived three equations which relate these parameters. Combining their derivatives can give us each of the six correlations in Fig. 5. We define the derivatives we need as follows.

From Eq. 10 with $b=\theta_{\mathrm{E}}$

[TABLE] 2. 2.

From Eq. 11

[TABLE] 3. 3.

From Eq. 14

[TABLE] 4. 4.

By assuming that $r_{\mathrm{eff}}\propto 1/\sqrt{\mu}$

[TABLE]

In each frame we then assume the correlation lies on a straight line with the gradient given by the derivative that combines those two parameters. The relevant gradients for each of the highlighted frames in Fig. 5 are then:

[TABLE]

Exact expressions for the derivatives can be found by symbolic computation software and numerically evaluated for the specific image plane configuration at hand.