Turning Gravitationally Lensed Supernovae into Cosmological Probes

TL;DR

This paper introduces a software tool for measuring time delays in gravitationally lensed supernovae, demonstrating that early detection improves accuracy and that Gaussian Process Regression effectively accounts for microlensing uncertainties.

Contribution

The paper presents an open source simulation and analysis package for lensed supernovae, enhancing time delay measurements and uncertainty characterization, especially regarding microlensing effects.

Findings

Early detection before peak improves measurement accuracy and precision.

Using Gaussian Process Regression effectively estimates microlensing uncertainties.

Simulations suggest improved time delay measurements with pre-peak detection.

Abstract

Recently, there have been two landmark discoveries of gravitationally lensed supernovae: the first multiply-imaged SN, "Refsdal", and the first Type Ia SN resolved into multiple images, SN iPTF16geu. Fitting the multiple light curves of such objects can deliver measurements of the lensing time delays, which are the difference in arrival times for the separate images. These measurements provide precise tests of lens models or constraints on the Hubble constant and other cosmological parameters that are independent of the local distance ladder. Over the next decade, accurate time delay measurements will be needed for the tens to hundreds of lensed SNe to be found by wide-field time-domain surveys such as LSST and WFIRST. We have developed an open source software package for simulations and time delay measurements of multiply-imaged SNe, including an improved characterization of the uncertainty caused by microlensing. We describe simulations using the package that suggest a before-peak detection of the leading image enables a more accurate and precise time delay measurement (by ~1 and ~2 days, respectively), when compared to an after-peak detection. We also conclude that fitting the effects of microlensing without an accurate prior often leads to biases in the time delay measurement and over-fitting to the data, but that employing a Gaussian Process Regression (GPR) technique is sufficient for determining the uncertainty due to microlensing.

Figures

Tables

Table 1. Table 1: The results of each time delay measurement technique for the before-peak and after-peak case studies. As expected, there is a decrease in accuracy and precision when you lose before-peak observations due to a lack of single inflection point and a decrease in the number of light curve points to fit.

\toprule   Simulation		Method	${\mu }_{\Delta t}$	${\sigma }_{\Delta t}$
Before Peak	Parallel (Section 3.2)	60.1	3.2
	Series (Section 3.3)	60.1	3.3
	Color Curves (Section 3.4)	60.1	2.6
After Peak	Parallel (Section 3.2)	58.0	5.0
	Series (Section 3.3)	59.1	5.2
	Color Curves (Section 3.4)	59.3	3.9

Table 2. Table A1: The SED models present in the SNCosmo package to describe SN evolution with wavelength and time.

Model Name	Type	Reference
nugent-sn1a	SN Ia	Nugent et al. (2002)
nugent-sn91t	SN Ia	Stern et al. (2004)
nugent-sn91bg	SN Ia	Nugent et al. (2002)
nugent-sn1bc	SN Ib/c	Levan et al. (2005)
nugent-hyper	SN Ib/c	Levan et al. (2005)
nugent-sn2p	SN IIP	Gilliland et al. (1999)
nugent-sn2l	SN IIL	Gilliland et al. (1999)
nugent-sn2n	SN IIn	Gilliland et al. (1999)
s11-2004hx	SN IIL/P	Sako et al. (2011)
s11-2005lc	SN IIP	Sako et al. (2011)
s11-2005hl	SN Ib	Sako et al. (2011)
s11-2005hm	SN Ib	Sako et al. (2011)
s11-2005gi	SN IIP	Sako et al. (2011)
s11-2006fo	SN Ic	Sako et al. (2011)
s11-2006jo	SN Ib	Sako et al. (2011)
s11-2006jl	SN IIP	Sako et al. (2011)
hsiao	SN Ia	Hsiao et al. (2007)
hsiao-subsampled	SN Ia	Hsiao et al. (2007)
salt2	SN Ia	Guy et al. (2010)
salt2	SN Ia	Guy et al. (2007)
salt2	SN Ia	Betoule et al. (2014)
salt2-extended	SN Ia	Pierel et al. (2018)
snf-2011fe	SN Ia	Pereira et al. (2013)
snana-2004fe	SN Ic	Kessler et al. (2010)
snana-2004gq	SN Ic	Kessler et al. (2010)
snana-sdss004012	SN Ic	Kessler et al. (2010)
snana-2006fo	SN Ic	Kessler et al. (2010)
snana-sdss014475	SN Ic	Kessler et al. (2010)
snana-2006lc	SN Ic	Kessler et al. (2010)
snana-2007ms	SN II-pec	Kessler et al. (2010)
snana-04d1la	SN Ic	Kessler et al. (2010)
snana-04d4jv	SN Ic	Kessler et al. (2010)
snana-2004gv	SN Ib	Kessler et al. (2010)
snana-2006ep	SN Ib	Kessler et al. (2010)
snana-2007y	SN Ib	Kessler et al. (2010)
snana-2004ib	SN Ib	Kessler et al. (2010)
snana-2005hm	SN Ib	Kessler et al. (2010)
snana-2006jo	SN Ib	Kessler et al. (2010)
snana-2007nc	SN Ib	Kessler et al. (2010)
snana-2004hx	SN IIP	Kessler et al. (2010)
snana-2005gi	SN IIP	Kessler et al. (2010)
snana-2006gq	SN IIP	Kessler et al. (2010)
snana-2006kn	SN IIP	Kessler et al. (2010)
snana-2006jl	SN IIP	Kessler et al. (2010)
snana-2006iw	SN IIP	Kessler et al. (2010)
snana-2006kv	SN IIP	Kessler et al. (2010)
snana-2006ns	SN IIP	Kessler et al. (2010)
snana-2007iz	SN IIP	Kessler et al. (2010)
snana-2007nr	SN IIP	Kessler et al. (2010)
snana-2007kw	SN IIP	Kessler et al. (2010)
snana-2007ky	SN IIP	Kessler et al. (2010)
snana-2007lj	SN IIP	Kessler et al. (2010)
snana-2007lb	SN IIP	Kessler et al. (2010)
snana-2007ll	SN IIP	Kessler et al. (2010)
snana-2007nw	SN IIP	Kessler et al. (2010)
snana-2007ld	SN IIP	Kessler et al. (2010)
snana-2007md	SN IIP	Kessler et al. (2010)
snana-2007lz	SN IIP	Kessler et al. (2010)
snana-2007lx	SN IIP	Kessler et al. (2010)
snana-2007og	SN IIP	Kessler et al. (2010)
snana-2007ny	SN IIP	Kessler et al. (2010)
snana-2007nv	SN IIP	Kessler et al. (2010)
snana-2007pg	SN IIP	Kessler et al. (2010)
snana-2006ez	SN IIn	Kessler et al. (2010)
snana-2006ix	SN IIn	Kessler et al. (2010)
whalen-z15b	PopIII	Whalen et al. (2013)
whalen-z15d	PopIII	Whalen et al. (2013)
whalen-z15g	PopIII	Whalen et al. (2013)
whalen-z25b	PopIII	Whalen et al. (2013)
whalen-z25d	PopIII	Whalen et al. (2013)
whalen-z25g	PopIII	Whalen et al. (2013)
whalen-z40b	PopIII	Whalen et al. (2013)
whalen-z40g	PopIII	Whalen et al. (2013)
mlcs2k2	SN Ia	Jha et al. (2007)

Table 3. Table A2: An example of the observation tables used to produce simulated light curves section 2.4 . Each SN is observed on a timescale of 0-100 Days, with a cadence of 3 days (Shu et al., 2018 ) for a total of 25 epochs.

Band	Time	Gain	Sky Noise	ZP	ZP System
	(Days)		(Counts)
F125W	0.0	50.0	0.008	26.3	AB
F125W	3.0	50.0	0.028	26.3	AB
F125W	6.0	50.0	0.008	26.3	AB
F125W	9.0	50.0	0.026	26.3	AB
F125W	12.0	50.0	0.011	26.3	AB
F125W	15.0	50.0	0.026	26.3	AB
F125W	18.0	50.0	0.013	26.3	AB
F125W	21.0	50.0	0.014	26.3	AB
F125W	24.0	50.0	0.029	26.3	AB
F125W	27.0	50.0	0.007	26.3	AB
F125W	30.0	50.0	0.029	26.3	AB
F125W	33.0	50.0	0.008	26.3	AB
F125W	36.0	50.0	0.024	26.3	AB
F125W	39.0	50.0	0.025	26.3	AB
F125W	42.0	50.0	0.021	26.3	AB
F125W	45.0	50.0	0.008	26.3	AB
F125W	48.0	50.0	0.010	26.3	AB
F125W	51.0	50.0	0.006	26.3	AB
F125W	54.0	50.0	0.013	26.3	AB
F125W	57.0	50.0	0.025	26.3	AB
F125W	60.0	50.0	0.012	26.3	AB
F125W	63.0	50.0	0.024	26.3	AB
F125W	66.0	50.0	0.007	26.3	AB
F125W	69.0	50.0	0.019	26.3	AB
F125W	72.0	50.0	0.025	26.3	AB
F125W	75.0	50.0	0.024	26.3	AB
F160W	0.0	50.0	0.009	26.0	AB
F160W	3.0	50.0	0.028	26.0	AB
F160W	6.0	50.0	0.027	26.0	AB
F160W	9.0	50.0	0.009	26.0	AB
F160W	12.0	50.0	0.008	26.0	AB
F160W	15.0	50.0	0.011	26.0	AB
F160W	18.0	50.0	0.025	26.0	AB
F160W	21.0	50.0	0.014	26.0	AB
F160W	24.0	50.0	0.027	26.0	AB
F160W	27.0	50.0	0.024	26.0	AB
F160W	30.0	50.0	0.025	26.0	AB
F160W	33.0	50.0	0.019	26.0	AB
F160W	36.0	50.0	0.015	26.0	AB
F160W	39.0	50.0	0.026	26.0	AB
F160W	42.0	50.0	0.022	26.0	AB
F160W	45.0	50.0	0.026	26.0	AB
F160W	48.0	50.0	0.029	26.0	AB
F160W	51.0	50.0	0.025	26.0	AB
F160W	54.0	50.0	0.017	26.0	AB
F160W	57.0	50.0	0.008	26.0	AB
F160W	60.0	50.0	0.028	26.0	AB
F160W	63.0	50.0	0.010	26.0	AB
F160W	66.0	50.0	0.006	26.0	AB
F160W	69.0	50.0	0.016	26.0	AB
F160W	72.0	50.0	0.027	26.0	AB
F160W	75.0	50.0	0.019	26.0	AB

Table 4. Table A3: The Gaussian luminosity functions for each SN type, presented in table 3 of Graur et al. ( 2014 ) .

Type	$M_R$	$\sigma$	Source
Ia	-19.37	0.47	(Wang et al., 2006)
Ib	-17.90	0.90	(Drout et al., 2011)
Ic	-18.30	0.60	(Drout et al., 2011)
IcBL	-19.00	1.10	(Drout et al., 2011)
II-P	-16.56	0.80	(Li et al., 2011)
II-L	-17.66	0.42	(Li et al., 2011)
IIn	-18.25	1.00	(Kiewe et al., 2012)

Table 5. Table A4: The Gaussian luminosity functions for each SN type, adapted from table of Richardson et al. ( 2014 ) .

Type	$M_B$	$\sigma$
Ia	-19.25$\pm$0.20	0.50
Ib	-17.45$\pm$0.33	1.12
Ic	-17.66$\pm$0.40	1.18
IIb	-16.99$\pm$0.45	0.92
II-L	-17.98$\pm$0.34	0.86
II-P	-16.75$\pm$0.37	0.98
IIn	-18.53$\pm$0.32	1.36

Equations

$$F(t)=A\frac{e^{-(t-t_0)/{\tau }_{\mathrm{fall}}}}{1+e^{-(t-t_0)/{\tau }_{\mathrm{rise}}}}+B,$$

$$t_{\rm{max}}=t_{0}+\tau_{\rm{rise}}\ln\Big{(}\frac{\tau_{\rm{fall}}}{\tau_{\rm{rise}}}\Big{)}$$

$$f(t_{\max })=A{x}^x(1-x{)}^{1-x}+B;x=\frac{{\tau }_{\mathrm{rise}}}{{\tau }_{\mathrm{fall}}}$$

$$R_E=\sqrt{\frac{4GM}{c^2}\frac{D_{ls}{D}_s}{D_l}},$$

Full text

1 Introduction

The theory to enable the use of a gravitationally lensed supernova (SN) resolved into multiple images as a cosmological tool was developed in the seminal work of Refsdal (1964). As the light for each of the multiple images follows a different path through the expanding universe and through the lensing potential, the SN images appear delayed by hours (for galaxy-scale lenses) or years (for cluster-scale lenses). These time delays are sensitive to various cosmological parameters, enabling new measurements of the Hubble Constant $H_{0}$ and the dark energy equation of state (Linder, 2011; Treu & Marshall, 2016). Time delay cosmography has been employed extremely successfully for decades using multiply-imaged quasars (e.g., Vuissoz et al., 2008; Suyu et al., 2010; Tewes et al., 2013b; Bonvin et al., 2017; Birrer et al., 2018), but using gravitationally lensed SNe with multiple images (hereafter glSNe) can extend and enhance this method for several reasons.

glSNe have a strong luminosity peak and occur on short time-scales, enabling relatively simple and accurate time delay measurements (e.g. Woosley et al., 2007; Sanders et al., 2015). In addition, the intrinsic luminosities of Type Ia and certain Type II SNe can be inferred independently of lensing (Phillips, 1993; Hamuy & Pinto, 2002; Poznanski et al., 2009; Kasen & Woosley, 2009), which helps to determine the absolute lensing magnifications. This piece of information, not available for lensed quasars, provides an independent check on the lens model and the characterization of the line of sight, thus helping to break the mass-sheet degeneracy (e.g., Holz, 2001; Rodney et al., 2015), which will otherwise introduce an additional uncertainty on $H_{0}$ (Kolatt & Bartelmann, 1998; Xu et al., 2016).

The first multiply-imaged core-collapse (SN Refsdal: Kelly et al., 2015) and Type Ia (SN iPTF16geu: Goobar et al., 2016) SNe have been discovered in just the past few years. However, it is anticipated that the next generation of telescopes, particularly the Large Synoptic Survey Telescope (LSST) and the Wide-Field Infrared Survey Telescope (WFIRST), will provide hundreds to thousands of Type Ia and CC glSNe observations over the next decade (Oguri & Marshall, 2010; Goldstein & Nugent, 2017). With this enormous volume of observed glSNe it will become impractical to employ such individual treatments of time delay measurements seen in Rodney et al. (2016) and Goobar et al. (2016). Incidentally, the time delay measurement techniques employed in these cases ignored the effects of “microlensing”, whereby each SN image is separately varied by lensing effects from stars in the lens plane (e.g., Dobler & Keeton, 2006; Bagherpour et al., 2006; Foxley-Marrable et al., 2018). This stellar-scale lensing operates at micro-arcsecond scales and can introduce uncertainty in the time delay measurement of $\sim 4$%, or $\sim 2.4$ days for a time delay of 60 days. (Goldstein et al., 2018b).

In this work we have developed an open source software package, written in Python, called Supernova Time Delays (SNTD)111SNTD. The SNTD package is capable of making accurate time delay measurements for glSNe of all types while including treatments of microlensing, and can immediately produce accurate simulations for wide-field time-domain surveys such as LSST and WFIRST. Section 2 describes the simulation toolkit in SNTD including updated light curve models, the treatment of microlensing, and examples of SNTD-generated light curves using hypothetical lensed SN case studies. Section 3 describes the time delay measurement capabilities of SNTD. This section uses the case-study lensing systems from section 2 to assess the accuracy and precision of these measurements, and includes a comparison of time delay measurement accuracy for before-peak and after-peak SN detections.

2 The SNTD Simulation Toolkit

The simulation toolkit within SNTD allows a user to simulate realistic glSN observations with a variety of SN and lensing properties for various SN classifications. This capability can be used for design and optimization of surveys and follow-up campaigns. The software used inside SNTD for generating light curves and the ways in which it has been modified for use with glSNe are described in sections 2.1 and 2.2. Microlensing is described in section 2.3, both how it can be simulated and its inclusion in light curve simulations. In section 2.4 we use a representative lensed system to create simulated glSN observations, with microlensing included, to be analyzed for time delays in section 3.

2.1 SNCosmo and Current Light Curve Templates

As discussed in section 1, SNe are an advantageous tool for time delay cosmography due to their light curve shape and time scales. Unlike the stochastic and heterogeneous light curves of AGNs, we have well-defined spectrophotometric time-series to describe a variety of SN types, with coverage from the ultraviolet (UV) to the near-infrared (near-IR). Within SNTD, these spectrophotometric time-series are used to fit observed glSN light curves, using the SNCosmo222SNCosmo Version 1.6 python toolkit (Barbary, 2014). SNCosmo has dozens of empirically defined time-series to describe the evolution of SNe Ia, Ib, Ic, II-P, II-L, and IIn, as well as the parametric SALT2 model for SNe Ia (see full table in the apendix). These SNCosmo models have coverage in the UV to optical wavelength range, and were extended by Pierel et al. (2018) to include more of the UV and near-IR.

SNCosmo is a powerful tool for supernova cosmology, enabling users to fit existing models to observed data at any redshift while considering the effects dust extinction from any point along the line of sight. The toolkit is also capable of using any of the spectrophotometric time-series and fitting capabilities to produce realistic simulations of observed light curves. However, the toolkit is not designed for the analysis of lensed SNe, having no capabilities for handling macrolensing, microlensing, or glSNe. SNCosmo is also not equipped for the specific case of SN Refsdal, the first glSN discovered, which happened to be an unusual SN that cannot be fully described by existing SN classification templates (e.g. Arnett et al., 1989; Kelly et al., 2016).

2.2 Extending SNCosmo

As SNCosmo is already widely used by the SN research community, we have opted to extend the package to provide optimized glSN capabilities that work within SNTD, as opposed to creating an entirely new framework. The updates made are as follows:

Parametric light curve model: while not uniquely useful to studies of glSNe, an important addition to the capabilities of SNCosmo is the inclusion of a parameterized model that operates within the framework of the existing package. This parameterized model allows a user to fit photometric data without knowledge of the SN classification or redshift, enabling time delay measurements in the cases where there is no spectroscopic follow-up or the data are not well fit with common light curve templates (e.g. Rodney et al., 2016).

A common practice for fitting photometric data with no underlying light curve template involves the use of splines or other flexible functions (i.e. Tewes et al., 2013a), but the simplicity of SN light curves enables the use of a parameterized model that has fewer free parameters. The parameterized light curve model defined by Bazin et al. (2009) is implemented in SNCosmo for use in SNTD, which calculates flux as a function of time in such a way as to be general enough to fit any light curve shape. While useful for fitting light curves in certain situations, this parameterized model describes only the shape of a SN light curve; it does not have any information about the SN SED, so it cannot leverage information about the color of an observed SN. The Bazin et al. (2009) model is defined according to:

[TABLE]

where $\tau_{rise}$ and $\tau_{fall}$ characterize the relative rise and decline times of the light curve, while $A$ and $B$ are constants. By taking the derivative, it is apparent that the time of peak and peak flux are:

[TABLE]

[TABLE]

While this model is extremely flexible, it still mimics the smooth rise and fall of a general SN light curve, which helps to keep a light curve fit relatively insensitive to the effects of sharp microlensing features.

Microlensing: we describe the phenomenon of microlensing, as well as SNTD’s method for simulating it, in section 2.3. Once the simulation is complete however, the effect can be simply described by a time-evolved “magnification curve”. This curve is added to a SN light curve as an SNCosmo “propagation effect,” in the same way as dust extinction is handled. These effects are simply models of how intervening structures affect a spectrum, as well as the associated light curve. The microlensing propagation effect causes time-dependent flux-scaling when applied to the simulated light curve points. Each image of a glSN will be along a different line of sight, subjecting it to unique microlensing effects. Therefore, during simulations, a separate microlensing effect is applied to each light curve as it is generated (See section 2.3 and 2.4).

Multiply-Imaged Supernovae: knowing that observations of glSNe stem from a single light curve is extremely valuable in measuring time delays (see section 3). As SNCosmo has no structures for handling this phenomenon, SNTD contains new objects capable of maintaining maximum information about the SN and each of its images for use in time delay measurements. All components of the new objects operate within the current SNCosmo framework so that, for example, each SN image has a usable light curve within existing SNCosmo functions.

2.3 Microlensing

The phenomenon of microlensing occurs when light rays from the expanding photosphere of a SN pass within the lensing potential of a set of stars in the lens plane, which have Einstein radii on the order of micro-arcseconds (Dobler & Keeton, 2006). This causes fluctuations in the light curve of $\sim 0.2$ to $>0.5$ magnitudes on timescales of weeks to months, which can make it difficult to obtain an accurate time delay measurement without accounting for microlensing in certain lens configurations (Dobler & Keeton, 2006; Goldstein et al., 2018b; Foxley-Marrable et al., 2018; Bonvin et al., 2019).

To simulate the effects of microlensing, we employ the microlens inverse ray-tracing code (Wambsganss, 1999). The microlens code defines random realizations of stars modeled as point-masses in the lens plane according to two parameters: $f_{*}\ \rm{and}\ q$. The parameter $f_{*}$ defines the relative fraction of convergence ($\kappa$) from discrete sources such as stars compared to smooth sources such as dark matter. The parameter $q$ is the stellar mass ratio ($m_{min}/m_{max}$) of the stellar mass function. The microlensing magnification distributions produced by microlens require a further two parameters, namely the shear due to the Singular Isothermal Ellipsoid (SIE) and external potentials ($\gamma$), and the local convergence ($\kappa$). The convergence causes a focusing of light rays, leading to an isotropic magnification of a source, while the shear introduces an anisotropy that distorts the shape of the magnified source (Narayan & Bartelmann, 1996).

Once the microcaustic field is generated by the microlens code for a particular set of parameters, it must be convolved with a model for the expanding photosphere of a SN. By default, we create a model for an expanding SN photosphere as flat-disk achromatic Gaussian brightness distribution. This photospheric model is a simplifying assumption also made in previous microlensing analyses (e.g. Foxley-Marrable et al., 2018), which other work has shown to be sufficiently accurate up to an ignored “microlensing time delay” on the order of 0.1 days (Bonvin et al., 2019). The microcaustic field defines a grid of positive and negative magnifications in the source plane. SNTD can produce a magnification curve—describing the effect of any given microlensing realization with time—by propagating the expanding photospheric model through the microcaustic field (Figure 1). As this default model is achromatic, the resulting magnification curve will be a function of only time. In principle SNTD is capable of replacing this simple model with a 2-D projected specific intensity profile, like that of Goldstein et al. (2018b). This method would better represent a SN atmosphere, which emits light differently as a function of radius from the center (e.g. Kasen, 2006) in such a way that a time delay effect can also be observed (e.g. Bonvin et al., 2018). This substitution could be made by a user, or may perhaps be the default in a future iteration of SNTD. The SN image can be placed randomly at any position in the microcaustic field, so that one can explore the variety of microlensing effects that may arise for a given set of lensing parameters ($\kappa$, $\gamma$, $f_{*}$, $q$). The size of each microcaustic field projected onto the source plane is $10R_{E}$ on a side, where $R_{E}$ is the Einstein radius of a typical deflector of mass $M$ in the lens plane:

[TABLE]

where $D_{s},\ D_{l},\ \rm{and}\ D_{ls}$ define the angular diameter distances to the source, lens, and between the lens and source, respectively.

2.4 SNTD Simulation Examples

In order to demonstrate the simulation and time delay measurement capabilities of SNTD, we created 2,000 simulated glSNe using an imagined galaxy-scale strong-lensing system that is similar to the lenses described by Shu et al. (2018). It should be noted that a true SN detection in such a galaxy would require excellent follow-up observations from space and/or large ground-based observatories to obtain sufficient data quality for time delay cosmography. This mock lensing galaxy is placed at redshift ($z_{\rm{lens}}$) of 0.5, and the source galaxy is placed at a redshift ($z_{\rm{source}}$) of 1.4. The overall magnification factor ($\mu$) for the leading image of a SN located at the light peak of the lensed galaxy is set to 5, and the relative time delay ($\Delta t$) between the two images is set to 60 days. The measurement of time delays for these mock lensed SNe is described in section 3.

Each simulated SN has an absolute magnitude in the B bandpass defined by a Gaussian luminosity function (see appendix for table), a redshift equal to the host galaxy’s $z_{source}$, and host galaxy extinction defined by the dust distributions used by Rodney et al. (2014). Magnification and time delay factors are applied to each SN image, defined by the lensing galaxy’s $\mu$ and $\Delta t$ described above, respectively.

A microcaustic field is randomly generated at a redshift equal to the lensing galaxy’s $z_{lens}$ described above, following the methodology described in section 2.3 using lensing parameters representative of the lenses presented by Shu et al. (2018). The resulting microlensing magnification curve is applied to each simulated light curve as an added propagation effect, causing positive and negative flux variations with time (Figure 1).

To simulate realistic light curves the observational strategies presented by Shu et al. (2018) are adopted, which describe a survey strategy that targets lensed galaxies with high star formation rates and prior knowledge of the lensing parameters. The photometric data are generated with observational parameters for the HST F125W and F160W filters and a fixed S/N curve based on the HST exposure time calculator. We simulate 1000 glSNe for which the SN is first detected before peak brightness, and another 1000 glSNe for which the first detection lands after peak brightness. With these two sets of simulations we can ask how much the precision of a SN time delay measurement is degraded when there is no detection before peak brightness. This will be an important consideration for the choice of cadence and depth in future surveys, and also for the allocation of resources for follow-up observations.

An example of an observations table used by SNCosmo to realize a set of light curves is given in the appendix. SNTD simulates the multiple images of a single SN by using the same absolute magnitude and host galaxy dust parameters discussed above for each light curve as they are inherent to the SN, then scaling the simulated flux measurements and shifting the time axis of the trailing image by the prescribed time delay and magnification parameters. Finally, the microlensing magnification curve described above and observational noise are applied to the simulated observations, causing further variations in the flux with time.

Before Peak Detections: We first generated light curves under the assumption that we were able to “observe” each leading SN image before its epoch of peak brightness. It is ideal to obtain this before-peak observation, so that the peak of the light curve is well-defined and more easily fit. A SN is simulated using the methods outlined above in this section, then a random epoch is chosen between the epoch at which the brightness of the SN reaches the detection threshold of the telescope, and the epoch of its peak brightness. This epoch defines the first observation of the leading image.

Once the SN is “observed”, we assume a set of non-disruptive follow-up observations are made by the HST. This necessitates a 2-week gap between the first epoch and the second, and thereafter a 5-day cadence is used to “observe” the SN until the trailing image is no longer detectable (Figure 2A). This process is repeated to create a mock catalog of 1000 realistic glSN light curves.

After Peak Detections: The process of generating a realistic catalog of 1000 glSN light curves is repeated under the assumption that we only detect each leading SN image after its epoch of peak brightness. It is more difficult to identify the peak of a light curve, and therefore the relative time delay, without any before-peak observations due to the lack of a primary inflection point.

The light curves for the glSNe detected after peak are generated by taking the first observation to be a random epoch between the epoch of peak brightness and the epoch at which the brightness of the SN falls below the detection threshold of the telescope. This detection epoch simultaneously defines the first observation of the trailing image, which is nearly always before peak in this case due to the relative time delay of 60 days. Then we again assume a set of non-disruptive follow-up observations are made by the HST. This results in a 2-week gap between the first epoch and the second, and thereafter a 5-day cadence is used to “observe” the SN until the trailing image is no longer detectable (Figure 2B).

3 The SNTD Time Delay Measurement Toolkit

We describe SNTD’s treatment of microlensing in section 3.1, and then the various tools available to measure the time delays and relative magnifications of glSNe. The methods include parallel analysis (Section 3.2), series analysis (Section 3.3), and color curve analysis (Section 3.4). Next, the simulated samples of glSNe created in section 2.4 are used to compare these three time delay measurement techniques.

3.1 Accounting for the Effects of Microlensing

Identifying and correcting for microlensing in a light curve is an extremely difficult problem, driven mainly by the stochasticity of the effect and the degeneracies present between the extrinsic and intrinsic model parameters. A pathological but plausible case of microlensing is capable of leaving the shape of a light curve relatively untouched, while shifting the peak in time (Figure 3). The effect of microlensing from a modeling perspective would be classified as low in this case, as the model performs well in fitting the data, but there is an added systematic error unaccounted for in the measurement of the time of peak that will propagate through to an error in the time delay.

Therefore, we seek to better characterize the systematic error in the time delay measurement due to microlensing, which if unaccounted for would cause an underestimation of the time delay measurement uncertainty. When not in the microlensing regime described in figure 3, fitting light curves for time delays without considering microlensing remains quite effective, due to the existence of templates for various SN classifications that reduce the flexibility of models. In these cases, the main deviations from the model should be due to microlensing effects, which are then visible in the residuals. However, simply fitting a curve to the residuals and attributing all of those deviations to microlensing is inadequate, because it will create a bias in your time delay measurement if your initial fit is poor, or may incorrectly attribute observational scatter to a microlensing effect (e.g. Tewes et al., 2013a). Therefore, we avoid directly fitting for a microlensing spline and instead use monte carlo simulations of realistic microlensing magnification curves to characterize the additional error on the time delay measurement:

A best fit model for the intrinsic SN light curve is found by fitting observed light curve data with a model that has no microlensing included (see section 3.2-3.4). This fit defines a measured time delay $\Delta t$. 2. 2.

Trends in the residuals of the best fit model are identified using Gaussian Process Regression (GPR), which produces a posterior distribution of possible microlensing functions. 3. 3.

A representative sample of N (typically 100) microlensing curves is then extracted from the GPR posterior (Figure 4). 4. 4.

The photometric data is then refit once for each GPR posterior sample, with each microlensing curve applied to the best fit model as a microlensing propagation effect (section 2.3). The flux calculated by the model is affected by the microlensing propagation effect, and otherwise the same model parameters that are measured in step 1 are allowed to vary.

The result of these steps is a time delay measurement error for each microlensing scenario, $\Delta t-\Delta t_{i}^{\mu}$, the combination of which defines a distribution of time delay errors due to microlensing. This distribution of microlensing errors characterizes an additional uncertainty on the time delay measurement, which is reported separately and then combined in quadrature with the model uncertainty. Including this measurement of uncertainty due to microlensing significantly improves the overall time delay uncertainty characterization, which is otherwise consistently underestimated when microlensing is ignored (Figure 5). While including the uncertainty due to microlensing certainly improves the overall error characterization, its deviations from Gaussian at $|\Delta_{\rm{True}}-\Delta_{\rm{Measured}}|/\sigma\gtrsim 1.0$ arise from extreme cases of microlensing. The agreement could likely be improved by increasing the number of draws from the microlensing GPR posterior, which was limited due to computation time to 100 per SN, or by combining the fitting methods described in sections 3.2 and 3.4 (see section 4).

3.2 Parallel Analysis

If the separate light curves of a multiply-imaged SN are each reasonably well-sampled, then SNTD is capable of making time delay measurements by fitting each of the light curves in parallel. The parameters for a given model are separated into those that are inherent to the SN and therefore constant across each image (i.e. redshift, light curve shape parameters, host galaxy extinction, etc.), and those that are impacted by lensing and therefore will be unique to each image (i.e. amplitude, time of peak, etc.). These parameter sets are denoted $\bf\theta_{SN}$ and $\bf\theta_{L}$, respectively.

We use the method of nested sampling to efficiently sample the full parameter space for $\bf\theta_{SN}$ and $\bf\theta_{L}$. Nested Sampling is effectively a Markov Chain Monte Carlo (MCMC) simulation that provides an estimate of the posterior probability distribution, but nested sampling additionally allows one to calculate the integral of the distribution and is more robust for finding global maxima (Skilling, 2004). SNCosmo contains a modified nested sampling algorithm, which SNTD uses to identify a best fit model for each separate image. The default is to employ uniform prior distributions,but any form for the prior probability could be used. The algorithm takes the best fit models identified separately, obtains a joint posterior for parameters in $\bf\theta_{SN}$, and marginalizes over the parameters in $\bf\theta_{L}$. The result is a model for each image, where the parameters in $\bf\theta_{SN}$ are the same for each image and the parameters in $\bf\theta_{L}$ are likely different. Particularly when fitting light curves with little information, this gives some constraints on the $\bf\theta_{SN}$ parameters by using the fact that the underlying light curve for each image must be the same, so any observed differences are only from extrinsic lensing effects like magnification, microlensing, and a relative time delay.

Once the best-fit light curve parameters are identified, each image is fit once more with the $\bf\theta_{SN}$ parameters varying in a tightly bound range around the values measured above, while the $\bf\theta_{L}$ parameters remain unrestricted. This results in fine-tuned fits to each individual image (Figure 6).

3.3 Series Analysis

SNTD contains a second method for determining time delays that is particularly useful in the cases of relatively sparse light curves. This method differs from that of parallel fitting in section 3.2 in that the intrinsic light curve shape is identified first by fitting all or a subset of the light curves, which is then used to measure the time delays of the multiple images simultaneously. In both cases, we utilize the knowledge that the underlying light curve is the same for every image, but in this case we apply that leverage more directly by employing only a single set of intrinsic SN light curve parameters $\theta_{SN}$ (see section 3.2) for the light curve model, in order to get simultaneous constraints on the relative time delays and magnifications for every image.

To perform the series time delay measurement, SNTD effectively employs a double layered MCMC simulation that instead uses nested sampling to estimate the posterior probability distributions (see section 3.2). The algorithm varies the $\bf\theta_{SN}$ parameters in the outer layer and the relative magnifications and time delays in the inner layer, attempting to simultaneously shift (time) and scale (flux) the photometric data from all images and fit it with a single varying light curve model (Figure 7). As in the parallel approach, the default for this customized “double-nested” sampling is to use uniform priors for all parameters, though in principle any informative priors could be used.

3.4 Color Curve Analysis

It has been proposed that the effect of microlensing on (at least) SNe Ia contains an achromatic period, meaning that it is not wavelength-dependent during this time (Goldstein et al., 2018b). Therefore color curves, which define the observed broad-band photometric color as a function of time, might be an effective way of measuring time delays and minimizing the complications of microlensing, as each filter should be affected identically. However, this method is only useful if the SN images are observed with multi-band photometry, the color curves are not featureless, and the achromatic phase is a true effect.

In cases where the SN classification is well-known, existing templates can be used to obtain expected color curves, which can then be used to constrain the functions being used to fit the data (Figure 8). Otherwise, we must resort to flexible functions such as splines or Chebyshev polynomials to fit the data. In either case, we can still leverage the fact that the color curve for each image of the SN should be the same, as in section 3.3. Therefore, the same “double-nested” sampling algorithm is employed to obtain time delays from color curves (see section 3.3). The algorithm once again attempts to shift (time) and scale (flux) the color curve data for all the SN images, while fitting the combined data with a single varying color curve model. SNTD allows the color curve flux data to be individually scaled in order to account for chromatic parameters in $\bf\theta_{L}$ (see section 3.2) that will affect each image differently, such as dust extinction or non-achromatic microlensing.

4 SNTD: Case Studies

We measure the relative time delays of the 2,000 simulated glSNe described in section 2.4, exploiting our knowledge of the SN classification by using SED templates to perform the light curve fitting. The simulations are separated into SNe “discovered” before peak and after peak brightness. We obtained three different sets of results for each case by using the methods described in sections 3.2-3.4 to measure the time delays, which are summarized in table 4 and figure 9. All three fitting methods measured with nearly identical accuracy and for the before-peak case, which is to be expected for these simulations as the light curves were sufficiently well-sampled to use any approach. Still, the fact that achromatic microlensing (used for these simulations, see section 2.3) does not affect the color curve method led to a higher precision. The series fitting routine did perform slightly better in the after-peak scenario when compared to the parallel routine, reflecting that it can be more robust in cases where one or more images has a lower quality light curve. Nevertheless, the accuracy decreased in the after-peak case study for both of these techniques, the systematic offset stemming from the lack of an inflection point to identify peak brightness, leading to an underestimation of the time delay. Likewise the precision for each technique is reduced for the after-peak case.

The difference between the two cases for the color curve routine is much less pronounced, corresponding to the relative independence of a color curve on peak brightness. While the lack of before-peak data does not directly harm the color curve fitting, it still leads to a decrease in accuracy as most of the simulated leading-image light curves only contained 3 or 4 data points before falling below the detection threshold (Figure 2). The performance of the color curve fitting method is somewhat idealized here however, as the simulated microlensing is achromatic and therefore completely absent in each color curve. It’s worth noting that the color curve routine’s relative insensitivity to peak brightness and comparatively small number of catastrophic outliers may suggest that combining the color curve fitting routine with either the parallel or series fitting methods could enable more accurate measurements. Such a reduction in the tails of figure 9A-B should also lead to an improvement in the uncertainty characterization described in section 3.1 and figure 5.