SN 2013ai: A Link between Hydrogen-rich and Hydrogen-poor Core-collapse Supernovae

Core-collapse supernovae (CCSNe) are the result of the explosion of massive stars (>8 M_⊙). These explosions are characterized by two main classes: those that show no hydrogen and those that do, Type I (SNe I) and Type II supernovae (SNe II), respectively (Minkowski 1941).

Type I CCSNe are produced when a massive star is stripped of most of its outer hydrogen and possibly its helium layers before explosion; thus. they are referred to as stripped-envelope SNe (SESNe, Clocchiatti & Wheeler 1997). It is believed that the outer layers of the SESNe progenitors are primarily removed either through strong winds (Woosley et al. 1993) or binary interaction (Nomoto et al. 1995; Podsiadlowski et al. 2004). These SESNe are further divided into two general spectroscopic groups: SNe Ib, which show strong helium lines in their optical spectra, and SNe Ic, which do not. SNe Ib and Ic form spectroscopically homogeneous groups with minimal variation (Modjaz et al. 2014). The light curves of SNe Ib and Ic are similar and can be characterized by a slow ∼20 day rise to maximum followed by a postmaximum decline of ∼0.75 magnitudes per 15 days in the V-band (e.g., Drout et al. 2011; Taddia et al. 2018; Prentice et al. 2019).

When a massive star retains most of its hydrogen envelope, it explodes as an SN II. SNe II were historically divided into groups based on the shape of their light curves. Those with a slow decline—a plateau—after maximum were classified as IIP and those with no plateau after maximum were classified as IIL (Barbon et al. 1979). However, recent studies have shown that there is likely a continuum in the decline rates and suggest that slow (IIP) and fast decliners (IIL) are not distinct groups (Anderson et al. 2014; Sanders et al. 2015; Galbany et al. 2016; Valenti et al. 2016; Rubin & Gal-Yam 2016; Pessi et al. 2019; de Jaeger et al. 2019). The different decline rates were proposed to correlate with the amount of hydrogen retained by the progenitor (e.g., Popov 1993; Anderson et al. 2014; Faran et al. 2014; Gutiérrez et al. 2014; Moriya et al. 2016; Hillier & Dessart 2019); however, see Morozova et al. (2017) for an alternative explanation based on the interaction with dense circumstellar material. On the contrary, the near-infrared (NIR) shows that there are spectroscopic distinctions between slow and fast decliners in the 1.06 μm region, which may be due to progenitor differences (Davis et al. 2019).

Preexplosion images of SNe IIP suggest that most of their progenitors are red supergiants (RSGs; e.g., Smartt et al. 2004; Maund et al. 2005; Smartt et al. 2009, 2015; Van Dyk 2017; Van Dyk et al. 2019). However, blue supergiants called SN 1987A–like can produce SNe II with longer rising light curves after the famous SN 1987A (see Arnett et al. 1989 for a review). Theory suggests that if an RSG progenitor has most of its hydrogen envelope stripped preexplosion, it can give rise to an SN IIb; however, the progenitor likely plays a role as it may not be that of a normal SN II (see Bersten et al. 2012). SNe IIb show hydrogen at early times and helium in their later spectra that quickly becomes the dominant feature (Filippenko et al. 1993). The required strong mass loss is believed to be due to a binary progenitor (e.g., Podsiadlowski et al. 1993; Aldering et al. 1994).

There is an open question as to whether or not there is a continuum between fast declining SNe II and SNe IIb. Arcavi et al. (2012) and Pessi et al. (2019) suggest that there is a discrete change between SNe II and SNe IIb. However, Pessi et al. (2019) discussed a peculiar SN II, SN 2013ai, that presented an unusually large rise time and was systematically "misclasified" as an SN IIb when performing clustering analysis. This, along with the spectral peculiarities discussed later in the current work, encouraged further analysis of SN 2013ai in order to elucidate if the object is in fact some kind of transitional event, in which case it can provide important clues about the stellar evolutionary pathways leading to CCSNe.

In this paper, we present optical and NIR observations and analysis of SN 2013ai. Section 2 outlines the observing techniques and reduction procedures. In Section 3, we describe the processes used to determine the reddening of SN 2013ai. In Sections 4 and 5, we present the photometric and spectroscopic properties of SN 2013ai and make comparisons to other CCSNe. The analysis of the preexplosion HST images and the possible progenitor scenarios are described in Section 6. In Section 7, the models and their outcomes are reviewed. The discussion of results and conclusions are in Sections 8 and 9, respectively.

SN 2013ai was discovered at R.A.(2000) = 6^h16^m18^s.35, decl. (2000) with the Zadko 1 m telescope at Gingin Observatory, Australia, on 2013 March 1 UT at an apparent R-band magnitude of 14.4 mag (Conseil et al. 2013). However, the SN was previously detected in prediscovery images as early as 2013 February 26 UT and not detected in images taken on 2013 February 21 UT, with a limiting magnitude of r = 20.7 from the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS; Chambers et al. 2016). Thus, an explosion date of 2013 February 24 UT (56347 MJD) ±2.7 days is adopted for the remainder of this work. The error listed on the explosion date is a flat distribution without any obvious bias toward the last nondetection or discovery.

The host of SN 2013ai, NGC 2207, is in the process of colliding with another galaxy, IC 2163, which makes any distance measurement based on galaxian properties uncertain. The host is a luminous infrared galaxy with a high star formation rate that produces SNe frequently, e.g., SNe 1975A, 1999ec, 2003H, 2018lab, and AT 2019eez (Kirshner et al. 1976; Jha et al. 1999; Graham et al. 2003; Sand et al. 2018; Stanek 2019). At a redshift of 0.009, NGC 2207 is too close for a reliable redshift distance as it is not in the Hubble flow.

NGC 2207 was host to SN 1975A, an SN Ia that provides the only redshift-independent distance measurements to the galaxy. Using the GELATO classification tool (Harutyunyan et al. 2008), the best match of the 1975 January 20 spectrum of SN 1975A published by Kirshner et al. (1976) was found to be SN 2003du, which had a B-band decline rate of Δm₁₅(B) = 1.07 ± 0.06 (Blondin et al. 2012). Pseudo equivalent width (pEW) measurements of the Si ii λ λ5972 and 6355 lines from this spectrum place SN 1975A in the middle of core-normal SNe in the Branch diagram (Branch et al. 2006, 2009). In addition, the pEW of the 5972 line is consistent with Δm₁₅(B) ∼ 1.1 (see Figure 17(b) of Folatelli et al. (2013)). Following the procedures of Riess et al. (1998) and using the reproduced photometry from Kirshner et al. (1976), the SN 1975A light curves were analyzed using SNooPy (Burns et al. 2011) to estimate a distance. Fixing the B-band decline rate in the range Δm₁₅(B) = 1.1–1.5 yielded a 0.1 mag difference in the distance moduli and a fit error of 0.15 mag from SNooPy. This procedure gives a final distance modulus of 33.5 ± 0.2 for NGC 2207, which we adopt in this paper.

SN 2013ai was classified as an SN II by the Public ESO Spectroscopic Survey of Transient Objects (PESSTO; Smartt et al. 2015) due to broad Hα emission in its optical spectrum (Klotz et al. 2013). Archival HST-WFPC2 images of the field are available, obtained using the F336W, F439W, F555W, and F814W filters, which contained the site of SN 2013ai on 1996 May 25 UT. However, no clear single progenitor could be detected (Milisavljevic et al. 2013). These data are reanalyzed in Section 6.

Optical photometric follow-up commenced soon after SN 2013ai was discovered, using ESO + EFOSC2, the 1.3 m SMARTS telescope ²⁹ + ANDICAM, and the 1 m Swope telescope + SITe3. Swope + SITe3 images were taken as a part of the Carnegie Supernova Project II (CSP-II; Phillips et al. 2019; Hsiao et al. 2019). In the NIR, imaging was obtained with NTT + SOFI and the 2.5 m du Pont telescope + RetroCam. A Swope V-band image taken at ∼10 days past explosion is shown in Figure 1.

Zoom In Zoom Out Reset image size

ANDICAM data were reduced with the dedicated pipeline that subtracts the overscan and bias from each image, followed by flat-fielding correction using dome flats. For EFOSC2 data, the reduction was performed using the PESSTO pipeline, as described in Smartt et al. (2015), which trims, debiases, and flat-fields images (using twilight sky flats). In addition, a fringe frame is used to correct the I-band data. The SN magnitude in each frame was measured via point-spread function (PSF) fitting photometry using the SNOoPy package ³⁰ , with errors estimated from artificial star tests. The photometric zero point for each image was determined from aperture photometry of local sequence stars, which in turn were calibrated to Landolt fields (Landolt 1992) observed on multiple photometric nights. Color terms were applied to the EFOCS2 and ANDICAM magnitudes using the values listed in Smartt et al. (2013) for the former and values listed on the SMARTS consortium webpages for the latter. The reduction of Swope photometry was performed as described in Phillips et al. (2019). Swope PSF photometry was performed in the same process described for ANDICAM and EFOSC2 data with standard fields observed on the same nights as the SN observations when possible. The Swope data have been converted to the standard system using the process outlined in Phillips et al. (2019). The light curves are presented without template subtraction as the SN is well separated from the host. A log of the optical photometry is given in Table 1.

The NIR light curves from Dupont + RetroCam were reduced in the standard manner following Phillips et al. (2019). The NIR NTT + SOFI data were corrected for crosstalk, flat-fielded (using dome flats), and scattered light, performed within the PESSTO pipeline. Separate off-target sky frames were taken as part of the sequence of observations, and these were used to subtract the sky background from the on-target images. The multiple exposures taken in each filter were combined to produce a deep image on which photometry was performed. The log of NIR photometry is given in Table 2. The photometry of SN 2013ai previously published by Valenti et al. (2016) is also included. The light curves in each band, during the photospheric phase, are presented in Figure 2.

Zoom In Zoom Out Reset image size

Through the PESSTO collaboration, an optical spectroscopic time series was obtained with EFOSC2 on the ESO-NTT (Buzzoni et al. 1984). Complementary, spectra were taken with the Wide-Field Spectrograph (WiFeS; Dopita et al. 2010) on the Australian National University (ANU) 2.3 m Telescope as part of the ANU WiFeS Supernova Programme (Childress et al. 2016). The EFOSC2 spectra were reduced, extracted, and calibrated with the PESSTO pipeline using standard techniques (Smartt et al. 2015). All spectra were bias-subtracted and divided by a normalized lamp flat. For the Gr#16 EFOSC2 spectra, flat fields were taken immediately after the spectrum, and at the same position on the sky, to allow for the removal of the strong fringing seen at longer wavelengths. Wavelength calibration was performed with respect to a HeAr arc lamp. After extracting each spectrum, the wavelength of the strong night sky emission lines were checked, and if necessary, a linear shift was applied to the dispersion solution. The spectra were flux calibrated using a sensitivity curve constructed from observations of spectrophotometric standard stars, and then telluric absorptions were removed using a scaled synthetic model of atmospheric transmission (Patat et al. 2011). The WiFeS spectra were reduced using the PyWiFeS software (Childress et al. 2014). The late-time spectrum taken with Gran Telescopio Canarias (GTC) on 2013 October 30 was reduced and extracted using standard IRAF (Tody 1986) routines and was flux-calibrated using a spectrophotometric standard star observed on the same night. A log of the spectroscopic observations is presented in Table 3, and the optical spectra during the photospheric phase are presented in Figure 3.

Zoom In Zoom Out Reset image size

An NIR spectrum of SN 2013ai was obtained using the Folded-port Infrared Echellette (FIRE; Simcoe et al. 2013) on the Magellan Baade telescope as part of CSP-II (Phillips et al. 2019; Hsiao et al. 2019). The spectrum was reduced and telluric-corrected following the procedure outlined in Hsiao et al. (2019) and has been previously published in Davis et al. (2019). All observations will be made public via WISeREP (Yaron & Gal-Yam 2012).

Further observations taken with Swift-XRT 25 days after explosion detected an X-ray source from the SN position at a level of ∼2 × 10⁻¹⁴ erg s⁻¹ cm⁻² (Margutti et al. 2013). Radio observations were taken using the Combined Array for Research in Millimeter-Wave Astronomy (CARMA) 10 and 11 days past explosion at 85 GHz; however, no radio source was detected at either epoch at a limit of 0.6 mJy (Zauderer et al. 2013).

Milky Way extinction toward SN 2013ai is taken from Schlafly & Finkbeiner (2011). In order to estimate the extinction of SN 2013ai, the equivalent width (EW) of the Na i absorption lines at the redshift of the host galaxy NGC 2207 were measured in all available WiFeS spectra as they have the highest resolution of R ∼ 7000. Assuming Gaussian profiles for both components, we obtain an average of D₁ = 0.68 ± 0.15 Å, D₂ = 0.91 ± 0.07 Å, and D₁ + D₂ = 1.52 ± 0.07 Å. Using Equation (9) from Poznanski et al. (2012), we obtain an E(B − V) = 0.84 ± 0.64 mag. However, according to Phillips et al. (2013), the error associated with this method can be much larger than predicted by relations derived for typical Milky Way dust and gas. As a check for consistency, we used the Balmer decrement to determine the dust extinction by computing the flux ratio of Hα and Hβ lines using the same WiFes spectra. The Balmer decrement was measured using the narrow host galaxy lines in the spectrum. The obtained color excess was E(B − V) = 0.70 ± 0.34 mag. Similar results were obtained using PESSTO low-resolution spectra from which we extracted a region near the SN without subtracting the background. The Balmer decrement was measured from a flux-calibrated Milky Way extinction-corrected spectrum. Errors were determined by a series of Monte Carlo realizations, assuming the uncertainty is dominated by the photometric errors and therefore flux calibration. The Balmer decrement was measured for each realization, with the error taken as the standard deviation. Finally, the diffuse interstellar band (DIB) at 5780 Å was analyzed. This band appears to be very weak or nonexistent in the available spectra. Thus, for this work, we adopt a color excess estimate of E(B − V)_host ∼ 0.8 ± 0.5 mag from the Balmer decrement and Na I D measurements.

The optical light curves of SN 2013ai have been previously studied by Valenti et al. (2016) and Pessi et al. (2019). The V-band light curve was characterized by Davis et al. (2019) using the parameters defined by Anderson et al. (2014). Here we combine the public data of SN 2013ai along with a previously unpublished data set to further analyze SN 2013ai.

The new data are in agreement with previous results that positioned SN 2013ai among typical fast declining SNe II with a decline rate per 100 days, s₂ (Anderson et al. 2014), of ∼2.0 (historically SNe IIL; see Barbon et al. 1979), and an absolute V-band maximum of −18.7 ± 0.06 mag. The only photometric difference SN 2013ai displays with respect to normal SNe II is that it presents an atypically long rise time (e.g., Valenti et al. 2016). From the Pessi et al. (2019) sample, the typical rise times of SNe II are 8.3 ± 2.0 days, 12.8 ± 2.4 days, and 16.0 ± 3.6 days in the B, V, and r bands, respectively. For SN 2013ai, the rise times are 13.0 ± 5.3 days, 18.9 ± 3.7 days, and 24.4 ± 4.0 days, in the B, V, and R bands, respectively. On average, the rise times of SN 2013ai are 2.4σ away from the mean values calculated in Pessi et al. (2019). Note that in Pessi et al. (2019) the computed rise time average for the SN II sample is longer than for other higher cadence samples (e.g., González-Gaitán et al. 2015; Gall et al. 2015; Valenti et al. 2016; Rubin & Gal-Yam 2016). Rubin & Gal-Yam (2016) studied the R-band light curves of 44 SNe II and found an average rise time of 7.8 ± 2.8 days, with the slowest rising object, PTF12hsx, taking 16.2 days to reach maximum. SN 2013ai rises significantly more slowly than any SN in the Rubin & Gal-Yam (2016) sample. The rise time of SN 2013ai is significantly longer than that of a typical SNe II.

Late-time photometric data may be used to estimate the ⁵⁶Ni mass using the procedure of Hamuy (2003) and the bolometric correction from Bersten & Hamuy (2009). However, the scarce late-time data do not allow for a decay slope to be measured, so only a lower limit of the ⁵⁶Ni mass can be given as the extent of the γ-ray trapping can only be estimated. For 100% trapping, we get M_⊙ with the error bars including the photometry error but not the reddening and distance uncertainties. There are large uncertainties associated with this lower limit as the light-curve decay is highly dependent on the amount of γ-ray trapping. For the case of SN 1993J (see Hoflich et al. 1993), which is a good approximation for SN 2013ai at late times (see Section 7), 79% of γ-rays had escaped at 350 days past explosion. However, ⁵⁶Co decays in two channels, electron capture and positron capture. Approximately 4% of the decays are positron capture and will remain trapped, giving a total trapping of ∼25%. This suggests that the ⁵⁶Ni mass lower limit from observations could be underestimated by up to a factor of 4, translating to a ⁵⁶Ni mass of 0.4–0.8 M_⊙. Similarly, assuming that the light curve is powered by ⁵⁶Co decay past 50 days from explosion, i.e., constant slope, the magnitude difference of pure ⁵⁶Co decay from the observed light curve gives a similar result.

SN 2013ai may be compared with two other SNe II that were found in the literature: ASASSN-14 kg (Valenti et al. 2016) and SN 2017it (Afsariardchi et al. 2019). ASASSN-14 kg was chosen as it was found to have the most similar light curve to SN 2013ai. SN 2017it was published as a high-⁵⁶Ni-mass SN II with a long rise, so despite its slow decline postmaximum, bears a comparison to SN 2013ai. Photometry from Valenti et al. (2016) shows a V-band rise time for ASASSN-14 kg of ∼18 days, while Afsariardchi et al. (2019) report a V-band rise time of ∼20 days for SN 2017it. In order to recalculate the rise times of these objects in a uniform fashion, the explosion epoch of the comparison SNe is taken as the midpoint between the last nondetection and the first detection while considering the error to be half the difference between said epochs. The light curves were then interpolated via a Gaussian process method using the Python library GPy. ³¹ Finally, the maximum date was obtained from the interpolated light curve, and the rise times were calculated. Results are presented in Table 4. It can be seen in Figure 4 that the behavior of SN 2013ai around the peak is similar to ASASSN-14 kg. SN 2017it, despite the long rise, exhibits a postmaximum decline similar to a normal SN II. For further comparison to the SN IIb subclass, SN 1993J is added as an example of a typical SN IIb light curve. The rise time of SN 2013ai is similar to an SN IIb, however, the postmaximum decline is much slower.

Zoom In Zoom Out Reset image size

To calculate the bolometric light curve, the Milky Way and NGC 2207 dust-reddening-corrected photometry was first converted into AB magnitudes using the spectral energy distribution of Vega. Monochromatic fluxes for each available date were numerically integrated in wavelength, obtaining the pseudo-bolometric flux (F_pbol). For the UV flux (F_UV), a linear extrapolation was used and integrated from the observed flux at the B band to zero flux at 2000 Å(see Bersten & Hamuy 2009; Folatelli et al. 2014 for further details). To estimate the IR flux (F_IR), the spectral energy distribution was fit by a grid of temperatures and angular sizes that define a blackbody for each point of the grid. Using the resultant χ² grid from each blackbody, a probability distribution to determine the best parameters for each epoch is calculated. Afterward, the flux was integrated from the effective wavelength of the I filter to 100,000 Å, where the flux is considered negligible.

The bolometric flux was calculated as F_bol = F_pbol + F_UV + F_IR. This bolometric light curve was compared to a subsample of the bolometric light curves presented by Faran et al. (2018). Their sample includes only objects with high-quality multiband data and rather good time coverage, which makes it useful for our goal of comparing SN 2013ai to a large sample of bolometric light curves computed in a uniform manner. Note that the bolometric light curves of Faran et al. (2018) were computed differently from the SN 2013ai bolometric light curves as they consider the effects in the bluer and redder wavelengths, e.g., line blanketing. However, given the uncertainties involved in calculating bolometric light curves, the comparisons presented would not significantly change if the bolometrics were recalculated. The subsample was selected such that only objects that show a clear maximum in the light curve, meaning there are data points before maximum, are considered. From the 29 objects presented in Faran et al. (2018), only 8 meet this requirement, although SN 2003hf is much more luminous than the others so it was excluded, giving 7 comparison objects. The bolometric light curve of SN 1993J (Richmond et al. 1994) is also presented for comparison. Figure 5 shows that SN 2013ai is the fastest decliner, except for the SN IIb SN 1993J. Using the explosion dates published by Faran et al. (2018), we obtain a mean bolometric rise time of the subsample of 8.6 ± 2.0 days while the bolometric rise time of SN 2013ai is 14.3 ± 4.1 days.

Zoom In Zoom Out Reset image size

The optical and NIR spectra of SN 2013ai were analyzed quantitatively by measuring expansion velocities and pEWs for every feature present. These techniques have been implemented successfully before to study large samples of SNe II and their diversity (e.g., Gutiérrez et al. 2017; Davis et al. 2019; de Jaeger et al. 2019).

5.1. Time Evolution

Figure 3 shows the optical spectra of SN 2013ai. In the optical, at all epochs, we see narrow optical emission in Hα, Hβ, [O ii] λ3727, [O iii] λ5007, and [S ii] λ λ 6717, 6731 at the redshift of NGC 2207. The earliest spectra, taken 11 days before V maximum, show only a broad Hα P Cygni profile with narrow Na i D λ λ5890 and 5896 likely from the interstellar medium (ISM). Hβ is not present until later times, first seen around 3.2 days before V maximum. A montage of the optical spectra is presented in Figure 3.

Starting with the spectrum taken 3 days after V maximum, the Hα absorption becomes split, likely due to Si ii λ6355 at these early times (Gutiérrez et al. 2017). Fe-group lines are also present during this epoch, primarily Fe ii λ λ 5169, 5267, and 5363. We also see the emergence of the Ca ii NIR triplet as a broad and blended P Cygni line around 8500 Å.

From 11 to 20 days after maximum, SN 2013ai became even redder, as more flux is lost to absorption by Cr, Sc, Ba, and Fe in the blue. The usual Fe-group lines seen in SNe II are present between 3500–5500 Å; however, these lines are weaker than what are seen in a normal SN II (see Gutiérrez et al. 2017; de Jaeger et al. 2019). Hα has a strong absorption component in its P Cygni profile, with a minimum at a velocity of 11,000 km s⁻¹, while the second absorption, which was tentatively associated with Si ii, is still present as a "notch" in the feature. There is a lack of absorption from Sc ii and Ba ii on the blue side of the Hα absorption line.

By 30 days after maximum, the notch in Hα is no longer present in the spectrum of SN 2013ai, suggesting that the notch was not due to high velocity (HV) H i (see Gutiérrez et al. 2017 for more information on this feature). The spectrum shows relatively little evolution over the following period, from 30 to 73 days postexplosion, the continuum continues to become redder, while the most notable change in the lines is an increase in the strength of the Na i D absorption.

The late-time optical spectrum obtained 230 days after V maximum with GTC+OSIRIS is shown in Figure 6. This spectrum was taken 250 days from the estimated explosion epoch and is shown in comparison to the SN IIP SNe 1999em and the SN IIb SN 1993J at similar times. The spectrum of SN 2013ai contains very little flux from the SN. Nonetheless, several lines characteristic of CCSNe in their nebular phase are detected, namely broad Hα emission, [Ca ii] λ7291, and possibly [O i] λ6300. The reddening-corrected [Ca ii] to [O i] line ratio appears approximately similar to that in SN 2004et (see Jerkstrand et al. 2012), suggesting that the core mass of the star that exploded was probably quite similar and not exceptionally massive. However, the weakness of these lines makes it difficult to accurately measure a flux ratio and suggests that the SN is not sufficiently nebular for an accurate progenitor to be assumed using the method of Jerkstrand et al. (2012). The Hα profile at this epoch is dominated by emission and shows a profile with a wide base and narrow top with no absorption. Together with the early X-ray observations, as noted in Section 2, this suggests interaction with circumstellar material (CSM).

Zoom In Zoom Out Reset image size

The NIR spectrum of SN 2013ai is shown in Figure 7 and exhibits mostly hydrogen and helium features, such as P_δ and B_γ and the He i/P_γ blend. In the NIR, SN 2013ai does not exhibit any lack of features like the absence of Sc ii and Ba ii in the optical.

Zoom In Zoom Out Reset image size

5.2. Comparison to Other Core-collapse SNe

Figure 8 shows optical spectra of SN 2013ai compared to other SNe at similar phases. Included are examples of well-studied SNe IIP, IIL, and IIb. This comparison highlights the uniqueness of SN 2013ai in the features seen, their strengths, and expansion velocities. As previously noted, there is a lack of features on the blue side of Hα and the Na/He absorption around 5990 Å is significantly shallower than that seen in a normal SN II/IIb.

Zoom In Zoom Out Reset image size

Figure 7 shows the 14.5 day NIR spectrum of SN 2013ai compared to other SNe II at similar epochs. The SN does not show any signs of being an NIR-weak SN II (Davis et al. 2019); however, at such early times, it is not always possible to determine an NIR spectroscopic subclass. The He i 1.083 μm line is stronger than that of SN 2013hj at a similar phase. Both SNe show a boxy P_β emission profile. Overall, the NIR features of SN 2013ai are much more normal than in the optical. However, similar to the optical, the features in the NIR have much higher velocities than other SNe IIP/L at similar times.

The top row of Figure 9 compares the pEW of all optical features present in the spectra of SN 2013ai with the mean from Gutiérrez et al. (2017) and all data from the de Jaeger et al. (2019) and Liu et al. (2016) samples. The metal lines, from 4500-5500 Å, of SN 2013ai evolve like those of a typical SN II. The 5990 Å He/Na blend is particularly interesting due to its weakness at later times. The pEWs of all features seen in the optical lie among normal SNe II values. Similarly, the top row of Figure 10 compares the pEW of all NIR spectral features of SN 2013ai with the sample from Davis et al. (2019). The NIR pEWs, like the optical, lie within the large comparative samples.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The bottom row of Figure 9 shows the absorption velocities of SN 2013ai compared to the samples of Liu et al. (2016), Gutiérrez et al. (2017), and de Jaeger et al. (2019). SN 2013ai exhibits high velocities in all features observed, for an SN II. The 5990 Å He/Na blend is particularly interesting due to its high velocities if assumed to be Na i λ5993, giving velocities over 13,000 km s⁻¹ at early times. While the Fe ii velocities are closer to normal SN II values, they are still noticeably higher than the comparative SN II samples, with the velocities more similar to SESNe. Despite the high velocities, the evolution is decreasing with time, like that of a normal SN II. Figure 10 shows the NIR absorption velocities of SN 2013ai compared to that of the data from Davis et al. (2019). Much like the optical, the velocities are significantly higher than the sample, ∼2σ, with the Pβ and Pδ features around 12,000–13,000 km s⁻¹ and He i λ1.083 μm around 11,500 km s⁻¹ at 22 days past maximum. The offset of velocities from typical SN II values is similar in the optical and NIR. The normal pEWs and high velocities of SN 2013ai suggest that the features are not broad, but are fast, like those of an SN IIb. However, the spectral features present are typical of an SN II.

The optical absorption velocities at 50 days past explosion, the middle of the plateau for a normal SN II, were also compared to the Gutiérrez et al. (2017) and de Jaeger et al. (2019) samples in order to see if any other SNe II have been found with velocities as high as SN 2013ai. No SNe with velocities as high as SN 2013ai were found, except for SN 2007ld from the Gutiérrez et al. (2017) sample, which has a 50 day Hα velocity over 11,000 km s⁻¹. SN 2007ld has a normal photometric and spectroscopic evolution; however, it shows some contamination in the Hα emission profile, which could be due to He i, a common sign of an SN IIb. Unfortunately, SN 2007ld has no spectra toward the end of its plateau, when an SN IIb can be easily identified. SN 2007ld does not exhibit a long rising light curve or a quick decline after maximum. Superluminous SNe (e.g., SN 2013fc, SN 2016gsd; Kangas et al. 2016; Reynolds et al. 2020) are not included for comparison as SN 2013ai is not superluminous. However, see SN 2016gsd (Reynolds et al. 2020) for a superluminous SN II that exhibits high velocities similar to SN 2013ai.

The spectra of SNe with similar rise times to SN 2013ai (see Section 4) were also examined. The classification spectrum of ASASSN-14 kg is plotted with the earliest spectrum of SN 2013ai in Figure 11. The spectra look similar at these early times. There are no later spectra taken of ASASSN-14 kg. ASASSN-14 kg has a well-defined explosion date, 59672.9 ± 0.5 MJD (Nicolas et al. 2014). Given the similarity of the light curves and early spectra of ASASSN-14g and SN 2013ai, it is likely that SN 2013ai is around the same phase from the explosion in the spectra seen in Figure 11. Afsariardchi et al. (2019) presented the Hα velocity of SN 2017it at 93 days past explosion, which is higher than that of a normal SN II. However, it is still ∼1000 km s⁻¹ slower than SN 2013ai. We were unable to find any SN II with expansion velocities comparable to SN 2013ai during the photospheric phase. However, the similarity of the early-time spectra and light curves of ASASSN-14 kg and SN 2013ai suggests that they could be alike.

Zoom In Zoom Out Reset image size

6.1. Preexplosion Observations

Preexplosion HST observations were available for SN 2013ai, consisting of HST+Wide-Field and Planetary Camera 2 (WFPC2) images taken on 1996 May 25; see Figure 12. The SN location lies on the WF2 chip, which has a pixel scale of 01 pixel⁻¹. Four exposures, consisting of two cr-split pairs, were taken in each of the F336W, F439W, F555W, and F814W filters, with total exposure times of 2000 s, 2000 s, 660 s, and 720 s, respectively. The pipeline reduced _c0f files were downloaded from the Mikulski Archive for Space Telescopes (MAST). ³² The images were first masked with their associated _c1f files, and the cr-split pairs were combined to reject cosmic rays using the crmask task within IRAF. Finally, the dithered images in each filter were aligned and coadded.

Zoom In Zoom Out Reset image size

Unfortunately, the SN position lies very close to the charge transfer trap 2-337 on the WF2 chip in all images (Whitmore & Wiggs 1995), which gives rise to a bad column artifact. While it is possible to attempt to reconstruct the flux lost due to the trap, as dithered exposures were available, we instead ensured that the bad column was masked in both cr-split pairs for each filter before coadding. As the cr-split pairs in each filter were offset by ∼3 pixels in the x direction, the defect (when masked) does not contribute to the final image when they are shifted and combined.

We also attempted to use the drizzle algorithm (Fruchter & Hook 1997) as implemented in astrodrizzle within the drizzlepac package to improve the spatial resolution of the F814W image. Before drizzling, the bad column close to the position of SN 2013ai was flagged in the data quality image. A series of drizzled images was produced to test the effect of varying the output pixel scale and the size of the drop. However, in all cases, a faint artifact was still visible at the position of SN 2013ai in the final image. It is likely that the artifact in the output image results from the interpolation over the masked bad column. Thus, the drizzled images were not considered any further.

To identify the position of SN 2013ai in our preexplosion images, we obtained Ks filter imaging of SN 2013ai with the Nasmyth Adaptive Optics System and Near-Infrared Imager and Spectrograph (NAOS-CONICA; NaCo) on the Very Large Telescope UT4 over 2013 March 3-4. SN 2013ai itself was used as a guide star for NAOS, while the S54 camera was employed on CONICA, yielding a pixel scale of 0.0543''/pixel, over a field of view of 56''×56''. Daytime calibration data consisting of flat fields and long and short exposure darks were reduced with the NaCo pipeline (version 4.3.1) to give a master flat field and a map of aberrant pixels on the detector. Using these, the individual science frames were masked and flat-fielded. All images of SN 2013ai were dithered on source, facilitating the subtraction of the sky background using the IRAF xdimsum package. After sky subtraction, the individual frames were aligned and combined to yield a single, deep image with a total exposure time of 5040 s.

The NaCo postexplosion image was aligned to the preexplosion WFPC2 F814W image using a geometric transformation derived from the pixel coordinates of sources common to both images. Two separate sets of sources were used for the alignment: a "good" set of 21 sources which were clearly point-like, detected with good S/N, and which appear to lie within NGC 2207, and a larger superset of 31 sources which we term the "complete" set, where the signal-to-noise ratio (S/N) of sources were lower, and which also included foreground and slightly extended objects. For both the "good" and "complete" sets, the pixel coordinates of each source were measured in the WFPC2 and NaCo frames, and three separate geometric transformations were then derived using the lists of matched coordinates, one allowing for translation, rotation, and scaling in x and y, and two that also included second- and third-order polynomial terms, respectively.

The pixel coordinates of SN 2013ai were measured in the NaCo image using three different algorithms (centroid, Gauss, and filter within IRAF), which had a standard deviation of only 0.05 pixels, or 3 mas. The average of the three algorithms was taken as the position of SN 2013ai in the NaCo image. This position was then transformed to the pixel coordinates of the preexplosion image using each of the transformations derived using both the "good" and "complete" sources. The standard deviation of the transformed positions was 16 mas, and we take this to be indicative of the uncertainty in SN position depending on the geometry or sources used for the alignment. We add this in quadrature with the uncertainty in the SN position in the NaCo observation (3 mas), and the average of the rms error found when fitting each transformation (58 mas), to obtain the total uncertainty of 60 mas on the location of SN 2013ai in the WFPC2 image.

6.2. Progenitor Age Estimate

The SN falls close to, but is not coincident, with a source in the WFPC2 image, which we designate "Source A". The center of Source A (see Figure 12), as measured using three centering algorithms within IRAF phot, lies 193 mas (1.93 pixels) from the transformed SN position. This offset is a factor of ∼3 greater than the uncertainty in the transformed position of SN 2013ai (60 mas), and so we discount the possibility that this is a single stellar progenitor of SN 2013ai.

By eye, Source A appears to be somewhat extended. To test this, 14 isolated, point-like sources were fit with Moffat profiles to determine their FWHM. The average FWHM of these sources was 1.6 ± 0.2 pixels. Source A is significantly broader, with an FWHM of 2.5 pixels. On these grounds, it appears that the source may be an unresolved cluster or complex (2.5 WF pixels corresponds to a physical scale of 46 pc at the distance of NGC 2207).

While the position of SN 2013ai lies outside the FWHM of Source A, it still appears to fall within the wings of the PSF. To quantify the degree to which SN 2013ai lies within the position of Source A, the latter was fit with a Moffat profile. Then, we integrated under this profile to determine the flux within an aperture with a radius corresponding to the distance from Source A to SN 2013ai. Assuming that SN 2013ai is associated with a cluster, it was found that % of the flux of Source A is in the position of SN 2013ai. In the remainder of this section, we consider the implications for the progenitor of SN 2013ai both if it was associated with Source A and if it was unrelated.

Photometry was performed on Source A using hstphot, (Dolphin 2000a) a stand-alone photometry package for use with WFPC2 data. hstphot includes corrections for charge transfer efficiency (CTE) losses, aperture corrections, and zero points appropriate for WFPC2 (Dolphin 2000b). While Source A is somewhat broader than most point sources in the field, it is still relatively well fit by a PSF. We measure PSF-fitted magnitudes for Source A in the HST flight system (Holtzman et al. 1995), with updated zero points per Andrew Dolphin's webpages ³³ of F336W = 22.68 ± 0.12 mag, F439W = 23.90 ± 0.15 mag, F555W = 23.40 ± 0.12 mag, and F814W = 22.60 ± 0.13 mag. Transforming to the standard UBVRI system, these correspond to U = 22.68 ± 0.12 mag, B = 23.89 ± 0.15 mag, V = 23.38 ± 0.12 mag, and I = 22.56 ± 0.13 mag. As a test of the PSF-fitting results, aperture photometry was also performed using a small (2 WF pixel) aperture. For all four filters, the difference between PSF fitting and aperture photometry was <0.1 mag, i.e., less than the photometric error.

To estimate an age for Source A, the chorizos SED-fitting package was used (Maíz-Apellániz 2004). A grid of Starburst99 models was fit (Leitherer et al. 1999) to the measured WFPC2 photometry of the source, constraining the reddening-law parameter R₅₄₉₅ to be 3.1 and the metallicity to the solar value, while allowing the cluster age and extinction to vary. The results of the fitting procedure are shown in Figure 13. The best-fitting model (shown in the left panel) has an age of 7.9 Myr and a relatively low extinction of E(4405-5495) = 0.12 mag. Comparing to the STARS models ³⁴ (Eggleton et al. 2011), such a cluster age would imply a progenitor mass of ∼21 M_⊙; however, we note that there are also reasonable fits for 11.2 Myr, giving a lower progenitor mass of ∼17 M_⊙ and higher extinction. The low-mass best fit of ∼17 M_⊙ is at the high end of SN II progenitor masses (Smartt et al. 2009; Sukhbold & Adams 2020).

Zoom In Zoom Out Reset image size

The extinction toward the cluster inferred by the fit is considerably lower than that seen toward SN 2013ai itself. There is a tail of solutions in the χ² plane that accommodates a reddening of E(4405-5495) = 0.7 mag, which is consistent with that seen toward the SN. However, these solutions require the cluster to be younger than ∼3 Myr, which would in turn lead to a progenitor that is extremely massive and probably inconsistent with the maximum mass of a star that can maintain its H envelope to the point of core collapse.

Further analysis was done using the Bayesian fitting technique of Maund (2018), which addresses cluster masses, distances, and filter transformations. Both the WFPC2 photometry and WFPC2 photometry converted to Johnson-Cousins photometry were included as separate runs. Given the luminosity of the cluster, the mass is constrained to lie in the range . At such masses, individual, luminous stars may have a disproportionate influence on the color of the cluster that might cause deviations from the expected color sequence predicted by STARBURST99 (for much more massive clusters of 10⁶ M_⊙).

The progenitor of SN 2013ai could be a single, isolated red supergiant as found for other SNe II. The relatively bright cluster nearby makes the identification of such a progenitor difficult, as the flux in the wings of the cluster may mask that of a fainter progenitor close by. The most convincing demonstration of such a progenitor would rely on image subtraction, where a flux deficit in new HST images taken after SN 2013ai has faded would indicate the absence of progenitor flux which was present in preexplosion images. As demonstrated by Maund et al. (2014) and Folatelli et al. (2015), difference imaging using late-time data can substantially improve the precision and sensitivity of limits placed on SN progenitors. While such a test will have to wait for new observational data, it is still possible to estimate the sensitivity of preexplosion data to a single red supergiant progenitor. From the photometry of point sources close to the SN position, we estimate a 5σ limiting magnitude of F814W ≳ 24.5 mag for the WFPC2 images. At the distance of NGC 2207, this would imply an absolute magnitude for the progenitor of F814W > −9 mag, which unfortunately only allows us to constrain the luminosity of the progenitor to log L/L_⊙ < 5.5. For single-star models, this luminosity corresponds to a progenitor zero-age main-sequence mass of <24 M_⊙(Eldridge & Tout 2004).

7.1. Modeling Background

7.2. Modeling Results

For the progenitor evolution of SN 2013ai, we use the stellar evolution code MESA (Paxton et al. 2011). To trigger the explosion, the explosion energy (E_exp) is put in the inner region as a thermal bomb, which causes the ejection of the envelope, leaving behind a neutron star or a black hole. For computational efficiency and to cover a sufficiently wide range of parameters, we first used the code developed by Bersten et al. (2011) to model SN 2013ai. To verify and further constrain these models, we then modeled SN 2013ai with the non-LTE hydrodynamical-radiation code HYDRA (Höflich 2003, 2009). A wide range of M_MS, E_exp, and mixing, including pure C/O cores, was evaluated. We will omit the "failed" attempts.

The comparison between a frequency-independent model from the code of Bersten et al. (2011) with the photometry is shown in Figure 14 using depth-dependent mean opacities based on SN 1993J and originating from a 20 M_⊙ star. The resulting models have 0.55 M_⊙ of ⁵⁶Ni. Such a high ⁵⁶Ni mass is problematic for the explosion of a massive star both from observations and theory (see Thielemann et al. 2018 for a review). The bolometric light curve, in principle, is a sensitive measure of the physical parameters (e.g., Suntzeff et al. 1992; Bersten et al. 2011). Despite the advantage of the bolometric light curve, three problem zones may be identified: (1) the reconstruction of the bolometric light curve, (2) distance uncertainties, and (3) asymmetries. Moreover, using the bolometric light curve reduces the information for finding model parameters. Due to these uncertainties, we use monochromatic light curves based on detailed non-LTE models.

Zoom In Zoom Out Reset image size

In Figures 15 and 16, we show the results of HYDRA starting with the progenitor of SN 1993J but modifying the explosion energy, the mass of the H-rich layers, the explosion energy, and the ⁵⁶Ni mass. The H-rich layer mass was adjusted to 0.2 M_⊙, and the ⁵⁶Ni mass to 0.3 M_⊙. These models are based on the results from the previously mentioned models, hereafter referred to as "LC models". The results between best-fit models of our simulations, presented as (LC, HYDRA), give consistent results within the uncertainties: M_MS = (20, 30) M_⊙, foe ³⁵ , and M_Ni = (0.55, 0.3) M_⊙.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Given the model best fits, the effect of uncertainties on the M_Ni, E_exp, C/O core mass and distance modulus must be discussed. Linear polarizations give evidence for a symmetry axis and bipolar explosions with P ≈ 1% indicating axis ratios; a measure of how nonspherical an object is of 1/2–1/4 (Wang et al. 1998; Maund et al. 2007a, 2007b; Stevance et al. 2016; Reilly et al. 2017). To estimate the uncertainty for the light curves, we use SN 1998bw/GRB980425 as a proxy for the explosion of a massive stripped-envelope SN (Höflich et al. 1999). Detailed analysis of polarization spectra resulted in an axis ratio of 0.7 for SN 1998bw seen from an angle of 30^o (Höflich 1995) but with significant uncertainties as discussed in Stevance et al. (2020). The possible effect of asymmetry on the V-band light curve of SN 2013ai is shown in Figure 17. Asymmetry may change the shape of the light curve by varying the maximum brightness up to ±1 mag, placing SN 2013ai near superluminous SNe (see Figure 17), and can shift the rise time by about 10 days.

Zoom In Zoom Out Reset image size

Even for the same model, the light curve rises more slowly and shows a lower maximum brightness that directly translates into uncertainties in E_exp and the C/O core mass. The error for the asymmetry and the shape may be regarded as an upper limit because even for purely jet-driven explosions of CCSNe, an extended outer envelope tends to become more spherical with time (Khokhlov et al. 1999; Höflich et al. 2004; Couch et al. 2009). The colors are found to be less sensitive as they are a measure of the physical conditions at the photosphere, namely, the flux density rather than the total flux. In either case, the tail of the light curve is less affected by polarization because the luminosity becomes isotropic within ∼0.2 mag 1–2 months after the explosion, which translates directly into an error for M_Ni. For the ⁵⁶Ni mass, both 0.3 and 0.55 M_⊙ fit the data; however, as previously mentioned, 0.3–0.55 M_⊙ of ⁵⁶Ni is higher than expected for a massive star explosion.

Though beyond the physics included in our simulations, note that the theoretical V-band light curve at 250 days is fainter by ∼1 mag than observed. This may indicate a significant contribution by the interaction between the SN ejecta and the ISM or a preexplosion stellar wind. This theory would also be consistent with the multitiered Hα emission profile seen in the nebular spectrum of SN 2013ai, which indicates interaction.

Given the similar early-time light-curve timing, expansion velocities, and density structure to SESNe, we suggest that the core of SN 2013ai is closer in physical properties to an SESNe than a normal SN II. This would give SN 2013ai a dense ∼10 M_⊙ C/O core and little hydrogen (0.2 M_⊙ from models) in its outer layers, consistent with the quick decline postmaximum. In addition, the explosion energies of SNe 2013ai and 1993J are comparable, 2.5 foe and 2.0 foe, respectively. In this scenario, the progenitor of SN 2013ai would have most of its hydrogen envelope stripped before explosion, much like an SN IIb.

This brings us to the differences between SN 2013ai and SESNe. The first difference is seen in the shape of the early time, ∼2 weeks premaximum, light curve. SN 1993J, for example, shows an initial maximum and subsequent decline in the light curve lasting about one week (see Figure 4) that is powered by the energy deposited by the shock in the very outer layers. This structure of the theoretical SN 2013ai light curve is very similar to that observed in SN 1993J but this phase has not been observed in SN 2013ai. In the SN 2013ai HYDRA models, the early peak is much less pronounced because of the lower H-shell mass and the longer diffusion times with limited mixing, as predicted by Bersten et al. (2012). Second, the ⁵⁶Ni mass is significantly larger in SN 2013ai than in SN 1993J, which could affect the lack of strong He i seen in the optical. However, in SN 2013ai, the much stronger He i line in the NIR λ1.083 μm is present and shows the presence of a He layer. To boost the optical He i lines in SN 1993J, models with ⁵⁶Ni mixed in the outer layers of the ejecta were needed. This is also consistent with the rapid drop in velocity seen at about 1 week after maximum because the continuum opacity is low. Finally, SN 2013ai shows a much slower postmaximum decline rate over 30 days (∼0.5 mag) compared to the rapid decline in SN 1993J, (∼2 mag). The difference can be attributed to strong ⁵⁶Ni mixing in SN 1993J that results in a steeply increasing γ-ray escape, whereas SN 2013ai has almost full trapping. While the models based on Bersten et al. (2011) show general agreement with the observations, the strong mixing needed is not consistent with little He i seen in the optical.

From observations, SN 2013ai seems to be rather unusual but, overall, it fits well within the picture of an energetic explosion of a massive star that lost most (but not all) of its hydrogen envelope, which puts it at the edge of becoming a SN Ib or SN Ic, physically. However, using the observational classification scheme of SNe, SN 2013ai lies between an SN II and an SN IIb as its spectra are those of an SN II and its light curve more similar to SESNe. Explosions of massive stars may appear either as SNe II, IIb, Ib, or Ic depending on the details of the mass loss and the amount of mixing. The open question remains whether the mass loss in SN 2013ai is line driven or the result of a binary evolution. However, given its similarities to an SN IIb, we cannot rule out that mass loss preexplosion occurred in a binary system for SN 2013ai. As such, the current classification scheme may mask diversity and similarity as SN 2013ai has the spectra of an SN II, but after analysis of the spectra and light curves, SN 2013ai lies between SNe II and SNe IIb.

SN 2013ai is a rare SN that exhibits a long light-curve rise, sustained high expansion velocities, and a light-curve decline similar to an SN IIL. Only one SN II was found to have comparable light curves and spectra, ASASSN-14 kg. However, the lack of ASASSN-14 kg data prevents a more detailed comparison between the two SNe. Some similarities were seen in the light-curve shape and expansion velocities when comparing SN 2013ai with SNe IIb. The data presented provide a link between different classes of SNe. In the current classification scheme, SN 2013ai lies between SNe II and SESNe with possible signs of interaction seen in the late-time Hα profile and early X-ray observations. The mass loss during the stellar evolution may be a decisive factor. However, similar data sets for a large number of SNe II are needed to better describe details of the explosion physics and to search for the diversity and address the question of the origin of this mass loss. We demonstrated the importance of the combination of optical and NIR spectra with light curves but, in lack of spectropolarimetry, we had to rely on similarity arguments and were unable to constrain some errors present in the early light-curve data due to asymmetry.

We thank the staff at the different observatories for performing the observations. Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Based on observations made with the 6.5 meter Magellan Telescopes at Las Campanas Observatory, Chile. Based on observations made with the Gran Telescopio Canarias (GTC), installed at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias, in the island of La Palma. Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere, Chile as part of PESSTO, (the Public ESO Spectroscopic Survey for Transient Objects Survey) ESO program ID 188.D-3003. Based on observations made with the Very Large Telescope (VLT) and Nasmyth Adaptive Optics System (NAOS) Near-Infrared Imager and Spectrograph (CONICA) program ID 090.D-0329.

The work of the CSP-II has been generously supported by the National Science Foundation under grants AST-1008343, AST-1613426, AST-1613455, and AST1613472. The CSP-II was also supported in part by the Danish Agency for Science and Technology and Innovation through a Sapere Aude Level 2 grant. M.F. is supported by a Royal Society—Science Foundation Ireland University Research Fellowship. P.H. acknowledges support by grants of the NSF AST-1008962 and NASA's ATP-1909476. L.G. was funded by the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 839090. This work has been partially supported by the Spanish grant PGC2018-095317-B-C21 within the European Funds for Regional Development (FEDER). T.W.C. acknowledges the EU Funding under Marie Skłodowska-Curie grant agreement No 842471. M.S. is supported by generous grants from Villum FONDEN (13261, 28021) and by a project grant (8021-00170B) awarded by the Independent Research Fund Denmark.

Facilities: Magellan:Baade - Magellan I Walter Baade Telescope, Swope - , MAST - , GTC - , VLT - , HST - , TAROT - , NTT - , CTIO:1.3m - .

Software: PyWiFeS (Childress et al. 2014), IRAF (Tody 1986), SNooPy (Burns et al. 2011), Astrodrizzle (Hack et al. 2012), chorizos (Maíz-Apellániz 2004)), MESA (Maíz-Apellániz 2004), HYDRA (Höflich 2003, 2009), STARS (Eggleton et al. 2011).