THE HIGH CADENCE TRANSIENT SURVEY (HITS). I. SURVEY DESIGN AND SUPERNOVA SHOCK BREAKOUT CONSTRAINTS

The core collapse of a massive star leads to the formation of a compact object and a shock wave that can unbind its outer layers in a core collapse SN. The shock can traverse the entire progenitor at very high velocities, v_shock ≈ 0.01 c, until it eventually emerges at the surface of the star in an SBO event. The shock's energy density will be dominated by radiation and its main source of opacity will be Compton scattering by free electrons (Weaver 1976). The shock thickness, ${d}_{\mathrm{shock}}$, can be estimated by equating the photon diffusion timescale, the typical time taken by photons to randomly walk out of the shock via electron scattering, and the advection timescale, the time taken by the shock to traverse one shock thickness. This results in a shock thickness of approximately ${d}_{\mathrm{shock}}\sim \tfrac{c\,{l}_{s}}{{v}_{\mathrm{shock}}}$, where l_s is the electron scattering mean-free path and v_shock is the shock speed, or an optical depth of ${\tau }_{\mathrm{shock}}\sim \tfrac{c}{{v}_{\mathrm{shock}}}$ (ignoring form factors that may be significant in some cases, see, e.g., Schawinski et al. 2008).

As the shock front approaches the stellar surface it will encounter a decreasing optical depth. When τ ≤ τ_shock the radiation dominated shock will leak its energy out of the stellar surface as a radiative precursor until the shock breaks out completely and becomes rarified. The SBO timescale, t_shock, will be given by the time that it takes for the radiative precursor to leak out into the stellar surface, ${t}_{\mathrm{shock}}\sim \tfrac{{d}_{\mathrm{shock}}}{{v}_{\mathrm{shock}}}$, but for an external observer it will also be affected by the light crossing time, ${t}_{\mathrm{light}}\sim \tfrac{{R}_{\star }}{c}$, where R_⋆ is the star's radius. Shock breakout timescales are expected to be typically of about an hour for RSGs, several minutes for blue supergiants (BSGs), and several seconds for stripped-envelope stars (Kistler et al. 2013, c.f. Schawinski et al. 2008).

The SBO properties will depend on many structure and composition parameters, but most strongly on the radius of the progenitor star (Calzavara & Matzner 2004). Stripped-envelope, BSG and RSG stars have very different radii: of the order of 5–10 R_⊙, 25–50 R_⊙, and 500–1000 R_⊙, respectively. Even for similar masses, the envelope structure can vary significantly depending on whether the envelope energy transport is radiative or convective. Because BSG stars have radiative envelopes and RSG stars, convective envelopes, they have very different effective polytropic indices (n = 3 in BSGs and n = 1.5 in RSGs), which also leads to different shock acceleration dynamics. Additionally, the presence of pre-SN winds can strongly change the shock propagation physics and the observable properties of the breakout event (Balberg & Loeb 2011; Moriya et al. 2011; Svirski et al. 2012; Svirski & Nakar 2014a).

The characteristic temperature of radiation during breakout depends mainly on the radius of the progenitor stars, being approximately proportional to the inverse of the progenitor radius squared (Calzavara & Matzner 2004), but the typical energy of the photons leaving the star will also depend on whether thermal equilibrium between the escaping radiation deep inside the star and the intervening gas is achieved and if so, the depth at which this occurs (Nakar & Sari 2010). Most photons leaving the star during SBO will be in X-ray and UV wavelengths, making the detection of these events very difficult from the ground. In fact, a few SBO candidates have been detected from space: three SNe II observed in the UV at early times with GALEX (Gezari et al. 2008; Schawinski et al. 2008) and a recent SN II observed in the optical using Kepler (Garnavich et al. 2016) have been associated with shocks breaking out of the envelopes of their RSG progenitors. Other SNe II appear to have shocks breaking into high density circumstellar material (CSM; Gezari et al. 2010, 2015; González-Gaitán et al. 2015; Garnavich et al. 2016; Khazov et al. 2016; Tanaka et al. 2016), which is supported by models of RSG winds (Mackey et al. 2014). A Swift X-ray transient has been associated with an stripped-envelope SNe Ibc, from what appears to be a Wolf Rayet star (Soderberg et al. 2008) and/or its surrounding wind (Svirski & Nakar 2014b). No direct detections of SBOs from BSG stars have been made, but there is indirect evidence for an SBO in SN 1987A (Lundqvist & Fransson 1996).

Although UV and X-ray detectors are better suited for SBO detections because of the typical temperatures encountered at the shock's front, Tominaga et al. (2011) suggested that optical detectors can also be used to find these events. In fact, their models suggest that it may be easier to detect SBOs in a systematic way with a new generation of large-field-of-view optical CCD cameras (see, e.g., Morokuma et al. 2014; Garnavich et al. 2016; Tominaga et al. 2016).

After SBO, the outer layers of the star will enter an adiabatic expansion and cooling phase, the so-called shock cooling phase. Nakar & Sari (2010) have divided the adiabatic expansion evolution after breakout into two distinct phases: planar and spherical expansion. During planar expansion the dominant contribution to the luminosity comes from the breakout shell, which would evolve adiabatically due to expansion at almost constant radius and with a luminosity ∝t^−4/3. During spherical expansion the radius cannot be considered constant and energy from the inner shells also contributes to the light curve, with a slower luminosity evolution ∝t^−0.17 to t^−0.35 for polytropic indices 1.5 ≤ n ≤ 3. For RSGs, the transition between planar and spherical expansion is expected at about 14 hr after explosion.

The evolution of the near-ultraviolet (NUV) and optical luminosity during adiabatic expansion can differ significantly from that of the total luminosity. While the total luminosity generally declines monotonically after SBO, the NUV and optical luminosity can decline and then rise again for several days up to maximum optical light (e.g., Schawinski et al. 2008; Tominaga et al. 2011).

The adiabatic approximation will be valid preferentially in the denser inner layers, as the outer layers become dominated by radiative diffusion with time. The radius where this transition occurs is known as the diffusion wave radius, which can be estimated by equating the star's radiation diffusion time with the time since explosion. In those stars, where radiative diffusion is negligible during this phase, with initial radii of R_⋆ ≤ 100 R_⊙, the product of ${R}_{\star }\,{E}_{\mathrm{in}}^{0.91}{M}_{\mathrm{ej}}^{-0.40}$ can be constrained by the luminosity evolution, where R_⋆ is the stellar radius, E_in is the explosion energy and M_ej is the mass of the ejecta (Chevalier 1992).

As radiative diffusion becomes important, the ionized He and H envelope can recombine in a wave that sweeps over the star inwards in mass coordinates, dominating the luminosity evolution in the plateau phase. The time when this transition occurs depends on the envelope mass and its structure, but it will usually start a few days after explosion (Chevalier 1992). The luminosity during this phase can evolve relatively slowly for a few months (see, e.g., Arcavi et al. 2012; Anderson et al. 2014a), depending on the explosion energy, ejected mass, initial radius, and composition (see, e.g., Popov 1993; Kasen & Woosley 2009), and will be followed by an exponentially decaying radioactive tail phase of evolution explained by the presence of newly synthesized ⁵⁶Ni. A schematic representation of the previous phases of evolution is shown in Figure 1.

Zoom In Zoom Out Reset image size

The most relevant variable characterizing the discovery potential of a survey telescope is the etendue—the product of the telescope's collecting area and the camera's field of view. Similarly, the data analysis challenges of a real-time survey are best characterized by the data acquisition rate, given by the ratio between the camera's image size and the typical exposure time plus overheads. A summary of some large-field-of-view cameras mounted in large collecting area telescopes is shown in Table 1. The large etendue of DECam is what enabled HiTS, but its large pixel number combined with the required real-time analysis became its main technical challenges.

As mentioned before, the most important science driver for HiTS was the detection of SBO events in the optical in a systematic fashion. In an initial stage, the strategy consisted of detecting SBO events in single band observations, giving preference to rapid monochromatic sampling over multiwavelength characterization, and attempting rapid follow-up of any of the detected events using our real-time analysis and Target of Opportunity follow-up capabilities. If the SBO detection is not achieved in real-time, the position of the associated rising SNe could be used to look for the SBO event in a post-processing, forced-photometry analysis (in which we fix both the position and shape of the point-spread function, PSF). Then, if the number of detected events was found to be significant, we would attempt the detection and spectral index characterization in a second phase of multiwavelength, high cadence, real-time analysis.

SBO characteristic timescales range from seconds to a few hours (Kistler et al. 2013), but we first focused on the longest timescales among SBO events, i.e., those originating from RSG stars and with characteristic timescales of an hour or more, in order not to compromise the survey effective volume. For this family of explosions, we use the RSG explosion models from Tominaga et al. (2011) and Nakar & Sari (2010), shown in Figures 2 and 3 for the same progenitor parameters, with zero-age main-sequence masses (M_ZAMS) of 13, 15, 20, and 25 M_⊙ and explosion energies of 1 foe⁹ for the 13, 15, and 20 M_⊙ models and 3 foe for the 25 M_⊙ model, in order to mimic the mass-energy relation found by Hamuy (2003) for observed SNe II and with some theoretical support from Müller et al. (2016) for core collapse SNe in general. To our knowledge, the models from Tominaga et al. (2011) are the only set of optical SBO models available in the literature that use exact solutions for realistic progenitor models and with a range of progenitor properties. We also use the analytic solutions from Nakar & Sari (2010) with the same parameters for comparison throughout the analysis.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We consider cadences between 1.3 and 2 hr (4, 5, or 6 epochs per field per night) in what follows, which can be considered to be on the long side of the SBO timescale distribution, but that would allow us to put constraints on the number of RSG SBO events that could be detected and characterized by future optical surveys with sub-hour cadences. For comparison, the Kepler ExtraGalactic survey (Olling et al. 2015) has a cadence of 30 minutes and the Dark Energy Survey time supernova survey (Dark Energy Survey Collaboration et al. 2016) has a typical cadence of about a week.

We focused on the detection of the initial peak seen in Figure 2, which we define as any detection with a signal-to-noise ratio (S/N) of at least 5 before 90% of the magnitude change between the optical peak maximum and subsequent minimum has occurred for the Tominaga et al. (2011) models (thick lines in Figure 2), or before the beginning of the spherical phase of evolution for the Nakar & Sari (2010) models (break in the dashed lines in Figure 3). To confirm that these restrictions would lead to recognizable SBOs we simulated 2400 light curves for each family of models with different extinctions, redshifts and progenitor properties using the empirical cadence and limiting magnitudes of the survey and found that we are typically 97% and 95% efficient in visually recognizing an SBO detection with the previous restrictions under the models of Tominaga et al. (2011) and Nakar & Sari (2010), respectively, with a purity of 99% for both families of models.

We considered whether to select the survey fields by their intrinsic properties or not, i.e., to perform a targeted or blind survey, respectively. Here, targeted refers to targeting clusters of galaxies instead of individual galaxies, given DECam's field of view. To choose the best strategy, we compared the outcomes of simulations of SBO events using the volumetric SN rate from Strolger et al. (2015) and the cluster SN rate from Graham et al. (2012) with typical galaxy cluster density profiles from Voigt & Fabian (2006). We found that targeting clusters of galaxies with DECam, we could obtain an estimated 10% increase in the number of events compared to an untargeted survey. Given that there is a relatively small number of known clusters at relatively low redshifts (Chon et al. 2013), of which only a small fraction remains visible for the entire night, we estimate that the time spent slewing the telescope from target to target would be more than 10% of the photon collecting time of a blind survey. Therefore, we decided to perform a blind search, targeting cluster, or supercluster fields only when they were adjacent to some of our blind fields, which would also allow us to obtain volumetric SN rates with fewer selection biases.

Having selected a blind field strategy, we find a reference right ascension (R.A.) that makes our targets achieve the lowest possible airmasses during all epochs, which is determined by the time of the year when the observations were initially allocated. Then we find a reference declination (decl.) that minimizes the average total extinction. We use an approximately rectangular grid of closely separated fields in the sky around this position, always switching to adjacent fields to minimize slew times¹⁰ , and balancing the exposure times with the rectangle dimensions to achieve an approximately constant airmass during each epoch of our observations. For example, we can observe along an approximately constant R.A. arc until the fields move by an hour angle equal to DECam's angular diameter, ${{\rm{\Theta }}}_{\mathrm{DECam}}$, which means switching to a new set of fields after approximately nine minutes.¹¹ Thus, the exposure time and number of approximately constant R.A. fields are related, as well as the exposure times plus overheads with the total number of fields and epochs, by the following approximate equations:

$$\begin{eqnarray}&&{N}_{\mathrm{fields}}=\displaystyle \frac{{T}_{\mathrm{night}}}{{N}_{\mathrm{epochs}}\,({T}_{\exp }+{T}_{\mathrm{overhead}})},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{N}_{\mathrm{decl}.}=\displaystyle \frac{{{\rm{\Theta }}}_{\mathrm{DECam}}}{15\,\deg /\mathrm{hr}\,({T}_{\exp }+{T}_{\mathrm{overhead}})},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{N}_{{\rm{R}}.{\rm{A}}.}=\displaystyle \frac{{N}_{\mathrm{fields}}}{{N}_{\mathrm{decl}.}},\end{eqnarray} \tag{ 3 }$$

where N_fields is the total number of fields to be visited, T_night is the duration of the night, N_epochs is the number of epochs per field, T_exp is the exposure time, T_overhead is the overhead time (maximum between slew and readout time), N_R.A. is the number of fields in the R.A. direction in our grid of fields to observe, and N_decl. is the number of fields in the decl. direction in our grid of fields to observe. The previous equations are only approximate and for a more realistic selection of the optimal combination of fields to observe, we simulated the exact evolution of airmass and extinction with time for different field configurations. Also note that because slew times are given by the largest R.A. or decl. difference in the Blanco equatorial mount, we try to slew the telescopes diagonally in a R.A.–decl. grid.

To estimate the decl. that would minimize the combined extinction from the Milky Way and the atmosphere, we used Milky Way extinction maps and the relation between atmospheric extinction and airmass. This requires having defined the exact number of fields to observe, which depend on the choices of the number of epochs per field per night, N_epochs, and the exposure time, T_exp (see Equation (1)), which are given by the solution of the survey design optimization problem.

In this analysis, we use two figures of merit: (1) the number of SNe that would be detected at least once within the initial optical peak caused by the shock breaking into the envelope of the RSG star and (2) the number of SNe that would be detected at least twice within the first day in rest-frame time after shock emergence. These two quantities are related, but their comparison will give a sense of the challenges associated with detecting the SBO optical peaks. Note that for a transient candidate to be considered stationary (as opposed to moving, e.g., asteroids) and to be detected in our online pipeline, we require a new event to be detected at least twice in image differences in the same position in the sky.

In order to simulate how many SNe would be detected with these restrictions, we use the previously described light curves computed by Tominaga et al. (2011) and Nakar & Sari (2010; see Figures 2 and 3), assuming a Λ–CDM cosmology with Ω_Λ = 0.73, Ω_M = 0.27, and H₀ = 71 km s⁻¹ Mpc⁻¹ and the DECam total efficiency. For the core collapse SN rate, we use the cosmic star formation rate density (SFRD) of Madau & Dickinson (2014), assuming a conversion efficiency between stellar mass and SNe, η, of 0.0091 ${M}_{\odot }^{-1}$ as measured by Strolger et al. (2015). Note that this conversion efficiency does not consider a large fraction of core collapse SNe missed by optical surveys (Mattila et al. 2012), so it is a conservative value for our purposes. We also tried the SFR of Horiuchi et al. (2011), a parametrization of Hopkins & Beacom (2006), which resulted in approximately 20% more events. Within a sufficiently large time interval T, we simulate a sampling function, e.g., 4, 5, or 6 epochs per night during several consecutive nights, and simulate a large set of explosions to estimate the likelihood of the event being detected at least once during the initial peak of a few hours seen in Figure 2 and the likelihood of the event being detected at least twice during the first rest-frame day after shock emergence. We record the distribution of detection ages to estimate more precisely the likelihood of detecting an event with a given age.

We weight these results taking into account the initial mass function (IMF) and normalizing their sum to 0.524 in order to reproduce the observed fraction of SNe II among core collapse SN explosions found by Li et al. (2011). In particular, we use a model weight proportional to a Salpeter-like IMF integrated in the vicinity of the model M_ZAMS:

$$\begin{eqnarray}&&w(M)\propto {\int }_{{a}_{M}}^{{b}_{M}}{m}^{-2.3}{dm},\end{eqnarray} \tag{ 4 }$$

where w(M) is the model weight associated with a zero-age main-sequence mass M and where a_M and b_M define the integration mass interval for a given model, chosen to be either at the low-mass limit of 8.5 M_⊙, equidistant between the masses of the available models, or at either 16.5 or 30 M_⊙ at the high-mass limit. The low-mass limit follows the low-mass constraint from Smartt et al. (2009) and the high-mass limit uses the high-mass constraint from the same work, but also considers the possibility that RSG stars are surrounded by CSM before explosion (see, e.g., Walmswell & Eldridge 2012; Gezari et al. 2015; González-Gaitán et al. 2015; Garnavich et al. 2016; Khazov et al. 2016; Tanaka et al. 2016), possibly becoming hidden from progenitor searches and allowing their masses to extend up to 30 M_⊙, which is consistent with the 95% confidence interval of Walmswell & Eldridge (2012). These Salpeter-like IMF distributions will be simply referred to as M16.5 and M30 distributions for the upper mass limits of 16.5 and 30 M_⊙, respectively.

Additionally, we have calculated model weights that would be representative of the estimated masses of observed RSG stars in the Milky Way (Levesque et al. 2005) and the Magellanic clouds (Levesque et al. 2006), assuming that these RSG stars are at the base of the giant branch and counting the number of stars closer to the mass of a given model to compute the weigths. The model properties as well as the integral intervals for Salpeter-like IMF distributions and associated weights for all SN II IMF distributions are shown in Table 2.

Table 2.  Explosion Model Properties from Tominaga et al. (2011) and Model Weights for the M16.5 and M30 Distributions Described in Section 3.4, as Well as those Representative of RSG Stars in the Milky Way (RSG_MW) and Magellanic Clouds (RSG_MC)

MZAMS	RpreSN	Z	Eexp	${a}_{M}^{{\rm{M}}16.5}$–${b}_{M}^{{\rm{M}}16.5}$	${w}^{{\rm{M}}16.5}$	${a}_{M}^{{\rm{M}}30}$–${b}_{M}^{{\rm{M}}30}$	wM30	${w}^{{\mathrm{RSG}}_{\mathrm{MW}}}$	${w}^{{\mathrm{RSG}}_{\mathrm{MC}}}$
13 M⊙	564 R⊙	0.02	1 foe	8.5–14.0	0.433	8.5–14.0	0.310	0.126	0
15 M⊙	507 R⊙	0.02	1 foe	14.0–16.5	0.091	14.0–17.5	0.086	0.154	0.091
20 M⊙	795 R⊙	0.02	1 foe	⋯	⋯	17.5–22.5	0.071	0.181	0.370
25 M⊙	1200 R⊙	0.02	3 foe	⋯	⋯	22.5–30.0	0.057	0.063	0.063

Note. The zero-age main-sequence (ZAMS), pre-sn radius, metallicity, and explosion energy are shown, as well as the integration intervals in M_⊙ and model weights from Equation (4). Note that the model weight's sum is 0.524, the observed fraction of SNe II in Li et al. (2011).

Download table as:  ASCIITypeset image

Assuming that all SNe II follow one of these explosion models, we estimate the number of events per unit redshift per field multiplying the cosmic SFRD, SFR(z); the conversion efficiency between stars and SNe, η; the total time of our simulation, T; the comoving volume per unit redshift bin per unit solid angle, $\tfrac{{dV}}{{dzd}{\rm{\Omega }}}(z);$ the field of view of DECam, ${\rm{\Delta }}{\rm{\Omega }};$ and the detection probabilities in the first day, ${P}_{M}^{{\rm{d}}}(z)$, or in the initial peak, ${P}_{M}^{{\rm{p}}}(z)$, obtained for one explosion model of mass M:

$$\begin{eqnarray}&&\displaystyle \frac{{{dn}}_{M}^{{\rm{d}}| {\rm{p}}}}{{dz}}=\eta \,\mathrm{SFR}(z)\,T\,\displaystyle \frac{{dV}}{{dzd}{\rm{\Omega }}}(z)\,{\rm{\Delta }}{\rm{\Omega }}\,{P}_{M}^{{\rm{d}}| {\rm{p}}}(z).\end{eqnarray} \tag{ 5 }$$

The detection probabilities are obtained simulating large sets of explosions in redshift bins of width 0.01 observed with the survey cadence and depth. We assume an exponential host galaxy extinction distribution A_V with a characteristic scale of λ_V = 0.187 and R_V = 4.5, as done in Strolger et al. (2015). We count the number of events satisfying the detection conditions and multiply this number by the efficiency with which we can visually recognize SBO light curves as described in Section 3.1.

For our survey design simulations, we assumed a median seeing at CTIO of 075 in r-band (Els et al. 2009), converted to full width at half maximum (FWHM) values for different frequencies using the relations found in the DECam exposure time calculator (ETC). Initially, we computed the survey depth at the zenith using the available DECam exposure time calculator (ETC) and then subtracted the excess extinction with respect to the zenith at a given airmass using the relations from Stone & Baldwin (1983). After our survey was completed, we found that this significantly underestimated the effect of airmass and modified the available ETCs to also take into account the effect of airmass on the FWHM, the atmospheric extinction and the sky emission. We also used a 063 DECam instrumental FWHM¹² and used a realistic evolution of the sky brightness¹³ to compute 50% completeness magnitudes (m₅₀).

The number of detected SNe per field up to a given redshift will then be the SN II IMF-weighted sum of the number of events per unit redshift per field for the different explosion models of mass M:

$$\begin{eqnarray}&&{N}^{{\rm{d}}| {\rm{p}}}(z)=\displaystyle \sum _{M}\,w(M){\int }_{0}^{z}\displaystyle \frac{{{dn}}_{M}^{{\rm{d}}| {\rm{p}}}}{{dz}}\,{dz}.\end{eqnarray} \tag{ 6 }$$

Different combinations of numbers of fields and numbers of epochs per night per field are associated with different exposure times, which in turn are associated with different limiting magnitudes and numbers of expected detections per night. In order to choose the best compromise between these quantities, keeping a cadence between 1.3 and 2 hr, we computed the total number of SNe that would be detected during SBO or that would be detected twice during first rest-frame day for different numbers of fields assuming 4, 5, and 6 epochs per night per field in different bands. We will show the results of these simulations using empirically derived limiting magnitudes in Section 5.

For our 2013A pilot campaign, we chose to observe in the u band based on the large temperatures during the SBO predicted by our models. However, in 2014A and 2015A, we switched to the g band after finding that this band offered the best compromise between increased detector efficiency and reduced SBO emission at larger wavelengths. Using these simulations, we chose 40 fields and 4 epochs per night per field for both our pilot 2013A u band campaign and real-time 2014A g band campaign. In 2015A, we changed to 50 fields and 5 epochs per night per field in order to explore faster cadences and to reduce the typical limiting magnitude of the survey, which facilitates follow-up observations with other telescopes. We will refer to the 2013, 2014, and 2015 HiTS survey campaigns as 13A, 14A, and 15A, respectively. In total, we were allocated four full nights in 13A, five full nights in 14A, and six full nights in 15A, always near the new moon. In 15A, we obtained three additional half nights at later times for multiwavelength follow up. The final coordinates of the fields observed in 13A, 14A, and 15A are shown in Figure 4 and their exact values are included in the Appendix Section in Tables 5–7. Note that the high cadence phase of observations were performed near the new moon: from March 12th to March 15th in 2013A (4 epochs per field per night), from March 1st to March 5th in 2014A (4 epochs per field per night), and from February 17th to February 22nd in 2015A (5 epochs per field per night), which explains a variation in optimal central R.A. by about two hours from 2013A to 2015A. Also note that, in order to minimize atmospheric extinction, a decl. close to the latitude of CTIO (−30°) is favoured, but in order to minimize galactic extinction large galactic latitudes are favoured.

Zoom In Zoom Out Reset image size

The data reduction pipeline consists of a series of sequential steps triggered in a distributed memory type cluster immediately after every observation, satisfying the following restrictions: (1) non-parallelizable steps cannot take longer than the exposure time plus overheads, e.g., data transfer, and (2) parallelizable steps cannot take longer than the exposure time plus overheads times the ratio between the total number of available cores and the product of the number of parallel tasks and the number of cores per task.

The general structure of the pipeline is the following. Data is first transferred to the distributed memory cluster to undergo an initial pre-processing phase with data stored in the different node's shared memory, including mosaic expansion, calibrations, cosmic-ray removal, and indexing of the processed data. Then, the image subtraction pipeline is triggered, again using the different node's shared memory, and including registration, convolution, difference, photometry, candidate classification, and indexing of the processed data. Finally, web visualization tools are updated with new candidates for human inspection.

The pipeline was written mostly in Python¹⁴ , using Sun Grid Engine (in 14A) and SLURM¹⁵ (in 15A) for job distribution within the NLHPC¹⁶ supercomputer. We use bash¹⁷ for many I/O operations, C for mosaic expansion (step 2 below), external software for pre-processing, cosmic-ray removal, and catalog generation (steps 3, 4, and 5 below), and Fortran 90¹⁸ subroutines parallelized via OpenMP¹⁹ and integrated with Python using F2PY²⁰ (Peterson 2009) for projection, convolution, difference object detection, and photometry (steps 7, 9, and 10 below). For the registration, convolution, and difference object detection, we use a similar approach to existing methods (e.g., Alard & Lupton 1999; Becker 2015), though with important differences in the representation of the convolution kernel and with our own Fortran routines for most of the computationally demanding steps. The following sequential steps are performed continuously during the night as data is acquired in the telescope.

• 1.  
  Data transfer: a bash script is left running on the telescope to upload raw DECam images from CTIO to the NLHPC (in 14A) or from La Serena to the NLHPC (in 15A) as soon as they arrive at the telescope's observer or La Serena pre-processing computer, respectively. The script looks for files that do not change in size after a lag of a few seconds and then transfers them to the NLHPC using rsync ²¹ in less than 10 s.
• 2.  
  Image expansion: the standard pre-processing DECam pipeline expands mosaic data into their individual CCDs serially, which can add a significant lag to our real-time pipeline. Thus, raw DECam mosaic files were expanded into their individual CCD images using a parallel custom made program, pimcopy, written in C and based on the CFITSIO library (Pence 1999). This allows us to perform the image expansion in less than three seconds for one mosaic image, about 60 times faster than a serial expansion.
• 3.  
  Pre-processing: with the expanded files, we run a pre-processing pipeline. In the 14A campaign, we used a custom pre-processing pipeline, which subtracted bias frames from the raw images, divided the resulting image by a flat-field frame, and computed inverse variance maps based on the Gaussian and Poisson noise components of the reduced images. The bias and flat-field frames were obtained for every CCD combining several bias and flat-field images observed during calibration time via median filtering and normalization. We also computed bad pixel masks for every CCD, applying thresholds over the individual flat-field images. However, we did not correct for CCD cross-talk effects, which we found later to be important because of the appearance of ghost images of bright stars. Therefore, in 15A, we used a modified version of the DECam community pipeline (DCP) for the pre-processing stage, including electronic bias calibration, cross-talk correction, saturation masking, bad pixel masking and interpolation, bias calibration, linearity correction, flat-field gain calibration, fringe pattern subtraction, bleed trail and edge bleed masking, and interpolation.
• 4.  
  Cosmic-ray removal: we remove cosmic rays from the pre-reduced images using a modified version of the public code CRBLASTER (Mighell 2010), using four cores per CCD via MPI.²² This code removes cosmic rays using Laplacian filter information, as in van Dokkum (2001). It takes about 20 s per CCD to remove cosmic rays, using 4 cores per CCD or 248 cores in total.
• 5.  
  Catalog generation and reference image: with the cosmic-ray-removed DECam images and previously obtained inverse variance maps, we use SExtractor (Bertin & Arnouts 1996) to derive source catalogs of the images in pixel coordinates. Given that we do not have previous reference images of the same area of the sky and that we are dealing with a variability of hours, we use the first image of the sequence that has a relatively low airmass and good atmospheric conditions as our reference image for the relative and absolute astrometric solutions, to define a projection grid, for image difference and for the relative and absolute magnitude calibrations of the image differences.
• 6.  
  Astrometric solution: with the previously generated catalogs, we match the brightest sources between the reference image catalog and the new image catalogs in pixel coordinates. Given that we chose the same pointing between different epochs and that the 4 m Blanco telescope has an equatorial mount, there is only a small offset with almost no rotation between epochs. Then we solve for a quadratic transformation between the pixel coordinates of both images using least square minimization, which was found to give consistent results to using the astrometric solution to apply the projection and convolution. To obtain the astrometric solution of the reference image, where all the projections are made, we solve for the astrometric solution using the positions of known stars from the USNO catalog (Monet et al. 2003) around our observations.²³
• 7.  
  Image projection: based on the previously derived transformations between the reference and new image in pixel coordinates we project the new image into the pixel grid of the reference image using a Lanzcos windowed sinc kernel with a kernel size of two. This kernel has been found to be a good compromise in terms of reduction of aliasing, sharpness, and minimal ringing (Turkowski & Gabriel 1990). This takes about 10 s using four cores per CCD.
• 8.  
  Kernel estimation: in order to match the image point-spread functions (PSFs) and scale of two images, we use pixel based kernels with a non-uniform pixel size distribution and an approximately circular shape. This is a different kernel model than that used in Becker (2015), which uses an analytic basis to describe the kernel accross the CCD. Our model is pixel based (see Masci et al. 2016), but with radially growing pixel sizes to reduce the number of free parameters per kernel from 625 to 81, which has a regularization effect over the values of the outermost pixels. An example of the previously determined convolution kernels for a given CCD is shown in Figure 5. DECam is composed of 62 CCDs, each one of 4k × 2k pixels. We divide every CCD in 18 (6 × 3) different regions with independently computed kernels. The size of the independent kernel regions was chosen by training kernels with only one bright star and then measuring the typical image subtraction residual change with distance for other test stars, such that the residual change at the kernel separation was comparable to the typical scatter of image subtraction residuals among different relatively bright stars. The kernels are derived via least square minimization of multiple pairs of stars taken from the same region of the CCD, selected for being isolated, bright, and within the linear response regime of the CCD. The kernels are defined so that the smaller full width at half maximum (FWHM) image is convolved into the larger FWHM image, in order to avoid effective deconvolutions.
• 9.  
  Image convolution and difference: we use the previously derived kernels to convolve the corresponding image and region of the CCD where it was defined. This changes both the shape and scale of the PSF of the convolved image. We then subtract the science and reference images, where one of them would have been convolved, to obtain difference and difference variance images in units of the image that has not been convolved. This takes less than 30 s using four cores per CCD.
• 10.  
  Image difference photometry: using the difference and difference variance images, we then perform an optimal photometry calculation (Naylor 1998) over the entire image assuming that the PSF is centered in every pixel. This is achieved in less than 10 s using four cores per CCD, producing a PSF integrated photometry image centered in every pixel and an associated PSF integrated photometry variance image, which are scaled up or down into the units of the reference image. We then select those pixels with S/Ns greater than five to obtain candidate image stamps, using the optimal photometry flux and variance, which are fed into a machine learning classifier. Computing S/Ns for measured fluxes centered in every pixel of the image allows us to perform normality tests over the S/N distribution for real-time quality control.
• 11.  
  Candidate filtering: every candidate is used to create candidate image stamps for which the dimensionality is reduced with custom-designed features that are input for a random forest classifier (e.g., Breiman 2001; Bailey et al. 2007; Goldstein et al. 2015) to be classified as either real or bogus. Because the training set was based on the 13A campaign (in the u band) to be used in the 14A/15A campaigns (in the g band), a key principle during the feature engineering process was to use as many dimensionless and scaleless quantities as possible. In order to not be dominated by unknown moving objects, e.g., asteroids, we consider that true non-moving transients are those that appear twice within a distance consistent with the astrometric errors and when at least one of the difference fluxes is positive with respect to the reference frame. A repetition of the candidate rules out moving transients, e.g., asteroids (to be presented in a separate publication), which are dominant among the candidates we detect, but only if it has at least one positive difference flux with respect to the reference frame. This is because a moving transient in the reference image should be present in all difference images with a negative difference flux with respect to the reference.
• 12.  
  Light-curve generation: once a candidate has passed all of the previous tests, a light curve is generated. By default, the light curve contains only those difference fluxes of candidates that were classified as true transients by the random forest classifier. Optionally, the light curve can contain the difference fluxes from all the available epochs, even when the candidate was not selected by its S/N or if it was not selected by the machine learning classifier. This can be computed revisiting the difference images and doing a forced photometry (fixing the PSF shape and position) using the optimal photometry method, which may take significantly longer times than all the previous steps because of the many more I/O operations required. For this reason, we revisit all the epochs only for a few visually selected candidates during the night, and for all the candidates classified as real non-moving transients during the day.
• 13.  
  Webpage visualization: as soon as a candidate is classified as a true candidate, it will be linked to a website that displays in visual form the location of all the candidates within the CCD. Candidates that are repeated and that have at least one positive difference flux with respect to the reference frame are marked with different colors to aid the visual inspection (using a combination of PHP²⁴ , JavaScript²⁵ and Highcharts²⁶ ). The positions of known stars, galaxies, moving objects, and variable stars are also shown with different colors. We mark those candidates that are found to be consistent with having periodic light curves, using a Lomb-Scargle periodogram analysis (Lomb 1976; Scargle 1982), to separate fast periodic stars from fast transients. All fields for a given campaign are displayed in a single map, with about three new candidates appearing at the location of the last field visited after every observation, giving higher visibility to those candidates that have not been inspected in order to aid with the process. The visual inspection of all fields and CCDs can be done in about one hour by a single person, but it can be easily accomplished in real-time with these visualization tools.
• 14.  
  Follow up: if a candidate is designated as a true non-moving transient by the previous tests, it will be followed-up only after going through the visual inspection test. We then decide whether to change DECam's observation plan, to obtain more observations in the same or other bands, or to trigger photometric or spectroscopic observations with other telescopes. Although we had the capability for triggering spectroscopic follow-up in a time comparable to about twice the cadence, no clear SBO candidates were found. Instead, during the last night of the 14A/15A high cadence run, we changed our observation to include multiple wavelengths and/or triggered spectroscopic observations toward some of the rising SNe.

Zoom In Zoom Out Reset image size

In total, we processed more than 10¹² pixels at a maximum rate of 4.5 Mpix s⁻¹, or 40 Mbps, generating about 10⁸ significant image subtraction candidates (with an S/N greater than five), of which about 10⁶ were classified as not bogus. Less than 10⁴ of these candidates were selected by the classifier twice in the same position in the sky, which were then visually inspected in a stream of about three new candidates per minute. The total lag time between shutter closing and end of candidate visualization for a given field was an approximately constant time of about five to six minutes taking into account all the previous steps. A visualization of the real-time pipeline steps is shown in Figure 6 for the 15A campaign, which had the shortest expsure times of the three available campaigns. In this figure, the yellow arrow labeled Image diff. pipeline contains steps 5–12. This figure shows why the visual inspection should be done as fast as the telescope slew plus exposure and that the DCP is normally run in parallel with the image difference pipeline from a previous exposure.

Zoom In Zoom Out Reset image size

In this section, we discuss the empirically determined depth of the survey, its relation to airmass and other observational variables, and how it compares to our survey design assumptions.

The depth of the survey in each of these campaigns was measured as an efficiency (probability of detection) versus magnitude for every observed epoch separately. This can be used to better estimate the expected number of events that should have been detected in the entire series of observations of a given field. The probability of detection at different magnitudes for a given set of representative observations is computed by building deep image stacks, co-adding selected epochs of a series of observations of the same field, and then measuring the fraction of real objects that were detected in single epoch observations relative to those in the deep stacks for a given magnitude range. Assuming that in the deep stacks all the true sources near the limiting magnitude of the individual images are detected, we measure the efficiency of detections in single epoch observations as a function of magnitude, taking into account all available CCDs simultaneously. We then fit the observed fraction with the following function.

$$\begin{eqnarray}&&P(m)=\displaystyle \frac{1}{2}\left[1+\mathrm{erf}\left(-\displaystyle \frac{m-{m}_{50}}{{\rm{\Delta }}{m}_{50}}\right)\right],\end{eqnarray} \tag{ 7 }$$

where P(m) is the probability of detection of a stellar-like source at a given magnitude, m; erf is the error function²⁷ ; m₅₀ is the best-fitting 50% completeness magnitude; and Δm₅₀ is a scale parameter approximately equal to half the width of the transition region to low efficiencies. We tried several analytic expressions, but we found Equation (7), using the error function, to better resemble the observed transition in efficiency. The best-fitting analytic approximation of the efficiencies for all the CCDs of one of the fields observed in 13A, 14A, and 15A is shown in Figure 7. The large spread in 50% completeness magnitudes is due to the airmass variations required in our observing strategy, shown in Figure 8, but also due to bad observing conditions in some epochs. Δm₅₀ best-fitting values were typically 0.6 ± 0.1 mag.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 9, we compare the best-fitting 50% completeness magnitudes (m₅₀) and the 5σ limiting magnitudes reported by the DECam ETCs versions 5 and 6 (ETC v5 and v6, respectively) modified to include the effect of airmass. Note that the ETC v5 is included because it was the latest available ETC at the time when the 15A semester observations were performed. We show the 5σ limiting magnitudes using the modified ETC v5 with the predicted FWHM and sky brightness (ETC v5 + pred. cond.), the modified ETC v5 with the observed FWHM and sky brightness (ETC v5 + obs. cond.), the modified ETC v6 with the predicted FWHM and sky brightness (ETC v6 + pred. cond.), the modified ETC v6 with the observed FWHM and sky brightness (ETC v6 + obs. cond.) as well as the empirical 50% completeness magnitude (Observed). Note that in the 15A campaign the fourth and fifth nights were strongly affected by cloud cover and resulted in only three epochs in total (observation numbers 14–16) and that the 15A campaign contained follow-up observations starting at observation number 22. In all plots, the measured 50% completeness magnitude matches better the 5σ limiting magnitude produced with the modified ETC v6 using the empirical FWHM. We found that there is a difference of approximately 0.6 mag between the ETC v5 and v6, but also that there were FWHM variations beyond what is expected from airmass variations, which are also required to explain some of the m₅₀ variations.

Zoom In Zoom Out Reset image size

As mentioned before, we initially used the relation between airmass and extinction from Stone & Baldwin (1983) to calculate an excess extinction with respect to the limiting magnitude at the zenith, which is given by the original ETC, and use this information to derive the limiting magnitude at a given airmass. However, this led to a significant overerestimation of the limiting magnitudes. Instead, in the modified version of the ETCs, we assume that the FWHM scales as $\cos {(z)}^{-3/5}$, where z is the zenithal angle, that the sky emission scales linearly with airmass and that the atmospheric extinction increases exponentially with airmass to derive the limiting magnitude at a given airmass. A Python version of our modified ETC can be found at https://github.com/fforster/HiTS-public.

The relation between m₅₀, FWHM, and airmass is shown in Figure 10. The expected relation between m₅₀ and FWHM for constant seeings of 075 (dotted–dashed gray line) and 10 (continuous gray line) in the r band are shown, scaled to the u or g bands and using a range of airmasses (between 1.0 and 1.6 for 2013A and between 1.1 and 1.9 for 14A and 15A) to derive the FWHMs, where we also assume an additional 063 DECam instrumental seeing. We use a fixed sky brightness of 22.6 mag arcsec⁻² in the u band for 13A and 22 mag arcsec⁻² in the g band for 14A and 15A, scaled linearly with airmass (in physical units) with the modified ETC v6. With the modified ETCs, we expect that for the typical range of airmasses the limiting magnitudes can vary as much as 1.1 mag. Thus, the upper-left areas of Figure 10 can be explained by an airmass effect, but the bottom-right areas are mostly due to changing observing conditions, including the appearance of clouds in some of the observing nights.

Zoom In Zoom Out Reset image size

The previous efficiencies are based on single epoch measurements on full DECam mosaics, but our detections are based on image differences. To account for the additional loss of efficiency due to image differencing, we injected artificial stars and studied their S/Ns at injection and after recovery, concluding that the effect was typically a decrease in the S/N of up to $\sqrt{2}$, which would be expected when subtracting two images with the same noise level, or a shift in the efficiency versus magnitude relation of up to $2.5\,\,{\mathrm{log}}_{10}\,(\sqrt{2})\approx 0.38$ mag. This is because we did not use deep templates for image subtraction, but used the first observation at a given band that approached the lowest airmass values achievable (about 1.05 in 13A and 1.2 in 14A/15A, see Figure 8) of the first night, which had photometric conditions (second observation in the main band for the three campaigns), effectively amplifying the noise by up to $\sqrt{2}$ since the images are dominated by sky shot noise. Therefore, we shifted our efficiency-magnitude relations by 0.38 mag to account for the effect of image differences, but also used the original unshifted relation as the best-case scenario in what follows.

Overall, there are four relevant factors that made the survey depth shallower than expected: (1) the ETC available when designing the survey (v5) overestimated the limiting magnitudes by approximately 0.6 mag with respect to the latest ETC version (v6), which matches our empirical 50% completeness magnitudes well; (2) the effect of airmass was larger than expected, with a variation of up to 1.1 mag in the 50% completeness magnitudes instead of less than 0.3 mag using the relations from Stone & Baldwin (1983); (3) the observational conditions were not always good, with nights that had significantly larger than expected seeing values or even clouds; and (4) the difference imaging could lead to a loss in depth of up to 0.38 mag with respect to the individual epoch measured depths. The combination of all the previous effects can be seen in the large magnitude differences seen in Figure 9.

A practical description of the main observational parameters for the 13A, 14A, and 15A observational campaigns is shown in Table 3. Note that, for the 15A campaign, we had three DECam half nights separated by 2, 5, and 20 days from the end of the high cadence phase, where we also observed selected fields in the r and i bands.

The lower than expected survey depths obtained in Section 5.1 (at least 0.6 mag of difference) should have an important effect on the predicted number of SBO events, which means that our initial assumptions during survey design had to be corrected. Thus, we have recomputed the expected number of events for different observational strategies using different assumptions about the survey depth. We use both the modified ETCs v5 and v6 for a different number of fields per night, assuming no airmass evolution or a realistic airmass evolution given a typical decl. and 4 or 5 epochs per field per night for 13A and 14A or 15A, respectively. We also test for the effect of the image difference process by assuming that it results in a loss of depth of 0.38 mag or assuming that it results in no loss in depth.

We found that using the ETC v5 instead of the ETC v6 led to an overestimation of the number of SNe that would be detected during SBO peak or twice during the first rest-frame day after shock emergence of approximately 2.5 in both cases. Using the relation between airmass and extinction from Stone & Baldwin (1983) to derive an excess extinction at a given airmass led to an overestimation of the number of SNe detected at SBO peak or twice within the same rest-frame day after shock emergence of approximately 1.8 and 1.4, respectively. Moreover, the loss in depth due to the difference imaging process can lead to an overestimation of up to approximately 1.8 in both quantities.

Using the ETC v6 and a realistic airmass evolution, we show the predicted number of SNe that would be detected at SBO optical peak or at least twice within the first rest-frame day after shock emergence as a function of the number of fields per night per epoch in Figure 11. We considered 4 continuous nights with 4 epochs per night in the u band, 5 continuous nights with 4 epochs per field per night in the g band, and 6 continuous nights with 5 epochs per field per night in the g band, as we did in 13A, 14A, and 15A, respectively. We also consider the two M_ZAMS distributions discussed in Section 3.4, M16.5 and M30, which differ in their upper mass cutoff limit.

Zoom In Zoom Out Reset image size

Figure 11 shows the predicted numbers of detections using the modified DECam ETC v6, for the M16.5 and M30 averaged ensemble of models as a function of the number of fields observed per night for the 13A, 14A, and 15A survey campaigns, as described in Section 3.4. We assume 4 epochs per field per night in the u band during 4 continuous nights for 13A, 4 epochs per field per night in the g band during 5 continuous nights for 14A or 5 epochs per field per night in the g band during 6 continuous nights for 15A. Shorter exposure times have smaller m₅₀, which favors more nearby SNe that are easier to follow up, but which requires a faster real-time pipeline. We found no strong deviations as a function of the number of fields within the parameter range under consideration, i.e., wide/shallow surveys result in a similar number of events than narrow/deep surveys, with less than 20% and 30% variations with respect to the mean for the M16.5 and M30 distributions, respectively. The actual observational campaigns were optimized using the ETC v5, which resulted in maxima at a larger number of fields per epoch than the ETC v6. Wider/shallower strategies enable more follow-up capabilities, favouring larger numbers of fields per epoch or shorter exposure times depending on the available resources.

In Figure 12, we compare the predicted number of SNe computed using the empirical observational conditions (airmass, FWHM, sky brightness, and transparency) versus using the expected observational conditions (FWHM, sky brightness, and transparency), in order to estimate the effect of cloud cover and bad seeing in the 13A, 14A, and 15A campaigns. From this, we conclude that in the 13A campaign the atmospheric conditions led to an overestimation factor of approximately 1.1 for the number of SNe that would have been discovered at the SBO optical peak and 1.2 for those discovered twice within the first day after shock emergence. In the 14A campaign, the overestimation factors were approximately 1.4 and 1.6, respectively, and in the 15A campaign, approximately 1.7 and 1.6, respectively.

Zoom In Zoom Out Reset image size

In total, we estimate that the expected limiting magnitudes (using the unmodified ETC v5, the relation between extinction and airmass from Stone & Baldwin 1983, expected atmospheric conditions, and negligible reference image noise) led to an SBO number overestimation factor of up to 14x compared to using the empirical limiting magnitudes.

The M_ZAMS distribution of the SNe found at SBO optical peak or twice within the first rest-frame day after shock emergence should differ considerably from the original SN II progenitor IMF distribution. More massive progenitor scenarios that become observed SNe II appear to have larger radii and more energetic explosions (Hamuy 2003; Müller et al. 2016), which would make them brighter after shock emergence. This should bias our discoveries toward larger M_ZAMS values in the M30 distribution. This can be seen in Figure 13, where we combine the 13A, 14A, and 15A campaigns and separate the predicted cumulative number of detections by associated progenitor model M_ZAMS for the M30 distribution (we lack enough model resolution for the M16.5 distribution) and the Tominaga et al. (2011) models. In the M30 distribution, about 70%–80% of the expected SBO optical peak detections should come from the weighted models with M_ZAMS of 20 or 25 M_⊙, and about 60%–70% in the case of the double detections at day 1 after shock emergence for each survey campaign. This means that the search for SBO or very young SNe will be strongly biased toward the more massive SN II explosions.

Zoom In Zoom Out Reset image size

The number of predicted events is sensitive to the upper mass limit of the Type II supernova (SN II) IMF distribution, with a larger upper mass limit leading to more detections because they are generally easier to detect than lower mass, generally smaller radii explosions. As shown in Figure 13, an upper mass limit of 30 M_⊙ would lead to about four times more SBO detections than an upper mass limit of 16.5 M_⊙, and about twice the number of SNe younger than one rest-frame day after shock emergence with at least two detections. We also tried a steeper IMF, with a γ of 2.3 instead of 1.3, and we expect approximately 50% fewer SBO detections for the M30 distribution.

We have also explored the sensitivity to explosion energy. More energetic models have brighter, narrower SBO optical peaks, and have significantly brighter shock cooling optical light curves. The 25 M_⊙ M_ZAMS model from Tominaga et al. (2011), which is the dominant model in the M30 distribution, would produce 30% fewer or more optical SBO peak detections if instead of the 3 foe model we had used the 1 and 9 foe models, respectively. However, it would produce 2.7 and 2.6 times fewer or more SNe younger than one day after shock emergence than the 3 foe model, respectively. Thus, the explosion energy is not a very important parameter for the predicted number of SBO detections, but it is a relevant parameter for the number of very young SNe detected in their shock cooling phase.

We have also explored analytic models for the SBO (Nakar & Sari 2010) and early days of evolution of SNe II (Nakar & Sari 2010; Rabinak & Waxman 2011) and have found that they are effectively a two parameter family in the optical, where the two parameters are the progenitor radius and the ratio between explosion energy and progenitor mass. Larger radii, but also larger energy to mass ratio explosions produce brighter/slower light curves in these models. We have found that the models from Nakar & Sari (2010) are generally in better agreement with those from Tominaga et al. (2011) during the rise to maximum light, but there are significant differences among all these models, which suggest that the physical interpretation of early SN light curves should be treated with caution.

In addition to the model uncertainties, it is not clear how the SBO optical peak should be observationally defined. We define the end of the SBO optical peak for the Tominaga et al. (2011) models as the time when 90% of the magnitude change between the SBO peak and subsequent minimum (maximum in magnitudes) has occured, which is shown with thick lines in Figure 2. When trying a more conservative definition, as the time when only 50% of the previous magnitude change has occurred, we found a decrease of only 1% in the number of predicted SBOs. If we define the SBO optical peak end as the time when only 10% of the previous magnitude change has occured, we see a decrease in the number of predicted SBOs of only 40%. This low sensitivity to the definition of the end of SBO optical peak means that most predicted detections are expected to occur only once and very close to the SBO optical peak maximum (minimum in magnitudes). In fact, the number of SBO optical peak double detections is expected to be very small, only 18% of the SBO detections in the 90% magnitude difference definition, 3% of the SBO detections in the 50% magnitude change definition, and none in our simulations for the 10% magnitude difference definition.

In Figure 14, we plot the cumulative distribution of expected SBO detections as a function of redshift using the Tominaga et al. (2011) models with the M16.5 and M30 distributions for the combined 13A, 14A, and 15A campaigns. A dotted horizontal gray line indicates the number of predicted events from which a non-detection in the three survey campaigns has a probability of less than 5%. We conclude that, although the search for SBO optical peaks should be sensitive to the highest mass limit of the SN II IMF distribution, the shallower than expected depth of the survey implies that we cannot rule out either the M16.5 or M30 distributions for this family of models.

Zoom In Zoom Out Reset image size

We then compare the previous results to those obtained using the brighter and longer-lived Nakar & Sari (2010) models and considering different model weigths, which resembles observed RSG stars discussed in Section 3.4. This is seen in Figure 15, where we show the total number of events that should have been detected with different assumptions about the SBO models, quality of the image difference and IMF weighting scheme. The corresponding non-detection probabilities are shown in Table 4, where we have found that the probability of having k detections in a Poisson distribution with mean λ is ${\lambda }^{k}{e}^{-\lambda }/k!$, which for k = 0 is e^−λ. The brighter and longer-lived models from Nakar & Sari (2010) models can be marginally rejected under all reasonable IMF distributions, especially the IMF distributions that favor larger progenitor masses. We cannot reject the Tominaga et al. (2011) models under any reasonable IMF distribution, but we marginally reject that most SBOs resemble those produced by their more massive 20 and 25 M_⊙ at ZAMS models shown in Figures 2 and 3.

Zoom In Zoom Out Reset image size

Although no good SBO candidates were detected, several very young SN candidates were detected in the three HiTS campaigns, a few of them apparently less than a day after shock emergence, and more than a hundred SN candidates of all ages in total, whose positions were made available to the community (Förster et al. 2014a, 2014b, 2015a, 2015b, 2015c, 2015d, 2015e, 2015f, 2015g, 2015h, 2015i, 2015j, 2015k, 2015l, 2015m). We could spectroscopically classify about a dozen of these (see, e.g., Anderson et al. 2014b, 2015a, 2015b, 2015c; Walton et al. 2014a, 2014b; Förster et al. 2015n), confirming their SN nature in all cases. We have also used the models from Tominaga et al. (2011) to study the number of SNe II expected to be detected twice during their rise (because our template images were taken at the beginning of each run, we should only detect rising SNe) and found that we expect about 80 and 130 SNe II for the M16.5 and M30 distributions, respectively. Because these numbers are derived using an SN efficiency computed from actual observations (Strolger et al. 2015), this suggests that the total number of events we detect is roughly consistent with the observed rate of SNe within the relevant redshift range. Considering a contribution from other SN types (which will be the subject of future publications) and model uncertainties, this also suggests that the presence of other types of transients is not required to explain the number of detected events in our sample, and may even suggest a preference for the M16.5 distribution. We will present a catalog of calibrated light curves that will include follow-up resources other than DECam in a future publication, including host galaxy redshifts and post-SN images for the most interesting candidates. Interestingly, the SN candidates apparently younger than one day after shock emergence appear to show initially faster than expected rise times (see González-Gaitán et al. 2015; Garnavich et al. 2016; Tanaka et al. 2016), which suggests that the evolution of SN II light curves during the first days may not be explained by standard assumptions about their CSM.