# The first Hubble diagram and cosmological constraints using superluminous supernovae

C. Inserra, M. Sullivan, C.R. Angus, E. Macaulay, R.C. Nichol, M. Smith, C.

Frohmaier, C.P. Gutiérrez, A. Möller, D. Brout, et al.

► To cite this version:

C. Inserra, M. Sullivan, C.R. Angus, E. Macaulay, R.C. Nichol, et al. The first Hubble diagram and cosmological constraints using superluminous supernovae. Mon.Not.Roy.Astron.Soc., 2021, 504 (2), pp.2535-2549. 10.1093/mnras/stab978. hal-02572608

HAL Id: hal-02572608

https://hal.science/hal-02572608

Submitted on 12 Aug 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The first Hubble diagram and cosmological constraints using superluminous supernovae

C. Inserra,¹★ M. Sullivan,² C. R. Angus,³ E. Macaulay,⁴ R. C. Nichol,⁵ M. Smith,²,⁶
C. Frohmaier,⁵ C. P. Gutiérrez,²,⁷,⁸ M. Vicenzi,⁵ A. Möller,⁹ D. Brout,¹⁰ P. J. Brown,¹¹
T. M. Davis,¹² C. B. D'Andrea,¹⁰ L. Galbany,¹³ R. Kessler,¹⁴,¹⁵ A. G. Kim,¹⁶ Y.-C. Pan,¹⁷
M. Pursiainen,² D. Scolnic,¹⁵ B. P. Thomas,⁵ P. Wiseman,² T. M. C. Abbott,¹⁸ J. Annis,¹⁹
S. Avila,²⁰ E. Bertin,²¹,²² D. Brooks,²³ D. L. Burke,²⁴,²⁵ A. Carnero Rosell,²⁶,²⁷
M. Carrasco Kind,²⁸,²⁹ J. Carretero,³⁰ F. J. Castander,³¹,³² R. Cawthon,³³ S. Desai,³⁴ H. T. Diehl,¹⁹
T. F. Eifler,³⁵,³⁶ D. A. Finley,¹⁹ B. Flaugher,¹⁹ P. Fosalba,³¹,³² J. Frieman,¹⁵,¹⁹ J. Garcia-Bellido,²⁰
E. Gaztanaga,³¹,³² D. W. Gerdes,³⁷,³⁸ T. Giannantonio,³⁹,⁴⁰ D. Gruen,²⁴,²⁵,⁴¹ R. A. Gruendl,²⁸,²⁹
J. Gschwend,²⁷,⁴² G. Gutierrez,¹⁹ D. L. Hollowood,⁴³ K. Honscheid,⁴⁴ D. J. James,⁴⁵ E. Krause,³⁵
K. Kuehn,⁴⁶ N. Kuropatkin,¹⁹ T. S. Li,¹⁴,¹⁹ C. Lidman,⁴⁷ M. Lima,²⁷,⁴⁸ M. A. G. Maia,²⁷,⁴²
J. L. Marshall,¹⁰ P. Martini,⁴⁴ F. Menanteau,²⁸,²⁹ R. Miquel,³⁰,⁴⁹ A. A. Plazas Malagón,⁵⁰ A. K. Romer,⁵¹
A. Roodman,²⁴,²⁵ M. Sako,¹⁰ E. Sanchez,²⁶ V. Scarpine,¹⁹ M. Schubnell,³⁸ S. Serrano,³¹,³²
I. Sevilla-Noarbe,²⁶ M. Soares-Santos,⁵² F. Sobreira,²⁶,⁵³ E. Suchyta,⁵⁴ M. E. C. Swanson,²⁹
G. Tarle,³⁸ D. Thomas,⁵ D. L. Tucker,¹⁹ V. Vikram,⁵⁵ A. R. Walker,¹⁸ Y. Zhang,¹⁹ J. Asorey,⁵⁶
J. Calcino,⁵⁷ D. Carollo,⁵⁸ K. Glazebrook,⁵⁹ S. R. Hinton,⁵⁷ J. K. Hoormann,⁵⁷ G. F. Lewis,⁶⁰
R. Sharp,⁴⁷ E. Swann,² and B. E. Tucker⁴⁷ (DES Collaboration)

*Affiliations are listed at the end of the paper*

Accepted 2021 April 1. Received 2021 March 9; in original form 2020 April 25

## ABSTRACT

We present the first Hubble diagram of superluminous supernovae (SLSNe) out to a redshift of two, together with constraints on the matter density, $\Omega_M$, and the dark energy equation-of-state parameter, $w(\equiv p/\rho)$. We build a sample of 20 cosmologically useful SLSNe I based on light curve and spectroscopy quality cuts. We confirm the robustness of the peak-decline SLSN I standardization relation with a larger data set and improved fitting techniques than previous works. We then solve the SLSN model based on the above standardization via minimization of the $\chi^2$ computed from a covariance matrix that includes statistical and systematic uncertainties. For a spatially flat $\Lambda$ cold dark matter ($\Lambda$CDM) cosmological model, we find $\Omega_M = 0.38^{+0.24}_{-0.19}$, with an rms of 0.27 mag for the residuals of the distance moduli. For a $w_0w_a$CDM cosmological model, the addition of SLSNe I to a 'baseline' measurement consisting of Planck temperature together with Type Ia supernovae, results in a small improvement in the constraints of $w_0$ and $w_a$ of 4 per cent. We present simulations of future surveys with 868 and 492 SLSNe I (depending on the configuration used) and show that such a sample can deliver cosmological constraints in a flat $\Lambda$CDM model with the same precision (considering only statistical uncertainties) as current surveys that use Type Ia supernovae, while providing a factor of 2-3 improvement in the precision of the constraints on the time variation of dark energy, $w_0$ and $w_a$. This paper represents the proof of concept for superluminous supernova cosmology, and demonstrates they can provide an independent test of cosmology in the high-redshift ($z > 1$) universe.

**Key words:** transients: supernovae – cosmology: dark matter – cosmology: cosmological parameters.

## 1 INTRODUCTION

Twenty years have passed since observations of Type Ia supernovae (SNe Ia) provided the first direct evidence for cosmic acceleration

(Riess et al. 1998; Perlmutter et al. 1999). The physical origin of this acceleration is unknown, but is often described by a phenomenon called ‘dark energy’. Combining SN Ia observations with measurements of large-scale structure (e.g. Percival et al. 2007; Anderson et al. 2014) and the cosmic microwave background (CMB; e.g. Hinshaw et al. 2013; Planck Collaboration XIII 2016) shows that

<footnote>* E-mail: inserrac@cardiff.ac.uk</footnote>
dark energy is the major component (≈70 per cent) of the energy density of the Universe at the present epoch. SNe Ia present a direct and mature method of probing this dark energy via its equation-of-state parameter $w$. Current SN-only measurements provide a precision of 20 per cent, dropping to 4–5 per cent when combined with measurements of the CMB (Scolnic et al. 2018; Abbott et al. 2019). However, SNe Ia at $z \gtrsim 1.2$ are exceptionally challenging to observe from the ground, and thus assembling large samples at these high redshifts is very time-consuming (Riess et al. 2018) due to both the faintness of SNe Ia and line blanketing in their ultraviolet (UV) spectra.

Hydrogen-free superluminous supernovae (SLSNe I; Quimby et al. 2011; Gal-Yam 2012) are significantly more luminous, do not suffer the same degree of line blanketing as SNe Ia and have been observed at higher redshifts than SNe Ia, photometrically out to $z \sim 4$ (Cooke et al. 2012) and spectroscopically out to $z \sim 2$ (Smith et al. 2018). These objects are characterized by a distinctive spectroscopic evolution linking them with massive stars (Pastorello et al. 2010), and show remarkable peak luminosities $\bar{M} < -21$ mag (De Cia et al. 2018; Inserra et al. 2018c; Lunnan et al. 2018; Angus et al. 2019). Their light-curve decline rates and colour evolution are similar, suggesting these events may be standardizable (Inserra & Smartt 2014) via a peak-decline relation in a synthetic band centred at 400 nm.

# 2 SUPERLUMINOUS SUPERNOVA DATA SAMPLE

## 2.1 The superluminous supernova definition and subtypes

The challenge in using SLSNe I as standardizable candles is to find a robust definition of the class that does not simply depend on their luminosity and, ideally, an association with a common explosion mechanism and progenitor scenario to decrease contamination.

In the previous work about SLSNe I standardization (Inserra & Smartt 2014) two observational subclasses of SLSNe I were used and, at the time, it was not immediately clear if these were distinct or if there was a continuum of properties bridging the gap between them. However, this distinction is important if they are to be utilized as standardizable candles, since the bulk of the population (and those showing the strongest correlation parameters) is SLSNe I with light-curve evolution similar to SN 2010gx (Pastorello et al. 2010, hereafter referred to as Fast). The other subtype, encompassing objects similar to SN 2007bi (Gal-Yam et al. 2009; Young et al. 2010, hereafter referred to as Slow), instead increases the scatter on the proposed correlations. This increase in the scatter may be due to the presence of interaction in these objects, which is observed in light curves and spectra (e.g. Yan et al. 2015; Nicholl et al. 2016; Inserra et al. 2017). More recent works have shown that a distinct division can indeed be made between these two classes (Inserra et al. 2018c; Quimby et al. 2018; Gal-Yam 2019a; Inserra 2019), Fast (F) and Slow (S). This distinction is possible via light curve, from peak to +30 d, and spectra information, at roughly +10 d and up to +30 d. This classification, based on K-means partitional cluster analysis, also requires photospheric velocity information derived from the Fe\II $\lambda5169$ line. Those evolving more slowly frequently show signatures of an interaction with a circumstellar medium (Nicholl et al. 2016; Inserra et al. 2017; Yan et al. 2017b; Inserra 2019), perhaps pointing to a different progenitor scenario. The first step in building a homogeneous sample is then to define what is a cosmologically useful SLSN I based on its spectrophotometric behaviour.

## 2.2 The sample construction

We begin to build our sample with all spectroscopically confirmed SLSNe I available in the literature, starting with the 40 SLSNe I from the compilation of Inserra et al. (2018c), and adding nine from the Panoramic Survey Telescope and Rapid Response System 1 (PanSTARRS-1; Lunnan et al. 2018), 15 from the Palomar Transient Factory (PTF) and intermediate PTF (iPTF; De Cia et al. 2018), and 17 from the Dark Energy Survey (DES; Angus et al. 2019).

We apply light-curve quality cuts to our sample in order to assemble a subsample with adequate photometric coverage in a synthetic rest-frame filter centred at 400 nm, which was previously used to test their standardization (Inserra & Smartt 2014). To fulfil our quality cuts, the objects need a light curve covering −15 to +30 d in the rest frame and without multiple peaks to remove ambiguity in identifying the main peak and thus measuring phases (first quality cut). These requirements do not exclude the presence of early time ‘bumps’ (Nicholl et al. 2015b; Smith et al. 2016). Furthermore, a spectrum taken between −15 and +30 d rest frame must also be available (second quality cut).

The literature sample has 23 such SLSNe I. We apply the same selection criteria to the DES SLSN I candidates (Angus et al. 2019). Of 17 events, 10 passed the light-curve quality criteria, and all have at least one spectrum in the required phase range. This retention fraction of 58 per cent is somewhat higher than what seen in the literature sample and can be explained by the DES cadence during the 6 months observing season and higher redshift of several objects (Diehl et al. 2016, 2018; Angus et al. 2019). Of the nine additional, and previously unpublished, events within the Pan-STARRS1 (PS1) Medium Deep survey (Lunnan et al. 2018) only three passed our quality cuts. We also examined events from the PTF/iPTF sample, but none had a sufficient sampling in the rest-frame 400 nm.

However, these published SLSN samples are very heterogeneous in their target selection and therefore we apply a homogeneous method to select our objects. To do this, our third quality cut is based on a statistical approach to identify SLSNe I from their multiband photometric behaviour and the distribution of the candidates on the hypersurface defined by four photometric variables (Four Observables Parameter Space – 4OPS; Inserra et al. 2018c). This parameter space uses the peak luminosity in the 400 nm filter, the decline in magnitudes in the 400 nm filter over the 30 d following peak brightness, the 400–520 colour at peak and the 400–520 colour at +30 d. This hypersurface provides information on the overall SLSN I population evolution, and it is a valuable alternative to identifying SLSNe I when only a single spectrum bearing resemblance to other SLSNe I is available, as is the case for several SLSN I candidates. However, one of the key relationships describing this hypersurface is the standardization relation. To remove any potential bias in the way we select cosmologically useful SLSNe I, we decide to use a slightly different hyperplane than the original. This is needed because if we would have used the original hypersurface, we would have implicitly selected SLSNe fitting the standardization relation and then creating a circular argument. The alterations made to the hypersurface here used with respect to the original one (Inserra et al. 2018c) are the following: (1) we replace panel A (Inserra et al. 2018c) with a different decline relation ($M(400)_{20}$ versus $\Delta M(400)_{30}$) to avoid introducing any biases over the fact that the peak-decline relation ($M(400)_{0}$ versus $\Delta M(400)_{30}$) is a consequence of our selection criteria; (2) we present the decline panel with the colour panel, which is panel D of the original 4OPS. A version with $M(400)_{0}$ rather than $M(400)_{20}$ as y-axis of the left-hand panel can be found in Appendix A.
We will refer to core SLSNe I as those that lie within 2.2σ of the hyperplane here constructed (see Fig. 1). The choice of 2.2σ is driven by the Chauvenet’s criterion value for a sample of this size. Chauvenet’s criterion is an outlier detection method used in experimental physics, but it has also been applied to supernova cosmology (e.g. Conley et al. 2011; Betoule et al. 2014), and can be used for Gaussian-distributed data sets such as the one presented here (Inserra & Smartt 2014). Moreover, each slice of the hypersurface described by the 2.2σ band contour spans ~3 mag in the luminosity scale, supporting the fact that we are not excluding a priori the luminosity extremes of the population. In some cases, like DES14X2byo, the supernova does not make the core population due to a different colour evolution with respect to the SLSN I core, i.e. missing the 2.2σ region in the colour panel (see right-hand panel of Fig. 1). We also note that all SLSNe I candidates not making into the 4OPS bands show an early ‘bump’ or undulation in the light curve after 30 d. Hence, we would have the same data set if the third quality cut was linked to the morphology of the whole light curve instead of the statistical description. This further inspection supports the fact that our selection criteria are not driven by any underlying relations. All explored relations, fit parameters, and statistical results of our sample are reported in Table 3.

The literature sample has 20 objects belonging to the core SLSN I population. Of the 10 DES objects only four reside in the hypersurface (of which one event is a Slow, see Fig. 1), while the others show a similar trend in their luminosity evolution but at lower luminosities, down to roughly M(400) ~ -19.3 mag. This suggests a population of transients similar to SLSNe I with peak magnitudes down to those of normal core-collapse SNe (Angus et al. 2019). Only two of the PS1 objects lay in the core distribution. The events lying outside the core population have a similar photometric behaviour to the DES SLSNe I that do not meet the selection criteria (see Fig. 1).

Hence, the final sample fulfilling all criteria comprises 26 SLSNe I (see Table 1), of which 15 belong to the Fast subclass and six to the Slow subclass; the remaining five, due to their high redshift that prevent measuring the Fe II line, have insufficient spectroscopic data to assign a subclass and hence will be labelled as ‘No Subclass’ (NS). Indeed, it is impossible to observe the Fe II line at $z \gtrsim 0.8$ with optical spectroscopy, although infrared spectroscopy and/or an analysis based on the ejecta velocity of UV lines (see Gal-Yam 2019b) might extend this spectroscopic division out to $z \sim 3.7$.

Because of the presence of interaction in the Slow subclass (Nicholl et al. 2016; Inserra et al. 2017; Yan et al. 2017b; Inserra 2019), which add an additional source of scatter, we use the Fast and those with NS in our analysis, for a final sample of 20 SLSNe I in total. There may, of course, be Slow subtypes present in the NS sample that may bias our results, but such contamination is unavoidable with the current data set if we want to probe the high-redshift region ($z > 1.5$) unexplored with SNe Ia. The impact of this contamination is beyond the scope of this paper, which is intended as a proof of concept. Nevertheless, their addition will provide a less prominent contamination than previous studies in which both Fast and Slow events were included in the analysis.

## 2.3 The reddening assumption

Before attempting any cosmological analysis or standardization relation, it is important that our assumption of negligible host galaxy reddening (or local reddening at the supernova location) is accurate. This assumption is supported by multiple factors. First, the colour distribution around peak for SLSNe I is quite narrow in the optical and UV irrespective of redshift; 0.46 mag in the

optical (Inserra & Smartt 2014; Inserra et al. 2018c) and 0.53 mag in the UV (Smith et al. 2018). This scatter would be significantly increased in the case of environmental reddening. Secondly, SLSNe I UV peak light distribution exhibits a small scatter ($\Delta M_{UV} < 1$ mag) regardless of the redshift (Smith et al. 2018), host reddening would have strongly affected the UV distribution causing a scatter larger than currently observed. Thirdly, SLSNe I spectra around peak epoch show the temperature sensitive O II lines ($12 000 < T(K) < 16 000$; Quimby et al. 2013; Mazzali et al. 2016), hence reddened sources would apparently lie outside of this temperature range. This is not observed. Finally, SLSNe I explode in dwarf, metal-poor galaxies similar to those hosting long gamma-ray bursts (e.g. Lunnan et al. 2014) that have low dust content as shown by high-redshift galaxies hosting gamma-ray bursts (Wiseman et al. 2017). Only a confirmed SLSN I and a candidate one, SN 2017egm (Chen et al. 2017b) and SN 2018don (Lunnan et al. 2020), both Slow events, show a significant host reddening ($E(B - V) > 0.2$). We also exclude any dust formation in the material surrounding SLSNe I (i.e. circumstellar material) because circumstellar material has only been indirectly observed (i.e. light-curve undulations) in the Slow type. Since SLSNe I, intrinsically, are stripped envelope SNe boosted in luminosity (Pastorello et al. 2010; Inserra et al. 2013), the distance and physical conditions of the material expelled by the star before undergoing its final demise are not favourable for early dust production causing a reddening excess of $E(B - V) > 0.02$ at the time-scale needed for our analysis (i.e. during the photospheric phase).

## 3 SUPERNOVA LIGHT CURVES

To estimate and model the light curves of SLSNe I around peak ($-15 \le$ phase (d) $\le$ 30), we first k-correct the observed magnitudes to the 400 and 520 nm synthetic filters with the SNAKE<fnref n="1" /> software package (Inserra et al. 2018a), which also estimates the uncertainties on the k-corrections. The apparent peak magnitude in rest-frame 400 nm band ($m(400)$) is given by

$$m(400) = m(X) - k_{X \to 400}, \qquad (1)$$

where $m(X)$ is the observed apparent magnitude in passband X and the passband is chosen from the observed filters available for each SLSN to be closest in wavelength to 400 nm after accounting for the cosmological redshift ($1 + z$). This process is known as cross-filter k-correction (Kim, Goobar & Perlmutter 1996). $k_{X \to 400}$ is the k-correction from this passband to the 400 nm passband. An analogous relation is used for the 520 nm passband. When observed spectra for a given SLSN I are not available, we use an average SLSN I time series spectral energy distribution (SED) to compute the k-correction (Prajs et al. 2017). We correct all our observed photometry for Milky Way extinction (Schlafly & Finkbeiner 2011), but make no corrections for extinction in the SLSN I host galaxies, which is assumed to be small (Leloudas et al. 2015; Nicholl et al. 2015a; Inserra 2019) as suggested by the small scatter in the colour distribution of SLSNe in the optical (Inserra et al. 2018c; Inserra 2019) and UV (Smith et al. 2018) (see Section 2.3). We then use Gaussian processes (GPs) regression (Bishop 2006; Rasmussen & Williams 2006) to fit the light curves, using the PYTHON package GEORGE and a Matern-3/2 kernel to perform our GP regression of SLSN I light curves and

<footnote n="1">¹https://github.com/cinserra/S3</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Our third criterion for SLSN I selection. A reduced and modified version of the Four Observables Parameter Space (4OPS) for SLSNe I. Data are taken from DES, literature, and other surveys' sample papers that made our first two quality cuts. The left-hand panel shows the magnitude at 20 d post-peak versus the decline rate over 30 d past peak. We have used this relationship to replace the peak–decline relation ($M(400)_0$ versus $\Delta M(400)_{30}$), which had been used in the original 4OPS paper. The right-hand panel shows the colour at peak versus the colour at 30 d post-peak. The literature objects, both Fast (circles) and Slow (open squares), are shown together to their best fit of the weighted linear regression (dashed, black line) and with the 2.2$\sigma$ confidence bands defined by the Chauvenet's criterion. We also include DES16C2nm, which was not presented before in the literature sample (Inserra et al. 2018c).</figcaption></figure>

<table><caption>**Table 1.** SLSN I complete sample. Literature SLSNe I have been broken down in terms of single object papers, for which the contribution from each survey to the total is reported between parentheses, and big sample papers of previously unpublished objects.</caption><thead><tr><th>Sample source</th><th>SLSN I candidates</th><th>Light curve and spectra (quality cuts 1 and 2)</th><th>4OPS (quality cut 3)</th><th>Reference</th></tr></thead><tbody><tr><td>Literature (4 DES – 6 PS1 – 11 PTF/iPTF)</td><td>40</td><td>23</td><td>20</td><td>Inserra et al. (2018c)</td></tr><tr><td>– PS1 (Medium Deep Survey)</td><td>9</td><td>3</td><td>2</td><td>Lunnan et al. (2018)</td></tr><tr><td>– PTF/iPTF</td><td>15</td><td>0</td><td>0</td><td>De Cia et al. (2018)</td></tr><tr><td>DES</td><td>17</td><td>9</td><td>3</td><td>Angus et al. (2019)</td></tr></tbody></table>

derive the uncertainties (see Inserra et al. 2018c, for a more in-depth description).<fnref n="2" />

## 4 SLSN I STANDARDIZATION

We next confirm that the previous observed relationships between peak luminosity ($M(400)_0$) and decline rate in magnitudes over 30 d ($\Delta M(400)_{30}$), here referred to as peak–decline, still hold. To do so, we first convert the rest-frame apparent magnitudes into absolute magnitudes (see Table 2) using the same cosmology of previous studies ($H_0 = 72$ km s$^{-1}$, $\Omega_{\text{matter}} = 0.27$, $\Omega_{\Lambda} = 0.73$) and employ a Bayesian approach to evaluate a weighted linear regression of these parameters, allowing for the uncertainties in both the x and y variables and intrinsic scatter (Kelly 2007). This process uses Bayesian inference that returns random draws from the posterior. Convergence to the posterior is performed using a Markov chain Monte Carlo with $10^5$ iterations. A weighted regression provides a standard deviation bigger than the unweighted one by a factor of roughly $\sum_{i=1}^n 1/\sigma_i$, where n is the sample size. For the peak–decline relation and a sample of 20 objects we retrieve a standard deviation similar to that of the previous unweighted study ($\sigma = 0.33$;

Inserra & Smartt 2014), suggesting that with a bigger sample we have decreased our scatter (see Table 3). This confirms that such relation is quantitatively useful to reduce the intrinsic scatter in the uncorrected peak magnitudes and hence provide a solid proof of concept that SLSNe I may be used as cosmological standardizable candles. The standard deviation of $\sigma = 0.33$ mag decreases, as expected, to $\sigma = 0.26$ mag if we use only the Fast subclass. Including the Slow subclass events substantially increases the dispersion to $\sigma = 0.74$ mag. Such a large dispersion further supports their exclusion.

We also retrieve a similar standard deviation to the previously published $M(400)_0$ versus $\Delta M(400)_{30} - M(520)_{30}$ relation (Inserra & Smartt 2014). We also explore other possible correlations to check if an equally strong relation as those above mentioned can be found at a shorter time-scale (phase < 30 d). We do not find any strong correlation (see Table 3), but the $M(400)_0$ versus $M(400)_{30} - M(520)_{30}$ (peak – colour) relation provides the lowest $\bar{\chi}^2$ ($\bar{\chi}^2 = 0.90$) and $\sigma = 0.19$ mag for the F+NS sample. We also consider correlating both decline and colour information with luminosity ($M(400)_0$ versus ($M(400)_{30} - M(520)_{30}$) $\times \Delta M(400)_{30}$), further reducing the scatter (see Table 3). However, the disadvantage of using such promising correlations is that they need a second, redder band (520 nm). Hence the size and redshift coverage of the sample is smaller than that defined by the peak–decline relation. Nevertheless, when the sample size becomes bigger than the current

<footnote n="2">²https://github.com/cinserra/Gaussian-Processes-GP-</footnote>
<table><caption>**Table 2.** SLSNe I and their peak–decline relation values, together with the quantities used in the third criterion approach.</caption><thead><tr><th>ID</th><th>z</th><th>SLSN I type</th><th>M(400)<sub>0</sub></th><th>ΔM(400)<sub>30</sub></th><th>M(400)<sub>20</sub></th><th>M(400)<sub>0</sub> - M(520)<sub>0</sub></th><th>M(400)<sub>30</sub> - M(520)<sub>30</sub></th></tr></thead><tbody><tr><td>Gaia16apd</td><td>0.102</td><td>F</td><td>-21.87 (0.04)</td><td>0.69 (0.06)</td><td>-21.18 (0.04)</td><td>-0.18 (0.07)</td><td>0.28 (0.07)</td></tr><tr><td>SN2011ke</td><td>0.143</td><td>F</td><td>-21.23 (0.09)</td><td>0.89 (0.09)</td><td>-20.34 (0.02)</td><td>0.04 (0.13)</td><td>0.59 (0.03)</td></tr><tr><td>SN2012il</td><td>0.175</td><td>F</td><td>-21.54 (0.10)</td><td>1.39 (0.17)</td><td>-20.15 (0.14)</td><td>-0.02 (0.11)</td><td>0.48 (0.13)</td></tr><tr><td>PTF11rks</td><td>0.190</td><td>F</td><td>-20.61 (0.05)</td><td>0.87 (0.07)</td><td>-19.74 (0.05)</td><td>0.20 (0.06)</td><td>1.16 (0.15)</td></tr><tr><td>SN2010gx</td><td>0.230</td><td>F</td><td>-21.73 (0.02)</td><td>0.76 (0.03)</td><td>-20.97 (0.02)</td><td>-0.11 (0.02)</td><td>0.53 (0.03)</td></tr><tr><td>SN2011kf</td><td>0.245</td><td>F</td><td>-21.74 (0.15)</td><td>0.52 (0.18)</td><td>-21.22 (0.02)</td><td>...</td><td>...</td></tr><tr><td>LSQ12dif</td><td>0.255</td><td>F</td><td>-21.52 (0.03)</td><td>0.76 (0.04)</td><td>-20.76 (0.03)</td><td>0.05 (0.03)</td><td>0.57 (0.10)</td></tr><tr><td>LSQ14mo</td><td>0.256</td><td>F</td><td>-21.04 (0.05)</td><td>1.30 (0.14)</td><td>-19.74 (0.13)</td><td>-0.08 (0.04)</td><td>0.61 (0.02)</td></tr><tr><td>PTF09cnd</td><td>0.258</td><td>F</td><td>-22.16 (0.08)</td><td>0.71 (0.14)</td><td>-21.45 (0.12)</td><td>...</td><td>...</td></tr><tr><td>SN2013dg</td><td>0.265</td><td>F</td><td>-21.35 (0.05)</td><td>1.03 (0.06)</td><td>-20.32 (0.03)</td><td>-0.26 (0.08)</td><td>0.56 (0.10)</td></tr><tr><td>PS1-10bzj</td><td>0.650</td><td>F</td><td>-21.03 (0.06)</td><td>1.23 (0.32)</td><td>-19.08 (0.31)</td><td>0.15 (0.11)</td><td>0.94 (0.25)</td></tr><tr><td>iPTF13ajg</td><td>0.740</td><td>F</td><td>-22.42 (0.07)</td><td>0.19 (0.10)</td><td>-22.23 (0.07)</td><td>-0.29 (0.09)</td><td>-0.11 (0.09)</td></tr><tr><td>DES15X3hm</td><td>0.860</td><td>F</td><td>-21.94 (0.06)</td><td>1.44 (0.07)</td><td>-21.19 (0.04)</td><td>-0.21 (0.04)</td><td>0.11 (0.06)</td></tr><tr><td>DES17X1amf</td><td>0.920</td><td>NS</td><td>-21.97 (0.07)</td><td>0.26 (0.15)</td><td>-21.78 (0.09)</td><td>-0.31 (0.10)</td><td>0.08 (0.10)</td></tr><tr><td>PS1-10ky</td><td>0.956</td><td>F</td><td>-22.05 (0.06)</td><td>0.61 (0.07)</td><td>-21.44 (0.04)</td><td>-0.06 (0.07)</td><td>0.25 (0.06)</td></tr><tr><td>PS1-11aib</td><td>0.997</td><td>NS</td><td>-22.05 (0.07)</td><td>0.31 (0.17)</td><td>-21.78 (0.40)</td><td>-0.45 (0.07)</td><td>-0.15 (0.07)</td></tr><tr><td>SCP-06F6</td><td>1.189</td><td>F</td><td>-22.19 (0.03)</td><td>0.57 (0.15)</td><td>-21.62 (0.15)</td><td>...</td><td>...</td></tr><tr><td>PS1-11tt</td><td>1.283</td><td>NS</td><td>-21.89 (0.16)</td><td>0.15 (0.20)</td><td>-21.83 (0.12)</td><td>...</td><td>...</td></tr><tr><td>PS1-11bam</td><td>1.565</td><td>NS</td><td>-22.45 (0.10)</td><td>0.36 (0.14)</td><td>-22.09 (0.10)</td><td>...</td><td>...</td></tr><tr><td>DES16C2nm</td><td>1.998</td><td>NS</td><td>-22.52 (0.10)</td><td>0.67 (0.11)</td><td>-22.13 (0.09)</td><td>...</td><td>...</td></tr></tbody></table>

<table><caption>**Table 3.** Fit parameters and statistical results of our sample.</caption><thead><tr><th>x</th><th>y</th><th>SLSN I type</th><th>N (objects)</th><th>β</th><th>α</th><th>σ</th><th>x̄<sup>2</sup></th></tr></thead><tbody><tr><td rowspan="3">ΔM(400)<sub>30</sub></td><td rowspan="3">M(400)<sub>0</sub></td><td>F</td><td>15</td><td>-23.09 ± 0.28</td><td>1.01 ± 0.17</td><td>0.26 ± 0.22</td><td>1.23</td></tr><tr><td>F+NS</td><td>20</td><td>-22.62 ± 0.19</td><td>0.72 ± 0.14</td><td>0.33 ± 0.22</td><td>1.73</td></tr><tr><td>F+NS+S</td><td>25</td><td>-22.31 ± 0.15</td><td>0.53 ± 0.13</td><td>0.36 ± 0.22</td><td>2.31</td></tr><tr><td rowspan="3">Δ(M(400)<sub>30</sub> - M(520)<sub>30</sub>)</td><td rowspan="3">M(400)<sub>0</sub></td><td>F</td><td>12</td><td>-22.76 ± 0.35</td><td>2.18 ± 0.62</td><td>0.29 ± 0.30</td><td>1.45</td></tr><tr><td>F+NS</td><td>14</td><td>-22.79 ± 0.30</td><td>2.22 ± 0.55</td><td>0.25 ± 0.25</td><td>1.20</td></tr><tr><td>F+NS+S</td><td>18</td><td>-22.31 ± 0.20</td><td>1.55 ± 0.41</td><td>0.32 ± 0.24</td><td>1.91</td></tr><tr><td rowspan="3">M(400)<sub>20</sub> - M(520)<sub>20</sub></td><td rowspan="3">M(400)<sub>0</sub></td><td>F</td><td>12</td><td>-21.97 ± 0.14</td><td>1.97 ± 0.47</td><td>0.31 ± 0.27</td><td>1.58</td></tr><tr><td>F+NS</td><td>14</td><td>-21.96 ± 0.12</td><td>1.89 ± 0.41</td><td>0.28 ± 0.24</td><td>1.41</td></tr><tr><td>F+NS+S</td><td>18</td><td>-21.77 ± 0.09</td><td>1.32 ± 0.33</td><td>0.33 ± 0.24</td><td>1.74</td></tr><tr><td rowspan="3">M(400)<sub>30</sub> - M(520)<sub>30</sub></td><td rowspan="3">M(400)<sub>0</sub></td><td>F</td><td>12</td><td>-22.27 ± 0.14</td><td>1.51 ± 0.27</td><td>0.21 ± 0.19</td><td>0.87</td></tr><tr><td>F+NS</td><td>14</td><td>-22.21 ± 0.11</td><td>1.41 ± 0.23</td><td>0.19 ± 0.17</td><td>0.95</td></tr><tr><td>F+NS+S</td><td>18</td><td>-21.95 ± 0.09</td><td>1.02 ± 0.21</td><td>0.28 ± 0.20</td><td>1.52</td></tr><tr><td rowspan="3">ΔM(400)<sub>30</sub></td><td rowspan="3">M(400)<sub>20</sub></td><td>F</td><td>15</td><td>-23.10 ± 0.31</td><td>1.55 ± 0.20</td><td>0.27 ± 0.23</td><td>1.82</td></tr><tr><td>F+NS</td><td>20</td><td>-22.67 ± 0.20</td><td>1.29 ± 0.15</td><td>0.31 ± 0.23</td><td>1.89</td></tr><tr><td>F+NS+S</td><td>25</td><td>-22.31 ± 0.16</td><td>1.07 ± 0.13</td><td>0.35 ± 0.23</td><td>2.34</td></tr><tr><td rowspan="3">M(400)<sub>30</sub> - M(520)<sub>30</sub> × ΔM(400)<sub>30</sub></td><td rowspan="3">M(400)<sub>0</sub></td><td>F</td><td>12</td><td>-22.21 ± 0.12</td><td>0.82 ± 0.14</td><td>0.19 ± 0.19</td><td>0.86</td></tr><tr><td>F+NS</td><td>14</td><td>-22.12 ± 0.09</td><td>0.73 ± 0.11</td><td>0.18 ± 0.16</td><td>0.61</td></tr><tr><td>F+NS+S</td><td>18</td><td>-21.95 ± 0.08</td><td>0.60 ± 0.11</td><td>0.25 ± 0.18</td><td>1.03</td></tr></tbody></table>

Note. Least-squares fits for a Bayesian weighted linear regression with weighted errors both in x and y of the form $\eta = \beta + \alpha \times x' + \epsilon$, where $x = x' + x_{\text{err}}$ and $y = \eta + y_{\text{err}}$. The $\sigma$ is the standard deviation of this fit. The last column gives the reduced $\chi^2$.

one ($\gtrsim 100$), these two relations including the 520 nm band might be more effective than the peak–decline.

## 5 SUPERNOVA MODEL

We therefore use the peak–decline standardization method in our analysis, i.e. our distance estimator assumes that SLSNe I with an identical light-curve decline rate have the same average intrinsic luminosity at all redshifts. The standardized distance modulus, $\mu_{\text{obs}}$, is then given by

$$ \mu_{\text{obs}} = m(400) - M(400) + \gamma \Delta M(400)_{30}, \quad (2) $$

where $m(400)$ is the peak apparent magnitude in rest-frame 400 nm band, and $M(400)$ (the peak absolute magnitude) and $\gamma$ are nuisance parameters in the distance estimate. This is compared to the model distance modulus, $\mu_{\text{model}}$, of $5 \log_{10}(d_L/10 \text{pc})$, where $d_L$ is the luminosity distance. The fit then minimizes the $\chi^2$ according to

$$ \chi^2 = \Delta\mu^T \cdot C^{-1} \cdot \Delta\mu, \quad (3) $$

where $C$ is the covariance matrix and $\Delta\mu$ is the vector of residuals $\Delta\mu = \mu_{\text{obs}} - \mu_{\text{model}}$. Note that the Hubble constant, $H_0$, enters in both $M(400)$ and $d_L$, and thus does not affect (and is not constrained by) the cosmological fit. We assume an unperturbed Friedmann–Lemaître–Robertson–Walker metric and a flat universe, i.e. $\Omega_M +$
$\Omega_\Delta = 1$, and the free parameters in the fit are therefore $\Omega_M$ and the two nuisance parameters $M(400)$ and $\gamma$.

The covariance matrix is defined as the sum of the statistical and systematic parts:

$$C = C_{\text{stat}} + C_{\text{sys}}. \qquad (4)$$

The statistical covariance matrix is diagonal, while the systematic covariance matrix can contain off-diagonal terms capturing the covariance between different events. Here we define systematic uncertainties as those terms whose effect on our final uncertainty budget could not be reduced by increasing the size of the SLSN I sample (i.e. reddening, surveys zero-points, Malmquist bias, and the light-curve fitting method). We note that, due to the limited size of our sample, our analysis is dominated by statistical uncertainties (i.e. uncertainties on fitted light-curve parameters).

In our minimization technique of equation (3) we use the following priors: the Joint Light-Curve Analysis (JLA) result (SN analysis only) as a prior on $\Omega_M$; the $M(400)$ and $\gamma$ values previously retrieved (Inserra & Smartt 2014) for the $M(400)_{0}$ versus $\Delta M(400)_{30}$ as a prior of the other two nuisance parameters ($M(400)$, $\gamma$). We also assume a Gaussian distribution for the form of the priors. To minimize equation (3) we use IMINUIT,<fnref n="3" /> a minimization technique based on MINUIT (James & Roos 1975) that runs over $10^5$ iterations, and a Markov chain Monte Carlo Ensemble sampler (Foreman-Mackey et al. 2013) with $5 \times 10^3$ iterations. Both provide non-Gaussian distributed uncertainties and similar results, although we note that IMINUIT uncertainties are always a factor of ~1.5 bigger than those from EMCEE and this is likely due to the small data set. This is confirmed by the analysis executed with a bigger (847 objects) simulated data set for which both algorithms give similar uncertainties (see Section 10).

## 5.1 Statistical uncertainties

The statistical uncertainties part is diagonal and includes uncertainties as follows:

$$C_{\text{stat},ii} = \sigma_{m_{400},i}^2 + \gamma^2 \sigma_{\Delta M(400)_{30}}^2 + \left( \frac{5(1+z_i)}{z_i(1+z_i/2)\log(10)} \right)^2 \times \sigma_{z,i}^2 + \sigma_{\text{lensing}}^2. \qquad (5)$$

Here $\sigma_{m_{400},i}$ and $\sigma_{\Delta M(400)_{30}}$ are the uncertainties on the fitted light-curve parameters. The third term is associated with our choice of an empty-universe approximation for the relation between redshift uncertainty and the associated magnitude uncertainty (Davis et al. 2011). The last term is a random uncorrelated scatter due to lensing, $\sigma_{\text{lensing}} = 0.055z$, following the prescription used for SNe Ia (Conley et al. 2011). The lensing dispersion for point sources depends on the line-of-sight density distribution, not the source properties, so this is appropriate even though our SLSN I population differs from the SNe Ia population. Further studies should address whether this functional form is appropriate for the high redshifts in our sample, however even at $z = 1.5$ the lensing dispersion is only $\sigma_{\text{lensing}} \sim 0.08$ mag. Since this is an order of magnitude lower than the dispersion in magnitudes (the first two terms) and the mean lensing magnification should be zero, we consider any possible lensing bias to be negligible.

## 5.2 Systematic uncertainties

The definition of systematic uncertainties is not always unambiguous and it depends on labels given by authors (e.g. Conley et al. 2011; Betoule et al. 2014; Scolnic et al. 2018). Here we interpret them as those terms whose effects on our final uncertainty budget could not be reduced by increasing the SLSN I sample. We also note that, due to the limited size of our sample, our analysis is dominated by statistical uncertainties and variations in the systematics have little leverage on our cosmological constraints.

There is no standard method for handling supernova systematic effects, but the most common approach is to initially fit the data set without any systematic effects (hence only with $C_{\text{stat},ii}$) and then marginalize over all the systematic terms by adding the systematic part of the covariance matrix to the statistical as in equation (4). The matrix is given by

$$C_{\text{sys},i,j} = \sigma_{\text{reddening},i,j}^2 + \sigma_{\text{ZP},i,j}^2 + \sigma_{\text{MalmquistBias},i,j}^2 + \sigma_{\text{model},i,j}^2. \qquad (6)$$

### 5.2.1 Reddening

The first term is related to reddening and we account only for Milky Way reddening along the line of sight, including an estimated 10 per cent random uncertainty for each SLSN I due to the conversion from dust column density to extinction (Schlegel, Finkbeiner & Davis 1998). The extinction is always $E(B-V) < 0.02$ mag and for this low value the extinction law is almost insensitive to the choice of $R_V$, hence the assumption of a Galactic value of $R_V = 3.1$ is appropriate (Cardelli, Clayton & Mathis 1989). At this stage we do not consider host galaxy reddening (see discussion in Section 2.3) since SLSNe I explode in dwarf galaxies with no reported host galaxy extinction for the majority of events (~85 per cent, data from Lunnan et al. 2014; Leloudas et al. 2015; Angus et al. 2016; Perley et al. 2016; Chen et al. 2017a; Schulze et al. 2018). When a host galaxy extinction value of $E(B-V) > 0.02$ is reported that is due to SED modelling (Schulze et al. 2018) rather than galaxy line analysis and hence exposed to larger uncertainties. However, in future analysis this term should be investigated more carefully and might be taken in consideration.

### 5.2.2 Zero-points

The second term is due to the uncertainties in the zero-points of each survey in each filter and each field. For example, in the case of DES we have four different filters and 10 SN deep fields, of which two are approximately 1 mag deeper in order to extend SN searches out to higher redshift. However, the difference between the zero-points of different fields is usually of the order of $10^{-3}-10^{-4}$ mag (Diehl et al. 2014) and so are their general uncertainties. These uncertainties are at a mmag level and hence we consider a general uncertainty for each filter in each survey. The general zero-points uncertainties are retrieved for DES (Diehl et al. 2014) and PS1 (Schlaflay et al. 2012; Tonry et al. 2012). For the rest of literature SLSNe I, which were not found by these surveys or did not have the majority of their data obtained from a single survey, we considered average<fnref n="4" /> zero-point uncertainties for each filter matching those reported by other large surveys (Tonry et al. 2012; Diehl et al. 2014; Smartt et al. 2015).

<footnote n="3">³https://github.com/iminuit/iminuit</footnote>

<footnote n="4">⁴Average evaluated from the uncertainties of DES and PS1.</footnote>
### 5.2.3 Malmquist bias

To model our search efficiency (i.e. how many SLSNe I are missed), for each SLSN I we simulated $10^4$ light curves using a Monte Carlo approach. Because of the relatively small sample of our analysis we treat it as a simple smoothed offset between nearby ($z<0.4$), medium redshift ($0.4 \le z \le 1.0$), and distant SNe ($z>1$). After correction for light-curve shape, we find the magnitude offset to be 0.02 mag up to $z=1.0$ and less than 0.05 mag for the high-redshift objects of our sample. This is in agreement with the fact that the majority of the $z>1$ SLSNe I were discovered before maximum light similar to those at lower redshift, suggesting that our overall sample is equally biased at all redshifts. However, a more in-depth treatment of the Malmquist bias might be needed with bigger samples. That should make use of the predicted number of SLSNe I for an unbiased survey to the number observed as a function of redshift, for which DES is an ideal testbed (Angus et al. 2019; Thomas et al., in preparation).

### 5.2.4 Light-curve fitting

The last term, $\sigma_{\text{model}}$, relates to our uncertainties in the interpolation used to fit our light curve, which means the kernel chosen for the GPs. We use a Matern-3/2 kernel that we find a suitable choice to avoid overfitting and retrieve a balanced precision/recall outcome (Inserra et al. 2018c). However, other kernels can be used, such as the Matern-5/2, and we have to take into account the differences in the light-curve outputs with different kernels. This term, in principle, could be reduced with more SLSNe I, but is correlated between different SLSNe I and therefore cannot be included in $C_{\text{stat,ii}}$.

### 5.2.5 Additional uncertainties

We do not introduce any additional term for contamination from other SN types since all our objects have spectroscopic confirmation and have passed the 4OPS criterion. We also did not include any systematic related to peculiar velocities or a discontinuous step in the local expansion (Hubble bubble; Jha, Riess & Kirshner 2007) since all our SLSNe I have $z > 0.1$. We also do not include any correction from host-galaxy properties since they all reside in similar galaxies at both low and high redshift (Leloudas et al. 2015) and no mass step function has been currently observed for SLSNe I as has been seen for Type Ia (Sullivan et al. 2010). A possible differential evolution with galaxy mass at high redshift was recently claimed in an analysis using rest-frame optical data of SLSN host galaxies (Schulze et al. 2018). However, this analysis does not hold if the rest frame is extended to include wavelengths bluer than the B band. Hence, we do not consider any additional source of uncertainties linked to galaxy evolution.

## 6 SLSN I HUBBLE DIAGRAM

Our sample size (20 F+NS SLSNe) is sufficient to provide a constraint on a single parameter driving the evolution of the expansion rate. In particular, in a flat universe with a cosmological constant (hereafter flat $\Lambda$CDM), SLSNe I alone can provide a measurement of the reduced matter density $\Omega_M$. The SLSN I Hubble diagram and the flat $\Lambda$CDM best fit, derived from the minimization above, are shown in Fig. 2. The fit parameters are given in the first row of Table 4. We find a best-fitting value of $\Omega_M = 0.38_{-0.19}^{+0.24}$ and rms = 0.27 mag for the residuals of the distance moduli, also shown in Fig. 2. For comparison, the dispersion in the JLA SN Ia sample is rms = 0.17 mag for 740 SNe Ia over $0.01 < z < 1.2$ (Betoule et al.

2014). This is also shown in Fig. 2, where the JLA sample and its residuals are overplotted. The redshift coverage makes it possible to assess the overall consistency of the SLSN I data with the flat $\Lambda$CDM model.

In Fig. 3, we plot the residuals of our sample without the decline-rate correction, as a function of the decline rate ($\Delta M(400)_{30}$). This shows the brighter-slower relationship of the standardization, with no apparent evolution in residuals across this relationship for our SLSN I sample. We search for any further significant trends between decline, colour, and Hubble residuals, but find none. Thus, if further parameters are capable of decreasing the scatter in the residuals of our cosmological fit, they are either related to quantities that we have not measured, or larger samples are required to investigate them.

## 7 $w_0 w_a$ CDM COSMOLOGY WITH SLSN I

We also explore how and to what extent this SLSN I sample can improve the constraints on the redshift-dependent equation-of-state of dark energy (Abbott et al. 2019), $w(a) = w_0 + (1-a)w_a$, where $a = (1+z)^{-1}$ is the cosmological scale factor. In our analysis, we include priors on the cosmological parameters from measurements of the CMB from Planck (Planck Collaboration XIII 2016). We assume a Gaussian distribution for the form of the prior, which we construct at the maximum likelihood value used by the Planck consortium (Planck Collaboration XIII 2016). We also include the JLA SNe Ia sample (Betoule et al. 2014). We adopt a flat $w_0 w_a$ CDM cosmological model, and a SLSN I likelihood of the form $\ln L \propto \chi^2$, where the $\chi^2$ is given by equation (3), and $\mu_{\text{model}}$ is now a function of $w_0$, $w_a$, $\Omega_M$, and $H_0$. To obtain convergence we use the CosmoMC tool (Lewis & Bridle 2002; Lewis 2013).

We find marginalized constraints on $w_0 = -0.904^{\pm 0.185}$ and $w_a = -0.594^{\pm 0.926}$, which are shown in Fig. 4. To measure the improvement on the constraints on $w_0$ and $w_a$, we evaluate the figure of merit (FoM) proposed by the Dark Energy Task Force (Albrecht et al. 2006), which is the area enclosed by the two standard deviation contours in the $w_0-w_a$ plane ($1/(\sigma_{w_0} \sigma_{w_a})$). Without SLSNe I, the FoM is 5.622, which increases by 4 per cent to 5.835 with the addition of SLSNe I. We then compare this improvement with constraints from Lyman $\alpha$ forest baryonic acoustic oscillations (BAO; de Sainte Agathe et al. 2019), which is also a high-redshift probe like SLSNe I. The FoM for the CMB+JLA+Lyman $\alpha$ is 5.9, suggesting that at present SLSNe are comparable to Lyman $\alpha$ BAO at such redshifts. That is also true of we analyse the uncertainties of the three data combinations (Base, Base+SLSNe, Base+Lyman $\alpha$) on $w_0$ and $w_a$ that are 0.187, 0.185, and 0.183 ($\sigma_{w_0}$) and 0.952, 0.926, and 0.928 ($\sigma_{w_a}$).

## 8 ENVIRONMENT

We investigate if there is an additional dependence on the global characteristics of SLSN I host galaxies. We retrieve SLSN I host galaxies information from Schulze et al. (2018), namely specific star formation rate (sSFR) and stellar mass ($M_s$), which used the SED fitting algorithm LEPHARE (Arnouts et al. 1999; Ilbert et al. 2006) and a Chabrier initial mass function (Chabrier 2003). This approach gives almost identical results to the MAGPHYS SED fitting (da Cunha, Charlot & Elbaz 2008) as shown by the SLSN I data set in literature (Chen et al. 2017a; Schulze et al. 2018). We then apply the same approach to the only DES SLSN I with host galaxy photometry that was not presented in the Schulze et al. (2018) sample (DES16C2nm; Smith et al. 2018). However, we note that the broad-band SED fitting approach used here is a relatively crude way to determine galaxy
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Upper panel: the Hubble diagram of our SLSN I sample (only F+NS subtypes) using the $\Delta M(400)_{30}$ standardization method. Overplotted (solid red line) is the best-fitting flat $\Lambda$CDM cosmology and its uncertainties (shaded red area) with $\Omega_M = 0.38_{-0.19}^{+0.24}$, measured only using SLSNe I. Lower panel: the residuals of each SLSN I from the best-fitting cosmology (red line with respect to red left label) as a function of redshift. The JLA SN Ia compilation from their best fit (dashed blue line) and residuals (grey dots versus blue line) is also shown as comparison. We chose the JLA sample because both studies use the same approach to derive the cosmological constraints, unlike the most recent Pantheon sample.</figcaption></figure>

<table><caption>**Table 4.** Best-fitting parameters for the flat $\Lambda$CDM model using SLSNe I alone and using the $M(400)_0$ versus $\Delta M(400)_{30}$ standardization relation.</caption><thead><tr><th>Sample</th><th>Ω<sub>M</sub></th><th>M (value)</th><th>γ</th><th>Residuals (rms)</th><th>Λ = 0 residuals (rms)</th><th>Redshift</th><th>χ<sup>2</sup>/degrees of freedom</th></tr></thead><tbody><tr><td>F+NS (stat+sys)</td><td>0.38<sup>+0.24</sup><sub>-0.19</sub></td><td>-22.40<sup>+0.38</sup><sub>-0.35</sub></td><td>0.62<sup>+0.19</sup><sub>-0.19</sub></td><td>0.27</td><td>0.31</td><td>~2.0</td><td>11.1/18</td></tr><tr><td>F+NS (stat)</td><td>0.37<sup>+0.31</sup><sub>-0.19</sub></td><td>-22.45<sup>+0.38</sup><sub>-0.35</sub></td><td>0.63<sup>+0.19</sup><sub>-0.19</sub></td><td>0.25</td><td>0.27</td><td>~2.0</td><td>11.0/18</td></tr><tr><td>F (stat+sys)</td><td>0.22<sup>+0.29</sup><sub>-0.21</sub></td><td>-23.20<sup>+0.56</sup><sub>-0.53</sub></td><td>0.97<sup>+0.25</sup><sub>-0.25</sub></td><td>0.26</td><td>0.40</td><td>~1.2</td><td>6.4/13</td></tr><tr><td>F (stat)</td><td>0.41<sup>+0.22</sup><sub>-0.20</sub></td><td>-22.83<sup>+0.43</sup><sub>-0.41</sub></td><td>0.88<sup>+0.22</sup><sub>-0.22</sub></td><td>0.21</td><td>0.26</td><td>~1.2</td><td>4.9/13</td></tr></tbody></table>

properties and hence this analysis is only an initial investigation into host galaxy dependencies.

In Fig. 5 (left-hand panels), we compare host galaxy properties with our residuals, but we do not find a mass step function (or any relation). We observe mild correlations between the host galaxy sSFR and the light-curve properties in terms of decline (Pearson $r = 0.45$) and colour (Pearson $r = 0.52$). This suggests that SLSN I in low-sSFR host galaxies are redder and faster decliners than those in high sSFR. There is no appreciable trend with the stellar mass of host galaxies. From this analysis it seems that there is not a systematic uncertainty in the residuals introduced by SLSN I host properties, but further

analysis with a bigger sample and a more precise estimates of global host properties and local to the SN environment is encouraged.

## 9 EXPLORING SLSNE I IN THE ULTRAVIOLET

Motivated by studies exploring the velocity and shape of lines of ionized elements in the UV (e.g. Gal-Yam 2019a), we measure the pseudo-equivalent width (pEW) of the C III/C II/Ti III and Mg II/C II blended lines at ~2200 and ~2800 Å, respectively, to look for a more quantitative method of distinguishing between Fast and Slow subtypes at high redshift. We check if, combining the UV line pEWs
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Left: the SLSN I residuals uncorrected for the decline parameter, plotted as a function of $\Delta M(400)_{30}$. This diagram computes distance modulus and returns the brighter-slower relationship. Right: histogram of the data distribution. The bin dimension has been chosen according to the Freedman-Diaconis estimator, which accounts for data variability and data size, and is optimized for smaller data sets.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Constraints on the dark energy equation-of-state parameters $w_0$ and $w_a$. We illustrate the one (filled) and two (unfilled) standard deviation contours, and the maximum likelihood values for the ‘Base’ configuration (JLA+CMB; blue dashed lines and square centroid) for estimates calculated including SLSNe I (dark red dot-dashed lines and circle centroid), or Lyman $\alpha$ BAO (orange dotted lines and triangle centroid). The horizontal and vertical dashed lines at $w_0 = -1$ and $w_a = 0$ correspond to a cosmological constant. The Base+SLSNe I results are $w_0 = -0.904^{+0.185}_{-0.926}$ and $w_a = -0.594^{+0.926}_{-0.904}$, which is a small improvement of 4 per cent with respect to the joint CMB+JLA data set, similar to that obtainable with the Base+Lyman $\alpha$ configuration.</figcaption></figure>

with light-curve evolution around peak ($-10 < \text{phase} < +30$ d), we can find similar clusters to those observed in the optical, using Fe II lines. The prefix ‘pseudo’ is used because the reference continuum level chosen does not represent the true underlying continuum level of the supernova. It defines the strength of the line with respect to the pseudo-continuum at any given time. We choose the Mg II/C II line since it is easy to sample at $0.1 \lesssim z \lesssim 3.0$, which covers this data set and the majority of future data sets (see Section 10), and it is a good proxy for the outermost layers of the carbon/oxygen-rich material (Mazzali et al. 2016). The C III/C II/Ti III line is also a good proxy, when available, and it is usually stronger than Mg II/C II. Such difference in strengths is likely due to a lower excitation potential of the lines at 4200 Å and a bigger contribution to the blending from

carbon lines. Nevertheless, it is harder to sample for our redshift baseline.

We collected 10 SLSNe I sampling these UV lines (Barbary et al. 2009; Chomiuk et al. 2011; McCrum et al. 2014; Vreeswijk et al. 2014; Lunnan et al. 2016; Yan et al. 2017a; Quimby et al. 2018; Smith et al. 2018; Angus et al. 2019). Eight of them also show the Fe II $\lambda$5169; six were identified as Fast, and two as Slow. The other two do not have Fe II lines sampled due to their higher redshift and are labelled as NS. For six of them (4F+2S) we have both pre- and post-peak spectra. We measure the pEW (see Gutiérrez et al. 2017, for further details in the methodology) and for each SN we do not observe any change in the pEW values from $-10$ to $+10$ d. Hence we group our measurements as ‘at peak epoch’ ($-10 < \text{phase} < +10$ d). We also measure the line velocities and find $18\,000 < v (\text{km s}^{-1}) < 23\,000$ in agreement with previously results (Chomiuk et al. 2011; Vreeswijk et al. 2014; Mazzali et al. 2016). However, due to the blending of several ions, we decide to only focus on pEWs, which are less sensitive to the signal-to-noise ratio and flux-calibration issues than velocities, flux ratios, and line depths (Folatelli et al. 2013).

We then compare the pEWs and their ratio ($R = \text{pEW}(\text{Mg II/C II})/\text{pEW}(\text{C III/C II/Ti III})$) with the light-curve decline ($\Delta M(400)_{30}$). In Fig. 6 no clear groups are observed comparing the pEW(Mg II/C II) or the ratio with the decline, panels (A) and (C), respectively. Instead, a mild correlation (Pearson $r = 0.76$) is shown in panel (A). When comparing the pEW(C III/C II/Ti III) with the decline (panel B) or with the pEW(Mg II/C II) in panel (D), we retrieve promising clusters. However, the data set is rather small and no further analysis can be done, although in future the pEW(Mg II/C II) versus pEW(C III/C II/Ti III) analysis could give useful information for their characterization at high redshift ($z \gtrsim 0.8$).

In the future we might use such UV information to add another nuisance parameter, e.g. pEW(Mg II/C II), to equation (2). This additional parameter would transform equation (2) into

$$ \mu = m_{400} - M(400) + \alpha \times (\Delta M(400)_{30} + \beta \times \text{pEW}(\text{Mg II/C II})), \quad (7) $$

with three nuisance parameters ($\alpha$, $\beta$, $M(400)$) and would allow SLSNe I Slow to be included in the standardization. Although this approach is appealing, we need a larger UV spectroscopic data set to confirm the findings described above.

## 10 SLSNE I COSMOLOGY: FUTURE AND IMPROVEMENTS

To understand the future potential of SLSNe I in cosmology, we consider SLSN I rates for the Euclid satellite (135 high-quality SNe in 5 yr; Inserra et al. 2018b), and SLSN I predictions for the deep drilling fields of the Legacy Survey of Space and Time (LSST) at the Vera Rubin Observatory. At the moment of writing this paper, the final cadence and number of deep drilling fields are still under evaluation. Hence, we assume two configurations. A generic set-up with 10 deep drilling fields, visited 180 d each year with a 5 d griz cadence. The single visit depths are 25.0, 24.7, 24.0, and 23.3 mag in griz, respectively (AB magnitudes for a $5\sigma$ point source; LSST Science Collaboration et al. 2009). The second configuration is that proposed in the white paper of the Dark Energy Science Consortium (DESC) at the Vera Rubin Observatory. This proposes five deep drilling fields, multiple visits per night in griz a 4 d cadence and depths of 25.3, 25.0, 24.8, and 24.5 mag. Following the methodology of the previous work of Prajs et al. (2017), we use an average peak
<figure><img src="image_7.png" /><figcaption>**Figure 5.** SLSNe I residuals from the best-fitting flat $\Lambda$CDM cosmology (left-hand panels), decline over 30 d (middle panels), and colour at 30 d (right-hand panels) as a function of host galaxy $M_*$ (upper row) and sSFR (bottom row). The error bars on the individual SLSNe I are taken from the SED fitting for the sSFR and $M_*$ axes and are the statistical errors propagated through the light-curve fitting for the residual and observables axes. Bayesian weighted linear regression (blue solid line) is also displayed with the exception of the residuals.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Spectroscopic, pEW(Mg II/C II) at 2800 Å, pEW(C III/C II/Ti III) at 2200 Å, and their ratio (R) versus photometric, $\Delta M(400)_{30}$, measurements. A mild trend is displayed in panel (A) (diamond markers), with a Pearson r coefficient of 0.76. There are promising relationships involving pEW(C III/C II/Ti III) in panels (B) and (D), although a larger UV sample is required to confirm this.</figcaption></figure>

SLSN I luminosity of $M(400)_0 = -21.756 \pm 0.495$ mag (where the average and standard deviations are determined from combing the F+NS subgroups), a model SED, and a spectral template for k-correction. This method is consistent with that previously used to predict the number of SLSN I in the LSST wide survey (Scovacricchi et al. 2016). Here, we only consider SLSNe I that have been detected four times in at least three filters that is, for consistency, the same that has been done for the *Euclid* SLSN I rates (Inserra et al. 2018b).

We then retrieve 929 and 441 SLSNe I in the range $0.25 < z < 3.95$ for the ‘generic’ and ‘DESC’ configuration, respectively. We note that even with an unfavourable LSST cadence to discover and monitor transients, such as normal supernovae, we would expect to recover SLSN I at $z > 1$ due to their intrinsic high luminosity and slow evolution that will be further exaggerated due to time dilation in the observer frame (Inserra et al. 2018b; Moriya et al. 2019). We note that our results are not as optimistic as those published in Villar,
Nicholl & Berger (2018). However, the methodology followed here is based on an observed luminosity distribution while that of Villar et al. (2018) is based on a prescription of a magnetar model that has been shown to provide non-physical results or discordant values with those of other prescriptions (see table A1 in Nicholl, Guillochon & Berger 2017).

We also make predictions on the number of suitable SLSNe I that will be observed by the Zwicky Transient Facility (ZTF; Bellm 2014) and the Asteroid Terrestrial-impact Last Alert System (ATLAS; Tonry et al. 2018) at $z < 0.25$. We assume the low-redshift SLSN I rates reported in literature (Quimby et al. 2013; Prajs et al. 2017, which include Slow events) and scale the star formation history (Li 2008) accordingly to construct a volumetric rate evolution. Based on assumptions of the ZTF and ATLAS observing strategy (Bellm 2014; Smith et al. 2020), and a $5\sigma$ depth of 20.5⁵ and 19.5 mag, respectively, we calculate the number of events ZTF and ATLAS would be capable of following to measure a decline of ~2 mag from peak. We retrieve a very conservative number of 21 SLSNe I per year out to $z \sim 0.25$,<fnref n="6" /> with information on the decline. Both surveys are expected to run for at least 3 yr, and hence we expect a total of at least 63 SLSNe I from them.

Considering the observed number of SLSNe I (Inserra 2019) and the relative fractions of Fast and Slow subgroups determined from the statistical analysis of the SLSN I population, in our simulated data set we envisage a division of 58 per cent Fast, 23 per cent Slow, and 19 per cent with NS. Hence our simulated samples have a total of 868 and 492 useful SLSNe I F+NS for our analysis. We run Monte Carlo simulations with 868 and 492 SLSNe I (LSST+Euclid+ZTF+ATLAS, see Table 5) following the redshift distribution of Fig. 7 (0.02 < z < 3.95). We randomly place them into the relation of Fig. 3 within $3\sigma$ from the best fit ($x_i, y_i$), and not within the $2.2\sigma$ discussed above (see Section 2), to account for increased uncertainties in the identification of SLSN I subclasses at high redshift. We associate random uncertainties in both $x$ and $y$ ($\lambda_{min}^{\text{err}} < x_{min}^{\text{err}} < x_{max}^{\text{err}}$ and $y_{min}^{\text{err}} < y_{min}^{\text{err}} < y_{max}^{\text{err}}$). We also assign a random distance for each redshift bin. Using this set-up, we run our cosmological fitter as previously done (see Section 6) and retrieve $\Omega_M = 0.262_{-0.018}^{+0.020}$, where the uncertainties are only statistical. With a similar size data set to those of current Type Ia cosmology (e.g. Pantheon and JLA), we also retrieved statistical uncertainties of the same order of those achieved using SNe Ia. This is promising since at our current stage it seems that statistical uncertainties are the major contributor to the SLSN I cosmology error budget. Estimating systematic sources of uncertainty that might occur at high redshift, such as dust evolution, is beyond the scope of this study but the occurrence of SLSNe I almost exclusively in low-metallicity environments might suggest of a typically low dust content for the vast majority of these events (Wiseman et al. 2017). With this set-up we also explored the $w_0w_a$ CDM cosmology and found that the FoM of CMB+JLA+SLSNe(ZTF+ATLAS+LSST-DESC+Euclid) and CMB+JLA+SLSNe(ZTF+ATLAS+LSST-generic + Euclid) are 11.5 and 15.5, respectively (see Fig. 8). While that of CMB+JLA is 5.6 suggesting that in the future, SLSNe I will help deliver more precise cosmological constraints (of a factor of 2–3) than this proof of concept. This analysis suggests that with a sample almost as large as that of current Type Ia cosmology (e.g. the Pantheon sample; Scolnic et al. 2018), we can retrieve a

similar statistical precision for $\Omega_M$ and can independently confirm Type Ia findings, as well as reach a redshift range that should be matter dominated but still unexplored with Type Ia SN cosmology.

A possible issue might be selecting 868 (or 492) cosmological useful SLSNe I among the total number of SLSNe I since optical spectroscopy can probe only Fe II line out to $z \sim 1$ and hence identifying their subtype would be challenging. In future, this may be solved with the European Extremely Large Telescope (E-ELT) and the High Angular Resolution Monolithic Optical and Near-infrared Integral field (HARMONI) spectrograph (first light in 2025) that is capable to probe the rest-frame region of interest out to $z \sim 3.7$ (Zieleniewski et al. 2015). However, it has recently been shown that UV absorption lines in SLSNe I Fast are generally sharper than Slow, in contrast to what happens to the absorption lines of less ionized elements in the optical. However, in SLSNe I Fast, all species at $\lambda > 1500$ Å show faster velocities than in Slow events (Gal-Yam 2019a). Such behaviour can also be appreciated in the spectroscopic UV comparison between low-redshift SLSNe I and the high-redshift SLSN I at $z \sim 2$ (Smith et al. 2018), where it has been demonstrated that SLSNe I Fast recede faster in the UV than Slow events, as well as displaying an irregular UV colour evolution. Such results have not yet been confirmed with a larger data set due to the paucity of UV observations. However, they are somewhat reassuring, as they may be used to pave the way to the identification of cosmological useful SLSNe I by means of optical facilities up to the highest redshift of our predictions. Thanks to their characteristic light-curve evolution (Smith et al. 2018; Inserra 2019; Inserra & Parrag, in preparation) and more accurate machine learning techniques (e.g. Ishida et al. 2019; Muthukrishna et al. 2019; Möller & de Boissière 2020), SLSN I identification might be possible without the need of spectroscopy. We note that considering and measuring an uncertainty parameter for high-redshift, photometrically identified SLSNe I is premature and beyond the scope of this work.

## 11 CONCLUSIONS

We examined a sample of 26 SLSNe I, 20 of which are useful for a cosmological analysis. We confirmed the previously established standardization relation of SLSNe I (Inserra & Smartt 2014) with a larger data set and improved light-curve fitting technique, and used the sample to make a measurement of the cosmological parameter $\Omega_M$. The resulting Hubble diagram contains the highest spectroscopically confirmed redshift SN to date ($z \sim 2$). From SLSN I data only, we find $\Omega_M = 0.38_{-0.19}^{+0.24}$ (stat+sys) and an rms = 0.27 mag. We also explored a $w_0w_a$ CDM cosmological model combining our SLSN I sample with the JLA sample and measurements from the CMB, finding that only a small improvement can be made in the constraints on $w_0$ and $w_a$ by 4 per cent in terms of their FoM. We have also simulated future data sets, and demonstrated their potential to reduce the current statistical uncertainties by a factor of 10 on $\Omega_M$, making them comparable to those found using current SN Ia samples. The FoM of the CMB+WMAP+JLA+SLSNe set-up will increase, providing an improvement of a factor 2–3 in the precision of cosmological constraints and also offering a longer redshift baseline. This represents a proof of concept of the current potential and future strengths of SLSN I in cosmology. The key output of this study is that it empowers the investigation of the behaviour of our Universe ($\Omega_M, w_0, w_a$) up to redshifts that cannot be explored using other SNe from the ground ($z > 1.5$).

<footnote n="5">⁵ZTF would need a depth of 22.66 mag to catch all SLSNe I out to $z = 0.5$.</footnote>

<footnote n="6">⁶From the ZTF/ATLAS ATels and Transient Name Server (TNS) AstroNotes we confirm that this number is conservative.</footnote>
<table><caption>**Table 5.** SLSN I future simulated data set out to redshift $z = 3.5$ from two configurations of the LSST deep drilling fields, the *Euclid* satellite and the low-redshift ZTF and ATLAS surveys. Results for the flat $\Lambda$CDM model (i.e. $\Omega_M$) and the $w_0w_a$CDM model figure of merit (FoM) are also reported. Precision on $\Omega_M$, due to statistical uncertainties only, has been increased as a consequence of adopting a larger sample.</caption><thead><tr><th>Surveys</th><th>No. of SLSNe I</th><th>Ω<sub>M</sub></th></tr></thead><tbody><tr><td>ZTF+ATLAS+LSST-generic+Euclid</td><td>868</td><td>0.278<sup>+0.018</sup><sub>-0.017</sub></td></tr><tr><td>ZTF+ATLAS+LSST-DESC+Euclid</td><td>492</td><td>0.269<sup>+0.022</sup><sub>-0.021</sub></td></tr><tr><td>Configuration</td><td colspan="2">Figure of merit (FoM)</td></tr><tr><td>CMB+JLA</td><td colspan="2">5.60</td></tr><tr><td>CMB+JLA+Lyman α</td><td colspan="2">5.90</td></tr><tr><td>CMB+JLA+SLSNe</td><td colspan="2">5.84</td></tr><tr><td>CMB+JLA+SLSNe(ZTF+ATLAS+LSST-generic+Euclid)</td><td colspan="2">15.50</td></tr><tr><td>CMB+JLA+SLSNe(ZTF+ATLAS+LSST-DESC+Euclid)</td><td colspan="2">11.50</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Predicted distribution of SLSNe I as a function of redshift. Bins are $\Delta z = 0.1$ for the LSST deep drilling fields SLSNe, while the *Euclid* rates, binned with a $\Delta z = 0.5$ (Inserra et al. 2018b), have been here resampled with a $\Delta z = 0.1$. A flat SLSN I distribution up to $z \sim 0.3$ from ZTF+ATLAS with a $\Delta z = 0.1$ bin size is also shown (see Section 10).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Constraints on the dark energy equation-of-state parameters $w_0$ and $w_a$ as for Fig. 4 where the Base configuration is given by CMB+JLA. We compare the Base+Lyman $\alpha$ configuration with the two future BASE+SLSNe configurations that both show the improvement in precision of the cosmological analysis (see the FoM in Table 5).</figcaption></figure>

## ACKNOWLEDGEMENTS

This paper has gone through internal review by the DES collaboration. It has Fermilab preprint number 19-115-AE and DES publication number 13387. We acknowledge support from EU/FP7-ERC grant 615929. RCN would like to acknowledge support from STFC grant ST/N000688/1 and the Faculty of Technology at the University of Portsmouth. LG was funded by the European Union's Horizon 2020 Framework Programme 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). Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft, and the Collaborating Institutions in the Dark Energy Survey. The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l'Espai (IEEC/CSIC), the Institut de Física d'Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, the National Optical Astronomy Observatory, the University of Nottingham, The Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, Texas A&M University, and the OzDES Membership Consortium. Based in part on observations at Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA)
under a cooperative agreement with the National Science Foun-
dation. The DES data management system is supported by the
National Science Foundation under grant numbers AST-1138766
and AST-1536171. The DES participants from Spanish institu-
tions are partially supported by MINECO under grants AYA2015-
71825, ESP2015-66861, FPA2015-68048, SEV-2016-0588, SEV-
2016-0597, and MDM-2015-0509, some of which include ERDF
funds from the European Union. IFAE is partially funded by the
CERCA program of the Generalitat de Catalunya. Research leading
to these results has received funding from the European Research
Council under the European Union Seventh Framework Programme
(FP7/2007-2013) including ERC grant agreements 240672, 291329,
and 306478. We acknowledge support from the Australian Research
Council Centre of Excellence for All-sky Astrophysics (CAASTRO),
through project number CE110001020, and the Brazilian Instituto
Nacional de Ciência e Tecnologia (INCT) e-Universe (CNPq grant
465376/2014-2).

This paper has been authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics. The United States Government retains and the publisher, by accepting the paper for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this paper, or allow others to do so, for United States Government purposes.

DATA AVAILABILITY

The data underlying this paper are available in the paper and in its
online supplementary material. The 4OPS code is available at the
following GitHub account https://github.com/cinserra/4OPS. Any
data not presented in the above table or references will be shared on
reasonable request to the corresponding author.

REFERENCES

Abbott T. M. C. et al., 2019, ApJ, 872, L30

Albrecht A. et al., 2006, preprint (arXiv:astro-ph/0609591)

Anderson L. et al., 2014, MNRAS, 441, 24

Angus C. R., Levan A. J., Perley D. A., Tanvir N. R., Lyman J. D., Stanway
E. R., Fruchter A. S., 2016, MNRAS, 458, 84

Angus C. R. et al., 2019, MNRAS, 487, 2215

Arnouts S., Cristiani S., Moscardini L., Matarrese S., Lucchin F., Fontana A.,
Giallongo E., 1999, MNRAS, 310, 540

Barbary K. et al., 2009, ApJ, 690, 1358

Bellm E., 2014, in Wozniak P. R., Graham M. J., Mahabal A. A., Seaman R., eds, The Third Hot-Wiring the Transient Universe Workshop. p. 27 (online at http://www.slac.stanford.edu/econf/C131113.1/)

Betoule M. et al., 2014, A&A, 568, A22

Bishop C., 2006, Pattern Recognition and Machine Learning. Springer-
Verlag, New York

Cardelli J. A., Clayton G. C., Mathis J. S., 1989, ApJ, 345, 245

Chabrier G., 2003, PASP, 115, 763

Chen T.-W., Smartt S. J., Yates R. M., Nicholl M., Krühler T., Schady P.,
Dennefeld M., Inserra C., 2017a, MNRAS, 470, 3566

Chen T.-W. et al., 2017b, ApJ, 849, L4

Chomiuk L. et al., 2011, ApJ, 743, 114

Conley A. et al., 2011, ApJS, 192, 1

Cooke J. et al., 2012, Nature, 491, 228

da Cunha E., Charlot S., Elbaz D., 2008, MNRAS, 388, 1595

Davis T. M. et al., 2011, ApJ, 741, 67

De Cia A. et al., 2018, ApJ, 860, 100

de Sainte Agathe V. et al., 2019, A&A, 629, A85

Diehl H. T. et al., 2014, in Peck A. B., Benn C. R., Seaman R. L., eds, Proc. SPIE Vol. 9149, Observatory Operations: Strategies, Processes, and Systems V. SPIE, Bellingham, p. 91490V

Diehl H. T. et al., 2016, in Peck A. B., Seaman R. L., Benn C. R., eds, Proc. SPIE Vol. 9910, Observatory Operations: Strategies, Processes, and Systems VI. SPIE, Bellingham, p. 99101D

Diehl H. T. et al., 2018, in Peck A. B., Seaman R. L., Benn C. R., eds, Proc. SPIE Vol. 10704, Observatory Operations: Strategies, Processes, and Systems VII. SPIE, Bellingham, p. 107040D

Folatelli G. et al., 2013, ApJ, 773, 53

Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

Gal-Yam A., 2012, Science, 337, 927

Gal-Yam A., 2019a, ARA&A, 57, 305

Gal-Yam A., 2019b, ApJ, 882, 102

Gal-Yam A. et al., 2009, Nature, 462, 624

Gutiérrez C. P. et al., 2017, ApJ, 850, 89

Hinshaw G. et al., 2013, ApJS, 208, 19

Ilbert O. et al., 2006, A&A, 457, 841

Inserra C., 2019, Nat. Astron., 3, 697

Inserra C., Smartt S. J., 2014, ApJ, 796, 87

Inserra C. et al., 2013, ApJ, 770, 128

Inserra C. et al., 2017, MNRAS, 468, 4642

Inserra C. et al., 2018a, MNRAS, 475, 1046

Inserra C. et al., 2018b, A&A, 609, A83

Inserra C., Prajs S., Gutierrez C. P., Angus C., Smith M., Sullivan M., 2018c, ApJ, 854, 175

Ishida E. E. O. et al., 2019, MNRAS, 483, 2

James F., Roos M., 1975, Comput. Phys. Commun., 10, 343

Jha S., Riess A. G., Kirshner R. P., 2007, ApJ, 659, 122

Kelly B. C., 2007, ApJ, 665, 1489

Kim A., Goobar A., Perlmutter S., 1996, PASP, 108, 190

Leloudas G. et al., 2015, MNRAS, 449, 917

Lewis A., 2013, Phys. Rev. D, 87, 103529

Lewis A., Bridle S., 2002, Phys. Rev. D, 66, 103511

Li L.-X., 2008, MNRAS, 388, 1487

LSST Science Collaboration et al., 2009, preprint (arXiv:0912.0201)

Lunnan R. et al., 2014, ApJ, 787, 138

Lunnan R. et al., 2016, ApJ, 831, 144

Lunnan R. et al., 2018, ApJ, 852, 81

Lunnan R. et al., 2020, ApJ, 901, 61

McCrum M. et al., 2014, MNRAS, 437, 656

Mazzali P. A., Sullivan M., Pian E., Greiner J., Kann D. A., 2016, MNRAS, 458, 3455

Möller A., de Boissière T., 2020, MNRAS, 491, 4277

Moriya T. J. et al., 2019, ApJS, 241, 16

Muthukrishna D., Narayan G., Mandel K. S., Biswas R., Hložek R., 2019, PASP, 131, 118002

Nicholl M. et al., 2015a, MNRAS, 452, 3869

Nicholl M. et al., 2015b, ApJ, 807, L18

Nicholl M. et al., 2016, ApJ, 826, 39

Nicholl M., Guillochon J., Berger E., 2017, ApJ, 850, 55

Pastorello A. et al., 2010, ApJ, 724, L16

Percival W. J., Cole S., Eisenstein D. J., Nichol R. C., Peacock J. A., Pope A. C., Szalay A. S., 2007, MNRAS, 381, 1053

Perley D. A. et al., 2016, ApJ, 830, 13

Perlmutter S. et al., 1999, ApJ, 517, 565

Planck Collaboration XIII, 2016, A&A, 594, A13

Prajs S. et al., 2017, MNRAS, 464, 3568

Quimby R. M. et al., 2011, Nature, 474, 487

Quimby R. M., Yuan F., Akerlof C., Wheeler J. C., 2013, MNRAS, 431, 912

Quimby R. M. et al., 2018, ApJ, 855, 2

Rasmussen C. E., Williams C. K. I., 2006, Gaussian Processes for Machine
Learning. The MIT Press, Cambridge, MA

Riess A. G. et al., 1998, AJ, 116, 1009

Riess A. G. et al., 2018, ApJ, 853, 126

Schlafly E. F., Finkbeiner D. P., 2011, ApJ, 737, 103

Schlafly E. F. et al., 2012, ApJ, 756, 158
Schlegel D. J., Finkbeiner D. P., Davis M., 1998, ApJ, 500, 525

Schulze S. et al., 2018, MNRAS, 473, 1258

Scolnic D. M. et al., 2018, ApJ, 859, 101

Scovacricchi D., Nichol R. C., Bacon D., Sullivan M., Prajs S., 2016, MNRAS, 456, 1700

Smartt S. J. et al., 2015, A&A, 579, A40

Smith M. et al., 2016, ApJ, 818, L8

Smith M. et al., 2018, ApJ, 854, 37

Smith K. W. et al., 2020, PASP, 132, 085002

Sullivan M. et al., 2010, MNRAS, 406, 782

Tonry J. L. et al., 2012, ApJ, 750, 99

Tonry J. L. et al., 2018, PASP, 130, 064505

Villar V. A., Nicholl M., Berger E., 2018, ApJ, 869, 166

Vreeswijk P. M. et al., 2014, ApJ, 797, 24

Wiseman P., Schady P., Bolmer J., Krühler T., Yates R. M., Greiner J., Fynbo J. P. U., 2017, A&A, 599, A24

Yan L. et al., 2015, ApJ, 814, 108

Yan L. et al., 2017a, ApJ, 840, 57

Yan L. et al., 2017b, ApJ, 848, 6

Young D. R. et al., 2010, A&A, 512, A70

Zieleniewski S., Thatte N., Kendrew S., Houghton R. C. W., Swinbank A. M., Tecza M., Clarke F., Fusco T., 2015, MNRAS, 453, 3754

APPENDIX A: 4OPS

<figure><img src="image_11.png" /><figcaption>**Figure A1.** A reduced version of the Four Observables Parameter Space (4OPS) for SLSNe I. Data are taken from DES, literature, and other surveys' sample papers that made our first two quality cuts. The left-hand panel shows the magnitude at peak versus the decline rate over 30 d post-peak. The right-hand panel shows the colour at peak versus the colour at 30 d post-peak. The literature objects, both Fast (circles) and Slow (open squares), are shown together to their best fit of the weighted linear regression (dashed, black line) and with the 2.2σ confidence bands defined by the Chauvenet's criterion. We also include DES16C2nm, which was not presented before in the literature sample (Inserra et al. 2018c).</figcaption></figure>
<footnote n="1">¹School of Physics and Astronomy, Cardiff University, Queens Buildings, The Parade, Cardiff CF24 3AA, UK</footnote>

<footnote n="2">²School of Physics and Astronomy, University of Southampton, Southampton SO17 1BJ, UK</footnote>

<footnote n="3">³DARK, Niels Bohr Institute, University of Copenhagen, Lyngbyvej 2, DK-2100 Copenhagen Ø, Denmark</footnote>

<footnote n="4">⁴Department of Physics and Astronomy, University of North Georgia, Dahlonega, GA 30597, USA</footnote>

<footnote n="5">⁵Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth PO1 3FX, UK</footnote>

<footnote n="6">⁶Institut de Physique des Deux Infinis, Université de Lyon 1, CNRS/IN2P3, F-69622 Villeurbanne Cedex, France</footnote>

<footnote n="7">⁷Finnish Centre for Astronomy with ESO (FINCA), FI-20014 University of Turku, Finland</footnote>

<footnote n="8">⁸Tuorla Observatory, Department of Physics and Astronomy, FI-20014 University of Turku, Finland</footnote>

<footnote n="9">⁹Universite Clermont Auvergne, CNRS/IN2P3, LPC, F-63000 Clermont-Ferrand, France</footnote>

<footnote n="10">¹⁰Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA</footnote>

<footnote n="11">¹¹George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, and Department of Physics and Astronomy, Texas A&M University, College Station, TX 77843, USA</footnote>

<footnote n="12">¹²School of Mathematics and Physics, University of Queensland, Brisbane, QLD 4072, Australia</footnote>

<footnote n="13">¹³Departamento de Física Teórica y del Cosmos, Universidad de Granada, E-18071 Granada, Spain</footnote>

<footnote n="14">¹⁴Department of Astronomy and Astrophysics, University of Chicago, Chicago, IL 60637, USA</footnote>

<footnote n="15">¹⁵Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA</footnote>

<footnote n="16">¹⁶Lawrence Berkeley National Lab, 1 Cyclotron Rd., Berkeley, CA 94720, USA</footnote>

<footnote n="17">¹⁷Division of Science, National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="18">¹⁸Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory, Casilla 603, La Serena, Chile</footnote>

<footnote n="19">¹⁹Fermi National Accelerator Laboratory, PO Box 500, Batavia, IL 60510, USA</footnote>

<footnote n="20">²⁰Instituto de Física Teórica UAM/CSIC, Universidad Autónoma de Madrid, E-28049 Madrid, Spain</footnote>

<footnote n="21">²¹Institut d'Astrophysique de Paris, CNRS, UMR 7095, F-75014 Paris, France</footnote>

<footnote n="22">²²Institut d'Astrophysique de Paris, Sorbonne Universités, UPMC Univ Paris 06, UMR 7095, F-75014 Paris, France</footnote>

<footnote n="23">²³Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="24">²⁴Kavli Institute for Particle Astrophysics and Cosmology, PO Box 2450, Stanford University, Stanford, CA 94305, USA</footnote>

<footnote n="25">²⁵SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA</footnote>

<footnote n="26">²⁶Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Madrid 28040, Spain</footnote>

<footnote n="27">²⁷Laboratório Interinstitucional de e-Astronomia - LIneA, Rua Gal. José Cristino 77, 20921-400 Rio de Janeiro, RJ, Brazil</footnote>

<footnote n="28">²⁸Department of Astronomy, University of Illinois at Urbana-Champaign, 1002 W. Green Street, Urbana, IL 61801, USA</footnote>

<footnote n="29">²⁹National Center for Supercomputing Applications, 1205 West Clark St, Urbana, IL 61801, USA</footnote>

<footnote n="30">³⁰Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, E-08193 Bellaterra (Barcelona), Spain</footnote>

<footnote n="31">³¹Institut d'Estudis Espacials de Catalunya (IEEC), E-08034 Barcelona, Spain</footnote>

<footnote n="32">³²Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s/n, E-08193 Barcelona, Spain</footnote>

<footnote n="33">³³Physics Department, University of Wisconsin-Madison, 2320 Chamberlin Hall, 1150 University Avenue, Madison, WI 53706-1390, USA</footnote>

<footnote n="34">³⁴Department of Physics, IIT Hyderabad, Kandi, Telangana 502285, India</footnote>

<footnote n="35">³⁵Department of Astronomy/Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721-0065, USA</footnote>

<footnote n="36">³⁶Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Dr., Pasadena, CA 91109, USA</footnote>

<footnote n="37">³⁷Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA</footnote>

<footnote n="38">³⁸Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA</footnote>

<footnote n="39">³⁹Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="40">⁴⁰Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="41">⁴¹Department of Physics, Stanford University, 382 Via Pueblo Mall, Stanford, CA 94305, USA</footnote>

<footnote n="42">⁴²Observatório Nacional, Rua Gal. José Cristino 77, 20921-400 Rio de Janeiro, RJ, Brazil</footnote>

<footnote n="43">⁴³Santa Cruz Institute for Particle Physics, Santa Cruz, CA 95064, USA</footnote>

<footnote n="44">⁴⁴Center for Cosmology and Astro-Particle Physics, The Ohio State University, Columbus, OH 43210, USA</footnote>

<footnote n="45">⁴⁵ASTRAVEO LLC, PO Box 1668, Gloucester, MA 01931, USA</footnote>

<footnote n="46">⁴⁶Australian Astronomical Optics, Macquarie University, North Ryde, NSW 2113, Australia</footnote>

<footnote n="47">⁴⁷The Research School of Astronomy and Astrophysics, Australian National University, Canberra, ACT 2611, Australia</footnote>

<footnote n="48">⁴⁸Departamento de Física Matemática, Instituto de Física, Universidade de São Paulo, CP 66318, 05314-970 São Paulo, SP, Brazil</footnote>

<footnote n="49">⁴⁹Institució Catalana de Recerca i Estudis Avançats, E-08010 Barcelona, Spain</footnote>

<footnote n="50">⁵⁰Department of Astrophysical Sciences, Princeton University, Peyton Hall, Princeton, NJ 08544, USA</footnote>

<footnote n="51">⁵¹Department of Physics and Astronomy, Pevensey Building, University of Sussex, Brighton BN1 9QH, UK</footnote>

<footnote n="52">⁵²Physics Department, Brandeis University, 415 South Street, Waltham, MA 02453, USA</footnote>

<footnote n="53">⁵³Instituto de Física Gleb Wataghin, Universidade Estadual de Campinas, 13083-859 Campinas, SP, Brazil</footnote>

<footnote n="54">⁵⁴Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA</footnote>

<footnote n="55">⁵⁵Argonne National Laboratory, 9700 South Cass Avenue, Lemont, IL 60439, USA</footnote>

<footnote n="56">⁵⁶Korea Astronomy and Space Science Institute, Yuseong-gu, Daejeon 305-348, Korea</footnote>

<footnote n="57">⁵⁷School of Mathematics and Physics, University of Queensland, Brisbane, QLD 4072, Australia</footnote>

<footnote n="58">⁵⁸INAF, Astrophysical Observatory of Turin, I-10025 Pino Torinese, Italy</footnote>

<footnote n="59">⁵⁹Centre for Astrophysics and Supercomputing, Swinburne University of Technology, Melbourne, VIC 3122, Australia</footnote>

<footnote n="60">⁶⁰Sydney Institute for Astronomy, School of Physics A28, The University of Sydney, Sydney, NSW 2006, Australia</footnote>

<footnote>This paper has been typeset from a TEX/LATEX file prepared by the author.</footnote>