The first Hubble diagram and cosmological constraints using superluminous supernovae

MNRAS 504, 2535–2549 (2021) doi:10.1093/mnras/stab978 Advance Access publication 2021 April 09 The first Hubble diagram and cosmological constraints using superluminous supernovae C. Inserra ,1‹M. Sullivan ,2C. R. Angus ,3E. Macaulay,4R. C. Nichol,5M. Smith ,2,6 C. Frohmaier ,5C. P. Guti´ errez ,2,7,8M. Vicenzi ,5A. M¨ oller,9D. Brout ,10 P. J. Brown,11 T. M. Davis ,12 C. B. D’Andrea,10 L. Galbany ,13 R. Kessler,14,15 A. G. Kim ,16 Y.-C. Pan,17 M. Pursiainen ,2D. Scolnic ,15 B. P. Thomas,5P. Wiseman ,2T. M. C. Abbott,18 J. Annis,19 S. Avila ,20 E. Bertin ,21,22 D. Brooks,23 D. L. Burke,24,25 A. Carnero Rosell ,26,27 M. Carrasco Kind ,28,29 J. Carretero,30 F. J. Castander,31,32 R. Cawthon ,33 S. Desai,34 H. T. Diehl,19 T. F. Eifler ,35,36 D. A. Finley,19 B. Flaugher,19 P. Fosalba ,31,32 J. Frieman,15,19 J. Garcia-Bellido ,20 E. Gaztanaga ,31,32 D. W. Gerdes,37,38 T. Giannantonio,39,40 D. Gruen ,24,25,41 R. A. Gruendl,28,29 J. Gschwend,27,42 G. Gutierrez ,19 D. L. Hollowood,43 K. Honscheid,44 D. J. James,45 E. Krause,35 K. Kuehn,46 N. Kuropatkin,19 T. S. Li ,14,19 C. Lidman ,47 M. Lima,27,48 M. A. G. Maia,27,42 J. L. Marshall,10 P. Martini,44 F. Menanteau,28,29 R. Miquel,30,49 A. A. Plazas Malag´ on ,50 A. K. Romer,51 A. Roodman,24,25 M. Sako,10 E. Sanchez,26 V. Scarpine,19 M. Schubnell,38 S. Serrano,31,32 I. Sevilla-Noarbe,26 M. Soares-Santos ,52 F. Sobreira,26,53 E. Suchyta ,54 M. E. C. Swanson,29 G. Tarle,38 D. Thomas ,5D. L. Tucker,19 V. Vikram,55 A. R. Walker,18 Y. Zhang,19 J. Asorey,56 J. Calcino,57 D. Carollo,58 K. Glazebrook,59 S. R. Hinton ,57 J. K. Hoormann,57 G. F. Lewis ,60 R. Sharp,47 E. Swann,2andB.E.Tucker 47 (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, 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 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 χ2computed from a covariance matrix that includes statistical and systematic uncertainties. For a spatially flat cold dark matter (CDM) cosmological model, we find M=0.38+0.24 −0.19, with an rms of 0.27 mag for the residuals of the distance moduli. For a w0waCDM 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 w0and waof 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 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, w0and wa. 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 E-mail: [email protected] (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 dark energy is the major component (≈70 per cent) of the energy C 2021 The Author(s) Published by Oxford University Press on behalf of Royal Astronomical Society Downloaded from https://academic.oup.com/mnras/article/504/2/2535/6219089 by Universidad de Granada - Biblioteca user on 30 June 2021 2536 C. Inserra et al. 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. 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 lightcurveevolutionsimilartoSN2010gx(Pastorelloetal.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 FeII λ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. 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 M(400)30) to avoid introducing any biases over the fact that the peak–decline relation (M(400)0versus 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)0rather than M(400)20 as y-axis of the left-hand panel can be found in Appendix A. MNRAS 504, 2535–2549 (2021) Downloaded from https://academic.oup.com/mnras/article/504/2/2535/6219089 by Universidad de Granada - Biblioteca user on 30 June 2021 SLSN I cosmology 2537 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 ∼3mag 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 FeII 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 FeII 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 thosewithNSinouranalysis,forafinalsampleof20SLSNeIintotal. 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 (MUV <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 OII lines (12 000 < T(K) <16000; 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 ≤phase(d) ≤30), we first k-correct the observed magnitudes to the 400 and 520 nm synthetic filters with the SNAKE1software 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)−kX→400,(1) where m(X) is the observed apparent magnitude in passband Xand 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). kX→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;Nicholletal.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)tofitthe light curves, using the PYTHON package GEORGE and a Matern-3/2 kernel to perform our GP regression of SLSN I light curves and 1https://github.com/cinserra/S3 MNRAS 504, 2535–2549 (2021) Downloaded from https://academic.oup.com/mnras/article/504/2/2535/6219089 by Universidad de Granada - Biblioteca user on 30 June 2021 2538 C. Inserra et al. 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)0versus 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). 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. Sample source SLSN I candidates Light curve and spectra (quality cuts 1 and 2) 4OPS (quality cut 3) Reference Literature (4 DES – 6 PS1 – 11 PTF/iPTF) 40 23 20 Inserra et al. (2018c) – PS1 (Medium Deep Survey) 9 3 2 Lunnan et al. (2018) – PTF/iPTF 15 0 0 De Cia et al. (2018) DES 17 9 3 Angus et al. (2019) derive the uncertainties (see Inserra et al. 2018c, for a more in-depth description).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 (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 (H0=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 xand yvariables 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 105iterations. A weighted regression provides a standard deviation bigger than the unweighted one by a factor of roughly n i=11/σi,wherenis 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 (σ=0.33; 2https://github.com/cinserra/Gaussian-Processes-GPInserra & 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.26magifwe use onlytheFastsubclass.IncludingtheSlowsubclass 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)0versus (M(400)30 −M(520)30) relation (Inserra &Smartt2014). 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)0versus 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)0versus (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 MNRAS 504, 2535–2549 (2021) Downloaded from https://academic.oup.com/mnras/article/504/2/2535/6219089 by Universidad de Granada - Biblioteca user on 30 June 2021 SLSN I cosmology 2539 Table 2. SLSNe I and their peak–decline relation values, together with the quantities used in the third criterion approach. ID zSLSN I type M(400)0M(400)30 M(400)20 M(400)0−M(520)0M(400)30 −M(520)30 Gaia16apd 0.102 F −21.87 (0.04) 0.69 (0.06) −21.18 (0.04) −0.18 (0.07) 0.28 (0.07) SN2011ke 0.143 F −21.23 (0.09) 0.89 (0.09) −20.34 (0.02) 0.04 (0.13) 0.59 (0.03) SN2012il 0.175 F −21.54 (0.10) 1.39 (0.17) −20.15 (0.14) −0.02 (0.11) 0.48 (0.13) PTF11rks 0.190 F −20.61 (0.05) 0.87 (0.07) −19.74 (0.05) 0.20 (0.06) 1.16 (0.15) SN2010gx 0.230 F −21.73 (0.02) 0.76 (0.03) −20.97 (0.02) −0.11 (0.02) 0.53 (0.03) SN2011kf 0.245 F −21.74 (0.15) 0.52 (0.18) −21.22 (0.02) ... ... LSQ12dlf 0.255 F −21.52 (0.03) 0.76 (0.04) −20.76 (0.03) 0.05 (0.03) 0.57 (0.10) LSQ14mo 0.256 F −21.04 (0.05) 1.30 (0.14) −19.74 (0.13) −0.08 (0.04) 0.61 (0.02) PTF09cnd 0.258 F −22.16 (0.08) 0.71 (0.14) −21.45 (0.12) ... ... SN2013dg 0.265 F −21.35 (0.05) 1.03 (0.06) −20.32 (0.03) −0.26 (0.08) 0.56 (0.10) PS1-10bzj 0.650 F −21.03 (0.06) 1.23 (0.32) −19.08 (0.31) 0.15 (0.11) 0.94 (0.25) iPTF13ajg 0.740 F −22.42 (0.07) 0.19 (0.10) −22.23 (0.07) −0.29 (0.09) −0.11 (0.09) DES15X3hm 0.860 F −21.94 (0.06) 1.44 (0.07) −21.19 (0.04) −0.21 (0.04) 0.11 (0.06) DES17X1amf 0.920 NS −21.97 (0.07) 0.26 (0.15) −21.78 (0.09) −0.31 (0.10) 0.08 (0.10) PS1-10ky 0.956 F −22.05 (0.06) 0.61 (0.07) −21.44 (0.04) −0.06 (0.07) 0.25 (0.06) PS1-11aib 0.997 NS −22.05 (0.07) 0.31 (0.17) −21.78 (0.40) −0.45 (0.07) −0.15 (0.07) SCP-06F6 1.189 F −22.19 (0.03) 0.57 (0.15) −21.62 (0.15) ... ... PS1-11tt 1.283 NS −21.89 (0.16) 0.15 (0.20) −21.83 (0.12) ... ... PS1-11bam 1.565 NS −22.45 (0.10) 0.36 (0.14) −22.09 (0.10) ... ... DES16C2nm 1.998 NS −22.52 (0.10) 0.67 (0.11) −22.13 (0.09) ... ... Table 3. Fit parameters and statistical results of our sample. xySLSN I type N(objects) βασ  χ2 M(400)30 M(400)0F15−23.09 ±0.28 1.01 ±0.17 0.26 ±0.22 1.23 F+NS 20 −22.62 ±0.19 0.72 ±0.14 0.33 ±0.22 1.73 F+NS+S25−22.31 ±0.15 0.53 ±0.13 0.36 ±0.22 2.31 (M(400)30 −M(520)30)M(400)0F12−22.76 ±0.35 2.18 ±0.62 0.29 ±0.30 1.45 F+NS 14 −22.79 ±0.30 2.22 ±0.55 0.25 ±0.25 1.20 F+NS+S18−22.31 ±0.20 1.55 ±0.41 0.32 ±0.24 1.91 M(400)20 −M(520)20 M(400)0F12−21.97 ±0.14 1.97 ±0.47 0.31 ±0.27 1.58 F+NS 14 −21.96 ±0.12 1.89 ±0.41 0.28 ±0.24 1.41 F+NS+S18−21.77 ±0.09 1.32 ±0.33 0.33 ±0.24 1.74 M(400)30 −M(520)30 M(400)0F12−22.27 ±0.14 1.51 ±0.27 0.21 ±0.19 0.87 F+NS 14 −22.21 ±0.11 1.41 ±0.23 0.19 ±0.17 0.95 F+NS+S18−21.95 ±0.09 1.02 ±0.21 0.28 ±0.20 1.52 M(400)30 M(400)20 F15−23.10 ±0.31 1.55 ±0.20 0.27 ±0.23 1.82 F+NS 20 −22.67 ±0.20 1.29 ±0.15 0.31 ±0.23 1.89 F+NS+S25−22.31 ±0.16 1.07 ±0.13 0.35 ±0.23 2.34 M(400)30 −M(520)30 ×M(400)30 M(400)0F12−22.21 ±0.12 0.82 ±0.14 0.19 ±0.19 0.86 F+NS 14 −22.12 ±0.09 0.73 ±0.11 0.18 ±0.16 0.61 F+NS+S18−21.95 ±0.08 0.60 ±0.11 0.25 ±0.18 1.03 Note. Least-squares fits for a Bayesian weighted linear regression with weighted errors both in xand yof the form η=β+α×x+,wherex=x+xerr and y=η+yerr.Theσis the standard deviation of this fit. The last column gives the reduced χ2. one (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, μobs, is then given by μobs =m(400) −M(400) +γM(400)30,(2) 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 log10(dL/10 pc), where dLis the luminosity distance. The fit then minimizes the χ2according to χ2=μT·C−1·μ,(3) where Cis the covariance matrix and μis the vector of residuals μ=μobs −μmodel. Note that the Hubble constant, H0,entersin both M(400) and dL, 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+ MNRAS 504, 2535–2549 (2021) Downloaded from https://academic.oup.com/mnras/article/504/2/2535/6219089 by Universidad de Granada - Biblioteca user on 30 June 2021 2540 C. Inserra et al. =1, and the free parameters in the fit are therefore Mand the two nuisance parameters M(400) and γ. The covariance matrix is defined as the sum of the statistical and systematic parts: C=Cstat +Csys.(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:theJointLight-CurveAnalysis(JLA)result(SNanalysisonly) as a prior on M;theM(400) and γvalues previously retrieved (Inserra & Smartt 2014) for the M(400)0versus 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,3a minimization technique based on MINUIT (James & Roos 1975) that runs over 105iterations, and a Markov chain Monte Carlo Ensemble sampler (ForemanMackey et al. 2013) with 5 ×103iterations. Both provide nonGaussian 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: Cstat,ii =σ2 m400,i +γ2σ2 M(400)30 +5(1+zi) zi(1 +zi/2)log(10)2 ×σ2 z,i +σ2 lensing.(5) Here σm400,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. 3https://github.com/iminuit/iminuit 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 Cstat,ii)andthen 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 Csys,i,j =σ2 reddening,i,j +σ2 ZP,i,j +σ2 MalmquistBias,i,j +σ2 model,i,j .(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 fromdustcolumndensity 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 RV, hence the assumption of a Galactic value of RV=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;Chenetal. 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−4mag (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 average4zero-point uncertainties for each filter matching those reported by other large surveys (Tonry et al. 2012; Diehl et al. 2014;Smarttetal.2015). 4Average evaluated from the uncertainties of DES and PS1. MNRAS 504, 2535–2549 (2021) Downloaded from https://academic.oup.com/mnras/article/504/2/2535/6219089 by Universidad de Granada - Biblioteca user on 30 June 2021 SLSN I cosmology 2541 5.2.3 Malmquist bias To model our search efficiency (i.e. how many SLSNe I are missed), for each SLSN I we simulated 104light curves using a Monte Carlo approach. Because of the relatively small sample of our analysis we treatit as a simple smoothedoffsetbetweennearby (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). 5.2.4 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 Cstat,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 Bband. Hence, we do not consider any additional source of uncertainties linked to galaxy evolution. 6 SLSNE 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 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. 7w0waCDM COSMOLOGY WITH SLSNE 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)=w0+(1 −a)wa,where a=(1 +z)−1is 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 w0waCDM cosmological model, and a SLSN I likelihood of the form lnL∝χ2, where the χ2is given by equation (3), and μmodel is now a function of w0,wa,M,andH0. To obtain convergence we use the CosmoMC tool (Lewis & Bridle 2002;Lewis2013) We find marginalized constraints on w0=−0.904±0.185 and wa=−0.594±0.926, which are shown in Fig. 4. To measure the improvement on the constraints on w0and wa, 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 w0–waplane (1/(σw0σwa)). Without SLSNe I, the FoMis5.622, which increases by 4 per cent to 5.835 withthe 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 α)onw0and wa that are 0.187, 0.185, and 0.183 (σw0) and 0.952, 0.926, and 0.928 (σwa). 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∗), 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 MNRAS 504, 2535–2549 (2021) Downloaded from https://academic.oup.com/mnras/article/504/2/2535/6219089 by Universidad de Granada - Biblioteca user on 30 June 2021 2542 C. Inserra et al. 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. Table 4. Best-fitting parameters for the flat CDM model using SLSNe I alone and using the M(400)0versus M(400)30 standardization relation. Sample MM(value) γResiduals (rms) =0 residuals (rms) Redshift χ2/degrees of freedom F+NS (stat+sys) 0.38+0.24 −0.19 −22.40+0.38 −0.35 0.62+0.19 −0.19 0.27 0.31 ∼2.0 11.1/18 F+NS (stat) 0.37+0.31 −0.19 −22.45+0.38 −0.35 0.63+0.19 −0.19 0.25 0.27 ∼2.0 11.0/18 F (stat+sys) 0.22+0.29 −0.21 −23.20+0.56 −0.53 0.97+0.25 −0.25 0.26 0.40 ∼1.2 6.4/13 F (stat) 0.41+0.22 −0.20 −22.83+0.43 −0.41 0.88+0.22 −0.22 0.21 0.26 ∼1.2 4.9/13 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 observemildcorrelations betweenthe hostgalaxy sSFR and the light-curve properties in terms of decline (Pearson r=0.45) andcolour(Pearson r=0.52).This suggeststhat SLSN Iin 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 MgII/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 MNRAS 504, 2535–2549 (2021) Downloaded from https://academic.oup.com/mnras/article/504/2/2535/6219089 by Universidad de Granada - Biblioteca user on 30 June 2021 SLSN I cosmology 2543 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 4. Constraints on the dark energy equation-of-state parameters w0 and wa. 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 includingSLSNeI(dark reddot–dashed linesand circlecentroid), orLyman α BAO (orange dotted lines and triangle centroid). The horizontal and vertical dashed lines at w0=−1andwa=0 correspond to a cosmological constant. The Base+SLSNe I results are w0=−0.904±0.185 and wa=−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 FeII 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 MgII/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 CIII/C II/Ti III line is also a good proxy, when available, and it is usually stronger than MgII/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;McCrumetal.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 FeII λ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 preand post-peak spectra. We measure the pEW (see Guti´ errez 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 18000 <v(kms−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(MgII/C II)/pEW(C III/C II/TiIII)) with the light-curve decline (M(400)30). In Fig. 6no clear groups are observed comparing the pEW(MgII/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(panelB)orwiththepEW(Mg II/CII) 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(MgII/C II) versus pEW(CIII/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 μ=m400 −M(400) +α ×(M(400)30 +β×pEW(Mg II/CII)),(7) 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. 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σ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 MNRAS 504, 2535–2549 (2021) Downloaded from https://academic.oup.com/mnras/article/504/2/2535/6219089 by Universidad de Granada - Biblioteca user on 30 June 2021