The Magnificent Five Images of Supernova Refsdal: Time Delay and Magnification Measurements

We processed all HST images using the DrizzlePac software suite. ³⁸ All images were first aligned to a common world coordinate system (WCS) using TweakReg. We adopted the WCS used for the CLASH, GLASS, and HFF images and catalogs, ³⁹ which is also the WCS generally used by published lens models for this cluster. We then co-added and resampled each image set to a scale of 003 pixel⁻¹ using AstroDrizzle (Fruchter et al. 2010).

In the WFC3 IR imaging of the MACS J1149 field, a bright foreground star (J = 14.0 mag Vega) offset by $\sim 6^{\prime\prime}$ from images S1–S4 produces four bright diffraction spikes, two of which can be seen near the bottom of Figure 1. We requested telescope roll angles for which the stellar diffraction spikes fell away from the SN images. However, for a significant fraction of observations no such angle was available.

The diffraction spikes require additional effort to obtain a measurement of the SN flux when they overlap the SN image, but we have found that it is possible to remove the diffraction spikes almost entirely. We performed all photometry on "difference images" that were constructed by subtracting a pre-SN "template" image from the co-added images of each visit containing the SNe. To create a template image without diffraction spikes, we constructed masks that enclose pixels contaminated by the diffraction spikes. We next constructed template images by combining all WFC3 IR images acquired before 2014 October 30 (preceding the first detection of the SN) and applying the pixel masks. The images that we use to construct the template include data acquired at multiple telescope roll angles, so the co-added template images have no regions that are completely masked.

As shown by Rodney et al. (2016), the similarity among the star's four diffraction spikes and rotation of the diffraction spikes make it possible to remove effectively if not entirely diffracted light at SN positions. We first co-add the single exposures taken during each visit and then subtract the template image from the resulting co-added image. As described above, the template images we created lack diffraction spikes. By subtracting the template from each single-epoch image, we removed the light from all galaxies and the sky background. Consequently, the difference images only show light collected during the follow-up epoch, and only from two sources: the multiply imaged SN and the diffraction spikes.

We next rotate the difference images by 90°, after setting all regions outside of those contaminated by the diffraction spikes to equal zero. Given the similarity among the four diffraction spikes, each rotated spike provides an excellent model for the adjacent spikes. After subtracting this masked and rotated image of the spikes, we are left with a difference image that has an increase in noise in a narrow region along the spike, but without substantial contamination.

In Table 2, we list the coordinates of images S1–SX measured in the CLASH astrometric system from co-additions of WFC3 IR F125W imaging acquired in 2015 and 2016.

Table 2. J2000 Positions of the Five Detected Images of SN Refsdal

Image	R.A.	Decl.
S1	11h49m35575	22°23'4425
S2	11h49m35454	22°23'4482
S3	11h49m35370	22°23'4392
S4	11h49m35474	22°23'4263
SX	11h49m36022	22°23'4810

Download table as:  ASCIITypeset image

Fitting a model of the point-spread function (PSF) to an unresolved point source such as S1–SX should, in principle, yield a measurement of the source's flux with an optimal signal-to-noise ratio (S/N). For a point source, pixels closest to the source's center have the greatest S/N. When fitting a PSF model to the data, one can assign greater weights to these central pixels. By contrast, when measuring instead an aperture magnitude, the flux within a fixed radius is summed and an aperture correction is applied, and all pixels in the aperture are assigned effectively the same weight.

We have adapted the DOLPHOT (Dolphin 2000) PSF-fitting package to be able to measure photometry from WFC3 IR difference imaging. When applying DOLPHOT, we use the default sky-fitting setting (FitSky=1) and model images S1–SX as point sources (Force1=1). When fitting images S1–SX, we instruct DOLPHOT to assign additional weight to central pixels (PSFPhot=2). Finally, we enforce the additional constraint that the flux of S1–SX acquired does not change between exposures taken within each imaging epoch, except for variation due to noise (ForceSameMag=0).

We estimate the uncertainty of the measured fluxes by using DOLPHOT to inject fake point sources into individual exposures. We next measure the fluxes of the fake sources and estimate their uncertainty by computing the standard deviation of measured fluxes.

For the purpose of comparison, we use PythonPhot ⁴⁰ (Jones et al. 2015) to measure the light curves of images S1–SX of SN Refsdal from co-added imaging. The PythonPhot package includes an implementation of PSF-fitting photometry measurement based on the DAOPHOT package (Stetson 1987). We have modified the code to inject fake stars at locations at which the local background is similar to the locations of images S1–SX.

The measured photometry for all five images is shown in Figure 2, along with a piecewise polynomial model fit (which will be described in Section 3.6). A full table of the measured photometry is available in the online version of the journal. Table 3 provides an excerpt.

Zoom In Zoom Out Reset image size

To assess the performance of the photometry pipelines, we selected a set of stars from a co-added image of the field acquired with the ACS WFC F606W broadband filter and required that the stars' F125W and F160W fluxes were similar to those of S1–SX over the evolution of their observed light curves. We next ran DOLPHOT in an image section around the location of each star and determined the star's position in each image. Using these positions, we then fit for the position, proper motion, and annual parallax. Finally, we ran DOLPHOT and PythonPhot on each of the comparison stars and calculated the standard deviation and reduced χ-squared statistics of the flux measurements.

We estimate the photometric uncertainties by injecting 100 fake point sources using DOLPHOT and PythonPhot. In the case of DOLPHOT, we select areas near each star without significant contamination. In Figures A1–A5, we plot the photometry of the stars, showing a comparison between the results of the photometry codes. From the DOLPHOT tests we measure a reduced χ² of 0.59−1.82 with a median of 1.06, which affirms the accuracy of our uncertainty estimates for our DOLPHOT photometry of the SN images. These test cases show that the DOLPHOT photometry provides more reliable photometry and uncertainty estimates. We therefore adopt the DOLPHOT photometry for all SN images in the following analysis.

To estimate the uncertainties and to evaluate and combine estimates of the time delay fitting methods that we will apply to the SN Refsdal data (Section 5), we create simulated photometric data that approximate the real light-curve data. To achieve this, we adopt a piecewise polynomial model as our "fiducial" light-curve model.

Spectroscopic follow-up observations and the peculiar light curve of SN Refsdal have shown that it was an SN 1987A–like SN with a peak B-band luminosity inferred to be a factor of 5–10 brighter than that of SN 1987A, given magnification favored by current cluster models (Kelly et al. 2016a). As shown in Figure 3, a piecewise polynomial function with separate, low-order polynomial functions for the photospheric and nebular phases can describe the light curves of both SN Refsdal and nearby SN 1987A–like SNe whose light curves have broadly similar shapes. Indeed, Rodney et al. (2016) found that a low-order polynomial could describe the SN Refsdal data taken through the photospheric phase.

Zoom In Zoom Out Reset image size

Table 5 describes fit statistics, including Akaike information criterion (AIC; Akaike 1974) values, for possible low-order models for each phase. While a fourth-order or fifth-order model for the photospheric phase and a first-order model for the nebular phase are preferred according to the AIC, we select a third-order model for the photospheric phase and a second-order model for the nebular phase. We choose the lower-order model for the photospheric phase because time-variable microlensing magnification may introduce structure into the observed light curves that we do not want to include as part of the fiducial model. Moreover, as Figure 3 demonstrates, nearby SN 1987A–like SN light curves (which lack microlensing) in similar rest-frame bands are well described by polynomials with this order. We select the higher, second-order model for the nebular phase, given that its evolution with time is expected to follow the nonlinear, exponential rate of energy deposition from the decay of radioactive ⁵⁶Co.

Table 5. Piecewise Model Comparison

Phot.	Neb.	χ2	npar	AIC	${\chi }_{\nu }^{2}$	ΔtX1
2	1	993.1	18	1029.1	1.66	372.254953
3	1	984.4	20	1024.4	1.65	374.454781
4	1	958.0	22	1002.0	1.61	379.730981
5	1	952.9	24	1000.9	1.61	375.789274
6	1	1037.9	26	1089.9	1.76	365.728932
2	2	993.1	20	1033.1	1.66	372.251261
3	2	1009.7	22	1053.7	1.70	374.956267
4	2	980.9	24	1028.9	1.65	381.557306
5	2	1036.8	26	1088.8	1.75	367.115460
6	2	993.1	28	1049.1	1.69	380.966961

Note. Results from fitting piecewise models with photospheric and nebular components, each of varying order.

Download table as:  ASCIITypeset image

We adopt the GLAFIC (Oguri 2010) estimates for the lensing parameter from a model that was developed for the MACS J1149 field and updated to predict SN Refsdal's reappearance (Kawamata et al. 2016). This model was among the most accurate in matching the observed time delays and magnifications for the images S1–S4 (Rodney et al. 2016). After unblinding, we found that it best matches the full set of observations of SN Refsdal in our accompanying paper (Kelly et al. 2023). Many models have been developed of the cluster lens and could also have been used for this purpose. Table 6 reports the relevant lens parameters from this model, used in the subsequent principal set of simulations. To investigate the importance of these parameters, we repeat our entire analysis using alternative lens parameters listed in Table 7.

The stars along the line of sight to each image of SN Refsdal comprise multiple components. These can include (1) the stellar populations in the early-type galaxy lens responsible for the Einstein-cross images S1–S4, (2) the MACS J1149 brightest cluster galaxy (BCG), and (3) the intracluster light (ICL). We model the BCG and ICL together as a single underlying stellar population with the same properties (age, initial mass function (IMF), etc.). Here we describe our models for these stellar populations, and we then describe the ray-tracing calculations that we use to compute simulated magnification maps in the source plane at z = 1.49. Finally, we briefly describe our simulations of millilensing by expected dark matter subhalos.

4.2.1. Models of the Early-type Lens, BCG, and ICL

To model the light distribution of the early-type lens and the BCG+ICL, we first extracted a $33^{\prime\prime} \times 39^{\prime\prime}$ (east–west × north–south, respectively) subsection of the HFF co-added images. Using SExtractor, we next identified a set of pixels associated with each detected object (a "segmentation" image). We assembled a merged mask of all sources except for the BCG and used GALFIT (Peng et al. 2002) to fit a Sérsic model centered on the BCG to the WFC3 F125W image subsection using the mask. During the fit, the sky is fixed to be zero, since the HFF imaging background is already removed. Morishita et al. (2017) found that the background estimate used to construct the HFF co-added images was accurate for the MACS J1149 field.

The best-fitting Sérsic model for the single profile fit to the light distribution of the BCG, which transitions to that of the ICL, in the F125W bandpass has a Sérsic index of 4.86, an effective radius of 196, and an F125W magnitude of 16.18 AB. Holding all parameters except for the total flux fixed, we next perform fits to subsections of the ACS WFC F435W, F606W, and F814W and WFC3 IR F105W, F125W, F140W, and F160W HFF co-additions. We subtract the best-fitting model for the BCG+ICL in each bandpass from the 33'' × 39'' image subsection. The top panel of Figure 4 shows the region around the SN images S1–S4 after subtracting this BCG+ICL model.

Zoom In Zoom Out Reset image size

The next step is to fit the early-type galaxy lens in the BCG+ICL-subtracted F125W image. We create a mask that excludes all pixels except those inside of an elliptical aperture with the same position angle, center, and ellipticity as the early-type galaxy. The aperture's semimajor axis is $1.8^{\prime\prime}$. Circular masks with radii of 015 are also placed over the locations of images of S1–S4 of SN Refsdal. Our GALFIT model consists of a Sérsic model for the early-type galaxy lens and a background level that accounts for flux from the superposed lensed host galaxy.

From the F125W image we find a Sérsic index of 2.22, an effective radius of 026, and a total magnitude of F125W = 20.16 AB for the early-type lens. The bottom panel of Figure 4 shows the S1–S4 region after this early-type galaxy has been subtracted.

We next repeat the fit for the BCG+ICL-subtracted co-added image sections in all other available HST bandpasses from 4000 Å to 1.6 μm, holding all parameters fixed except for the galaxy's total flux and the local surface brightness. Figure 5 shows the resulting measured SEDs for the BCG+ICL, the early-type galaxy lens, and the host galaxy of SN Refsdal. The inferred SED of the BCG+ICL matches that of the early-type galaxy.

Zoom In Zoom Out Reset image size

4.2.2. SED Fitting to Estimate the Stellar Mass Density at the SN Image Positions

Using the estimated shear and convergence along the lines of sight to each image, we generated caustic magnification patterns in the source plane with an adapted version of the MULE inverse ray-tracing code (Dobler et al. 2015). However, we obtain similar results instead using the microlens code (Wambsganss 1990, 1999). We next use a spherically symmetric model for SN 1987A (Dessart & Hillier 2019) at phases of 10.9, 21.3, 31.1, 41.5, 50.1, 66.8, 80.8, 97.8, 117.0, 177.0, 220.0, and 300.0 days to generate the microlensing simulations. We rescale the radius of the photosphere models linearly with time and interpolate between the rescaled model to approximate the effect of the expansion of the SN during its photospheric phase.

In our simulation, each pixel has a length of 0.002 Einstein radii for a solar-mass lens. For the photospheric expansion velocity of −8356 ± 1105 km s⁻¹ at a phase of −47 ± 8 days relative to maximum brightness measured by Kelly et al. (2016a), the SN photosphere would expand in radius by a single pixel each 1.25 days in the rest frame at z = 1.49. Each ray-tracing simulation creates a map of the magnification across a grid of 1024 × 1024 pixels, and we place the simulated SN photosphere at the center of the pixel grid. Figure 6 plots example source-plane magnification maps at the positions of S1–SX. The fact that we use a spherically symmetric explosion model may present a limitation of the simulation. A more sophisticated treatment would use a two-dimensional or even three-dimensional explosion model (e.g., Goldstein et al. 2018; Bonvin et al. 2019). To approximate an asymmetric explosion, we perturb the circular isophotes of the symmetric model to become elliptical with a ratio of the major axes of b/a = 0.8. The isophotal brightness of each elliptical isophote matches that of the circular isophote whose diameter is equal to the length of the ellipse's major axis. As shown in Table 9, the elliptical photosphere yields almost identical constraints on the SX–S1 time delay.

Zoom In Zoom Out Reset image size

In Figure 7, we show the simulated impact of microlensing on the F125W–F160W color. This indicates that stellar microlensing of the photosphere of an SN 1987A–like SN affects its apparent F125W and F160W magnitudes equally within ≲0.05 mag. Furthermore, Figure 8 shows that microlensing is not expected to strongly affect the shape of the light curve of SN 1987A–like SNe during the photospheric phase. Microlensing can also impact the overall magnification of each of the images, causing either increases or decreases to the apparent total magnification of individual SN images. Figure 9 shows the distributions of total magnifications inferred from our simulations, including effects from both "macro" lensing and microlensing.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Using the ray-tracing technique presented by Gilman et al. (2019), we compute the expected millilensing due to dark matter subhalos for a projected mass density in substructure bracketing the range of expected values based on theoretical and observational arguments (Gilman et al. 2020). In practice, we explore 10⁷ M_⊙ kpc⁻² in the main lens plane, corresponding to 0.5% substructure mass fraction, as well as 4 × 10⁷ M_⊙ kpc⁻², or a 2% mass fraction in substructure. We model the line-of-sight structure using a Sheth–Tormen halo mass function (Sheth et al. 2001). In these simulations, halos and subhalos are rendered in a range of 10⁶–10¹⁰ M_⊙. Masses outside of this range are either massive enough to contain a visible counterpart or too small to matter statistically (Gilman et al. 2020). Table 10 lists the 16% and 84% quantiles of the fractional change in magnification due to 0.5% and 2% substructure fractions, and Figure 10 plots the distributions we measure for the case where the substructure fraction is 0.5%. These show that magnification due to millilensing ranges from a ∼4% to ∼18% effect among the five SN images.

Zoom In Zoom Out Reset image size

The light curve of any SN will not, of course, be described exactly as a piecewise polynomial. Baklanov et al. (2021) modeled the light curves of SN Refsdal measured through 2015 presented by Rodney et al. (2016) and Kelly et al. (2016c) using a radiation hydrodynamics code. We first fit the Baklanov et al. (2021) F125W and F160W model using the same piecewise polynomial model that we apply to the data. We next model these residuals as a realization of a Gaussian process (GP; see Rasmussen & Williams 2006) with a squared exponential covariance function and estimate correlated deviations from the piecewise polynomial with a timescale of 10 days. We then model residuals of the piecewise model from the actual data as a combination of a correlated component common to all images and an uncorrelated ("white-noise") component whose variance differs between images. The correlated component is modeled as a realization of a GP with a squared exponential kernel and fixed timescale of 10 days. For each filter, we estimate the standard deviations of the correlated and uncorrelated perturbations, which are listed in Table 11. We use our GP model to generate realizations of perturbations to the best-fitting piecewise polynomial model that are informed by both the Baklanov et al. (2021) model and the residual variation seen in the real data.

We now have all the components we believe are needed to assemble simulated light curves that mimic the SN Refsdal data. These simulations will then be used to evaluate our light-curve fitting algorithms (Section 5) and define weights when combining the inferences from all models together (Section 5.4).

Figure 11 shows the steps to produce a single simulated data set, which was chosen at random from the 1000 primary simulated light curves. We use the fiducial light-curve model (a piecewise polynomial, described in Section 3.6) to generate 1000 simulated light curves with the same cadence and photometric uncertainties as for SN Refsdal, as well as similar apparent magnitudes to those of SN Refsdal. We apply additional perturbations from our GP model of the expected variation around the fiducial model. We next apply our simulations of the microlensing magnifications derived from the inferred stellar mass densities, shear, and convergence along the lines of sight to each image (as described in Section 4.3). Since the uncertainties in our measurements are background dominated, we do not change photometric uncertainties when rescaling by microlensing fluctuations. We include the effects of millilensing (Section 4.4) until after fitting each microlensed light curve, since its effect is only a constant change in magnification.

Zoom In Zoom Out Reset image size

In Figure 12, for the F125W and F160W light curves of SX, we compare the piecewise polynomial model that best fits the simulated light curves in Figure 11, against the same model shifted to reproduce the actual relative time delay. In the case of this simulated light curve, the best-fitting time delay is within 3 days of the true delay. Indeed, the estimated time delay can only be substantially affected by microlensing if the apparent shape of the light curve is altered, which requires a significant change in magnification with time. Given the low values of κ_⋆ for the images of SN Refsdal, the optical depth for microlensing is small, so these events should be comparatively rare.

Zoom In Zoom Out Reset image size

Four separate teams attempted to recover the relative time delays and magnification ratios from the 1000 sets of simulated light curves, using the methods described in Section 5. The fitting was fully blind, meaning that the true values of the simulated time delays and magnification ratios were hidden from the light-curve fitting teams. Section 5.4 presents the analysis of these simulation tests and their application when we combine their results.

5.1.1. Piecewise Photospheric and Nebular Polynomial Model

5.1.2. Bayesian Gaussian Process Model

5.1.3.  Supernova Time Delays Parametric Model

We fit the light-curve data using version 2 of the Supernova Time Delays (SNTD) software package (Pierel & Rodney 2019). Here we model SN Refsdal's flux as a function of time using

$$\begin{eqnarray}&&F(t)=A\displaystyle \frac{{e}^{-(t-{t}_{0})/{\tau }_{\mathrm{fall}}}}{1+{e}^{-(t-{t}_{0})/{\tau }_{\mathrm{rise}}}}+B,\end{eqnarray} \tag{ 1 }$$

presented by Bazin et al. (2009) as a generic model for SN light curves, where τ_rise and τ_fall characterize the relative rise and decline times of the light curve, respectively, and A and B set its normalization. While this model is flexible, its simple shape minimizes sensitivity to microlensing features in the light curves of S1–SX.

In order to measure the time delays of the simulated light curves described in Section 4.6, we employed SNTD's parallel method according to which the light curves of S1–SX are first fit separately before deriving joint likelihood estimates for the light-curve model parameters (see Pierel & Rodney 2019). We note that SNTD was designed to use spectrophotometric models of SNe, to leverage, when available, prior knowledge of the intrinsic color and shape of well-studied SN species including Type Ia SNe or common types of core-collapse SNe. SN Refsdal's light curve, however, was peculiar and similar to that of SN 1987A (Kelly et al. 2016c), the explosion of a blue supergiant, and a useful spectrophotometric model was unfortunately not available.

We applied the parallel algorithm to the F125W and F160W light curves separately, since Equation (1) only describes the flux in a single band. For each filter, we did two iterations. We first identified the time τ_pk as the epoch with the maximum observed flux and carried out a "coarse" fit subject to the following bounds on the fitting parameters:

$$\begin{eqnarray*}&&\begin{array}{l}{\tau }_{\mathrm{pk}}-150\lt {t}_{0}\lt {\tau }_{\mathrm{pk}}+150\\ 20\lt {t}_{\mathrm{rise}}\lt 80\\ 60\lt {t}_{\mathrm{fall}}\lt 300\\ 0.001\lt A\lt 4\\ 0\lt B\lt 2.\end{array}\end{eqnarray*}$$

From this coarse run we identified the set of parameters that returned the best χ² as (t_0,c, t_rise,c, t_fall,c, A_c, B_c). We adopted the first three of those parameter values as estimates for the time of peak, rise time, and fall time, respectively.

Our next step was to carry out a "fine" fit using the same bounds for A and B as in the coarse run while applying the following constraints:

$$\begin{eqnarray*}&&\begin{array}{l}{t}_{0,{\rm{c}}}-50\lt {t}_{0}\lt {t}_{0,{\rm{c}}}+50\\ {t}_{\mathrm{rise},{\rm{c}}}-3\lt {t}_{\mathrm{rise}}\lt {t}_{\mathrm{rise},{\rm{c}}}+3\\ {t}_{\mathrm{fall},{\rm{c}}}-10\lt {t}_{\mathrm{fall}}\lt {t}_{\mathrm{fall},{\rm{c}}}+10.\end{array}\end{eqnarray*}$$

After completing this procedure for the F125W and F160W bands separately, we defined the final time delay and magnification estimates from the band preferred by the Bayesian evidence for the fine-run fitting results.

During the testing phase with simulated light curves (i.e., before applying any fitting method to the actual Refsdal data), we experimented with more aggressive constraints on the parameter values. Specifically, we experimented with tighter bounds on the model parameters when fitting to the light curves of S2, S3, and S4 and then using these parameters to model the light curves of S1 and SX. However, this procedure resulted in a ∼7-day bias in the estimated time of peak for the S1 and SX light curves, and we discarded this approach.

5.1.4.  Python Curve Shifting3 Free-knot Spline Model

Each of the four methods described in Section 5, when given the light curves of images S1–SX, can yield a single estimate for eight lensing observables ${ \mathcal O }$. These observables consist of the four time delays between images S2, S3, S4, and SX relative to S1 (Δt_2,1,...,Δt_x,1) and the four respective magnification ratios (μ₂/μ₁,...,μ_X /μ₁). The choice of image S1 as the reference image is arbitrary. In order to combine the estimates from multiple methods with appropriate weighting for each method, we need a consistent way to define a probability distribution $p({ \mathcal O }| D;{M}_{k})$ for each method M_k given the light-curve data D. We address this by performing fits to the 1000 simulated light curves, which were described in Section 4.6.

For the ith simulated light curve, we calculate the residual of each input value ${{ \mathcal O }}_{\mathrm{input};i}$ from the estimated value ${{ \mathcal O }}_{\mathrm{est};i}$,

$$\begin{eqnarray}&&\delta {{ \mathcal O }}_{\mathrm{resid};i}={{ \mathcal O }}_{\mathrm{est};i}-{{ \mathcal O }}_{\mathrm{input};i}.\end{eqnarray} \tag{ 2 }$$

We use a histogram of these residuals to construct a probability distribution P_resid.

We next compare the ability of the four fitting methods to recover the input relative time delays. As shown in Figures 13, the polynomial method has the least bias, and the 16th, 50th, and 84th percentiles of its residual distribution are most closely separated. The rms scatter of the polynomial method, as listed in Table 12, is also smallest. In some cases, the performance of the GP and PyCS fitting methods approaches that of the polynomial method. The light curve of S4 has a lower S/N than those of the other images, however, and in that case, both methods are not as effective in recovering its relative time delay as the polynomial method.

Zoom In Zoom Out Reset image size

Considering the codes' ability to recover magnification ratios, Figure 13 and Table 12 demonstrate that all fitting codes yield comparably accurate and precise magnification ratios S2/S1 and S3/S1. In the case of the S4/S1 ratio, the GP and PyCS methods show a greater bias than the polynomial and SNTD methods. In Table 13, we report the fraction of simulated light curves for which the difference between measurement and truth ($\delta {{ \mathcal O }}_{\mathrm{resid};i}$) is more than 30 days (for time delays) or >1 (for magnification ratios).

The IMF of the stellar populations in the early-type galaxy lens (e.g., Treu 2010; Conroy & van Dokkum 2012) and the intracluster medium may have a Salpeter (1955) IMF instead of the Chabrier (2003) IMF that we have assumed for our simulations. In Tables A1 and A2, we report the same set of statistics but instead using microlensing simulations computed for a Salpeter (1955) IMF. In Figure A6, we plot the residual distributions given a Salpeter (1955) IMF. The residual distributions for simulations that adopt a Salpeter IMF show a very modest increase (≲10%) in the difference between the 16th and 84th percentiles in the case of the SX–S1 time delay and larger (∼30%) increases for the S2–S1, S3–S1, and S4–S1 images.

We next want to compute the probability distribution for a property of the light curve ${ \mathcal O }$, given a single-point estimate of the same observable ${{ \mathcal O }}_{\mathrm{est};{\rm{D}};k}$, derived from the actual light-curve data D using method M_k , to compute this probability distribution. To accomplish this, we use the distribution of residuals $P(-\delta {{ \mathcal O }}_{\mathrm{resid}})$ assembled from fits to the simulations. This probability distribution $P({ \mathcal O }| D;{M}_{k})$ is given by

$$\begin{eqnarray}&&P({ \mathcal O }| D;{M}_{k})=P({ \mathcal O }| {{ \mathcal O }}_{\mathrm{est};{\rm{D}};k})={P}_{\mathrm{resid}}({ \mathcal O }-{{ \mathcal O }}_{\mathrm{est};{\rm{D}};k}).\end{eqnarray} \tag{ 3 }$$

We computed the relative time delays and magnifications of the SN light curves after creating 1000 light curves after performing bootstrapping with replacement. Since the uncertainties we measured were smaller than those we computed from the simulated light curves, which included perturbations from the piecewise polynomial model, microlensing, and millilensing, we adopt the results from the simulations.

For each observable, we sought to identify an optimal combination of the set of four probability distributions $P({ \mathcal O }| D;{M}_{k})$ corresponding to the respective methods. We find the nonnegative linear combination of probability distributions $P({ \mathcal O }| D;{M}_{k})$ that maximizes the likelihood of the input observable values used to construct the simulations,

$$\begin{eqnarray}&&\mathop{\max }\limits_{w\in {{ \mathcal S }}_{1}^{K}}\displaystyle \frac{1}{n}\sum _{i=1}^{n}\mathrm{log}\sum _{k=1}^{K}{w}_{k}P({{ \mathcal O }}_{\mathrm{input};i}| {D}_{i};{M}_{k}),\end{eqnarray} \tag{ 4 }$$

subject to ${{ \mathcal S }}_{1}^{K}=\{w\in {[0,1]}^{K}:{\sum }_{k\,=\,1}^{K}{w}_{k}=1\}$. When computing these weights using fits to the simulated light curves, we use a modified probability function ${P}_{\mathrm{penal}}({ \mathcal O }| D;{M}_{k})$ that has the effect of penalizing bias. The penalized probability distribution is given by

$$\begin{eqnarray}&&{P}_{\mathrm{penal}}(x)={P}_{\mathrm{resid}}(x-\mathrm{median}({\rm{\Delta }}{{ \mathcal O }}_{\mathrm{resid}})).\end{eqnarray} \tag{ 5 }$$

For a distribution where the median residual, $\mathrm{median}(\delta {{ \mathcal O }}_{\mathrm{resid}})$, is equal to zero, P_penal = P_resid. Given best-fitting weights w_k , the probability distribution for the actual light curve becomes

$$\begin{eqnarray}&&P({ \mathcal O }| D)=\sum _{k=1}^{K}{w}_{k}P({ \mathcal O }| D;{M}_{k}).\end{eqnarray} \tag{ 6 }$$

We note that maximizing the likelihood of the relative time delays and magnifications used to construct the simulated light curve corresponds to minimizing the Kullback–Leibler divergence between a single-valued distribution where the residual $\delta {{ \mathcal O }}_{\mathrm{resid};i}=0$ and the stacked set of probability distributions $p({ \mathcal O }| D)$. The use of the Kullback–Leibler divergence as a means to identify an optimal combination instead of predictive models was developed by Yao et al. (2018).

The Kullback–Leibler divergence allows us to consider full probability distributions, as opposed to single-point estimates. For this analysis, we use the residual distribution for all simulated light curves to construct the probability for every light curve as described in Equation (2); the framework we describe could accommodate a separate probability distribution for each simulated light curve.

To calculate the likelihood of the input light-curve parameters in Equation (4), we need to choose how to construct a probability density function from the set of residuals $\delta {{ \mathcal O }}_{\mathrm{resid};i}$. The probability density must always be nonzero, with the bin widths sufficiently large to avoid small number statistics. We create bins such that each bin includes three residuals (i.e., each contains three values) in each of the <5% and >95% quantiles, 15 residuals in the 5%–30% and 70%–95% quantiles, and 30 residuals in the 30%–70% quantiles. Moreover, we adjust these bins so that the bin width is at least 0.5 days for the time delays, or 0.01 for the magnification ratios.

For each observable, we carry out an optimization using Equation (4) to identify the weights that maximize the likelihood of the parameters used to construct the light curves. Table 14 lists the weights that we obtain. The piecewise polynomial method receives the majority of the weight for all of the relative time delay measurements, while multiple codes receive the majority of the weight for the estimates of magnification ratios. Marginal probability distributions from the combination of four methods are shown as a "corner plot" in Figure 17.

In Figure A7, we show the separate constraints from all the methods.

We next apply the piecewise polynomial modeling code to the F125W and F160W light curves separately in order to evaluate whether the inferred SX–S1 time delays are consistent. To estimate the uncertainties associated with the measurement for each filter, we use the simulated sets of light curves in the respective filters. We find that the F125W light curve yields a constraint on the SX–S1 relative time delay of ${377.6}_{-5.04}^{+5.13}$ days, while the F160W light curve, whose measurements have, on average, lower S/Ns, yields ${370.6}_{-10.07}^{+9.78}$ days. The estimates from the F125W and F160W light curves are in 0.6σ agreement.

We explore the effects of alternative assumptions about the lens model, the mass density of stars making up the ICL, and the symmetry of the SN photosphere. In Table 7, we list the convergence and shear at the positions of images S1–SX according to the Grillo-g and Sharon-g models. Table 9 lists the constraints on the SX–S1 time delay when we instead adopt the Grillo-g and Sharon-g lensing parameters. We also obtain inferences on the SX–S1 time delay after constructing an elliptical photosphere with a ratio of major axes of b/a = 0.8 and, separately, after doubling the inferred value of κ_⋆, the density of the stars that account for the intracluster medium. Table 9 shows that adopting these alternative assumptions yields almost identical constraints on the relative time delay between images SX and S1.

We next derive a covariance matrix for Δt and μ/μ₁. The covariance matrix can be used to compute likelihoods for model predictions, although we develop a more precise approach with less sensitivity to outlying values in the accompanying paper. To compute the measurement covariance matrix, we use the 1000 simulations of the light curves of images S1–SX, which include the effects of microlensing and photometric noise.

We can compute the covariance matrix as

$$\begin{eqnarray}&&{{\rm{K}}}_{{X}_{i}{X}_{j}}=\mathrm{cov}[{X}_{i},{X}_{j}]={\rm{E}}[({X}_{i}-{\rm{E}}[{X}_{i}])({X}_{j}-{\rm{E}}[{X}_{j}])].\end{eqnarray} \tag{ 7 }$$

Here we have X and Y correspond to the distributions of residuals of the estimated value subtracted from the actual value for the simulation,

$$\begin{eqnarray}&&\delta {\rm{\Delta }}{t}_{i,1}^{\mathrm{est}}={\rm{\Delta }}{t}_{i,1}^{\mathrm{est}}-{\rm{\Delta }}{t}_{i,1}^{\mathrm{true}}.\end{eqnarray} \tag{ 8 }$$

As described in Section 5.4, we have identified the weights w_j,obs for each model M_j and observable such that the linear combination of the probability distributions,

$$\begin{eqnarray}&&p({ \mathcal O }| D)=\sum _{j=1}^{K}{w}_{j,\mathrm{obs}}\,p(\mathrm{observable}| D;{M}_{j}),\end{eqnarray} \tag{ 9 }$$

maximizes the likelihood of the value of the observable used to create the simulated light curve. Consequently, we have a set of weights w_j,obs where ∑_j w_j,obs = 1 for each relative delay Δt_1,j and magnification ratio μ_k /μ₁ between the jth image and S1.

In general, considering a discrete set of possible realizations (x_i , y_i ) of (X, Y) with probability p_i , the covariance is

$$\begin{eqnarray}&&\mathrm{cov}(X,Y)=\sum _{i=1}^{n}{p}_{i}({x}_{i}-E(X))({y}_{i}-E(Y)).\end{eqnarray} \tag{ 10 }$$

In our case, the probability p_i depends on the fitting method from which the x_i or y_i are drawn. The probability ${p}_{i}=(1/{n}_{\mathrm{sim}}){w}_{i,X}\times {w}_{j,Y}$.

In our case, to create our realizations of the stacked probability distribution, we need to draw from the probability distributions. The probability of each pair of realizations corresponds to the product of the weights w_i,X × w_j,Y that the respective fitting methods receive for the observable. We can write

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{cov}(X,Y) & = & \displaystyle \sum _{i=1}^{{n}_{\mathrm{sim}}}\displaystyle \sum _{j\,=\,1}^{{n}_{\mathrm{mod}}}\displaystyle \sum _{k\,=\,1}^{{n}_{\mathrm{mod}}}\displaystyle \frac{1}{{n}_{\mathrm{sim}}}{w}_{j,X}\times {w}_{k,Y}\\ & & \times \,({x}_{i,j}-E(X))({y}_{i,k}-E(Y)),\end{array}\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{cov}(X,Y) & = & \displaystyle \sum _{j=1}^{{n}_{\mathrm{mod}}}\displaystyle \sum _{k\,=\,1}^{{n}_{\mathrm{mod}}}{w}_{j,X}\times {w}_{k,Y}\\ & & \times \,\displaystyle \sum _{i=1}^{{n}_{\mathrm{sim}}}\displaystyle \frac{1}{{n}_{\mathrm{sim}}}({x}_{i,j}-E(X))({y}_{i,k}-E(Y)).\end{array}\end{eqnarray} \tag{ 12 }$$

While there may be model covariance between measured or predicted image separations and the time delay and magnification, we do not include it here.

Unlike the effects of microlensing, the effects of millilensing were not included when we created simulated light curves. Since the Einstein radii of stars in the foreground cluster can be comparable in size to the expanding SN photosphere, microlensing can cause changes in the shape of the light curve, as well as the mean magnification of each image (e.g., Dobler & Keeton 2006). Given the much greater masses of subhalo deflectors expected to be responsible for millilensing, their effect should be to change the mean magnification but not affect the shape of the light curve.

We are interested in the covariance of random draws x_i and y_i from the linear combinations of distributions X and Y corresponding to two observables. The probability of drawing from the distribution X_j corresponding to the jth model is w_j,X , and the probability of drawing from the distribution Y_k corresponding to the kth model is w_k,Y . Therefore, the probability of drawing from both is p = w_j,X × w_k,Y . We present this estimate for the covariance matrix in Table 17 and Figure 20. When computing the covariance matrix, we have excluded light curves for which a combined estimate from all fitting codes using the codes' respective weights represents a >30-day outlier for a time delay, or a 1.5 outlier for a magnification.

Zoom In Zoom Out Reset image size