The first Hubble diagram and cosmological constraints using superluminous supernovae

We present the ﬁrst Hubble diagram of superluminous supernovae (SLSNe) out to a redshift of two, together with constraints on the matter density, (cid:2) M , and the dark energy equation-of-state parameter, w ( ≡ p / ρ ). We build a sample of 20 cosmologically useful SLSNe I based on light curve and spectroscopy quality cuts. We conﬁrm the robustness of the peak–decline SLSN I standardization relation with a larger data set and improved ﬁtting techniques than previous works. We then solve the SLSN model based on the above standardization via minimization of the χ 2 computed from a covariance matrix that includes statistical and systematic uncertainties. For a spatially ﬂat (cid:5) cold dark matter ( (cid:5) CDM) cosmological model, we ﬁnd (cid:2) M = 0 . 38 + 0 . 24 − 0 . 19 , with an rms of 0.27 mag for the residuals of the distance moduli. For a w 0 w 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 conﬁguration used) and show that such a sample can deliver cosmological constraints in a ﬂat (cid:5) 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.

density of the Universe at the present epoch.SNe Ia present a direct and mature method of probing this dark energy via its equationof-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 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 ∼ 4 (Cooke et al. 2012) and spectroscopically out to z ∼ 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 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.

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 lightcurve 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 λ5169 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.

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(Diehl et al. , 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 M(400) 30 ) to avoid introducing any biases over the fact that the peak-decline relation (M(400) 0 versus 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 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 ∼ 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.

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 ( 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).

S U P E R N OVA L I G H T C U RV E S
To estimate and model the light curves of SLSNe I around peak (−15 ≤ phase (d) ≤ 30), we first k-correct the observed magnitudes to the 400 and 520 nm synthetic filters with the SNAKE1 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 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 → 400 is the kcorrection 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 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).
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.derive the uncertainties (see Inserra et al. 2018c, for a more in-depth description). 2

S L S N I S TA N DA R D I Z AT I O N
We next confirm that the previous observed relationships between peak luminosity (M(400) 0 ) and decline rate in magnitudes over 30 d ( 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 , matter = 0.27, = 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 n i=1 1/σ i , where n is the sample size.For the peakdecline relation and a sample of 20 objects we retrieve a standard deviation similar to that of the previous unweighted study (σ = 0.33; 2 https://github.com/cinserra/Gaussian-Processes-GP-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 σ = 0.33 mag decreases, as expected, to σ = 0.26 mag if we use only the Fast subclass.Including the Slow subclass events substantially increases the dispersion to σ = 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 (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 χ 2 ( χ 2 = 0.90) and σ = 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 ) × 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 one ( 100), these two relations including the 520 nm band might be more effective than the peak-decline.

S U P E R N OVA M O D E L
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, μ obs , is then given by where m(400) is the peak apparent magnitude in rest-frame 400 nm band, and M(400) (the peak absolute magnitude) and γ are nuisance parameters in the distance estimate.This is compared to the model distance modulus, μ model , of 5 log 10 (d L /10 pc), where d L is the luminosity distance.The fit then minimizes the χ 2 according to where C is the covariance matrix and μ is the vector of residuals μ = μ obs − μ 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.M + = 1, and the free parameters in the fit are therefore M and the two nuisance parameters M(400) and γ .
The covariance matrix is defined as the sum of the statistical and systematic parts: (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 M ; the M( 400) and γ values previously retrieved (Inserra & Smartt 2014) for the M(400) 0 versus M(400) 30 as a prior of the other two nuisance parameters (M(400), γ ).We also assume a Gaussian distribution for the form of the priors.To minimize equation (3) we use IMINUIT, 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 × 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).

Statistical uncertainties
The statistical uncertainties part is diagonal and includes uncertainties as follows: Here σ m 400 ,i and σ M(400) 30 are the uncertainties on the fitted lightcurve 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, σ 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 σ lensing ∼ 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.

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 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 sys,i,j = σ 2 reddening,i,j + σ 2 ZP,i,j + σ 2 MalmquistBias,i,j + σ 2 model,i,j . (6)

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.

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 (Schlafly 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 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).

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 ≤ z ≤ 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).

Light-curve fitting
The last term, σ 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 stat,ii .

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.

S L S N E I H U B B L E D I AG R A M
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 CDM), SLSNe I alone can provide a measurement of the reduced matter density M .The SLSN I Hubble diagram and the flat 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 M = 0.38 +0.24 −0.19 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 CDM model.
In Fig. 3, we plot the residuals of our sample without the declinerate correction, as a function of the decline rate ( 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.

w 0 w a C D M C O S M O L O G Y W I T H S L S N E 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 ∝ χ 2 , where the χ 2 is given by equation ( 3), and μ model is now a function of w 0 , w a , 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±0.185 and w a = −0.594±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/(σ w 0 σ wa )).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 α 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 α is 5.9, suggesting that at present SLSNe are comparable to Lyman α BAO at such redshifts.That is also true of we analyse the uncertainties of the three data combinations (Base, Base+SLSNe, Base+Lyman α) on w 0 and w a that are 0.187, 0.185, and 0.183 (σ w 0 ) and 0.952, 0.926, and 0.928 (σ wa ).

E N V I RO N M E N T
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 * ), 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  400) 30 standardization method.Overplotted (solid red line) is the best-fitting flat CDM cosmology and its uncertainties (shaded red area) with M = 0.38 +0.24  −0.19 , 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.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.

E X P L O R I N G S L S N E I I N T H E U LT R AV I O L E T
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  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.

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 α 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 and w a = −0.594±0.926 , 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 α configuration.
with light-curve evolution around peak (−10 < 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 ࣠ z ࣠ 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 λ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 < phase < +10 d).We also measure the line velocities and find 18 000 < v (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 = pEW(Mg II/C II)/pEW(C III/C II/Ti III)) with the light-curve decline ( 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 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 with three nuisance parameters (α, β, 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.

0 S L S N E I C O S M O L O G Y: F U T U R E A N D I M P ROV E M E N T S
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;   SLSN I luminosity of M(400) 0 = −21.756± 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σ depth of 20.55 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 ∼ 0.25,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σ from the best fit (x i , y i ), and not within the 2.2σ 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 (x err min < x err i < x err max and y err min < y err i < y err max ).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 M = 0.262 +0.020 −0.018 , 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 setup we also explored the w 0 w 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 Patheon sample; Scolnic et al. 2018), we can retrieve a similar statistical precision for 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 ∼ 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 ∼ 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 λ > 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 ∼ 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.

C O N C L U S I O N S
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 M .The resulting Hubble diagram contains the highest spectroscopically confirmed redshift SN to date (z ∼ 2).From SLSN I data only, we find M = 0.38 +0.24  −0.19 (stat+sys) and an rms = 0.27 mag.We also explored a w 0 w 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 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 ( M , w 0 , w a ) up to redshifts that cannot be explored using other SNe from the ground (z > 1.5).
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 CDM model (i.e.M ) and the w 0 w a CDM model figure of merit (FoM) are also reported.Precision on M , due to statistical uncertainties only, has been increased as a consequence of adopting a larger sample.

Surveys
No 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 α 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).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).

Figure 1 .
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 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σ 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).
Light curve and spectra (quality cuts 1 and 2) 4OPS (quality cut 3

Figure 2 .
Figure 2. Upper panel: the Hubble diagram of our SLSN I sample (only F+NS subtypes) using the M(400) 30 standardization method.Overplotted (solid red line) is the best-fitting flat CDM cosmology and its uncertainties (shaded red area) with M = 0.38+0.24−0.19 , 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.

Figure 3 .
Figure 3. Left: the SLSN I residuals uncorrected for the decline parameter, plotted as a function of 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.

Figure 5 .
Figure 5. SLSNe I residuals from the best-fitting flat 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.

Figure 6 .
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, 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.

Figure 7 .
Figure 7. Predicted distribution of SLSNe I as a function of redshift.Bins are z = 0.1 for the LSST deep drilling fields SLSNe, while the Euclid rates, binned with a z = 0.5 (Inserra et al. 2018b), have been here resampled with a z = 0.1.A flat SLSN I distribution up to z ∼ 0.3 from ZTF+ATLAS with a z = 0.1 bin size is also shown (see Section 10).

Figure A1 .
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).

Table 2 .
SLSNe I and their peak-decline relation values, together with the quantities used in the third criterion approach.

Table 3 .
Fit parameters and statistical results of our sample.
Note.Least-squares fits for a Bayesian weighted linear regression with weighted errors both in x and y of the form η = β + α × x + , where x = x + x err and y = η + y err .The σ is the standard deviation of this fit.The last column gives the reduced χ 2 .

Table 4 .
Best-fitting parameters for the flat CDM model using SLSNe I alone and using the M(400) 0 versus M(400) 30 standardization relation.
Prajs et al. (2017)b), 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σ point source; LSST ScienceCollaboration 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 ofPrajs et al. (2017), we use an average peak 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, Fundac ¸ã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 Inovac ¸ã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)