Lens modelling of the strongly lensed Type Ia supernova iPTF16geu

TL;DR

This paper refines lens models of the strongly lensed Type Ia supernova iPTF16geu, constraining the lens density profile and the Hubble constant through detailed observations and analysis of fluxes, time delays, and sub-structure effects.

Contribution

It provides a detailed lens model of iPTF16geu, incorporating sub-structure and microlensing, and offers independent constraints on the density slope and Hubble constant.

Findings

Smooth lens models cannot fully explain the fluxes, indicating sub-structure is needed.

The lens density slope is constrained to around 1.8, consistent with high-resolution imaging and velocity dispersion data.

A lower bound on the Hubble constant is established at approximately 40 km/s/Mpc.

Abstract

In 2016, the first strongly lensed Type Ia supernova, iPTF16geu at redshift $z=0.409$ with four resolved images arranged symmetrically around the lens galaxy at $z=0.2163$, was discovered. Here, refined observations of iPTF16geu, including the time delay between images, are used to decrease uncertainties in the lens model, including the the slope of the projected surface density of the lens galaxy, $\Sigma\propto r^{1-\eta}$, and to constrain the universal expansion rate $H_0$. Imaging with HST provides an upper limit on the slope $\eta$, in slight tension with the steeper density profiles indicated by imaging with Keck after iPTF16geu had faded, potentially due to dust extinction not corrected for in host galaxy imaging. Since smaller $\eta$ implies larger magnifications, we take advantage of the standard candle nature of Type Ia supernovae constraining the image magnifications, to obtain an independent constraint of the slope. We find that a smooth lens density fails to explain the iPTF16geu fluxes, regardless of the slope, and additional sub-structure lensing is needed. The total probability for the smooth halo model combined with star microlensing to explain the iPTF16geu image fluxes is maximized at $12\,\%$ for $\eta\sim 1.8$, in excellent agreement with Keck high spatial resolution data, and flatter than an isothermal halo. It also agrees perfectly with independent constraints on the slope from lens velocity dispersion measurements. Combining with the observed time delays between the images, we infer a lower bound on the Hubble constant, $H_0 \gtrsim 40\,{\rm km \ s^{-1} Mpc^{-1}}$ at $68.3\,\%$ confidence level.

Tables

Table 1. Table 1: Observed image positions, arrival times and magnifications of iPTF16geu images (Dhawan et al., 2020 ) . The radial coordinate is given in arcsec, the angle ϕ italic-ϕ \phi in radians North of East and the time delay in days.

Image	$r$ [”]	$\varphi$	$t_{\mathrm{obs}}$ [days]	${\mu }_{\mathrm{obs}}$
1	$0.251\pm 0.001$	$4.468\pm 0.002$	$0\pm 0.33$	$35.6\pm 1.1$
2	$0.324\pm 0.001$	$2.679\pm 0.003$	$-0.23\pm 0.99$	$15.7\pm 1.1$
3	$0.297\pm 0.002$	$1.013\pm 0.006$	$-1.43\pm 0.64$	$7.5\pm 1.1$
4	$0.276\pm 0.001$	$5.860\pm 0.005$	$1.36\pm 1.07$	$9.1\pm 1.1$

Table 2. Table 2: Observed time delay differences, Δ ​ t i ​ j ≡ t j − t i Δ subscript 𝑡 𝑖 𝑗 subscript 𝑡 𝑗 subscript 𝑡 𝑖 \Delta t_{ij}\equiv t_{j}-t_{i} , and magnification ratios, r i ​ j ≡ μ i / μ j subscript 𝑟 𝑖 𝑗 subscript 𝜇 𝑖 subscript 𝜇 𝑗 r_{ij}\equiv\mu_{i}/\mu_{j} , between iPTF16geu images.

Images ($ij$)	$\Delta t_{ij}$ [days]	$r_{ij}$
$12$	$0.23\pm 1.04$	$2.27\pm 0.17$
$13$	$1.43\pm 0.81$	$4.74\pm 0.74$
$14$	$-1.36\pm 1.12$	$3.91\pm 0.50$

Table 3. Table 3: Stellar mass fractions, f ∗ subscript 𝑓 f_{*} , at the iPTF16geu image positions as derived from HST F ​ 814 ​ W 𝐹 814 𝑊 F814W photometry approximated by the SDSS i 𝑖 i -band and lens modelling assuming a lens mass slope η = 1.8 𝜂 1.8 \eta=1.8 . The columns correspond 5th, 25th, 50th (median), 75th and 95th percentile mass fractions.

Image	$f_{\ast }$ (5 %)	$f_{\ast }$ (25 %)	$f_{\ast }$ (50 %)	$f_{\ast }$ (75 %)	$f_{\ast }$ (95 %)
1	$0.073$	$0.130$	$0.198$	$0.271$	$0.422$
2	$0.069$	$0.122$	$0.186$	$0.256$	$0.397$
3	$0.065$	$0.115$	$0.175$	$0.241$	$0.373$
4	$0.067$	$0.118$	$0.180$	$0.247$	$0.384$

Table 4. Table 4: For each image, the convergence, κ 𝜅 \kappa , the shear, γ 𝛾 \gamma , and the (signed) magnification predicted for a smooth halo, μ th subscript 𝜇 th \mu_{\rm th} , assuming a lens mass slope of η = 1.8 𝜂 1.8 \eta=1.8 are listed, together with the one-tailed p 𝑝 p -values, derived from the microlensing PDFs shown in figure 6 using stellar mass fractions, f ∗ subscript 𝑓 f_{*} , listed in table  3 . The total probability for the observed magnifications is p tot = 0.12 subscript 𝑝 tot 0.12 p_{\rm tot}=0.12 .

Image	$\kappa$	$\gamma$	${\mu }_{\mathrm{th}}$	$p$
1	$0.644\pm 0.004$	$0.444\pm 0.004$	$-{14.2}_{-1.3}^{+1.1}$	$0.08$
2	$0.546\pm 0.004$	$0.349\pm 0.004$	${11.9}_{-0.7}^{+0.9}$	$0.21$
3	$0.653\pm 0.005$	$0.451\pm 0.005$	$-{12.0}_{-1.2}^{+1.0}$	$0.43$
4	$0.566\pm 0.003$	$0.362\pm 0.003$	${17.4}_{-1.2}^{+1.5}$	$0.22$

Table 5. Table 5: Derived parameters for the lens mass distribution of the iPTF16geu system assuming a lens mass slope of η = 1.8 𝜂 1.8 \eta=1.8 from HST F ​ 814 ​ W 𝐹 814 𝑊 F814W -band data when iPTF16geu is active and Keck J 𝐽 J -band reference imaging. The inclination angle ϕ mass subscript italic-ϕ mass \phi_{\rm mass} of the major axis is defined in observers coordinates, East of North (i.e. with respect to the y 𝑦 y -axis in a Cartesian coordinate system).

Parameter	HST F814W	Keck $J$
${\theta }_E$ [”]	$0.292\pm 0.001$	$0.291\pm 0.001$
${\varphi }_{\mathrm{mass}}$ [∘]	$65.5\pm 0.3$	$65.0\pm 0.8$
$q$	$0.88\pm 0.01$	$0.77\pm 0.01$

Equations

$$h\equiv \frac{H_0}{100km/s/Mpc}.$$

$$\Delta t=\frac{D_{\mathrm{s}}{D}_{\mathrm{l}}}{D_{\mathrm{ls}}}(1+z_{\mathrm{l}})\tau ,\tau \equiv \frac{1}{2}(\vec{\theta }-\vec{\beta }{)}^2-\Psi (\vec{\theta }),$$

$$D_{\Delta t}\equiv \frac{D_{\mathrm{s}}{D}_{\mathrm{l}}}{D_{\mathrm{ls}}}(1+z_{\mathrm{l}}),$$

$$\Psi (\vec{\theta })=2\frac{D_{\mathrm{ls}}}{D_{\mathrm{s}}{D}_{\mathrm{l}}}\int \Phi dl,$$

$$\kappa (\vec{\theta })=4\pi \frac{D_{\mathrm{l}}{D}_{\mathrm{ls}}}{D_{\mathrm{s}}}\int \rho dl=\frac{\Sigma (\vec{\theta })}{{\Sigma }_{\mathrm{crit}}},$$

$$\Sigma (\vec{\theta })\equiv \int \rho dl\mathrm{and}{\Sigma }_{\mathrm{crit}}\equiv \frac{1}{4\pi }\frac{D_{\mathrm{s}}}{D_{\mathrm{l}}{D}_{\mathrm{ls}}}.$$

$$\vec{\beta }=\vec{\theta }-\vec{\alpha }(\vec{\theta }),$$

$$\mathcal{A}(\vec{\theta })=\frac{\delta \vec{\beta }}{\delta \vec{\theta }},$$

$$\mu (\vec{\theta })=\frac{1}{\det \mathcal{A}(\vec{\theta })}=\frac{1}{(1-\kappa {)}^2-{\gamma }^2},$$

$${\gamma }_1=\frac{1}{2}({\Psi }_{11}-{\Psi }_{22}),{\gamma }_2={\Psi }_{12}={\Psi }_{21}.$$

$${\Psi }^{\prime }=\xi \Psi +\frac{1-\xi }{2}{\theta }^2+\vec{s}\cdot \vec{\theta }+u,$$

$${\vec{\alpha }}^{\prime }=\xi \alpha +(1-\xi )\vec{\theta }+\vec{s},$$

$$1-{\kappa }^{\prime }=\xi (1-\kappa ).$$

$${\mu }^{\prime }=\frac{\mu }{{\xi }^2}\mathrm{and}\Delta t^{\prime }=\xi \Delta t.$$

$$\frac{\delta {\mu }^{\prime }}{{\mu }^{\prime }}$$

$$\frac{\delta h^{\prime }}{h^{\prime }}$$

$$\kappa =\frac{b}{2\sqrt{s^2+(1-ϵ){\theta }_1^2+(1+ϵ){\theta }_2^2}},$$

$$q^2=\frac{1-ϵ}{1+ϵ}.$$

$${ϵ}_e$$

$$e$$

$$\kappa =\frac{b^{\eta -1}}{2{\left[s^2+(1-ϵ){\theta }_1^2+(1+ϵ){\theta }_2^2\right]}^{(\eta -1)/2}},$$

$$\Delta \tau$$

$$=\frac{r_{\mathrm{I}}^2-r_{\mathrm{II}}^2}{2}\left[1+\frac{2({\vec{\theta }}_{\mathrm{II}}\cdot {\vec{\alpha }}_{\mathrm{II}}-{\vec{\theta }}_{\mathrm{I}}\cdot {\vec{\alpha }}_{\mathrm{I}}+{\Psi }_{\mathrm{I}}-{\Psi }_{\mathrm{II}})}{r_{\mathrm{I}}^2-r_{\mathrm{II}}^2}\right].$$

$$\Delta {\tau }_{\mathrm{SIS}}=\frac{r_{\mathrm{I}}^2-r_{\mathrm{II}}^2}{2}.$$

$$p\equiv \frac{\Delta r}{<r>}=\frac{2(r_{\mathrm{I}}-r_{\mathrm{II}})}{r_{\mathrm{I}}+r_{\mathrm{II}}}\text{ and }\delta \equiv {\varphi }_{\mathrm{I}}+\pi -{\varphi }_{\mathrm{II}}$$

$$\Delta \tau =\Delta {\tau }_{\mathrm{SIS}}(\eta -1)\left[1-\frac{(\eta -2{)}^2}{12}{p}^2-\frac{\eta -2}{4}{\delta }^2\right].$$

$$D_{\mathrm{l}}=745.8\mathrm{Mpc},D_{\mathrm{s}}=1157.0\mathrm{Mpc}\mathrm{and}{D}_{\mathrm{ls}}=513.2\mathrm{Mpc}.$$

$${\mu }_{\mathrm{central}}\approx \frac{1}{(1-\kappa {)}^2}.$$

$$\kappa >1+\frac{1}{\sqrt{{\mu }_{\max }}}\to {\Sigma }_{\mathrm{central}}>{\Sigma }_{\mathrm{crit}}\left(1+\frac{1}{\sqrt{{\mu }_{\max }}}\right),$$

$$I(R)=I_e\exp \left[-b_n{\left(\frac{R}{R_e}\right)}^{1/n}-1\right],$$

$$R\equiv \sqrt{x_1^2+x_2^2/q_S^2},$$

$$\frac{d}{dr}(\nu {\sigma }_r^2)+\frac{2\beta }{r}\nu {\sigma }_t^2=-\nu \frac{d\Phi }{dr}\equiv -\nu {\Phi }^{\prime },$$

Full text

††pubyear: 2019††pagerange: Lens modelling of the strongly lensed Type Ia supernova iPTF16geu–Lens modelling of the strongly lensed Type Ia supernova iPTF16geu

1 Introduction

The expansion history of the Universe can be constrained by measuring redshifts and distances of standard candles such as Type Ia supernovae (SNe Ia) (Goobar & Leibundgut, 2011). The redshift gives the growth since the time when the light was emitted. This time is measured, together with the spatial curvature of the Universe, by the distance to the SNe Ia as inferred from its apparent magnitude.

Weak gravitational lensing from inhomogeneities in the matter distribution will cause a scatter in the distance measurements, possibly degrading the accuracy of the cosmological parameters derived from the expansion history. The first tentative detection of the gravitational magnification of SNe Ia was made in Jönsson et al. (2007), see also Mörtsell et al. (2001b); Jönsson et al. (2006); Nordin et al. (2014); Rodney et al. (2015); Rubin et al. (2018). In principle, the effect can be corrected for by cross-correlating the SNe Ia observations with data on the foreground galaxies responsible for the lensing effect (Amanullah et al., 2003; Jönsson et al., 2006, 2008, 2009). The induced scatter can also be used to measure the masses of the foreground galaxies (Jönsson et al., 2010) and to constrain the fraction of matter inhomogeneities in compact objects (Rauch, 1991; Metcalf & Silk, 1999; Seljak & Holz, 1999; Goliath & Mörtsell, 2000; Mörtsell et al., 2001a; Mörtsell, 2002; Zumalacarregui & Seljak, 2018), see also Dhawan et al. (2018).

In this paper, we study the first resolved strongly lensed SNe Ia. It is well known that the time delay between images in strong gravitational lensing systems can be used to constrain the Hubble constant, $H_{0}$, (Refsdal, 1964), see also Goobar et al. (2002); Mörtsell et al. (2005); Mörtsell & Sunesson (2006). SNe Ia are especially useful in this respect since (Kolatt & Bartelmann, 1998; Oguri & Kawano, 2003)

the time delay can potentially be measured with high accuracy, 2. 2.

their standard candle nature can partly break the mass and source sheet degeneracy, 3. 3.

their transient nature allows for accurate reference imaging.

Given the persistent tension [currently at the $4.4\,\sigma$-level (Riess et al., 2019)] between $H_{0}$ inferred from local distance indicators and the Cosmic Microwave Background (CMB) (Aghanim et al., 2016), independent measurements of the current expansion rate will shed light on possible explanations (e.g., Mörtsell & Dhawan, 2018).

The first observed SN Ia with expected multiple images is PS1-10afx at redshift $z=1.388$ and a flux magnification $\mu\sim 30$ (Chornock et al., 2013; Quimby et al., 2013). However, the strong lensing nature of the system was not verified by high spatial resolution imaging.

Searches for high-redshift core-collapse supernovae, intrinsically fainter than SNe Ia, has profited from observations through lensing clusters of galaxies (Goobar et al., 2009; Amanullah et al., 2011; Petrushevska et al., 2016). Four images of the core-collapse supernova Refsdal at $z=1.49$ were also observed in the line of sight of a massive cluster (Kelly et al., 2016b). Lens models of the cluster predicted that a fifth image should reappear within 500 days, in agreement with the subsequent detection of the image after $310-380$ days (Kelly et al., 2016a). Unlike SNe Ia, lensed core-collapse supernovae can not be used to measure the lensing magnification directly.

The first strongly lensed SN Ia, iPTF16geu at redshift $z=0.409$, was identified by its high magnification $\mu>50$ (Goobar et al., 2017). Subsequent high spatial resolution imaging confirmed the multiple images of the SNe Ia. In Goobar et al. (2017), the positions of the SN images with respect to the lensing galaxy were used to construct a lensing model, an isothermal ellipsoid galaxy (Kassiola & Kovner, 1993; Kormann et al., 1994) with ellipticity $\epsilon_{e}=0.15\pm 0.07$ and mass $M=(1.70\pm 0.06)\cdot 10^{10}\,M_{\odot}$ within a radius of $\sim 1\,{\rm kpc}$. The total magnification of the SN images was not well constrained by the model, but the adopted smooth lens halo predicted brightness differences between the SN images in disagreement with observations, providing evidence for substructures in the lensing galaxy, possibly in the forms of stars. The time delays between images were predicted to be shorter than 35 hours at $99.9\,\%$ confidence level. The system was also studied in More et al. (2017) where the anomalous image flux ratios were confirmed and time delays between images were predicted to be less than a day. Also, the slope of the lens mass distribution was constrained to be $\eta=2.1\pm 0.1$, consistent with an isothermal profile for which $\eta=2$. In Yahalomi et al. (2017), anomalies between macro-model predictions and the observed iPTF16geu flux ratios were investigated, finding that the discrepancies are too large to be due to microlensing alone. Using the macrolens models from More et al. (2017), the probability for microlensing to explain the observed flux ratio was found to be $0.3\,\%$. Allowing for alternative macromodels, with varying contributions from external shear and ellipticity, increased the probability to $3\,\%$.

Substantial effort in obtaining follow-up data, including reference imaging after iPTF16geu had faded, has resulted in smaller uncertainties in SN image positions and lens galaxy properties, as well as the first observational constraints on the time delay between SN images, as presented in an accompanying paper (Dhawan et al., 2020), summarized in section 4.

In this paper, we take advantage of the better observational constraints on the system to refine the lens model and investigate the dependence on the expansion rate $H_{0}$. We also quantify to what degree lens galaxy stellar microlensing can explain the observed image flux anomalies taking the new data into account.

We use geometrized units for which $c=G=1$ and express the dimensionless Hubble constant, $h$, in units of 100 km/s/Mpc,

[TABLE]

2 Gravitational lensing

When light from a source at angular position $\vec{\beta}$ and redshift $z_{\rm s}$, passes a mass at $z_{\rm l}$ along its path, it is deflected and observed at an angle $\vec{\theta}$. The time delay, $\Delta t$, compared to an unlensed image is given by

[TABLE]

where the combination

[TABLE]

is denoted the time delay distance. Here, $D_{\rm l},D_{\rm s}$ and $D_{\rm ls}$ are angular distances to the lens, source and between the lens and source, respectively, and $\Psi(\vec{\theta})$ is the scaled, projected Newtonian potential, $\Phi$, of the lens

[TABLE]

integrated along the path of the light ray. It is related to the surface mass density of the lens, or convergence through the Poisson equation, $\kappa=\nabla^{2}\Psi/2$, giving

[TABLE]

where

[TABLE]

Using Fermat’s principle that light rays traverse paths of stationary time with respect to variations of the path, we obtain the lens equation

[TABLE]

where $\vec{\alpha}(\vec{\theta})=\nabla\Psi(\vec{\theta})$. Given the Jacobian matrix of the lens mapping,

[TABLE]

the magnification, $\mu(\vec{\theta})$, is given by

[TABLE]

where the shear $\gamma=\sqrt{\gamma_{1}^{2}+\gamma_{2}^{2}}$ and (here, subscripts denote partial derivatives with respect to angle components $\theta_{i}$)

[TABLE]

2.1 Mass sheet degeneracy

If we rescale $\Psi$ using the scalars $\xi$ and $u$, and the vector $\vec{s}$ as

[TABLE]

the deflection angle $\vec{\alpha}=\nabla\Psi$ rescales as

[TABLE]

and the convergence according to

[TABLE]

If we also rescale the (unobserved) source position $\vec{\beta}^{\prime}=\xi\vec{\beta}-\vec{s}$, the image positions will be unchanged. This is the mass sheet degeneracy; given only the observed image positions, we are free to rescale the projected mass111Possibly restricted by physical considerations, such as keeping the projected mass positive.. Magnification and time delay predictions, however, will change according to

[TABLE]

From equation 2, the inferred $h$ is proportional to the predicted $\Delta t^{\prime}$ from the lens model, $h^{\prime}=\xi h$. Error propagation gives

[TABLE]

where $\Delta m$ is the magnification expressed in magnitudes. For a SN Ia, $|\Delta m|\approx 0.1$ and the fractional uncertainty on $h$ from the mass sheet degeneracy is of order 5 %, ignoring additional substructure magnifications, see section 8.3.

3 Lens model

For a cored isothermal ellipsoid, the convergence is given by (Kassiola & Kovner, 1993; Kormann et al., 1994)

[TABLE]

where $b$ is a mass normalization corresponding to the Einstein radius $\theta_{\rm E}$ for $s=0$, $s$ is the core radius and $\epsilon$ is related to the minor and major axis ratio $q$ as

[TABLE]

In terms of the ellipticity, $\epsilon_{e}$, and eccentricity, $e$,

[TABLE]

We can generalize equation 16 as

[TABLE]

for which $\theta_{\rm E}=b\sqrt{(1+q^{2})/2q}$. In the core-less circularly symmetric case ($\epsilon=s=0$), equation 19 is a good approximation of the projection of a three-dimensional density profile $\rho\propto r^{-n}$, where $r^{2}=\theta_{1}^{2}+\theta_{2}^{2}$, for $\eta=n$. For a singular isothermal sphere (SIS), $\epsilon=s=0$ and $\eta=n=2$. We refer to its elliptic generalization $\epsilon\neq 0$ as a singular isothermal ellipsoid (SIE). For an SIE lens, the image magnification and convergence are related by $\mu_{\rm SIE}=(1-2\kappa_{\rm SIE})^{-1}$.

In terms of the dimensionless $\tau$, the time delay between two lensed images, I and II, of a single source is

[TABLE]

For a SIS halo, $\Psi=\vec{\theta}\cdot\vec{\alpha}$ and

[TABLE]

Converting to polar coordinates, for two images on opposite sides of the lens, we can gain some analytical insight by Taylor expanding in (Mörtsell et al., 2005)

[TABLE]

obtaining

[TABLE]

For symmetric systems, like iPTF16geu, for which $p$ and $\delta$ are small, the time delays are simply given by the difference of the squared image distances from the lens centre, times a factor $(\eta-1)/2$.

Note that the time delay only depends on the slope of the projected surface density at radii between the images (Falco et al., 1985). The main uncertainty affecting the precision of $\Delta\tau$, and thus $h$, is usually the slope of the surface mass density, $\eta$, in the annulus between the images.

4 Summary of observations

We here summarize the observations, and the derived characteristics of iPTF16geu, as described in more detail in Goobar et al. (2017) and Dhawan et al. (2020).

For the iPTF16geu lens system, the redshifts are $z_{\rm l}=0.2163$ and $z_{\rm s}=0.409$, corresponding to angular distances, assuming flat Planck parameter values of $h=0.678$ and $\Omega_{\rm M}=0.308$ (Ade et al., 2016),

[TABLE]

Expressing positions in arc seconds [arcsec, ”], the time delay distance in equation 3, is given by (again assuming flat Planck parameter) $D_{\Delta t}=57.27\,{\rm days}$.

The discovery of iPTF16geu prompted observations with the Hubble Space Telescope (HST) Wide Field Camera 3 using the ultra violet and near-infrared (NIR) channels, as well as with laser guided star adaptive optics (LGS-AO) at the Very Large Telescope (VLT) and Keck (Goobar et al., 2017). Subsequent data while iPTF16geu was active were collected until it disappeared behind the Sun, as well as after it had faded below the detection limit using HST and LGS-AO NIR observations at the Keck telescope in $J,H$ and $K_{S}$ bands. For $J$ and $H$, we obtained 9 exposures of 20 seconds each following a dithering pattern. For the $K_{S}$ band, 18 exposures of 65 seconds were acquired (Dhawan et al., 2020). Together with the data presented in Goobar et al. (2017), this resulted in 3 epochs for $J$ and 2 epochs for the $H$ and $K_{S}$ bands, respectively. The HST data are summarized in table 1 in Dhawan et al. (2020).

Multiband photometry for the resolved images of iPTF16geu is obtained using a combination of a forward modelling approach and template subtractions. The data described above is used to fit for the global lightcurve parameters, including lightcurve peak and shape, as well as color excess from intrinsic color variations or dimming by dust in the host galaxy. The light curves are also corrected for Milky Way dust. Additionally, the resolved images are used to fit for the time offsets between the four SN images, the SN image positions, as well as extinction in the lensing galaxy for each individual line of sight. Assuming a total-to-selective extinction ratio $R_{V}=A_{V}/E(B-V)=2$ in the host and lens galaxies, we find $E(B-V)_{\rm host}=0.29\pm 0.05$ mag (common for all images) and $E(B-V)_{\rm lens}$ for images 1, 2, 3 and 4 of $0.06\pm 0.08$, $0.17\pm 0.08$, $0.42\pm 0.09$ and $0.94\pm 0.07$ mag, respectively, showing that especially image 4 is heavily extinguished in the lens galaxy. Combining the intrinsic source luminosity standardized from the light curve shapes and colours, with the light curve peak heights and individual image dust extinctions, the magnification factor for each image can be computed.

The derived image positions, arrival times and magnifications are constrained to the values given in table 1 (Dhawan et al., 2020). The total magnification of iPTF16geu is $\mu_{\rm tot}=67.8^{+2.6}_{-2.9}$. The corresponding time delay differences, $\Delta t_{ij}\equiv t_{j}-t_{i}$, and magnification ratios, $r_{ij}\equiv\mu_{i}/\mu_{j}$, are presented in table 2. The fact that data is consistent with zero time delay between images means that we can only hope to be able to put a lower limit on $h$.

5 Macrolens image modelling

In this section, we constrain the parameters of the macrolens model given in equation 16 using imaging of the lens and host galaxy light distribution, as well as the iPTF16geu image positions. Here, the fluxes of the individual iPTF16geu image fluxes are not included in the fits. In subsequent sections, we also take advantage of the standard nature of the supernova source to make use of the iPTF16geu image magnifications in the modelling procedure.

5.1 Lens mass constraints from iPTF16geu image positions

In line with the analysis in Goobar et al. (2017), we first use only the SN image positions to constrain the parameters of the lens mass model, including $b$, $\epsilon$, the slope $\eta$, the core size $s$, the orientation of the major axis, the position of the centre of mass of the lensing galaxy, and the corresponding iPTF16geu source position, employing the lensmodel software package (Keeton, 2001a, b). Since the image fluxes can be subjected to systematics effects such as lensing by substructure and dust extinction not properly corrected for, this represents the most conservative approach to the problem222In fact, it is not possible to successfully reproduce the observed iPTF16geu fluxes with the employed lens model, even when allowing for varying slope and core size, $\eta$ and $s$..

The slope and core size are not well-constrained by the SN positions, since varying $\eta$ and/or $s$ can be compensated for by changing $b$ in a way that keeps the total mass within the images constant. When varying the slope of the lens, flux ratios vary little, but the total magnification a lot. We can thus get a wide range of magnifications, even without invoking the mass sheet degeneracy.

At least in principle, the central surface density can be constrained by the fact that we do not observe a central image of iPTF16geu. Close to the centre of the lens, the image magnification is given by

[TABLE]

If we can observationally constrain $\mu_{\rm central}<\mu_{\rm max}$, then (for $\kappa>1$)

[TABLE]

where $\Sigma_{\rm crit}\approx 10.50\,{\rm kg}/{\rm m}^{2}$ for iPTF16geu. Since we expect the central image to be subject to large dust extinction not accounted for in the analysis, equation 25 represents the most optimistic limit obtainable for the lens galaxy central surface density. Assuming a Hsiao template (Hsiao et al., 2007), with a Bessell B band magnitude of -19.3, reddened by $A_{V}=0.23$ magnitudes in the Milky Way and $A_{V}=0.58$ magnitudes in the host galaxy, the $3\,\sigma$ upper limit on a possible fifth central image corresponds to a demagnification of $1.2$ magnitudes, or $\mu_{\rm max}=0.33$, implying $\Sigma_{\rm central}>29\,{\rm kg}/{\rm m}^{2}$. At least in principle, the lack of a detected central image can be used also to constrain the slope of a parametrized lens mass model. However, in the case of iPTF16geu, also in the case of a very flat halo with $\eta=1.2$, the central image would be below our detection limit, even when ignoring dust extinction.

5.2 Lens mass constraints from imaging data

In order to constrain also the slope of the deflector, we include full imaging data using lenstronomy, an open-source python package using forward modelling to reconstruct strong lensing systems (Birrer & Amara, 2018). We simultaneously reconstruct the lens mass model, the SN images (when applicable) and the surface brightness distributions of the lens and the lensed iPTF16geu host galaxy. Compared to systems of strongly lensed quasars, we can take advantage of the fact that the transient nature of SNe allows for accurate reference imaging, i.e., we have a clear view of the lensed host when the SN has faded.

The lens mass is described by equation 19, using zero core size, $s=0$. For the lens galaxy and (unlensed) host galaxy light intensity, we assume elliptical Sérsic profiles

[TABLE]

where $I_{e}$ is the intensity at the half-light radius $R_{e}$, $b_{n}=1.9992n-0.3271$ (Birrer & Amara, 2018) and

[TABLE]

where $x_{1}$ and $x_{2}$ are coordinates aligned with the major and minor axes of the galaxy light distribution and $q_{S}$ is the axis ratio.

5.2.1 HST imaging data

We expect to get the most reliable and conservative results from the lens modeling from HST data, given the relatively well known and stable point spread functions (PSFs), determined empirically in Appendix B of Dhawan et al. (2020). However, due to their limited resolution, compared to the size of the lens system, HST data may not be very effective in constraining the slope $\eta$ of the lens mass distribution. Out of the available HST filters, we concentrate on $F814W$ data, being least subject to dust extinction and displaying the most evident separation between the SN, the lens galaxy and the host galaxy fluxes.

Minimizing the residual difference between the observed and the model image for a range of lens mass slopes in the interval $1.2<\eta<2.5$, we find that the $F814W$ data only provide an upper limit of $\eta\lesssim 1.7$ at $95\,\%$ confidence level (CL), in contrast to previous work in which the slope was constrained to $\eta=2.1\pm 0.1$ (More et al., 2017). The best fit model, for which $\eta=1.2$ has a reduced chi-square of $\chi_{\nu}^{2}=0.99735$. Low values for $\eta$ are also preferred by $F814W$ reference imaging after iPTF16geu had faded, although the overall quality of the fit is somewhat decreased compared to the imaging when iPTF16geu was active. We note that a possible reason for this, and also for the slightly discrepant value of $\eta$ compared to constraints from Keck imaging, velocity dispersion measurements and iPTF16geu image magnifications discussed in upcoming sections, is that the host images are affected by dust absorption as evident from the iPTF16geu images. Also, figure 3 in Johansson et al. (2020) show evidence of colour gradients in the host galaxy along the Einstein ring, possibly related to lens galaxy dust extinction. Given the longer integration time of the HST reference images, the noise level is decreased, and the sensitivity to not correcting the host galaxy images for the differential dust extinction (which is beyond the scope of this paper) is increased.

5.2.2 Keck imaging data

Since constraining the slope relies on high spatial resolution data, preferably at long wavelengths in order to minimize sensitivity to dust extinction, we also investigate to what extent Keck host imaging data can provide additional information about the lens mass distribution. The Keck LGS-AO PSF may vary substantially both in time and over the focal plane, potentially making the lens modelling less reliable than for HST data, despite the higher resolution of the Keck LGS-AO data. Since iPTF16geu and reference imaging in the $H$-band were obtained with different instruments, and conditions during $K_{S}$-observations were unstable, we will use data obtained in the $J$-band, having the highest quality out of the three Keck bands.

We approximate the PSF as a Moffat profile being constant in time and space. Since there are no isolated stars in the field, we determine its parameters from the SN images, after subtracting the reference image. We then model the lens mass distribution from the reference image of the lensed host galaxy, since we expect this to be less sensitive to the exact shape of the PSF, compared to imaging when iPTF16geu is active. The lens model parameters obtained are consistent with that from HST data, except that it favours a slightly higher slope of $\eta=1.8\pm 0.15$ at $95\,\%$ CL, see figure 1. Adding information on the iPTF16geu image positions, when fitting the $J$-band image slightly shifts the best fit value of the slope to $\eta\sim 1.7$.

Summarizing, there is a slight tension between the lens mass slope best fitting the HST and Keck data. One reason could be that the assumed mass and light distribution models are too restrictive, another that the tension is caused by incomplete descriptions of the image PSFs. Also, given the large extinction evident in the iPTF16geu images, we can not rule out the possibility that the HST host images are affected by dust extinction, not corrected for in the present analysis. In the next section we investigate how the slope can be independently constrained using the observed magnifications of the iPTF16geu images.

Figure 2 and 3 show the observed and reconstructed HST $F814W$-band and Keck $J$-band image for $\eta=1.8$ (for reasons that will soon become evident). We obtain very similar results for lower values of $\eta$, the only exception being that the intrinsic size and luminosity of the host galaxy becomes smaller since the magnification increases as $\eta$ decreases. From left to right, the upper panel shows the reconstructed lens galaxy light distribution, the unlensed host and the lensed host image. In the lower panel, the observed image is shown to the left and in the middle the total reconstructed image, i.e. the sum of the lens and the lensed host light intensity. The reconstructed images have been convolved with the PSF of the instrument, but do not include background or Poisson noise. To the right, normalized residuals after subtracting the reconstructed (model) image and the observed image (data) are shown.

6 Macrolens velocity dispersion modelling

An independent velocity dispersion measurement of the lensing galaxy can be used to break the mass sheet degeneracy, and to constrain the slope of the lens mass distribution. Johansson et al. (2020) present high-resolution spectroscopic observations of iPTF16geu, and derive a velocity dispersion of the lens galaxy, $\sigma_{\rm obs}=129\pm 4\,{\rm km/s}$ using stellar continuum fitting. This value is lower, but consistent with the value derived from emission lines in low-resolution spectra, $\sigma_{\rm obs}=163^{+41}_{-27}\,{\rm km/s}$ (Goobar et al., 2017).

The theoretical velocity dispersion for a spherically symmetric system can be derived from the equations of stellar hydrodynamics:

[TABLE]

where $\sigma_{t}$ and $\sigma_{r}$ are the velocity dispersions in the tangential and radial direction, respectively, $\beta=1-(\sigma_{t}/\sigma_{r})^{2}$ is the velocity anisotropy, $\nu$ is the three dimensional density of velocity dispersion tracers (luminous matter) and $\Phi$ is the total gravitational potential. The prime indicates differentiation with respect to $r$. Assuming constant $\beta$,

[TABLE]

To estimate the expected velocity dispersion of the iPTF16geu lens galaxy, we approximate the light and mass distributions as spherically symmetric and deproject the best-fit Sérsic profile for the lens luminosity distribution and the power-law projected total mass profile to obtain $\nu(r)$ and $\Phi(r)$ respectively333We have used codes developed by Matthew Auger and Alessandro Sonnenfeld available at https://github.com/astrosonnen/spherical_jeans..

The lens luminosity distribution is derived from HST $F475W$, expected to most closely trace the star distribution used when observationally constraining the velocity dispersion, see Johansson et al. (2020). In doing this, we keep the lens mass distribution obtained from $F814W$ data fixed. The model velocity dispersion, $\sigma_{\rm mod}$, is given by a line-of-sight luminosity weighted average over the effective spectroscopic aperture of the observations, here approximated as a circular aperture with radius $0.9$ arcsec. In figure 4, we show $\sigma_{\rm mod}$ as a function of lens mass slope $\eta$. Compared with the observed velocity dispersion of $\sigma_{\rm obs}=129\pm 4\,{\rm km/s}$ indicated by the horizontal band, we see that $1.6<\eta<2$ are prefered, in reassuring agreement with the values obtained from the iPTF16geu image fluxes in section 8.3. The main uncertainty is the anisotropy parameter, $\beta$, for which we conservatively assume $\beta=0.2^{+0.2}_{-0.5}$ (Gerhard et al., 2001; Chae et al., 2019), yielding $\sigma_{\rm mod}=132^{+4}_{-7}\,{\rm km/s}$ for $\eta=1.8$.

7 Macrolens image magnifications modelling

Since the inferred value of the Hubble constant depends directly on the slope of the projected lens mass $\eta$, and smaller $\eta$ implies larger image magnifications, in this section we take advantage of the standard candle nature of SNe Ia observationally constraining the image magnifications as listed in table 1, in order to resolve the tension between the inferred slope from HST and Keck data. In figure 5, the predicted magnifications for the individual iPTF16geu images as a function of the slope $\eta$, as well as the total magnification are plotted together with the observed values indicated by the horizontal bands. Two things are immediately apparent: Ignoring substructure lensing effects, the observed total magnification suggests that $\eta\sim 1.75$ (close to the value preferred by Keck data, see figure 1, and velocity dispersion observations, see figure 4): Also, regardless of the value of $\eta$, a smooth lens density fails to explain the individual iPTF16geu image fluxes, and additional sub-structure lensing is needed. This is the subject of the upcoming section.

8 Microlens modelling

In order to use the iPTF16geu image magnifications when modelling the macrolens mass distribution, we also have to take into the fact that substructure in the form of e.g. stars in the lens galaxy may give rise to additional (de)magnification of the iPTF16geu images.

8.1 Stellar microlensing

Lens galaxy stars can microlens strongly lensed compact sources, such as quasars or SNe, significantly altering their observed fluxes (Chang & Refsdal, 1979). Although the principles of microlensing are the same as those of galaxy-scale strong gravitational lensing (see section 2), there are some important differences between the two scenarios. First, the significantly lower deflector masses involved in microlensing ($m\sim M_{\odot}$) lead to microimage time delays and image separations that are too small to be observed with current optical instrumentation (typically of order 10 -6 seconds and 10*-6* arcsec, respectively; Moore & Hewitt, 1996). Second, individual deflectors are replaced by fields of lens stars, having complex caustic patterns that vary over spatial scales of microarcseconds. These patterns can magnify or demagnify compact sources by several magnitudes over the value expected from a smooth lens model (e.g., Kayser et al., 1986; Schneider & Weiss, 1987; Wambsganss, 1992). This effect has been observed in many lensed quasars (e.g., Irwin et al., 1989; Witt et al., 1995; Tewes et al., 2013), and simulations indicate that it should be ubiquitous in the multiple images of strongly SNe (Dobler & Keeton, 2006; Goldstein et al., 2018).

To determine if microlensing can explain the flux anomalies in iPTF16geu, we use a set of publicly available microlensing magnification patterns produced by the GERLUMPH project (Vernardos et al., 2014). The maps were generated by the inverse ray-shooting method (Wambsganss, 1999), in which a uniform surface density of rays is followed from the observer through a random field of point-mass deflectors in the lens galaxy, and collected in pixels on the source plane. The ray count in each pixel is proportional to the lensing magnification that a source at the position of the pixel would experience.

The statistical properties of the microlensing are determined by the macrolensing model convergence and shear, $\kappa$ and $\gamma$, as well as the fraction of matter in stars, $f_{*}$, at the location of each image. For our best fit lens mass models, images 1 and 3 correspond to saddle points of the time delay surface in equation 2 (for which $\det\mathcal{A}<0$) and images 2 and 4 to minima (for which $\det\mathcal{A}>0$ and $\operatorname{Tr}\mathcal{A}>0$). In the simulations, fields of stars were realized at the location of each image by assuming a uniform deflector mass $m=M_{\odot}$, (Vernardos et al., 2014). Studies of lens galaxy star fields have established that $m=0.3M_{\odot}$ is a more representative value (Poindexter & Kochanek, 2010), but other studies indicate that the deflector mass has a negligible effect on the microlensing magnification probability distributions (Wambsganss, 1992; Lewis & Irwin, 1995; Wyithe & Turner, 2001; Schechter et al., 2004).

8.2 Stellar mass fractions

With the knowledge of stellar mass-to-light ratios, $m/L$, we can estimate the stellar mass surface density at the iPTF16geu image positions from the lens light intensity. Combined with the total surface mass density as derived from the lens modelling, we obtain the stellar mass fraction, $f_{*}$, needed to determine the stellar microlensing probabilities for the SN images. In Kauffmann et al. (2003), $m/L$ as a function of the host galaxy magnitude are presented for 4 out of the 5 Sloan Digital Sky Survey (SDSS) band, $g,r,i$ and $z$. Since the scatter in $m/L$ decreases at longer wavelengths, we use lens galaxy imaging in the HST $F814W$-band, which overlaps to large extent with the SDSS $i$-band. The lens light distribution is derived using the parametric fit described in Dhawan et al. (2020).

The absolute magnitude of the lens galaxy is $k$-corrected to the $i$-band at a redshift of $z=0.1$ is $M_{i}=-20.4$. Using figure 14 in Kauffmann et al. (2003), gives a median $(m/L)_{i}=1.70$, in solar units for which the absolute magnitude is $M_{\odot,i}=4.58$ at $z=0.1$. The 5th, 25th, 75th and 95th percentile mass-to-light ratios are given by $(m/L)_{i}=[0.63,\,1.12,\,2.34,\,3.63]$. The corresponding stellar mass fractions, $f_{*}$, at the SN image positions for $\eta=1.8$ (with similar results for other values of $\eta$) are presented in table 3.

8.3 Stellar microlensing probabilities

In figure 6, we show the microlensing magnification probability density functions (PDFs) for each image, assuming a slope of $\eta=1.8$. Here, $\mu_{\rm tot}$ is the total magnification of each image and $\mu_{\rm th}$ the (absolute value of the) smooth lens model magnification. Thanks to the standard candle nature of SNe Ia, the required microlensing magnification can be estimated (vertical dotted lines) and compared to the microlensing PDFs, assuming different values of $\eta$. The probability for each image is obtained by taking a weighted average of the one-tailed $p$-value for the corresponding $f_{*}$ (listed in table 3 for the case of $\eta=1.8$). The total probability, $p_{\rm tot}$, for stellar microlensing to explain the observed iPTF16geu image fluxes is obtained by combining the individual image probabilities using the Fisher’s combined probability test (Fisher, 1925). As shown in figure 7, this probability is maximized for $\eta=1.8$, yielding $p_{\rm tot}=0.12$, primarily driven by the large microlensing magnification required for image 1. In table 4, we list the convergence $\kappa$, the shear $\gamma$, the (signed) smooth lens model magnification $\mu_{\rm th}$, the observed magnifications $\mu_{\rm obs}$, and the probability $p$ for stellar microlensing to accommodate the difference between $\mu_{\rm th}$ and $\mu_{\rm obs}$ for the case of $\eta=1.8$. We note that the probability derived is higher than that obtained in Yahalomi et al. (2017). This can be attributed to the difference in the assumed slope of the lens mass distributions, our updated magnifications factors taking dust extinction into account, and differences in the way the individual probabilities for each image are combined. We conclude that stellar microlensing, with an acceptable level of probability can explain the observed flux anomalies for iPTF16geu. Note that the use of microlensing to place constraints on the macrolens model is only possible since the source is a precisely calibrated standard candle.

9 Lens model and Hubble constant results

In this section, we summarize the results on the lens mass distribution, and derive limits on the Hubble constant using the observed time delay between the iPTF16geu images.

9.1 Lens model parameters

Given the constraint on the lens mass slope derived from the iPTF16geu image magnifications and velocity dispersion observations, strengthened by their agreement with LGS-AO Keck data, in the following we assume a value of $\eta=1.8$. In table 5, we present the resulting lens mass distribution parameters as derived from HST $F814W$ data, and Keck $J$ reference imaging. The two independent fits display a reassuring agreement, the only notable difference being that Keck data indicates a slightly higher ellipticity of the lens mass distribution.

Figure 8 shows the critical curves, corresponding to image positions with infinite magnification, i.e., $\det\mathcal{A}(\vec{\theta})=0$, and caustics, the corresponding source positions, for the best fit lens model with slope $\eta=1.8$ overlayed on Keck LGS-AO $J$-band imaging data.

Given the high degree of circular symmetry of the lens system, the ellipticity of the lensing galaxy is degenerate with a possible external shear component for light deflections close to the critical curve. If the external shear is caused by a SIS galaxy, its magnitude is given by

[TABLE]

where $b_{\rm ext}$ is the Einstein radius of the external perturber, $r$ the projected distance between the lensing galaxies and $\sigma_{v}$ the velocity dispersion of the perturber. The largest expected contribution from galaxies in a square field with side length 100 arcsec centred on the iPTF16geu system is from a galaxy with $r\approx 58$ arcsec at $\phi\approx 40^{\circ}$ and an SDSS photometric redshift estimate of $z=0.336\pm 0.0365$ (Blanton et al., 2017). From its $g$-band magnitude of $m_{g}=20.77$, we can estimate its velocity dispersion to $\sigma_{v}\approx 187$ km/s (Mitchell et al., 2005; Jönsson et al., 2008). This corresponds to a shear contribution of $\gamma_{\rm ext}=1.4\cdot 10^{-4}$ for iPTF16geu. Other individual galaxies in the field contribute at most $20\,\%$ of this value and we conclude that the major part of the shear is caused by the ellipticity of the lensing galaxy itself.

9.2 Time delays and $h$

The measured time delays in table 2 can be used to constrain the value of the time delay distance $D_{\Delta t}$, using equation 2, and finally the Hubble constant $h$. In practice, this is done using lenstronomy in a manner very similar to when constraining only the lens and source properties. The only difference is that we also include the time delays as observational constraints. In principle, including time delays may also change the derived lens properties. However, since in the case of iPTF16geu the time delay uncertainties are large, in practice, they will only constrain $D_{\Delta t}$.

Given the large uncertainties in $D_{\Delta t}$, we have assumed fixed (Planck) values for cosmological parameters other than $h$ when translating from $D_{\Delta t}$ to $h$. Figure 9, shows the confidence contours for $D_{\Delta t}$, together with the lens mass parameters $\theta_{\rm E}$, $\phi_{\rm mass}$ and $q$, derived from HST $F814W$-data assuming a lens mass slope of $\eta=1.8$. Since the observed time delays are consistent with being zero at $1\,\sigma$, we can only (weakly) constrain $D_{\Delta t}$ from above, and consequently $h$ from below to $h>0.5$ and $h>0.2$ at $68.3\,\%$ and $95\,\%$ CL, respectively. Since the derived Hubble constant scales with the observed time delays, $\Delta t_{ij}$, and the slope of the lens mass distribution as $h\propto(\eta-1)/\Delta t_{ij}$, conservatively assuming a very shallow lens profile with $\eta=1.6$, the constraints are slightly weakened to $h>0.4$ and $h>0.15$ at $68.3\,\%$ and $95\,\%$ CL.

As discussed in Dhawan et al. (2020), we expect to be able to decrease the uncertainty on $h$ considerably, given a strongly lensed SN Ia observed closer to its maximum light, and with larger time delays between the images.

10 Conclusions

With the aid of follow-up observations of iPTF16geu, we derive tighter constraints on the lens model, as well as the first constraint on the Hubble constant from observations of a strongly lensed SN Ia, following Refsdal’s original proposal (Refsdal, 1964). The projected mass distribution of the lens galaxy is determined, including the slope which is constrained by the iPTF16geu image magnifications and velocity dispersion measurements to $\eta\sim 1.8$, flatter than an isothermal profile for which $\eta=2$.

The fact that iPTF16geu exploded very close to the centre of the inner caustic of the lens, makes the predicted images for a smooth lens to be similar, in terms of radial position, magnification and arrival time. Since the majority of the data is post first maximum, the measurement of the time delays between images is challenging. Combined with the fact that the absolute time delays are small, the resulting fractional uncertainties are large. We expect constraints on the Hubble constant to be dominated by uncertainties in the arrival time of the images and obtain a relatively weak lower limit on the current expansion rate of $h\gtrsim 0.4$ at $68.3\,\%$ confidence level.

Observational limits on the image fluxes are better constrained, and show beyond doubt that substantial additional (de)magnification from substructure, possibly stars, in the lens has to take place. Compared to macrolens predictions, images 1 and 2 need an additional magnification of $\Delta m_{1}=-1.0$ and $\Delta m_{2}=-0.3$, whereas images 3 and 4 require a demagnification of $\Delta m_{3}=0.5$ and $\Delta m_{4}=0.7$, possibly from substructure lensing. Estimating the stellar mass fractions at the iPTF16geu image positions, we derive stellar microlensing magnification probabilities for each image. The total probability for stellar microlensing to explain the observed flux anomalies is $p_{\rm tot}\sim 12\,\%$. We conclude that the iPTF16geu flux "anomalies" are within stellar microlensing predictions. This is the first time that microlensing has been used in combination with a precisely calibrated standard candle to place constraints on the macrolens model in a strong lensing system.

Acknowledgements

We thank the anonymous referee for many insightful comments helping to improve the quality of the manuscript. Also, thanks to Simon Birrer, Rahul Gupta and Mattia Bulla for help running lenstronomy. AG acknowledges support from the Swedish National Space Agency and the Swedish Research Council.