Toward Rate Estimation for Transient Surveys. I. Assessing Transient Detectability and Volume Sensitivity for iPTF

The iPTF real-time image subtraction pipeline (henceforth ISP) was hosted at the National Energy Research Scientific Computing Center (NERSC). A complete exposure of 11 working CCDs was transferred to NERSC immediately after data acquisition to search for new candidates. The pipeline preprocessed the images to remove bias and correct for flat-fielding. It solved for astrometry and photometry, and performed image subtraction using the HOTPANTS algorithm (Becker 2015). New candidates were assigned a real–bogus classification score between 0 and 1 corresponding to bogus and real, respectively (Bloom et al. 2013). Additionally, candidates would be cross-matched to external catalogs to remove asteroids, active galactic nuclei, and variable stars.

We follow the clone stamping technique used by Frohmaier et al. (2017) for PTF to perform our fake point-source injections. The parameters describing these fake transients are single epoch—they represent the intrinsic properties of the object and observing conditions at a particular epoch. In other words, here we assess the detectability given the transient was in the field of view of the instrument.

The computational cost for performing injections into all iPTF images and running ISP on them is significant. Therefore, we carry out the process in a single iPTF field 100019. We choose this field since the distribution of the transient population in this field is an accurate representation of the transient population in the sky observed from Palomar (see Figure 1 of Frohmaier et al. 2017).

The fake injections are bright stars chosen from each original image. These are objects having the following properties:

$$\begin{eqnarray}\begin{array}{rcl}{m}_{* } & \in & [13.5,16];\quad {\mathtt{CLASS}}\_{\mathtt{STAR}}\in [0.5,1.0]\\ {\mathtt{FWHM}} & \in & [1.0,3.0];\quad {\mathtt{ELLIP}}\in [0.0,0.3].\end{array}\end{eqnarray} \tag{ 1 }$$

Here ${m}_{* }$ is the apparent magnitude, and ${\mathtt{CLASS}}\_{\mathtt{STAR}}$ is a quantity having a value between 0 (not star-like) and 1 (star like). ${\mathtt{FWHM}}$ is the full width at half maximum, in pixels. ${\mathtt{ELLIP}}$ is the ellipticity of the object. These quantities are reported after running SExtractor (Bertin & Arnouts 1996) on the original images. The reason we choose objects in this range is because we want the point-spread function (PSF) to be well estimated, which is the case for bright stars having a high signal-to-noise ratio ≳100 (${m}_{* }\leqslant 16$). At the same time we want to avoid pixel saturation and therefore select stars with ${m}_{* }\geqslant 13.5$. Objects falling in a 50 pixel wide edge boundary are left out since they could potentially be affected by image subtraction artifacts.

A square of side length ∼9'',⁵ centered around the star and local-background subtracted, constitutes a stamp. A stamp containing any other object apart from the source star is avoided. The local background refers to that reported by SExtractor. The stamp is scaled by an appropriate scaling factor to create a point-source transient of desired magnitude. Each transient is allocated a host galaxy.⁶ We follow Frohmaier et al. (2017) regarding the location in the host and place our stamp at a random pixel location within an elliptical radius⁷ of 3 pixels. This value contains a sufficient amount of the flux from the galaxy.

This procedure is performed on all the images in field 100019 of iPTF, 10-fold, with a total of ≈2.24 × 10⁶ injected transients. The transient magnitudes are chosen uniformly between the 15th and 22nd magnitudes with the constraint that the stamp is one magnitude fainter than the original star. We only rescale to fainter magnitudes because we do not want artifacts like noise residuals from the average background subtraction to be scaled up as noise spikes. Therefore, ${m}_{\mathrm{inj}}$ is as follows:

$$\begin{eqnarray}{m}_{\mathrm{inj}}\sim \left\{\begin{array}{ll}U(15,22); & {m}_{* }\in (13.5,14)\\ U({m}_{* }+1,22); & \mathrm{otherwise}\end{array}\right..\end{eqnarray} \tag{ 2 }$$

An example of an injected transient in a galaxy and the new object recovered by the ISP is shown in Figure 1.

Zoom In Zoom Out Reset image size

The recovery efficiency ε is defined as the ratio of the number of injections recovered in a part of the parameter space to the total number of injections in that part. Let our injections be described by parameters ${\boldsymbol{\lambda }}$, then

$$\begin{eqnarray}&&\varepsilon ({\boldsymbol{\lambda }})=\displaystyle \frac{{N}_{\mathrm{rec}}({\boldsymbol{\lambda }})d{\boldsymbol{\lambda }}}{{N}_{\mathrm{tot}}({\boldsymbol{\lambda }})d{\boldsymbol{\lambda }}}.\end{eqnarray} \tag{ 3 }$$

The quantities in the numerator and denominator are the number of recovered and total injections, respectively, $\in \left({\boldsymbol{\lambda }},{\boldsymbol{\lambda }}+d{\boldsymbol{\lambda }}\right)$. Here ${\boldsymbol{\lambda }}$ includes both intrinsic source properties of the transient and its environment along with the observing conditions. Examples of intrinsic properties include the magnitude of the transient and the surface brightness of the host galaxy, whereas those for observing conditions include airmass or sky brightness. While we control fake transient brightness, the observing conditions are those of the images themselves. Since images across the full survey time are used, the parameter space of the observing conditions is automatically spanned.

We determine recovery based on the spatial cross-matching of the injections with new objects reported after running the ISP. To determine the tolerance to be imposed during the cross-matching, we define ${{\rm{\Theta }}}_{\mathrm{IQ}}$ as

$$\begin{eqnarray}&&{{\rm{\Theta }}}_{\mathrm{IQ}}=\displaystyle \frac{\sqrt{{({x}_{\mathrm{inj}}-{x}_{\mathrm{rec}})}^{2}+{({y}_{\mathrm{inj}}-{y}_{\mathrm{rec}})}^{2}}}{{\rm{\Phi }}},\end{eqnarray} \tag{ 4 }$$

where ${{\rm{\Theta }}}_{\mathrm{IQ}}$ is the distance between the injected and the recovered sources in units of the seeing, Φ.

We choose the threshold of ${{\rm{\Theta }}}_{\mathrm{IQ}}$ such that 99% of the found injections lie within this threshold, which has a value of ${{\rm{\Theta }}}_{\mathrm{IQ}}^{99 \% }=0.48$ (see Figure 2). We also impose a real–bogus score threshold ${\mathtt{RB}}{\mathtt{2}}\geqslant 0.1$ on the new object. This threshold on ${\mathtt{RB}}{\mathtt{2}}$ is inspired from survey operation thresholds. Out of the ≈2.24 × 10⁶ injections, we recover ≈1.62 × 10⁶.

Zoom In Zoom Out Reset image size

The single-parameter efficiency is the marginalized version of Equation (3). Suppose our parameter of interest is θ and the other "nuisance" parameters are given by ${\boldsymbol{\gamma }}$, such that in Equation (3), ${\boldsymbol{\lambda }}=\{\theta ,{\boldsymbol{\gamma }}\}$. The single-parameter efficiency is

$$\begin{eqnarray}&&\varepsilon (\theta )=\displaystyle \frac{\left[{\int }_{{\boldsymbol{\gamma }}}{N}_{\mathrm{rec}}(\theta ,{\boldsymbol{\gamma }})d{\boldsymbol{\gamma }}\right]d\theta }{\left[{\int }_{{\boldsymbol{\gamma }}}{N}_{\mathrm{tot}}(\theta ,{\boldsymbol{\gamma }})d{\boldsymbol{\gamma }}\right]d\theta }.\end{eqnarray} \tag{ 5 }$$

In Figure 3 we show the single-parameter efficiencies. The expected trend of missing faint transients is seen in the plot for ${m}_{\mathrm{inj}}$. We find that the recovery efficiency starts to drop for transients by the 20th magnitude and sensitivity is almost nil by the 22nd magnitude.

Zoom In Zoom Out Reset image size

In this section, we make a selection of parameters from the full parameter set, ${\boldsymbol{\lambda }}$, to those on which the detectability depends strongly. In other words, the detectability is a multivariate function of all the possible parameters that influence the detection of a transient. We identify the minimal set that captures maximum variability. There can be correlations among a pair of parameters. For example, the sky brightness, F_sky, and the limiting magnitude, ${m}_{\mathrm{lim}}$, are correlated—a bright sky hinders the depth and results in a low value of limiting magnitude. The variation of the marginalized efficiencies shown in Figure 3 assist us with the choice of such a parameter set. Since the trend in the single-parameter efficiencies are similar to those from PTF, we select the parameters considered by Frohmaier et al. (2017) with a minor difference in the usage of the galaxy surface brightness directly, as used in Frohmaier et al. (2018), in place of the F_box⁸ parameter used in the former. This is justified because our fakes were injected in galaxies.

We choose, the following set to represent the dependence of detectability:

$$\begin{eqnarray}&&{\boldsymbol{\beta }}=\{m,{S}_{\mathrm{gal}},{F}_{\mathrm{sky}},{{\rm{\Phi }}}_{\mathrm{IQ}},{m}_{\mathrm{lim}}\}.\end{eqnarray} \tag{ 6 }$$

Here m is the apparent magnitude of the transient, ${S}_{\mathrm{gal}}$ is the host galaxy surface brightness, ${F}_{\mathrm{sky}}$ is the sky brightness, ${{\rm{\Phi }}}_{\mathrm{IQ}}$ is the ratio of the astronomical seeing to that of the reference image, and ${m}_{\mathrm{lim}}$ is the limiting magnitude. The quantities m and ${S}_{\mathrm{gal}}$ are natural in capturing detectability. Sky brightness affects the detectability in a strong way, as is apparent from Figure 3. The ${{\rm{\Phi }}}_{\mathrm{IQ}}$ parameter captures the variability of the atmosphere. Finally, the limiting magnitude, ${m}_{\mathrm{lim}}$, although correlated with ${F}_{\mathrm{sky}}$, captures longer exposure times and the status of instrument electronics.

With this set, we use the machinery of supervised learning provided by the scikit-learn library (Pedregosa et al. 2011) to train a binary classifier based on the results of the ISP. Once trained, the classifier outputs a probability of detection given arbitrary but physical values of ${\boldsymbol{\beta }}$. We denote this trained classifier by $\hat{\varepsilon }$:

$$\begin{eqnarray}&&\hat{\varepsilon }=\hat{\varepsilon }(m,{S}_{\mathrm{gal}},{F}_{\mathrm{sky}},{{\rm{\Phi }}}_{\mathrm{IQ}},{m}_{\mathrm{lim}}).\end{eqnarray} \tag{ 7 }$$

The scikit-learn library provides a suite of classifiers. We choose the nonparametric KNearestNeighbor classifier based on speed and accuracy given our large volume of training data. Our complete data set comprises ∼2.24 × 10⁶ fake point-source injections of which ∼1.62 × 10⁶ (∼6.2 × 10⁵) are found (missed) by the ISP. We train the classifier using 11 neighbors—twice the number of dimensions plus one to break ties. The observation of a fiducial transient is a point in this parameter space. To decide if that point is "missed" or "found," we use a majority vote from the nearest 11 neighbors. To cross-validate the performance, the data set is split into a training set containing 90% of the full data set, and a testing set containing the remaining 10%. We checked that increasing the number of neighbors does not significantly increase the correctness of predictions made by the classifier. We note that one could use a different threshold for this classification. For example, a different option could be to use more than three "found" neighbors to call the arbitrary point as found. However, it comes at a cost of misclassification. From the predictions of the classifier on the testing set, we find the systematic uncertainty of the classifier to be ≈6% i.e., 6 out of 100 predictions made by the classifier are expected to be either true-negative or false-positive cases. The result does not change much if the size of the training and testing set is varied (see Table 1). A comparison between the predictions made by the trained classifier and the original ISP efficiency with the transient magnitude is presented in Figure 4. We see that the behavior of the ISP is reproduced by feeding the classifier with only a few thousand points randomly chosen from the parameter space.

Zoom In Zoom Out Reset image size

We use SN Ia light curves from the SALT2 model (Guy et al. 2007). In particular, we use the Python implementation of SALT2 provided in the sncosmo library (Barbary 2014). This model is based on observations of SNe Ia by the SDSS and SNLS surveys. The free parameters of the model include the stretch (x₁) and color (C) parameters of the SN Ia. Regarding the range of these parameters, we follow the same range as Frohmaier et al. (2017, see Table 1 and Equation (4) therein). The ranges cover the possible light-curve morphologies of SNe Ia (Betoule et al. 2014). We show an example light curve, at a redshift of z = 0.01 with an intrinsic M_B = −19.05 in Figure 5. When propagating the flux, we also take into account the extinction due to host galaxy dust and the Milky Way (MW) dust. We use the MW dust map by Fitzpatrick (1999), which is a part of the sncosmo package. For the host galaxy extinction, we use the distribution of E(B − V) of SN Ia in their host galaxies (Hatano et al. 1998). Dust extinction plays a significant role in the detectability of light curves as the SNe can be dimmed by as much as 1–1.5 mag.

Zoom In Zoom Out Reset image size

We simulate ≈5 × 10⁶ SN Ia light curves uniformly in a comoving volume up to a redshift, ${z}_{\max }^{\mathrm{Ia}}=0.28$,⁹ uniform in peak time distribution in the observer frame. We assume a flat ΛCDM cosmology with Hubble constant H₀ =69.3 km s⁻¹ Mpc⁻¹ and matter to critical density Ω_m = 0.287 (Hinshaw et al. 2013).¹⁰ We associate a host galaxy surface brightness to each of these SNe using the distribution of surface brightness from iPTF data.

The epochs when the SN Ia is observed come from the iPTF observing schedule. At each observation, we obtain the transient magnitude at that epoch from the light curve and the observing conditions from the iPTF survey database. The single-epoch classifier then tells us the epochs when the transient was detected. An example is shown in Figure 5 where the vertical lines in the upper and lower panels, respectively, represent the observations and detections at each epoch.

To understand rates, one must have a good estimate of the survey sensitivity to particular transient types. Let ${{\rm{\Lambda }}}_{\mathrm{SNe}}$ be the expected count of SNe seen during survey time. Then, with R as the intrinsic rate we have

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Lambda }}}_{\mathrm{SNe}} & = & \int f(t;\mathop{\underbrace{{M}_{B},z,\,...}}\limits_{{\boldsymbol{\kappa }}})\mathop{\overbrace{\displaystyle \frac{{dN}}{{{dt}}_{e}{{dV}}_{c}}}}\limits^{R}\displaystyle \frac{1}{1+z}\displaystyle \frac{{{dV}}_{c}}{{dz}}{dzdtd}{\boldsymbol{\kappa }}\\ & = & R\int f(t;\mathop{\underbrace{{M}_{B},z,\,...}}\limits_{{\boldsymbol{\kappa }}})\displaystyle \frac{1}{1+z}\displaystyle \frac{{{dV}}_{c}}{{dz}}{dzdtd}{\boldsymbol{\kappa }}\\ & = & R\langle {VT}\rangle ,\end{array}\end{eqnarray} \tag{ 8 }$$

where the integral runs over time of observation and comoving volume up to ${z}_{\max }^{\mathrm{Ia}}=0.28$. The selection function, $f(...)\in \{0,1\}$, is to be interpreted as the weight assigned to regions in spacetime. The value of the selection function is a consequence of running a particular instance of SN Ia through the observing schedule and inferring detectability based on the single-epoch classifier in Equation (7). Therefore, the selection function depends on the observer time, t, which captures the duty cycle and cadence. Also, it depends on the intrinsic properties of the supernova like the absolute intrinsic magnitude, M_B, the redshift, z, at which it was simulated, the sky location, and so on. These are collectively represented by ${\boldsymbol{\kappa }}$ in Equation (8). Since we have distributed the supernovae uniformly in comoving volume, the integral is approximated in the Monte Carlo sense:

$$\begin{eqnarray}\begin{array}{rcl}\langle {VT}\rangle & = & \int f(t;\mathop{\underbrace{{M}_{B},z,\,...}}\limits_{{\boldsymbol{\kappa }}})\displaystyle \frac{1}{1+z}\displaystyle \frac{{{dV}}_{c}}{{dz}}{dzdtd}{\boldsymbol{\kappa }}\\ & \approx & \displaystyle \frac{{N}_{\mathrm{rec}}}{{N}_{\mathrm{tot}}}T\int \displaystyle \frac{1}{1+z}\displaystyle \frac{{{dV}}_{c}}{{dz}}{dz},\end{array}\end{eqnarray} \tag{ 9 }$$

where ${N}_{\mathrm{rec}}$ is the number of SNe recovered from this simulation campaign, ${N}_{\mathrm{tot}}$ is the total number simulated, and T is the four-year period of iPTF over which we performed the simulations.¹¹ We obtain the result

$$\begin{eqnarray}&&\langle {VT}{\rangle }_{\mathrm{Ia}}=(2.93\pm 0.21)\times {10}^{-2}\,{\mathrm{Gpc}}^{3}\,\mathrm{yr},\end{eqnarray} \tag{ 10 }$$

where the error includes the $\sim 1/\sqrt{N}$ statistical error from Monte Carlo integration and the 6% systematic error of the single-epoch detectability classifier computed in Section 4.2, the latter being the dominant source of error. The distribution of the detected SNe Ia in sky is shown in Figure 6 colored by redshift. Using the recovered SNe Ia, the median sensitive comoving volume is found to be 0.305 Gpc³. We report the redshift corresponding to this value as the median sensitive redshift to SNe Ia, ${z}_{\mathrm{median}}^{\mathrm{Ia}}=0.099$, shown in Figure 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In contrast to the well-defined Ia light curves with their typical timescales of several weeks, we also wanted to explore longer-timescale light curves as a limiting case. Therefore, we consider SNe IIp and compute their spacetime sensitive volume in similar lines as Section 5.2. In general, type II supernovae (SNe II) vary in light-curve morphology and are categorized in various subtypes (Li et al. 2011). Specifically, type IIp light curves have a distinct "plateau" feature after the rise lasting for about 100 days after explosion, as shown in Figure 8. The intrinsic brightness, M_B ∼ −16.75, is significantly lower than that of SNe Ia (Richardson et al. 2014). Hence, we expect the spacetime sensitive volume to be lower than that of the SNe Ia. When considering the Ia light curves in Section 5.1, the SALT2 model parameters were used to tune possible light-curve morphologies. Here we take a simpler approach and consider a time-series model from Gilliland et al. (1999; named nugent-sn2p in the sncosmo package) to compute the flux up to 100 days from the explosion time. Thus, while simulating the SNe IIp in spacetime, the only change to the light-curve shape is the "stretch" depending on the cosmological redshift.

Zoom In Zoom Out Reset image size

We simulate ∼9.1 × 10⁵ SN IIp light curves uniform in sky location, observer time, and comoving volume up to a redshift z = 0.1. Like the SNe Ia, each SN IIp is assigned a host galaxy surface brightness from the surface brightness distribution of galaxies in iPTF and an E(B − V) extinction value from IIp extinction distribution in Hatano et al. (1998). In this case, we use the criteria that the light curve must be recovered a minimum of five epochs, brighter than the 20th magnitude in a span of 3 weeks within the 100 days postexplosion. The iPTF observing schedule along with the single-epoch classifier is used to compute the detectability in each epoch. We obtain the result

$$\begin{eqnarray}&&\langle {VT}{\rangle }_{\mathrm{IIp}}=(7.80\pm 0.76)\times {10}^{-4}\,{\mathrm{Gpc}}^{3}\,\mathrm{yr},\end{eqnarray} \tag{ 11 }$$

where the error includes the statistical error from the Monte Carlo integration and the 6% systematic uncertainty from the single-epoch classifier (see Section 4.2). The median sensitive redshift is found to be ${z}_{\mathrm{median}}^{\mathrm{IIp}}=0.038$.

The computation of the rate posterior assumes the likelihood of observing N candidate events is an inhomogeneous Poisson process (Loredo & Wasserman 1995; Farr et al. 2015). Our search will filter the SN Ia population based on the model presented in Section 5 at the expense of some contamination from other transient types, potentially with similar light-curve morphologies. If the mean count of these impurities is Λ₀, the likelihood function is

$$\begin{eqnarray}\begin{array}{rcl}p\left(N| {{\rm{\Lambda }}}_{0},{{\rm{\Lambda }}}_{\mathrm{SNe}}\right) & \propto & {\left({{\rm{\Lambda }}}_{0}{p}_{0}+{{\rm{\Lambda }}}_{\mathrm{SNe}}{p}_{\mathrm{SNe}}\right)}^{N}\\ & & \times \,\exp \left(-{{\rm{\Lambda }}}_{0}-{{\rm{\Lambda }}}_{\mathrm{SNe}}\right),\end{array}\end{eqnarray} \tag{ 12 }$$

where ${p}_{\mathrm{SNe}}$ (p₀) is the a priori weight that a transient is (is not) an SN Ia after the filtering process. With a suitable choice of prior, we can use Bayes's theorem to obtain the posterior. Considering the Jeffreys's prior,

$$\begin{eqnarray}&&p\left({{\rm{\Lambda }}}_{0},{{\rm{\Lambda }}}_{\mathrm{SNe}}\right)=\displaystyle \frac{1}{\sqrt{{{\rm{\Lambda }}}_{0}}}\displaystyle \frac{1}{\sqrt{{{\rm{\Lambda }}}_{\mathrm{SNe}}}},\end{eqnarray} \tag{ 13 }$$

the posterior takes the form

$$\begin{eqnarray}\begin{array}{rcl}p\left({{\rm{\Lambda }}}_{0},{{\rm{\Lambda }}}_{\mathrm{SNe}}\right|N) & \propto & p\left(N| {{\rm{\Lambda }}}_{0},{{\rm{\Lambda }}}_{\mathrm{SNe}}\right)p\left({{\rm{\Lambda }}}_{0},{{\rm{\Lambda }}}_{\mathrm{SNe}}\right)\\ & \propto & \displaystyle \frac{{\left({{\rm{\Lambda }}}_{0}{p}_{0}+{{\rm{\Lambda }}}_{\mathrm{SNe}}{p}_{\mathrm{SNe}}\right)}^{N}}{\sqrt{{{\rm{\Lambda }}}_{0}{{\rm{\Lambda }}}_{\mathrm{SNe}}}}\\ & & \times \,\exp \left(-{{\rm{\Lambda }}}_{0}-{{\rm{\Lambda }}}_{\mathrm{SNe}}\right).\end{array}\end{eqnarray} \tag{ 14 }$$

Integrating out the nuisance parameter, Λ₀, we have the marginalized posterior on ${{\rm{\Lambda }}}_{\mathrm{SNe}}=R\langle {VT}\rangle$ or, equivalently, on R:

$$\begin{eqnarray}\begin{array}{rcl}p\left(R| N\right) & = & {\int }_{0}^{\infty }p\left({{\rm{\Lambda }}}_{0},{{\rm{\Lambda }}}_{\mathrm{SNe}}\right|N)d{{\rm{\Lambda }}}_{0}\\ & \propto & \displaystyle \frac{{e}^{-R\langle {VT}\rangle }}{\sqrt{R\langle {VT}\rangle }}\times \left[{\left(R\langle {VT}\rangle {p}_{\mathrm{SNe}}\right)}^{N}\right.\\ & & +\,\left.\displaystyle \frac{N}{2}{p}_{0}{\left(R\langle {VT}\rangle {p}_{\mathrm{SNe}}\right)}^{N-1}\right],\end{array}\end{eqnarray} \tag{ 15 }$$

where we expand Equation (14) and integrate, keeping terms up to linear order in p₀ since we expect that ${p}_{0}\ll {p}_{\mathrm{SNe}}$.

SN Ia rates have been studied earlier in the literature (Gal-Yam et al. 2007; Dilday et al. 2008; Brown et al. 2019). Deep field instruments have provided estimates of the Ia rate out to high redshift (Gal-Yam et al. 2007). The iPTF, being an all-sky survey has a comparatively lower sensitivity to SNe Ia at ${z}_{\mathrm{median}}^{\mathrm{Ia}}=0.099$, evaluated in Section 5. The SDSS-II supernova survey has estimated the volumetric SN Ia rate at z ≈ 0.1 to be ${R}_{\mathrm{SNIa}}^{\mathrm{SDSS} \mbox{-} \mathrm{II}}\sim {2.9}_{-0.75}^{+1.07}\times {10}^{-5}\,{\mathrm{Mpc}}^{-3}\,{\mathrm{yr}}^{-1}$ (Dilday et al. 2008). Using our estimate of the spacetime sensitive volume from Equation (10), an estimate of the count of SNe Ia in iPTF is 630–1160. This is consistent with 1035 objects tagged "SN Ia" during the survey time.

While the number of transients tagged as "SN Ia" by human scanners during the iPTF survey time seem consistent with our ballpark above, the systematic uncertainty of such a classification remains unquantified. The quantities p₀, ${p}_{\mathrm{SNe}}$, and N in Equation (15) require a systematic search into the iPTF archival data to retrieve the candidate count and systematic errors associated with such a classification. We defer this and the computation of SN Ia volumetric rate to a future work in the series.

The methodology developed here facilitates the computation of spacetime volume sensitivities of general transient types. Of particular interest are the fast transients in iPTF archival data as discussed in Ho et al. (2018). Also, the observation of the "kilonova" resulting from the binary neutron star merger, GW170817 (Abbott et al. 2017a, 2017b, 2017c), hints toward the association of transients to binary neutron star mergers. There is no evidence of detection of such a transient in the iPTF data, in which case rate upper limits could be placed due to nondetection.