Evidence for Late-stage Eruptive Mass Loss in the Progenitor to SN2018gep, a Broad-lined Ic Supernova: Pre-explosion Emission and a Rapidly Rising Luminous Transient

ZTF observing time is divided between several different surveys, conducted using a custom mosaic camera (Dekany et al. 2016) on the 48 inch Samuel Oschin Telescope (P48) at Palomar Observatory. See Bellm et al. (2019a) for an overview of the observing system, Bellm et al. (2019b) for a description of the surveys and scheduler, and Masci et al. (2019) for details of the image processing system.

Every 5σ point-source detection is saved as an "alert." Alerts are distributed in avro format (Patterson et al. 2019) and can be filtered based on a machine learning–based real-bogus metric (Duev et al. 2019; Mahabal et al. 2019), light-curve properties, and host characteristics (including a star-galaxy classifier; Tachibana & Miller 2018). The ZTF collaboration uses a web-based system called the GROWTH marshal (Kasliwal et al. 2019) to identify and keep track of transients of interest.

ZTF18abukavn was discovered in an image obtained at UT 2018 September 9 03:55:18 (start of exposure) as part of the ZTF extragalactic high-cadence partnership survey, which covers 1725 deg² in six visits (3g, 3r) per night (Bellm et al. 2019b). The discovery magnitude was $r=20.5\pm 0.3\,\mathrm{mag}$, and the source position was measured to be $\alpha ={16}^{{\rm{h}}}{43}^{{\rm{m}}}{48.22}^{{\rm{s}}}$, $\delta =+{41}^{{\rm{d}}}{02}^{{\rm{m}}}43\buildrel{\rm{s}}\over{.} 4$ (J2000), coincident with a compact galaxy (Figure 1) at $z={\rm{0.03154}}$ or $d\approx {\rm{143}}\,\mathrm{Mpc}$. As described in Section 2.3, the redshift was unknown at the time of discovery; we measured it from narrow galaxy emission lines in our follow-up spectra. The host redshift along with key observational properties of the transient are listed in Table 1.

Zoom In Zoom Out Reset image size

Table 1.  Key Observational Properties of SN2018gep and Its Host Galaxy

Parameter	Value	Notes
z	0.03154	From narrow host emission lines
${L}_{\mathrm{peak}}$	$\gtrsim 3\times {10}^{43}\,\mathrm{erg}$	Peak UVOIR bolometric luminosity
${t}_{\mathrm{rise}}$	0.5–3 days	Time from t0 to ${L}_{\mathrm{peak}}$
Erad	${10}^{50}\,\mathrm{erg}$	UVOIR output, ${\rm{\Delta }}t=0.5$–40 days
${M}_{r,\mathrm{prog}}$	−15	Peak luminosity of pre-explosion emission
${E}_{\gamma ,\mathrm{iso}}$	$\lt 4.9\times {10}^{48}$	Limit on prompt gamma-ray emission from Fermi/GBM
LX	$\lt 2.5\times {10}^{41}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$	X-ray upper limit from Swift/XRT at Δt = 0.4–14 days
	$\lt {10}^{40}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$	X-ray upper limit from Chandra at Δt = 15 and ${\rm{\Delta }}t=70\,\mathrm{days}$
$\nu {L}_{\nu }$	$\approx {10}^{37}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$	9 GHz radio luminosity from VLA at Δt = 5 and Δt = 16
${M}_{* ,\mathrm{host}}$	$1.3\times {10}^{8}\,{M}_{\odot }$	Host stellar mass
SFRhost	0.12 ${M}_{\odot }$ ${\mathrm{yr}}^{-1}$	Host star formation rate
Host metallicity	1/5 solar	Oxygen abundance on O3N2 scale

Download table as:  ASCIITypeset image

As shown in Figure 2, the source brightened by over two magnitudes within the first three hours. These early detections passed a filter written in the GROWTH marshal that was designed to find young SNe. We announced the discovery and fast rise via the Astronomer's Telegram (Ho et al. 2018c), and reported the object to the IAU Transient Server (TNS²⁸ ), where it received the designation SN2018gep.

Zoom In Zoom Out Reset image size

We triggered ultraviolet (UV) and optical observations with the UV/Optical Telescope (UVOT; Roming et al. 2005) aboard the Neil Gehrels Swift Observatory (Gehrels et al. 2004), and observations began 10.2 hr after the ZTF discovery (Schulze et al. 2018a). A search of IceCube data found no temporally coincident high-energy neutrinos (Blaufuss 2018).

Over the first two days, the source brightened by two additional magnitudes. A linear fit to the early g-band photometry gives a rise of 1.4 ± 0.1 mag ${\mathrm{hr}}^{-1}$. This rise rate is second only to the IIb SN 16 gkg (Bersten et al. 2018) but several orders of magnitude more luminous at discovery (${M}_{g,\mathrm{disc}}\approx -17\,\mathrm{mag}$).

To establish a reference epoch, we fit a second-order polynomial to the first three days of the g-band light curve in flux space, and define t₀ as the time at which the flux is zero. This gives t₀ as being 25 ± 2 minutes prior to the first detection, or ${t}_{0}\,\approx$ UT 2018 September 9 03:30. The physical interpretation of t₀ is not straightforward, since the light curve flattens out at early times (see Figures 2 and 3). We proceed using t₀ as a reference epoch but caution against assigning it physical meaning.

Zoom In Zoom Out Reset image size

From ${\rm{\Delta }}t\approx 1\,\mathrm{day}$ to ${\rm{\Delta }}t\approx 60\,\mathrm{days}$, we conducted a photometric follow-up campaign at UV and optical wavelengths using Swift/UVOT, the Spectral Energy Distribution Machine (SEDM; Blagorodnova et al. 2018) mounted on the automated 60 inch telescope at Palomar (P60; Cenko et al. 2006), the optical imager (IO:O) on the Liverpool Telescope (LT; Steele et al. 2004), and the Lulin 1 m Telescope (LOT).

Basic reductions for the LT IO:O imaging were performed by the LT pipeline.²⁹ Digital image subtraction and photometry for the SEDM, LT, and LOT imaging was performed using the Fremling Automated Pipeline (FPipe; Fremling et al. 2016). Fpipe performs calibration and host subtraction against Sloan Digital Sky Survey reference images and catalogs (SDSS; Ahn et al. 2014). SEDM spectra were reduced using pysedm (Rigault et al. 2019).

The UVOT data were retrieved from the NASA Swift Data Archive³⁰ and reduced using standard software distributed with HEAsoft version 6.19.³¹ Photometry was measured using uvotmaghist with a 3'' circular aperture. To remove the host contribution, we obtained a final epoch in all broadband filters on UT 2018 October 18 and built a host template using uvotimsum and uvotsource with the same aperture used for the transient.

Figure 3 shows the full set of light curves, with a cross denoting the peak of the r-band light curve for reference. The position of the cross is simply the time and magnitude of our brightest r-band measurement, which is a good estimate given our cadence. The photometry is listed in Table 5 in Appendix A. Note that despite the steep spectral energy distribution (SED) at early times, the K-correction is minimal. We estimate that the effect is roughly 0.03 mag, which is well within our uncertainties. In Figure 4 we compare the rise time and peak absolute magnitude to other rapidly evolving transients from the literature.

Zoom In Zoom Out Reset image size

The first spectrum was taken 0.7 day after discovery by the Spectrograph for the Rapid Acquisition of Transients (SPRAT; Piascik et al. 2014) on the LT. The spectrum showed a blue continuum with narrow galaxy emission lines, establishing this as a luminous transient (${M}_{{\rm{g}},\mathrm{peak}}=-19.7$). Twenty-three optical spectra were obtained from ${\rm{\Delta }}t=0.7$–61.1 days using SPRAT, the Andalusia Faint Object Spectrograph and Camera (ALFOSC) on the Nordic Optical Telescope (NOT), the Double Spectrograph (DBSP; Oke & Gunn 1982) on the 200 inch Hale telescope at Palomar Observatory, the Low Resolution Imaging Spectrometer (LRIS; Oke et al. 1995) on the Keck I 10 m telescope, and the Xinglong 2.16 m telescope (XLT+BFOSC) of NAOC, China (Wang et al. 2018). As discussed in Section 3.2, the early ${\rm{\Delta }}t\lt 5\,\mathrm{days}$ spectra show broad absorption features that evolve redward with time, which we attribute to carbon and oxygen. By ${\rm{\Delta }}t\sim 8\,\mathrm{days}$, the spectrum resembles a stripped-envelope SN, and the usual broad features of a Ic-BL emerge (Costantin et al. 2018).

We use the automated LT pipeline reduction and extraction for the LT spectra. LRIS spectra were reduced and extracted using Lpipe (Perley 2019). The NOT spectrum was obtained at parallactic angle using a 1'' slit, and was reduced in a standard way, including wavelength calibration against an arc lamp, and flux calibration using a spectrophotometric standard star. The XLT+BFOSC spectra were reduced using the standard IRAF routines, including corrections for bias, flat field, and removal of cosmic rays. The Fe/Ar and Fe/Ne arc lamp spectra obtained during the observation night are used to calibrate the wavelength of the spectra, and the standard stars observed on the same night at similar airmasses as the supernova were used to calibrate the flux of the spectra. The spectra were further corrected for continuum atmospheric extinction during flux calibration, using mean extinction curves obtained at Xinglong Observatory. Furthermore, telluric lines were removed from the data.

Swift obtained three UV-grism spectra between UT 2018 September 15 3:29 and 6:58 (${\rm{\Delta }}t\approx 6.4\,\mathrm{days}$) for a total exposure time of 3918 s. The data were processed using the calibration and software described by Kuin et al. (2015). During the observation, the source spectrum was centered on the detector, which is the default location for Swift/UVOT observations. Because of this, there is second-order contamination from a nearby star, which was reduced by using a narrow extraction width (13 instead of 25). The contamination renders the spectrum unreliable at wavelengths longer than 4100 Å, but is negligible in the range 2850–4100 Å due to absorption from the interstellar medium. Below 2200 Å, the spectrum overlaps with the spectrum from another star in the field of view.

The resulting spectrum (Figure 5) shows a single broad feature between 2200 Å and 3000 Å (rest frame). One possibility is that this is a blend of the UV features seen in superluminous supernovae (SLSNe). Line identifications for these features vary in the SLSN literature, but are typically blends of Ti iii, Si iii, C ii, C iii, and Mg ii (Quimby et al. 2011; Howell et al. 2013; Mazzali et al. 2016; Yan et al. 2017).

Zoom In Zoom Out Reset image size

The spectral log and a figure showing all the spectra are presented in Appendix B. In Section 3.2 we compare the early spectra to spectra at similar epochs in the literature. We model one of the early spectra, which shows a "W" feature that has been seen in SLSNe, to measure the density, density profile, and element composition of the ejecta. From the Ic-BL spectra, we measure the velocity evolution of the photosphere.

The nominal ZTF pipeline only generates detections above a 5σ threshold. To extend the light curve further back in time, we performed forced photometry at the position of SN2018gep on single-epoch difference images from the IPAC ZTF difference imaging pipeline. The ZTF forced photometry point spread functions (PSF)–fitting code will be described in Y. Yao et al. (2019, in preparation). As shown in Figure 2, forced photometry uncovered an earlier 3σ r-band detection.

Next, we searched for even fainter detections by constructing deeper reference images than those used by the nominal pipeline, and subtracting them from 1 to 3 day stacks of ZTF science images. The reference images were generated by performing an inverse-variance weighted coaddition of 298 R-band and 69 g-band images from PTF/iPTF taken between 2009 and 2016 using the CLIPPED combine strategy in SWarp (Bertin 2010; Gruen et al. 2014). PTF/iPTF images were used instead of ZTF images to build references as they were obtained years prior to the transient, and thus less likely to contain any transient flux. No cross-instrument corrections were applied to the references prior to subtraction. Pronounced regions of negative flux on the PTF/iPTF references caused by crosstalk from bright stars were masked out manually.

We stacked ZTF science images obtained between UT 2018 February 22 and 2018 August 31 in a rolling window (segregated by filter) with a width of 3 days and a period of 1 day, also using the CLIPPED technique in SWarp. Images taken between 2018 Sep 01 and t₀ were stacked in a window with a width of 1 day and a period of 1 day. Subtractions were obtained using the HOTPANTS (Becker 2015) implementation of the Alard & Lupton (1998) PSF matching algorithm. Many of the ZTF science images during this period were obtained under exceptional conditions, and the seeing on the ZTF science coadds was often significantly better than the seeing on the PTF/iPTF references. To correct for this effect, ZTF science coadds were convolved with their own PSFs, extracted using PSFEx, prior to subtraction. During subtraction, PSF matching and convolution were performed on the template and the resulting subtractions were normalized to the photometric system of the science images. We show two example subtractions in Figure 6.

Zoom In Zoom Out Reset image size

Using these newly constructed deep subtractions, PSF photometry was performed at the location of SN2018gep using the PSF of the science images. To estimate the uncertainty on the flux measurements made on these subtractions, we employed a Monte Carlo technique, in which thousands of PSF fluxes were measured at random locations on the image, and the PSF-flux uncertainty was taken to be the 1σ dispersion in these measurements. We loaded this photometry into a local instance of SkyPortal (Van der Walt et al. 2019), an open-source web application that interactively displays astronomical data sets for annotation, analysis, and discovery.

We detected significant flux excesses at the location of SN2018gep in both g and r bands in the weeks preceding t₀ (i.e., its first detection in single-epoch ZTF subtractions). The effective dates of these extended prediscovery detections are determined by taking an inverse-flux variance weighted average of the input image dates. The detections in the week leading up to explosion are ${m}_{g}\sim {m}_{r}\approx 22$, which is approximately the magnitude limit of the coadd subtractions. However, in an r-band stack of images from August 24–26 (inclusive), we detect emission at ${m}_{r}\sim 21.5$ at 5σ above the background.

Assuming that the rapid rise we detected was close to the time of explosion, this is the first definitive detection of pre-explosion emission in a Ic-BL SN. There was a tentative detection in another source, PTF 11qcj (Corsi et al. 2014), 1.5 and 2.5 yr prior to the SN. In Section 4 we discuss possible mechanisms for this emission, and conclude that it is likely related to a period of eruptive mass loss immediately prior to the explosion. We note that it is unlikely that this variability arises from active galactic nucleus (AGN) activity, due to the properties of the host galaxy (Section 3.3).

With forced photometry and faint detections from stacked images and deep references, we can construct a light curve that extends weeks prior to the rapid rise in the light curve, shown in Figure 7.

Zoom In Zoom Out Reset image size

We observed the field of SN2018gep with the Karl G. Jansky Very Large Array (VLA) on three epochs: on 2018 September 14 UT under the Program ID VLA/18A-242 (PI: D. Perley; Ho et al. 2018b), and on 2018 September 25 and 2018 November 23 UT under the Program ID VLA/18A-176 (PI: A. Corsi). We used 3C286 for flux calibration, and J1640+3946 for gain calibration. The observations were carried out in X- and Ku-band (nominal central frequencies of 9 GHz and 14 GHz, respectively) with a nominal bandwidth of 2 GHz. The data were calibrated using the automated VLA calibration pipeline available in the CASA package (McMullin et al. 2007) then inspected for further flagging. The CLEAN procedure (Högbom 1974) was used to form images in interactive mode. The image rms and the radio flux at the location of SN2018gep were measured using imstat in CASA. Specifically, we report the maximum flux within pixels contained in a circular region centered on the optical position of SN2018gep with radius comparable to the FWHM of the VLA synthesized beam at the appropriate frequency. The source was detected in the first two epochs, but not in the third (see Table 2). As we discuss in Section 4, the first two epochs were conducted in a different array configuration than the third epoch, and may have had a contribution from host galaxy light.

Table 2.  Radio Flux Density Measurements for SN2018gep

Start Time	Δt	Instrument	ν	fν	Lν	${\theta }_{\mathrm{FWHM}}$	Int. Time
(UTC)	(days)		(GHz)	(μJy)	(erg ${{\rm{s}}}^{-1}$ ${\mathrm{Hz}}^{-1}$)	''	(hr)
2018-09-12 17:54	3.6	AMI	15	<120	$\lt 2.9\times {10}^{27}$	$43.53\times 30.85$	4
2018-09-23 15:35	14.5	AMI	15	<120	$\lt 2.9\times {10}^{27}$	$39.3\times 29.29$	4
2018-10-20 14:01	41.4	AMI	15	<120	$\lt 2.9\times {10}^{27}$	$43.53\times 30.85$	4
2018-09-15 02:33	6.0	SMA	230	<590	$\lt 1.4\times {10}^{28}$	$4.828\times 3.920$	1.25
2018-09-14 01:14	4.9	VLA	9.7	34 ± 4	$8.3\times {10}^{26}$	$7.06\times 5.92$	0.5
2018-09-25 00:40	15.9	VLA	9	24.4 ± 6.8	$6.0\times {10}^{26}$	$7.91\times 6.89$	0.7
2018-09-25 00:40	15.9	VLA	14	26.8 ± 6.8	$6.6\times {10}^{26}$	$4.73\times 4.26$	0.5
2018-11-23 13:30	75.4	VLA	9	<16	$\lt 3.9\times {10}^{26}$	$3.52\times 2.08$	0.65
2018-11-23 13:30	75.4	VLA	14	<17	$\lt 4.2\times {10}^{26}$	$2.77\times 1.32$	0.65

Note. For VLA measurements: the quoted errors are calculated as the quadrature sums of the image rms, plus a 5% nominal absolute flux calibration uncertainty. When the peak flux density within the circular region is less than three times the rms, we report an upper limit equal to three times the rms of the image. For AMI measurements: nondetections are reported as 3σ upper limits. For SMA measurements: nondetections are reported as a 1σ upper limit.

Download table as:  ASCIITypeset image

We also obtained three epochs of observations with the AMI large array (AMI-LA; Zwart et al. 2008; Hickish et al. 2018), on UT 2018 September 12, 2018 September 23, and 2018 October 20. AMI-LA is a radio interferometer comprised of eight 12.8 m diameter antennas that extend from 18 m up to 110 m in length and operates with a 5 GHz bandwidth around a central frequency of 15.5 GHz.

We used a custom AMI data reduction software package reduce_dc (e.g., Perrott et al. 2013) to perform initial data reduction, flagging, and calibration of phase and flux. Phase calibration was conducted using short interleaved observations of J1646+4059, and for absolute flux calibration we used 3C286. Additional flagging and imaging were performed using CASA. All three observations resulted in null detections with 3σ upper limits of $\approx 120\,\mu$ Jy in the first two observations, and a 3σ upper limit of $\approx 120\,\mu$ Jy in the last observation.

Finally, we observed at higher frequencies using the Submillimeter Array (SMA; Ho et al. 2004) on UT 2018 September 15 under its target-of-opportunity program. The project ID was 2018A-S068. Observations were performed in the sub-compact configuration using seven antennas. The observations were performed using RxA and RxB receivers tuned to LO frequencies of 225.55 GHz and 233.55 GHz, respectively, providing 32 GHz of continuous bandwidth ranging from 213.55 to 245.55 GHz with a spectral resolution of 140.0 kHz per channel. The atmospheric opacity was around 0.16–0.19 with system temperatures around 100–200 K. The nearby quasars 1635+381 and 3C345 were used as the primary phase and amplitude gain calibrators with absolute flux calibration performed by comparison to Neptune. Passband calibration was derived using 3C454.3. Data calibration was performed using the MIR IDL package for the SMA, with subsequent analysis performed in MIRIAD (Sault et al. 1995). For the flux measurements, all spectral channels were averaged together into a single continuum channel and an rms of 0.6 mJy was achieved after 75 minutes on-source.

The full set of radio and submillimeter measurements are listed in Table 2.

We observed the position of SN2018gep with Swift/XRT from Δt ≈ 0.4–14 days. The source was not detected in any epoch. To measure upper limits, we used web-based tools developed by the Swift-XRT team (Evans et al. 2009). For the first epoch, the 3σ upper limit was 0.003 ct s⁻¹. To convert the upper limit from count rate to flux, we assumed³² a Galactic neutral hydrogen column density of $1.3\times {10}^{20}\,{\mathrm{cm}}^{-2}$, and a power-law spectrum with photon index ${\rm{\Gamma }}=2$. This gives³³ an unabsorbed 0.3–10 keV flux of $\lt 9.9\times {10}^{-14}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}$, and ${L}_{X}\lt 2.5\times {10}^{41}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$.

We obtained two epochs of observations with the Advanced CCD Imaging Spectrometer (Garmire et al. 2003) on the Chandra X-ray Observatory via our approved program (Proposal No. 19500451; PI: Corsi). The first epoch began at 9:25 UTC on 2018 10 October (${\rm{\Delta }}t\approx 15\,\mathrm{days}$) under ObsId 20319 (integration time 12.2 ks), and the second began at 21:31 UTC on 2018 December 4 (${\rm{\Delta }}t\approx 70\,\mathrm{days}$) under ObsId 20320 (integration time 12.1 ks). No X-ray emission is detected at the location of SN2018gep in either epoch, with 90% upper limits on the 0.5–7.0 keV count rate of $\approx 2.7\times {10}^{-4}\,\mathrm{ct}\,{{\rm{s}}}^{-1}$. Using the same values of hydrogen column density and power-law photon index as in our XRT measurements, we find upper limits on the unabsorbed 0.5–7 keV X-ray flux of $\lt 3.2\times {10}^{-15}$ erg cm⁻² s⁻¹, or (for a direct comparison to the XRT band) a 0.3–10 keV X-ray flux of $\lt 4.2\times {10}^{-15}$ erg cm⁻² s⁻¹. This corresponds to a 0.3–10 keV luminosity upper limit of ${L}_{X}\lt 1.0\times {10}^{40}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$.

We created a tool to search for prompt gamma-ray emission (GRBs) from Fermi-Gamma-ray Burst Monitor (GBM) (Gruber et al. 2014; von Kienlin et al. 2014; Narayana Bhat et al. 2016), the Swift Burst Alert Telescope (BAT; Barthelmy et al. 2005), and the IPN, which we have made available online.³⁴ We did not find any GRB consistent with the position and t₀ of SN2018gep.

Our tool also determines whether a given position was visible to BAT and GBM at a given time, using the spacecraft pointing history. We use existing code³⁵ to determine the BAT history. We find that the position of SN2018gep was in the BAT field of view from UTC 03:13:40 to 03:30:38, before Swift slewed to another location.

We also find that at t₀ SN2018gep was visible to the Fermi GBM (Meegan et al. 2009). We ran a targeted GRB search in 10–1000 keV Fermi/GBM data from three hours prior to t₀ to half an hour after t₀. We use the soft template, which is a smoothly broken power law with low-energy index −1.9 and high-energy index −2.7, and an SED peak at 70 keV. The search methodology (and parameters of the other templates) are described in Blackburn et al. (2015) and Goldstein et al. (2016). No signals with a consistent location were found. For the 100 s integration time, the fluence upper limit is $2\times {10}^{-6}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}$. This limit corresponds to a 10–1000 keV isotropic energy release of ${E}_{\gamma ,\mathrm{iso}}\lt 4.9\times {10}^{48}\,\mathrm{erg}$. Limits for different spectral templates and integration times are shown in Figure 8.

Zoom In Zoom Out Reset image size

We measure line fluxes using the Keck optical spectrum obtained at ${\rm{\Delta }}t\approx 61\,\mathrm{days}$ (Figure 25). We model the local continuum with a low-order polynomial and each emission line by a Gaussian profile of FHWM ∼5.3 Å. This is appropriate if Balmer absorption is negligible, which is generally the case for starburst galaxies. For the host of SN2018gep, the Balmer decrement between Hβ, Hγ, and Hδ does not show any excess with respect to the expected values in Osterbrock & Ferland (2006). The resulting line fluxes are listed in Table 7.

We retrieved archival images of the host galaxy from Galaxy Evolution Explorer (GALEX) Data Release (DR) 8/9 (Martin et al. 2005), SDSS DR9 (Ahn et al. 2012), Panoramic Survey Telescope And Rapid Response System (PanSTARRS, PS1) DR1 (Chambers et al. 2016), Two-Micron All Sky Survey (2MASS; Skrutskie et al. 2006), and Wide-Field Infrared Survey Explorer (WISE; Wright et al. 2010). We also used UVOT photometry from Swift, and NIR photometry from the Canada–France–Hawaii Telescope Legacy Survey (CFHTLS; Hudelot et al. 2012).

The images are characterized by different pixel scales (e.g., SDSS $0\buildrel{\prime\prime}\over{.} 40$/px, GALEX $1^{\prime\prime}$/px) and different point spread functions (e.g., SDSS/PS1 1''–2'', WISE/W2 $6\buildrel{\prime\prime}\over{.} 5$). To obtain accurate photometry, we use the matched-aperture photometry software package Lambda Adaptive Multi-Band Deblending Algorithm in R (LAMBDAR; Wright et al. 2016) that is based on a photometry software package developed by Bourne et al. (2012). To measure the total flux of the host galaxy, we defined an elliptical aperture that encircles the entire galaxy in the SDSS/$r^{\prime}$-band image. This aperture was then convolved in LAMBDAR with the PSF of a given image that we specified directly (GALEX and WISE data) or that we approximated by a two-dimensional Gaussian (2MASS, SDSS and PS1 images). After instrumental magnitudes were measured, we calibrated the photometry against instrument-specific zero-points (GALEX, SDSS, and PS1 data), or, as in the case of 2MASS and WISE images, against a local sequence of stars from the 2MASS Point Source Catalogue and the AllWISE catalog. The photometry from the UVOT images was extracted with the command uvotsource in HEAsoft and a circular aperture with a radius of 8''. The photometry of the CFHT/WIRCAM data was performed with the software tool presented in Schulze et al. (2018b).³⁶ To convert the 2MASS, UVOT, WIRCAM, and WISE photometry to the AB system, we applied the offsets reported in Blanton & Roweis (2007), Breeveld et al. (2011), and Cutri et al. (2013). The resulting photometry is summarized in Table 8.

By interpolating the UVOT and ground-based photometry, we construct multiband SEDs and fit a Planck function on each epoch to measure the evolution of luminosity, radius, and effective temperature. To estimate the uncertainties, we perform a Monte Carlo simulation with 600 trials, each time adding noise corresponding to a 15% systematic uncertainty on each data point, motivated by the need to obtain a combined ${\chi }^{2}$/dof ∼ 1 across all epochs. The uncertainties for each parameter are taken as the 16-to-84 percentile range from this simulation. The SED fits are shown in Appendix A, and the resulting evolution in bolometric luminosity, photospheric radius, and effective temperature is listed in Table 3. We plot the physical evolution in Figure 9, with a comparison to iPTF16asu and AT2018cow.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3.  Physical Evolution of AT2018gep from Blackbody Fits

${\rm{\Delta }}t$	$L({10}^{10}{L}_{\odot })$	R (au)	T (kK)
0.05	${0.04}_{-0.02}^{+0.04}$	${21}_{-6}^{+14}$	${13}_{-4}^{+5}$
0.48	${7.4}_{-4.1}^{+8.6}$	${22}_{-5}^{+7}$	${46}_{-13}^{+16}$
0.73	${4.5}_{-2.8}^{+5.5}$	${31}_{-6}^{+11}$	${35}_{-11}^{+12}$
1.0	${2.2}_{-1.2}^{+2.1}$	${46}_{-9}^{+18}$	${24}_{-6}^{+6}$
1.7	${3.5}_{-2.1}^{+4.2}$	${46}_{-10}^{+22}$	${27}_{-8}^{+9}$
2.7	${1.3}_{-0.4}^{+1.2}$	${78}_{-20}^{+22}$	${16}_{-3}^{+5}$
3.2	${3.5}_{-1.3}^{+2.2}$	${50}_{-8}^{+14}$	${26}_{-5}^{+6}$
3.8	${2.9}_{-0.8}^{+1.7}$	${56}_{-11}^{+11}$	${23}_{-3}^{+5}$
4.7	${1.7}_{-0.3}^{+0.7}$	${69}_{-14}^{+16}$	${18}_{-2}^{+3}$
5.9	${0.88}_{-0.08}^{+0.17}$	${100}_{-21}^{+14}$	${13}_{-0}^{+1}$
8.6	${0.46}_{-0.06}^{+0.08}$	${220}_{-39}^{+46}$	${7.4}_{-0.5}^{+0.6}$
9.6	${0.33}_{-0.03}^{+0.04}$	${200}_{-24}^{+33}$	${7.1}_{-0.4}^{+0.4}$
10.0	${0.31}_{-0.03}^{+0.04}$	${210}_{-28}^{+34}$	${6.9}_{-0.4}^{+0.4}$
11.0	${0.28}_{-0.03}^{+0.04}$	${220}_{-33}^{+35}$	${6.5}_{-0.3}^{+0.4}$
13.0	${0.25}_{-0.03}^{+0.04}$	${260}_{-42}^{+50}$	${5.8}_{-0.3}^{+0.3}$
14.0	${0.22}_{-0.03}^{+0.04}$	${270}_{-47}^{+60}$	${5.5}_{-0.3}^{+0.4}$
16.0	${0.17}_{-0.03}^{+0.04}$	${260}_{-58}^{+76}$	${5.3}_{-0.5}^{+0.5}$
18.0	${0.15}_{-0.02}^{+0.04}$	${300}_{-64}^{+77}$	${4.7}_{-0.4}^{+0.4}$
21.0	${0.11}_{-0.02}^{+0.03}$	${250}_{-58}^{+83}$	${4.7}_{-0.4}^{+0.4}$
25.0	${0.073}_{-0.013}^{+0.02}$	${240}_{-85}^{+95}$	${4.5}_{-0.5}^{+0.9}$
38.0	${0.034}_{-0.007}^{+0.012}$	${180}_{-55}^{+86}$	${4.2}_{-0.5}^{+0.6}$

A machine-readable version of the table is available.

Download table as:  DataTypeset image

The bolometric luminosity peaks between ${\rm{\Delta }}t=0.5\,\mathrm{day}$ and ${\rm{\Delta }}t=3\,\mathrm{days}$, at $\gt 3\times {10}^{44}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$. In Figure 10 we compare this peak luminosity and time to peak luminosity with several classes of stellar explosions. As in iPTF16asu, the bolometric luminosity falls as an exponential at late times ($t\gt 10\,\mathrm{days}$). The total integrated UV and optical (≈2000–9000 Å) blackbody energy output from Δt = 0.5–40 days is $\sim {10}^{50}\,\mathrm{erg}$, similar to that of iPTF16asu. The earliest photospheric radius we measure is ∼20 au, at ${\rm{\Delta }}t=0.05\,\mathrm{day}$. Until ${\rm{\Delta }}t\approx 17\,\mathrm{days}$ the radius expands over time with a very large inferred velocity of v ≈ 0.1c. After that, it remains flat, and even appears to recede. This possible recession corresponds to a flattening in the temperature at $\sim 5000\,{\rm{K}}$, which is the recombination temperature of carbon and oxygen. This effect was not seen in iPTF16asu, which remained hotter (and more luminous) for longer. Finally, the effective temperature rises before falling as ∼t⁻¹. We interpret these properties in the context of shock-cooling emission in Section 4.

3.2.1. Comparisons to Early Spectra in the Literature

We obtained nine spectra of SN2018gep in the first five days after discovery. These early spectra are shown in Figure 11, when the effective temperature declined from 50,000 K to 20,000 K. To our knowledge, our early spectra have no analogs in the literature, in that there has never been a spectrum of a stripped-envelope SN at such a high temperature (excluding spectra during the afterglow phase of GRBs).³⁷ Two of the earliest spectra in the literature, one at ${\rm{\Delta }}t=2\,\mathrm{days}$ for Type Ic SN PTF10vgv (Corsi et al. 2012) and one at ${\rm{\Delta }}t=3\,\mathrm{days}$ for Type Ic SN PTF12gzk (Ben-Ami et al. 2012) are redder and exhibit more features than the spectrum of SN2018gep. We show the comparison in Figure 11.

Zoom In Zoom Out Reset image size

At ${\rm{\Delta }}t\approx 4\,\mathrm{days}$, a "W" feature emerges in the rest-frame wavelength range 3800–4350 Å. In the second-from-bottom panel of Figure 11 we make a comparison to "W" features seen in SN 2008D (e.g., Modjaz et al. 2009), which was a Type Ib SN associated with an X-ray flash (Mazzali et al. 2008), and in a typical pre-max stripped-envelope SLSN (Type I SLSN; Moriya et al. 2018; Gal-Yam 2019b). The absorption lines are broadened much more than in PTF12dam (Nicholl et al. 2013) and probably more than in SN2008D as well. Finally, SN2018gep cooled more slowly than SN2008D: only after 4.25 days did it reach the temperature that SN 2008D reached after >2 days.

3.2.2. Origin of the "W" feature

The lack of comparison data at such early epochs (high temperatures) motivated us to model one of the early spectra in order to determine the composition and density profile of the ejecta. We used the spectral synthesis code JEKYLL (Ergon et al. 2018), configured to run in steady state using a full NLTE solution. An inner blackbody boundary was placed at a high continuum optical depth (∼50), and the temperature at this boundary was iteratively determined to reproduce the observed luminosity. The atomic data used is based on what was specified in Ergon et al. (2018), but has been extended as described in Appendix C. We explored models with C/O (mass fractions: 0.23/0.65) and O/Ne/Mg (mass fractions: 0.68/0.22/0.07) compositions taken from a model by Woosley & Heger (2007)³⁸ and a power-law density profile, where the density at the inner border was adjusted to fit the observed line velocities. Except for the density at the inner border, various power-law indices were also explored, but in the end an index of −9 worked out best.

Figures 12 and 13 show the model with the best overall agreement with the spectra and the SED (as listed in Table 6 the spectrum was obtained at high airmass, making it difficult to correct for telluric features). The model has a C/O composition, an inner border at 22,000 km ${{\rm{s}}}^{-1}$ (corresponding to an optical depth of ∼50), a density of 4 × 10⁻¹² g cm⁻³ at this border, and a density profile with a power-law index of −9. In Figure 12 we show that the model does a good job of reproducing both the spectrum and the SED of SN2018gep. In particular, it is interesting to note that the "W" feature seems to arise naturally in C/O material at the observed conditions. A similar conclusion was reached by Dessart (2019), whose magnetar-powered SLSN-I models, calculated using the NLTE code CMFGEN, show the "W" feature even when nonthermal processes were not included in the calculation (as in our case).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the model, the "W" feature mainly arises from the O ii 2p²(3P)3 s 4P $\leftrightarrow$ 2p²(3P)3p 4D° (4639–4676 Å), O ii 2p²(3P)3 s 4P $\leftrightarrow$ 2p²(3P)3p 2D (4649 Å), and O ii 2p²(3P)3 s 4P $\leftrightarrow$ 2p²(3P)3p 4P° (4317–4357 Å) transitions. The departure from LTE is modest in the line-forming region, and the departure coefficients for the O ii states are small. The spectrum redward of the "W" feature is shaped by carbon lines, and the features near 5700 and 6500 Å arise from the C ii 3 s 2S $\leftrightarrow$ 3p 2P° (6578, 6583 Å) and C iii 2s3p 1P° $\leftrightarrow$ 2s3d 1D (5696 Å) transitions, respectively. In the model, the C ii feature is too weak, suggesting that the ionization level is too high in the model. There is also a contribution from the C iii 2s3s 3S $\leftrightarrow$ 2s3p 3P° (4647–4651 Å) transition to the red part of the "W" feature, which could potentially be what is seen in the spectra from earlier epochs. In addition, there is a contribution from Si iv 4 s 2S $\leftrightarrow$ 4p 2P° (4090, 4117 Å) near the blue side of the "W" feature, which produces a distinct feature in models with lower velocities and which could explain the observed feature on the blue side of the "W" feature.

In spite of the overall good agreement, there are also some differences between the model and the observations. In particular the model spectrum is bluer and the velocities are higher. These two quantities are in tension and a better fit to one of them would result in a worse fit to the other. As mentioned above, the ionization level might be too high in the model, which suggests that the temperature might be too high as well. It should be noted that adding host extinction (which is assumed to be zero) or reducing the distance (within the error bars) would help in making the model redder (in the observer frame), and the latter would also help in reducing the temperature. The (modest) differences between the model and the observations could also be related to physics not included in the model, like a nonhomologous velocity field, departures from spherical asymmetry, and clumping.

The total luminosity of the model is 6.2 × 10⁴³ erg s⁻¹, the photosphere is located at ∼33,000 km ${{\rm{s}}}^{-1}$, and the temperature at the photosphere is ∼17,500 K, which is consistent with the values estimated from the blackbody fits (although the blackbody radius and temperature fits refer to the thermalization layer). As mentioned, we have also tried models with an O/Ne/Mg composition. However, these models failed to reproduce the carbon lines redwards of the "W" feature. We therefore conclude that the (outer) ejecta probably have a C/O-like composition, and that this composition in combination with a standard power-law density profile reproduces the spectrum of SN2018gep at the observed conditions (luminosity and velocity) 4.2 days after explosion.

In our model, the broad feature seen in our Swift UVOT grism spectrum is dominated by the strong Mg ii (2796, 2803 Å) resonance line. However, a direct comparison is not reliable because the ionization is probably lower at this epoch than what we consider for our model.

3.2.3. Photospheric Velocity from Ic-BL Spectra

At ${\rm{\Delta }}t\gtrsim 7.8\,\mathrm{days}$, the spectra of SN2018gep qualitatively resemble those of a stripped-envelope SN. We measure velocities using the method in Modjaz et al. (2016), which accommodates blending of the Fe ii λ5169 line (which has been shown to be a good tracer of photospheric velocity; Branch et al. 2002) with the nearby Fe ii λλ4924, 5018 lines.

At earlier times, when the spectra do not resemble typical Ic-BL SNe, we use our line identifications of ionized C and O to measure velocities. As shown in Figure 14, the velocity evolution we measure is comparable to that seen in Ic-BL SNe associated with GRBs (more precisely, low-luminosity GRBs; LLGRBs) which are systematically higher than those of Ic-BL SNe lacking GRBs (Modjaz et al. 2016). However, as discussed in Section 2.7, no GRB was detected.

Zoom In Zoom Out Reset image size

We infer a star formation rate of $0.09\pm 0.01\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ from the Hα emission line using the Kennicutt (1998) relation converted to use a Chabrier initial mass function (IMF); (Chabrier 2003; Madau & Dickinson 2014). We note that this is a lower limit as the slit of the Keck observation did not enclose the entire galaxy. We estimate a correction factor of 2–3: the slit diameter in the Keck spectra was 10, and the extraction radius was $\sim 1\buildrel{\prime\prime}\over{.} 75$ in the February observation and $\sim 1.21^{\prime\prime}$ in the March observation. The host diameter is roughly 4''.

We derive an electron temperature of $13,{100}_{-1000}^{+900}\,{\rm{K}}$ from the flux ratio between [O iii] λ4641 and [O iii] λ5007, using the software package PyNeb version 1.1.7 (Luridiana et al. 2015). In combination with the flux measurements of [O ii] $\lambda \lambda$ 3226, 3729, [O iii] λ4364, [O iii] λ4960, [O iii] λ5008, and Hβ, we infer a total oxygen abundance of ${8.01}_{-0.09}^{+0.10}$ (statistical error; using Equations (3) and (5) in Izotov et al. 2006). Assuming a solar abundance of 8.69 (Asplund et al. 2009), the metallicity of the host is ∼20% solar.

We also compute the oxygen abundance using the strong-line metallicity indicator O3N2 (Pettini & Pagel 2004) with the updated calibration reported in Marino et al. (2013). The oxygen abundance in the O3N2 scale is $8.05\pm 0.01(\mathrm{stat})\,\pm 0.10(\mathrm{sys})$.³⁹

We also estimate mass and star formation rate by modeling the host SED; see Appendix D for a table of measurements, and details on where we obtained them. We use the software package LePhare version 2.2 (Arnouts et al. 1999; Ilbert et al. 2006). We generated 3.9 × 10⁶ templates based on the Bruzual & Charlot (2003) stellar population synthesis models with the Chabrier IMF (Chabrier 2003). The star formation history (SFH) was approximated by a declining exponential function of the form $\exp \left(t/\tau \right)$, where t is the age of the stellar population and τ the e-folding timescale of the SFH (varied in nine steps between 0.1 and 30 Gyr). These templates were attenuated with the Calzetti attenuation curve (Calzetti et al. 2000) varied in 22 steps from $E(B-V)=0$ to 1 mag.

As shown in Figure 15, the SED is well characterized by a galaxy mass of $\mathrm{log}M/{M}_{\odot }={8.11}_{-0.08}^{+0.07}$ and an attenuation-corrected star formation rate of ${0.12}_{-0.05}^{+0.08}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. The derived star formation rate is comparable to the measurement inferred from Hα. The attenuation of the SED is marginal, with $E{(B-V)}_{\mathrm{star}}=0.05$, and consistent with the negligible Balmer decrement (Section 2.8).

Zoom In Zoom Out Reset image size

Figure 16 shows that the host galaxy of SN2018gep is even more low-mass and metal-poor than the typical host galaxies of Ic-BL SNe, which are low-mass and metal-poor compared to the overall core-collapse SN population. The figure uses data for 28 Ic-BL SNe from PTF and iPTF (Modjaz et al. 2019; Taddia et al. 2019) and a sample of 11 long-duration GRBs (including LLGRBs, all at z < 0.3). We measured the emission lines from the spectra presented in Taddia et al. (2019) and used line measurements reported in Modjaz et al. (2019) for objects with missing line fluxes. The photometry was taken from S. Schulze et al. (2019, in preparation). Photometry and spectroscopy were taken from a variety of sources.⁴⁰ The oxygen abundances were measured in the O3N2 scale like for SN2018gep and their SEDs were modeled with the same set of galaxy templates. For reference, the mass and SFR of the host of AT2018cow was $1.4\times {10}^{9}\,{M}_{\odot }$ and $0.22\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, respectively (Perley et al. 2019). The mass and SFR of the host of iPTF16asu was ${4.6}_{-2.3}^{+6.5}\times {10}^{8}\,{M}_{\odot }$ and $0.7\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, respectively (Whitesides et al. 2017).

Zoom In Zoom Out Reset image size

The majority of stripped-envelope SNe have light curves powered by the radioactive decay of ${}^{56}\mathrm{Ni}$. As discussed in Kasen (2017), this mechanism can be ruled out for light curves that rise rapidly to a high peak luminosity, because this would require the unphysical condition of a nickel mass that exceeds the total ejecta mass. With a peak luminosity exceeding ${10}^{44}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ and a rise to peak of a few days, SN2018gep clearly falls into the disallowed region (see Figure 1 in Kasen 2017). Thus, we rule out radioactive decay as the mechanism powering the peak of the light curve.

We now consider whether radioactive decay could dominate the light curve at late times ($t\gg {t}_{\mathrm{peak}}$). The left panel of Figure 17 shows the bolometric light curve of SN2018gep compared to several other Ic-BL SNe from the literature (Cano 2013), whose light curves are thought to be dominated by the radioactive decay of ${}^{56}\mathrm{Ni}$ (although see Moriya et al. 2017 for another possible interpretation). The luminosity of SN2018gep at $t\sim 20\,\mathrm{days}$ is about half that of SN1998bw, and double that of SN2010bh and SN2006aj. By modeling the light curves of the three Ic-BL SNe shown, Cano (2013) infers nickel masses of 0.42 ${M}_{\odot }$, 0.12 ${M}_{\odot }$, and 0.21 ${M}_{\odot }$, respectively. On this scale, SN2018gep has ${M}_{\mathrm{Ni}}\sim 0.1$–$0.2\,{M}_{\odot }$.

Zoom In Zoom Out Reset image size

The right panel of Figure 17 shows the light curve of SN2018gep compared to that of AT2018cow (Perley et al. 2019). To estimate the nickel mass of AT2018cow, Perley et al. (2019) compared the bolometric luminosity at $t\sim 20\,\mathrm{days}$ to that of SN2002ap (whose nickel mass was derived via late-time nebular spectroscopy; Foley et al. 2003) and found ${M}_{\mathrm{Ni}}\lt 0.05\,{M}_{\odot }$. On this scale, we would expect ${M}_{\mathrm{Ni}}\lesssim 0.05\,{M}_{\odot }$ for SN2018gep as well.

Finally, Katz et al. (2013) and Wygoda et al. (2019) present an analytical technique for testing whether a light curve is powered by radioactive decay. At late times, the bolometric luminosity is equal to the rate of energy deposition by radioactive decay Q(t), because the diffusion time is much shorter than the dynamical time: ${L}_{\mathrm{bol}}(t)=Q(t)$. At any given time, the energy deposition rate Q(t) is

$$\begin{eqnarray}&&Q(t)={Q}_{\gamma }(t)\left(1-{e}^{-{\left({t}_{0}/t\right)}^{2}}\right)+{Q}_{\mathrm{pos}}(t)\end{eqnarray} \tag{ 1 }$$

where ${Q}_{\gamma }(t)$ is the energy release rate of gamma-rays and t₀ is the time at which the ejecta becomes optically thin to gamma-rays. The expression for ${Q}_{\gamma }(t)$ is

$$\begin{eqnarray}&&\displaystyle \frac{{Q}_{\gamma }(t)}{{10}^{43}\,\mathrm{erg}\,{{\rm{s}}}^{-1}}=\displaystyle \frac{{M}_{\mathrm{Ni}}}{{M}_{\odot }}\left(6.45{e}^{-t/8.76\mathrm{days}}+1.38{e}^{-t/111.4\mathrm{days}}\right).\end{eqnarray} \tag{ 2 }$$

${Q}_{\mathrm{pos}}(t)$ is the energy deposition rate of positron kinetic energy, and the expression is

$$\begin{eqnarray}&&\displaystyle \frac{{Q}_{\mathrm{pos}}(t)}{{10}^{41}\,\mathrm{erg}\,{{\rm{s}}}^{-1}}=4.64\displaystyle \frac{{M}_{\mathrm{Ni}}}{{M}_{\odot }}\left(-{e}^{-t/8.76\mathrm{days}}+{e}^{-t/111.4\mathrm{days}}\right).\end{eqnarray} \tag{ 3 }$$

The dotted line in Figure 17 shows a model track with ${M}_{\mathrm{Ni}}=0.28\,{M}_{\odot }$ and ${t}_{0}=30\,\mathrm{days}$. Lower nickel masses produce tracks that are too low to reproduce the data, and larger values of t₀ produce tracks that drop off too rapidly. Thus on this scale it seems that ${M}_{\mathrm{Ni}}\sim 0.3\,{M}_{\odot }$, similar to other Ic-BL SNe (Lyman et al. 2016). Because the data have not yet converged to model tracks, we cannot solve directly for t₀ and ${M}_{\mathrm{Ni}}$ using the technique for Ia SNe in Wygoda et al. (2019).

We can also try to solve directly for t₀ and ${M}_{\mathrm{Ni}}$ using the technique for Ia SNe in Wygoda et al. (2019). The first step is to solve for t₀ using Equation (1) and a second equation resulting from the fact that the expansion is adiabatic,

$$\begin{eqnarray}&&{\int }_{0}^{t}Q(t^{\prime} )\,t^{\prime} \,{dt}^{\prime} ={\int }_{0}^{t}{L}_{\mathrm{bol}}(t^{\prime} )\,t^{\prime} \,{dt}^{\prime} .\end{eqnarray} \tag{ 4 }$$

The ratio of Equation (1) to Equation (4) removes the dependence on ${M}_{\mathrm{Ni}}$, and enables t₀ to be measured. However, as shown in Figure 18, the data have not yet converged to model tracks.

Zoom In Zoom Out Reset image size

One way to power a rapid and luminous light curve is to deposit energy into CSM at large radii (Nakar & Sari 2010; Nakar & Piro 2014; Piro 2015). Since this is a Ic-BL SN, we expect the progenitor to be stripped of its envelope and therefore compact ($R\sim 0.5\,{R}_{\odot }\sim {10}^{10}\,\mathrm{cm};$ Groh et al. 2013), although there have never been any direct progenitor detections for a Ic-BL SN.

With this expectation, extended material at larger radii would have to arise from mass loss. This would not be surprising, as massive stars are known to shed a significant fraction of their mass in winds and eruptive episodes; see Smith (2014) for a review.

First we perform an order-of-magnitude calculation to see whether the rise time and peak luminosity could be explained by a model in which shock interaction powers the light curve ("wind shock breakout"). Assuming that the progenitor ejected material with a velocity v_w at a time t prior to explosion, the radius of this material at any given time is

$$\begin{eqnarray}\begin{array}{rcl}{R}_{\mathrm{sh}} & = & {R}_{* }+{v}_{{\rm{w}}}t\\ & \approx & (8.64\times {10}^{12}\,\mathrm{cm})\left(\displaystyle \frac{{v}_{w}}{1000\,\mathrm{km}\,{{\rm{s}}}^{-1}}\right)\left(\displaystyle \frac{t}{\mathrm{day}}\right).\end{array}\end{eqnarray} \tag{ 5 }$$

For material ejected 15 days prior to explosion, traveling at 1000 km ${{\rm{s}}}^{-1}$, the radius would be ${R}_{\mathrm{CSM}}\sim {10}^{14}\,\mathrm{cm}$ at the time of explosion. The shock crossing timescale is t_cross:

$$\begin{eqnarray}&&{t}_{\mathrm{cross}}\sim {R}_{\mathrm{CSM}}/{v}_{s}\approx (0.4\,\mathrm{day})\left(\displaystyle \frac{R}{{10}^{14}\,\mathrm{cm}}\right){\left(\displaystyle \frac{{v}_{s}}{0.1c}\right)}^{-1}\end{eqnarray} \tag{ 6 }$$

where v_s is the velocity of the shock. The shock heats the CSM with an energy density that is roughly half of the kinetic energy of the shock, so ${e}_{s}\sim (1/2)(\rho {v}_{s}^{2}/2)$. The luminosity is the total energy deposited divided by t_cross,

$$\begin{eqnarray}\begin{array}{rcl}{L}_{\mathrm{BO}} & \sim & \displaystyle \frac{{E}_{\mathrm{BO}}}{{t}_{\mathrm{cross}}}\sim \displaystyle \frac{{v}_{s}^{3}}{4}\displaystyle \frac{{dM}}{{dR}}\\ & = & (8\times {10}^{44}\,\mathrm{erg}\,{{\rm{s}}}^{-1}){\left(\displaystyle \frac{{v}_{s}}{0.1c}\right)}^{3}\left(\displaystyle \frac{{dM}}{{M}_{\odot }}\right){\left(\displaystyle \frac{{dR}}{{10}^{14}\mathrm{cm}}\right)}^{-1}\end{array}\end{eqnarray} \tag{ 7 }$$

assuming a constant density. Thus, for shock velocities on the order of the observed photospheric radius expansion (0.1c), and a CSM radius on the order of the first photospheric radius that we measure ($3\times {10}^{14}\,\mathrm{cm}$), it is easy to explain the rise time and peak luminosity that we observe.

To test whether shock breakout (and subsequent post-shock cooling) can explain the evolution of the physical properties we measured in Section 3, we ran one-dimensional numerical radiation hydrodynamics simulations of an SN running into a circumstellar shell with CASTRO (Almgren et al. 2010; Zhang et al. 2011). We assume spherical symmetry and solve the coupled equations of radiation hydrodynamics using a gray flux-limited nonequilibrium diffusion approximation. The setup is similar to the models presented in Rest et al. (2018) but with parameters modified to fit SN2018gep.

The ejecta is assumed to be homologously expanding, characterized by a broken power-law density profile, an ejecta mass ${M}_{\mathrm{ej}}$, and energy ${E}_{\mathrm{ej}}$. The ejecta density profile has an inner power-law index of n = 0 (that is, $\rho (r)\propto {r}^{-n}$) then steepens to an index n = 10, as is appropriate for core collapse SN explosions (Matzner & McKee 1999). The circumstellar shell is assumed to be uniform in density with radius ${R}_{\mathrm{CSM}}$ and mass ${M}_{\mathrm{CSM}}$. We adopt a uniform opacity of $\kappa =0.2$ cm² g⁻¹, which is characteristic of hydrogen-poor electron scattering.

The best-fit model, shown in Figure 19, used the following parameters: ${M}_{\mathrm{ej}}=8\,{M}_{\odot }$, ${E}_{\mathrm{ej}}=2\times {10}^{52}\,\mathrm{erg}$, ${M}_{\mathrm{CSM}}=0.02\,{M}_{\odot }$, and ${R}_{\mathrm{CSM}}=3\times {10}^{14}\,\mathrm{cm}$. The inferred kinetic energy is consistent with typical values measured for Ic-BL SNe (e.g., Cano et al. 2017; Taddia et al. 2019), and R_CSM is similar in value to the first photospheric radius we measure (at ${\rm{\Delta }}t=0.05\,\mathrm{day}$; see Figure 9).

Zoom In Zoom Out Reset image size

The inferred values presented here are likely uncertain to within a factor of a few, given the degeneracies of the rise time and peak luminosity with the CSM mass and radius. Qualitatively, a larger CSM radius will result in a higher peak luminosity and longer rise time. The peak luminosity is relatively independent of the CSM mass, which instead affects the photospheric velocity and temperature (i.e., a larger CSM mass slows down the post-interaction velocity to a greater extent and increases the shock-heated temperature). A full discussion of the dependencies of the light curve and photospheric properties on the CSM parameters will be presented in an upcoming work (D. Khatami et al. 2019, in preparation.).

In this framework, the shockwave sweeps through the CSM prior to peak luminosity, so that at maximum luminosity the outer parts of the CSM have been swept into a dense shell moving at SN-like velocities (${v}_{\mathrm{post}-\mathrm{shock}}\approx 3{v}_{s}/4$). This scenario was laid out in Chevalier & Irwin (2011) and discussed in Kasen (2017). This explains the high velocities we measure at early times and the absence of narrow emission features in our spectra. For another discussion of the absence of narrow emission lines due to an abrupt cutoff in CSM density, see Moriya & Tominaga (2012). Following Chevalier & Irwin (2011), the rapid rise corresponds to shock breakout from the CSM, and begins at a time ${R}_{\mathrm{CSM}}/{v}_{\mathrm{sh}}$ after the explosion, where v_sh is the velocity of the shock. The time to peak luminosity (1.2 days) is longer than this delay time by a factor (${R}_{w}/{R}_{d}$). Given the best-fit ${R}_{w}=3\,\times {10}^{14}\,\mathrm{cm}$, and assuming ${R}_{d}\sim {R}_{w}$, we find ${v}_{\mathrm{sh}}=0.1c$, and an explosion time $\sim 1\,\mathrm{day}$ prior to t₀. This model also predicts an increasing temperature while the shock breaks out (i.e., during the rise to peak bolometric luminosity).

Other Ic SNe have shown early evidence for interaction in their light curves, but in other cases the emission has been attributed to post-shock cooling in expanding material rather than shock breakout itself. For example, the first peak observed in iPTF14gqr (De et al. 2018) was short-lived ($\lesssim 2\,\mathrm{days}$) and attributed to shock-cooling emission from material stripped by a compact companion. iPTF14gqr is different in a number of ways from SN2018gep: the spectra showed high-ionization emission lines, including He ii, and the explosion had a much smaller kinetic energy (${E}_{K}\approx {10}^{50}\,\mathrm{erg}$) and smaller velocities (10,000 km ${{\rm{s}}}^{-1}$). The main peak in iPTF16asu was also modeled as shock-cooling emission rather than shock breakout (Whitesides et al. 2017).

Under the assumption that the light curve represented post-shock cooling emission, De et al. (2018) and Whitesides et al. (2017) both used one-zone analytic models from Piro (2015) to estimate the properties of the explosion and the CSM. This approximation assumes that the emitting region is a uniformly heated expanding sphere. In iPTF14gqr the inferred properties of the extended material were ${M}_{e}\sim 8\times {10}^{-3}\,{M}_{\odot }$ at ${R}_{e}\sim 3\,\times {10}^{13}\,\mathrm{cm}$. In iPTF16asu the inferred properties of the extended material were ${M}_{e}\sim 0.45\,{M}_{\odot }$ at ${R}_{e}\sim 1.7\,\times {10}^{12}\,\mathrm{cm}$. The fit also required a more energetic explosion than iPTF14gqr ($4\times {10}^{51}\,\mathrm{erg}$). By applying the same framework to the decline of the bolometric light curve of SN2018gep, we arrive at similar values to those inferred for iPTF16asu, as shown in Figure 20.

Zoom In Zoom Out Reset image size

We model the main peak of SN2018gep as shock breakout rather than post-shock cooling emission. Our motivation for this choice is that the timescale over which we detect the precursor emission is more consistent with a large radius and lower shell mass. From the shell mass and radius, we can also estimate the mass-loss rate immediately prior to explosion,

$$\begin{eqnarray}&&\displaystyle \frac{\dot{M}}{{M}_{\odot }\,{\mathrm{yr}}^{-1}}\approx 32\left(\displaystyle \frac{{M}_{\mathrm{sh}}}{{M}_{\odot }}\right)\left(\displaystyle \frac{{v}_{w}}{1000\mathrm{km}\,{{\rm{s}}}^{-1}}\right){\left(\displaystyle \frac{{R}_{\mathrm{sh}}}{{10}^{14}\mathrm{cm}}\right)}^{-1}.\end{eqnarray} \tag{ 8 }$$

For our best-fit parameters ${M}_{\mathrm{sh}}=0.02\,{M}_{\odot }$ and ${R}_{\mathrm{sh}}=3\,\times {10}^{14}\,\mathrm{cm}$, and taking ${v}_{w}=1000\,\mathrm{km}\,{{\rm{s}}}^{-1}$, we find $\dot{M}\,\approx 0.6\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, 4–6 orders of magnitude higher than what is typically expected for Ic-BL SNe (Smith 2014).

In the shock breakout model, the shock sweeps through confined CSM and passes into lower density material. Thus, it is not surprising that we do not observe the X-ray or radio emission that would indicate interaction with high-density material. From our VLA observations of SN2018gep, the radio flux marginally decreased from ${\rm{\Delta }}t=5\,\mathrm{days}$ to ${\rm{\Delta }}t=75\,\mathrm{days}$. This could be astrophysical, but could also be instrumental (change in beam size due to change in VLA configuration). Using the relation of Murphy et al. (2011), the estimated contribution from the host galaxy (for an SFR of ${0.12}_{-0.05}^{+0.08}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$; see Section 3.3) is

$$\begin{eqnarray}\begin{array}{rcl}\left(\displaystyle \frac{{L}_{1.4\mathrm{GHz}}}{\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1}}\right) & \approx & 1.57\times {10}^{28}\left(\displaystyle \frac{{\mathrm{SFR}}_{\mathrm{radio}}}{{M}_{\odot }\,{\mathrm{yr}}^{-1}}\right)\\ & \approx & 1.9\times {10}^{27}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1}.\end{array}\end{eqnarray} \tag{ 9 }$$

Taking a spectral index of −0.7 (a synchrotron spectrum), the expected 9 GHz luminosity would be between $3.0\,\times {10}^{26}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1}$ and $8.6\times {10}^{26}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1}$. From Table 2, the measured spectral luminosity is $8.3\times {10}^{26}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1}$ (at 10 GHz) in the first epoch, and $6\times {10}^{26}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1}$ (at 9 GHz) in the second epoch. The slit covering fraction of our LRIS observations is again relevant here; as discussed in Section 3.3, the true SFR is likely a factor of a few higher than what we inferred from modeling the galaxy SED. So, it is plausible that the first two radio detections are entirely due to the host galaxy.

In the third epoch, the luminosity (at 9 GHz) is $\lt 3.9\times {10}^{26}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1}$, although the difference from the first two epochs may be due to the different array configuration. Taking the peak of the 9–10 GHz light curve to be $8.3\times {10}^{26}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1}$ at ${\rm{\Delta }}t\approx 5\,\mathrm{days}$, Figure 21 shows that SN2018gep would be an order of magnitude less luminous in radio emission than any other Ic-BL SN. If the luminosity truly decreased, then the implied mass-loss rate is $\dot{M}\sim 3\times {10}^{-6}$, consistent with the idea that the shock has passed from confined CSM into much lower density material.

Zoom In Zoom Out Reset image size

If the emission is constant and due entirely to the host galaxy, the point shown in Figure 21 is an upper limit in luminosity. Assuming that the peak of the SED of any radio emission from the SN is not substantially different from the frequencies we measure (i.e., that the spectrum is not self-absorbed at these frequencies), we have a limit on the 9 GHz radio luminosity of ${L}_{p}\lesssim {10}^{27}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{Hz}}^{-1}$ at ${\rm{\Delta }}t\approx 5$–15 days.

The shell mass and radius also give an estimate of the optical depth: $\tau \approx \kappa M/{r}^{2}\approx 100\gt \gt 1$, which means that the shell would be optically thick. The lack of detected X-ray emission is consistent with the expectation that any X-ray photons produced in the collision would be thermalized by the shell and reradiated as blackbody emission.

Finally, assuming that the rapid rise to peak is indeed caused by shock breakout, we examine whether our model is consistent with our detections in the weeks prior to explosion. Material ejected 10 days prior to the explosion at the escape velocity of a Wolf-Rayet star (${v}_{\mathrm{esc}}\sim 1000\,\mathrm{km}\,{{\rm{s}}}^{-1}$) would lie at $R\sim {10}^{14}\,\mathrm{cm}$, which is consistent with our model. Assuming that the emission mechanism is internal shocks between shells of ejected material traveling at different velocities, we can estimate the amount of mass required:

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}\epsilon {{Mv}}^{2}=L\tau \end{eqnarray} \tag{ 10 }$$

where $v\approx 1000\,\mathrm{km}\,{{\rm{s}}}^{-1}$, $\epsilon \approx 0.5$ is the efficiency of thermalizing the kinetic energy of the shells, M is the shell mass, $L\,\approx {10}^{39}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ is the luminosity we observe, and $\tau \approx 10\,\mathrm{days}$ is the timescale over which we observe the emission. We find $M\,\approx 0.02\,{M}_{\odot }$, again consistent with our model.

We conclude that the data are consistent with a scenario in which a compact Ic-BL progenitor underwent a period of eruptive mass loss shortly prior to explosion. In the terminal explosion, the light curve was initially dominated by shock breakout through (and post-shock cooling of) this recently ejected material.

Finally, we return to the question of the emission detected in the first few minutes, which showed an inflection point prior to the rapid rise to peak (Figure 2). Given the pre-explosion activity and inference of CSM interaction, it is not surprising that the rise is not well-modeled by a simple quadratic function. One possibility is that we are seeing ejecta already heated from earlier precursor activity. Another possibility is that we are seeing the effects of a finite light travel time. For a sphere of $R\sim 3\times {10}^{14}\,\mathrm{cm}$, the light-crossing time is ∼20 minutes. The slower rising phase could represent the time for photons to reach us across the extent of the emitting sphere.

In Table 4, we summarize the key properties inferred from Section 4.

Table 4.  Key Model Properties of SN2018gep

Parameter	Value	Notes
${t}_{\mathrm{rise}}$	1.2 days
ESN	$2\times {10}^{52}\,\mathrm{erg}$
${M}_{\mathrm{ej}}$	8 ${M}_{\odot }$
MCSM	0.02 ${M}_{\odot }$
RCSM	$3\times {10}^{14}\,\mathrm{cm}$
$\dot{M}$	$0.6\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$	Assuming ${v}_{w}=1000\,\mathrm{km}\,{{\rm{s}}}^{-1}$
${M}_{\mathrm{Ni}}$	<0.2–0.3 ${M}_{\odot }$

Download table as:  ASCIITypeset image