Detailed spectrophotometric analysis of the superluminous and fast evolving SN

SN 2019neq was a very fast evolving superluminous supernova. At a redshift z = 0.1059, its peak absolute magnitude was − 21.5 ± 0.2 mag in g band. In this work, we present data and analysis from an e xtensiv e spectrophotometric follow-up campaign using multiple observational facilities. Thanks to a nebular spectrum of SN 2019neq, we investigated some of the properties of the host galaxy at the location of SN 2019neq and found that its metallicity and speciﬁc star formation rate are in a good agreement with those usually measured for SLSNe-I hosts. We then discuss the plausibility of the magnetar and the circumstellar interaction scenarios to explain the observed light curves


INTRODUCTION
It is widely accepted that supernovae (SNe) are a possible final stage of the life of massive stars and white dwarfs in close binary systems; they are observationally classified as type-I (SNe I) or type-II (SNe II) depending on whether they are hydrogen-poor or hydrogen-rich, respectively.SNe are usually discovered by untargeted wide-field surveys and they can be identified in their local environment performing differential photometry (see also Sec. 2 and Kessler et al. 2015): in fact, their light curves (LCs) have a magnitude at maximum spanning the range −14 to −19 mag in the optical bands.SN light curves are normally interpreted as being powered by 56 Ni radioactive decay (e. g.Nadyozhin 1994) and by the energy deposited in the ejecta by the shock break-out.However, observations of the so-called superluminous SNe (SLSNe) have complicated this picture, as their absolute magnitude at maximum can be even brighter than −21 mag in optical bands (e.g.Gal-Yam 2012, 2019).Similar to their lowerluminosity siblings, these events are grouped in SLSNe I and SLSNe II based on the strength of their H features, however their luminosity cannot be comfortably explained by the energy released from 56 Ni decay if the classical neutrino-wind driven core-collapse scenario is assumed.In fact, the explosion of standard SN progenitors synthesizes ∼ 0.1 M ⊙ , while SLSNe require ≳ 1 − 10 M ⊙ of 56 Ni (e. g.Umeda & Nomoto 2008;Gal-Yam et al. 2009;Kasen et al. 2011;Dessart et al. 2012).In principle, the explosion of a pair-instability SN (e. g.Yoshida et al. 2016) is a channel to synthesize much more 56 Ni than standard core-collapse SNe, but this scenario is disfavoured by the observed spectra of SLSNe (see e.g.Kozyreva & Blinnikov 2015;Moriya et al. 2019;Mazzali et al. 2019).Two main alternatives have been considered to reproduce the observational properties of SLSNe.The first requires a shock where the SN ejecta collides with circumstellar material (CSM) lost by the progenitor star before its explosion (Woosley et al. 2007;Sorokina et al. 2016;Tolstov et al. 2017;Woosley 2017).In this scenario, the shock driven by their SN ejecta converts the kinetic energy in radiation via collisional excitation and ionisation processes.Alternatively, the deposition in the ejecta of the spindown radiation of a newly-born, highly-magnetized neutron star (a magnetar, e.g.Woosley et al. 2007;Kasen & Bildsten 2010) can boost SLSNe luminosities.However, the magnetar model does not naturally predict the bumps often seen in SLSNe-I LCs both before and after maximum (see e.g.Nicholl et al. 2016a;Smith et al. 2016;Vreeswĳk et al. 2017;Nicholl et al. 2017;Lunnan et al. 2020;Fiore et al. 2021;Gutiérrez et al. 2022;Chen et al. 2023a,b;West et al. 2023;Lin et al. 2023, but see Moriya et al. 2022) but is able to account for the huge variety of LC evolutionary timescales shown by SLSNe I (Inserra et al. 2013;Chatzopoulos et al. 2013;Nicholl et al. 2014Nicholl et al. , 2015aNicholl et al. , 2017)).
Observationally, SLSNe I have blue spectra (with a blackbody temperature spanning  BB = 10000 − 20000 K) during the premaximum/maximum phases, with prominent absorptions between 3000-5000 Å.These features are usually identified as O ii (e. g.Mazzali et al. 2016), although there is no general consensus on this (as discussed later).After 2-3 weeks from maximum, SLSNe I enter a new phase in which their spectra bear a resemblance to those of broad lined SNe Ic (SNe Ic BL) at their maximum luminosity (e.g.Pastorello et al. 2010).Photometrically, their LCs are very heterogeneous and may evolve over a wide range of timescales (see e.g.Fig. 5 of De Cia et al. 2018).Inserra et al. (2018) suggested to distinctly separate rapidly declining SLSNe I from slow ones, but larger SLSNe-I samples point towards a continuous distribution of timescales (Nicholl et al. 2015a;De Cia et al. 2018;Lunnan et al. 2018;Angus et al. 2019).Recently, Könyves-Tóth & Vinkó (2021) and Könyves-Tóth (2022) proposed a new SLSNe-I subclassification based on their features in pre-maximum/maximum spectra: type-W SLSNe I, whose absorptions are well fitted by O ii at reasonable velocities and type-15bn SLSNe I, whose early spectral features are not easily explained by O ii and show similarities with the coeval spectra of SN 2015bn (Nicholl et al. 2016a).
In the present work, we deal with the SLSN I SN 2019neq, which was classified by Könyves-Tóth & Vinkó (2021) in the type-W subgroup.At a redshift  ≃ 0.1059 (see Sec. 3), SN 2019neq is located at RA = 17 h 54 m 26.s 736, Dec = +47 • 15 ′ 40.′′ 62 and it is likely associated with the galaxy SDSS J175426.70+471542.3.It was discovered on 2019 August 9 by the Zwicky Transient Facility (ZTF, Bellm et al. 2019) with an apparent magnitude of  = 20.4mag (Perley et al. 2019), and named with the internal designation of ZTF19abpbopt.A few days later, a spectrum observed with the Palomar 60-inch+SED Machine revealed a hot continuum with some unidentified features (Perley et al. 2019).After ∼ 18 days, a spectrum of SN 2019neq was taken at the Liverpool Telescope (LT, Steele et al. 2004) (Roque de Los Muchachos Observatory, La Palma, Spain) equipped with the SPRAT (SPectrograph for the Rapid Acquisition of Transients, Piascik et al. 2014) instrument and classified as a SLSN I ( Könyves-Tóth et al. 2019).A  = 17.2 mag detection implied that it was still in the rising phase with a rest-frame rate of ∼ 0.2 mag/day.SN 2019neq was already included in several studies: Könyves-Tóth et al. (2020a) presented and discussed the ZTF photometry of SN 2019neq and three photospheric spectra.Furthermore, SN 2019neq was included in four sample papers (Hosseinzadeh et al. 2022;Chen et al. 2023a,b;Pursiainen et al. 2023): in particular, Hosseinzadeh et al. (2022) and Chen et al. (2023a,b) interpret postmaximum undulations in the -and -filter LCs of SN 2019neq as SLSNe-I bumps, i. e. they attributed the undulations either to a variable energy source or to an external one (e.g.CSM interaction).Pursiainen et al. (2023) analyzed linear-polarimetry data of 7 Hpoor SLSNe I and concluded that SLSNe I with oscillating LCs usually show an increase of the degree of polarimetry: therein, SN 2019neq was included in the non-oscillating SLSNe I subsample.In this work, we present a deep photometric and spectroscopic dataset of SN 2019neq and, based on new coeval photometry, disfavour that the undulations seen in its -and -filters LCs have an astrophysical origin.Our analysis is in agreement with the ejecta mass estimated by Könyves-Tóth et al. (2020b) in an independent way.
In detail, in Sec. 2 and 3 we introduce the photometric and the spectroscopic observations of SN 2019neq, respectively; in Sec. 4 we compare SN 2019neq LCs and spectra with other SLSNe I (Sec.4.1), we study its color and temperature evolution (Sec.4.2), and the evolution of its expansion velocity via the O ii spectral absorption features (Sec. 4.3).Finally, we discuss the magnetar and the CSMinteraction interpretations (Sec.4.4) using the nebular spectrum of SN 2019neq (Sec.4.4.1) and modelling its observed multicolor LCs (Sec.4.4.2).Throughout the paper, we will assume a flat Universe with Ω M = 0.31, Ω Λ = 0.69 and with a Hubble constant  0 = 71 ± 3 km s −1 Mpc −1 ; we took this value of  0 as an average among the estimates provided by Planck Collaboration et al. ( 2016), Khetan et al. (2021) and Riess et al. (2021).

Observations and data reductions
Ultraviolet/optical/near-infrared imaging data of SN 2019neq were obtained via different facilities.In detail, we used the NOT Unbiased Transients Survey 2 (NUTS2 1 , Mattila et al. 2016;Holmbo et al. 2019) at the 2.56-m Nordic Optical Telescope (NOT)+ALFOSC/NOTCam, the 2.0-m Liverpool Telescope (LT)+IO:O, La Palma observatory, Spain; the 1.82-m Copernico Telescope+AFOSC and Schmidt telescopes at the Asiago Astrophysical Observatory, Italy; the Tsinghua-NAOC (National Astronomical Observatories of China) Telescope (TNT)+BFOSC (Beĳing Faint Object Spectrograph and Camera), Xinglong Observatory, China (Wang et al. 2008;Huang et al. 2015).2, 2, 1, , ,  follow up was triggered with the Neil Gehrels Swift Obser-vatory+Ultraviolet/Optical Telescope (UVOT, Gehrels et al. 2004).Swift/UVOT 2, 2, 1, , , -filter frames were reduced 2 with the heasoft package (version 6.25 (HEASARC) 2014).The ground-based , , , , , ,  photometric frames were preliminary pre-processed, i.e. they were corrected for overscan, bias and flat field.Magnitude measurements were performed by means of the ecsnoopy pacakage 3 (Cappellaro 2014).Before measuring the SN flux, we accounted for the contribution of the background contaminating the SN light either with a polynomial interpolation or with the template-subtraction technique.The latter was performed with the hotpants tool (Becker 2015).When possible, we used template frames obtained with the same instrumental setting used for the scientific observations and observed after that the SN faded well below the detection limit.When deep template frames were not available, we estimated the background contribution interpolating a low-order polynomial to the area surrounding the SN position.This was interpreted as the background level and then subtracted from the photometric frames contaminated by the SN background.SN magnitudes were then measured fitting a point spread function to the SN (Stetson 1987).A detailed description of the image reduction and measurements procedures can be found in Fiore et al. (2021).Instrumental , , , , , ,  magnitudes were calibrated on a sequence of non-saturated field stars from the SDSS (Sloan Digital Sky Survey) and the Pan-STARRS (Panoramic Survey Telescope and Rapid Response System, Chambers et al. 2016) surveys, respectively.For , ,  filters, we converted the reference star magnitudes from Sloan to Johnson system following Chonis & Gaskell (2008)., ,  s magnitudes were calibrated with reference to field stars present in the 2MASS catalogue (Two-Micron All Sky Survey, Skrutskie et al. 2006).Given the very rapid evolution of SN 2019neq, deep , , , , , ,  template frames were observed on 2020 June 12 (MJD = 59012, only 281 rest-frame days after the maximum luminosity) by NUTS2.Since the SN was not detected in the template frames within a 3 detection limit, we assumed that the luminosity of SN 2019neq faded well below the galaxy luminosity.Finally, it was not possible to template subtract the , ,  s -filter frames since no suitable template was available for these filters.Therefore we estimated the background contribution fitting a low-order polynomial.We also include  and  ATLAS (Asteroid Terrestrial-impact Last Alert System, Tonry et al. 2018a) photometry, which we converted to standard  and  via Tonry et al. (2018b).The results of our magnitude measurements are listed in Tabs.C1, C2, C3, C4.Finally, we applied to these data the S-and K-corrections to bring back the many instrumental photometric systems to the corresponding standard and rest-frame ones (see Appendix A for details).
2 Their calibration was done using the updated version (November 2020) of the sensitivity corrections. 3ecsnoopy is a python package for SN photometry using PSF fitting and/or template subtraction developed by E. Cappellaro.A package description can be found at http://sngroup.oapd.inaf.it/ .

Observed and bolometric light curves
The S-corrected LCs are shown in Fig. 1.To investigate the presence of a possible pre-maximum bump feature in the LC, and to better constrain the rising phase towards the maximum luminosity, we added publicly-available ZTF -and -filter photometry.The latest non-detection fainter than the polynomial fit of the early  LC corresponds to 2019 August 7 (MJD = 58703.312, = 20.7 mag), and the first detection in  filter ( = 20.4mag) was about one day later (MJD = 58704.307).It is therefore possible to estimate that the explosion date occurred on MJD ≃ 58704 ± 1 (2019 August 9).Early -filter detection limits exclude the occurrence of a pre-maximum bump within ∼ 79 days before the estimated explosion date.To estimate the maximum-luminosity epoch, we fitted a high-order polynomial to the -filter LC and find that the maximum occurred on MJD MAX = 58731 ± 2 (2019 September 5) at a magnitude  MAX = 17.07 ± 0.10 mag (uncertainties were established by varying the order of the fitted polynomial).Given a luminosity distance  L = 482.2+21.3  −19.5 Mpc, corresponding to a distance modulus  = 38.416+0.094  −0.090 mag, and a Galactic extinction   = 0.104 mag (Schlafly & Finkbeiner 2011), the peak absolute magnitude of SN 2019neq is   = −21.5 ± 0.2 mag.As no narrow absorption lines from the Na iD doublet (Poznanski et al. 2012) are seen in the spectra, we assume no extinction from the host galaxy.In the following, we refer to the rest-frame time with respect to maximum as "phase" ().The assumed explosion epoch implies that the phase at the point of explosion is −24 days.
The pseudo-bolometric LC of SN 2019neq was computed by integrating the broad band photometry, adopting as reference the epochs of the -filter photometry (see Fig. 2, Tab.C10 for the Kcorrection and Tab.C11 for the pseudo-bolometric LC).SN2019neq rises towards the maximum luminosity very rapidly at a rate  rise ≈ −0.13 mag d −1 , it reaches a maximum luminosity of (2.04 ± 0.04) × 10 44 erg s −1 and fades at a rate  decline ≈ 0.05 mag d −1 , which is ∼ 2.6 times slower than the rising rate.This ratio is somewhat higher than the typical  decline / rise ≈ 2 ratio for SLSNe I (Nicholl et al. 2015a).

LC bumps in SN 2019neq
LCs of SLSNe I often display pre-maximum and/or post-maximum bumps, which are usually interpreted as signatures of CSMinteraction.Their occurrence seems to be more common for the slower-evolving events (see e. g.Yan et al. 2015Yan et al. , 2017;;Nicholl et al. 2015bNicholl et al. , 2016b)), but this connection is debated and might simply be due to the fact that slow SLSNe I can be observed for a longer time.Hosseinzadeh et al. (2022) and Chen et al. (2023b) noticed postand pre-maximum undulations in the SN 2019neq LCs respectively, and included it in sample studies of bumpy SLSNe I.In particular, Hosseinzadeh et al. (2022) used ZTF, Pan-STARRS and AT-LAS data, together with , ,  observations from the 1.2m tele-scope+KeplerCam at the Fred Lawrence Whipple Observatory for SN 2019neq (Szentgyorgyi et al. 2005).Their bumpy-SLSNe I sample consists of 34 objects whose LCs deviate from the best-fit magnetar-powered mosfit LC, which introduces a model dependency in the definition of a 'bump'.Their KeplerCam , ,  measurements are in good agreement with our coeval photometric data, which indeed seem to show a flattening at about 80 days after maximum (see Fig. 1, upper-right panel) and reveal ≲ 0.3 mag-amplitude undulations occurring on a timescale of ≲ 2 days.However, given the interruption of the follow-up at this phase, we are not able to reach  a firm conclusion on the veracity of physical origin of these LC features.Chen et al. (2023b) discussed ZTF photometry which we also use in this work.In the pre-maximum data (at a phase of ∼ −20), the  and possibly  LCs appear to display a hump (see Fig. 1, lower-right panel).This deviation from the general rising slope of the LCs is however encompassed by the errorbars, although the possible correspondence two bands (the only in which we have very early detections) might provide some additional support.However, as evidence for the presence of bumps is relatively weak, we do not consider data from Hosseinzadeh et al. (2022) in the following analysis.

Observations and data reductions
Optical spectra of SN 2019neq were obtained with the 1.82m Coper-nico+AFOSC telescope, with NOT+ALFOSC via NUTS2, with the Xinglong 2.16-m telescope+BFOSC and with the 10.4m Gran Telescopio CANARIAS (GTC)+OSIRIS (Optical System for Imaging and low-Intermediate-Resolution Integrated Spectroscopy, Cepa et al. 2000).The raw two-dimensional spectroscopic frames were preprocessed, wavelength-calibrated, extracted and flux-calibrated with standard iraf procedures called via the graphical user interface4 foscgui5 .BFOSC spectra were reduced with iraf procedures directly.
The flux calibration of each spectrum was then checked against the coeval optical photometry (which was interpolated in case of missing epochs).

The spectra
The spectral evolution of SN 2019neq is shown in Fig. 3 (see also Tab.C12).Pre-maximum/maximum spectra of SN 2019neq have a hot blue continuum with black-body temperatures of about  BB ≈ 16000 K.The hot continuum is almost featureless, with only a Wshaped O ii absorption feature and the O i  7774 feature visible in the red part.At bluer wavelengths, the Ca ii H&K doublet is also present.After 10 days from maximum, the continuum cools, the Wshaped features gradually disappear in favour of the Fe ii and Mg ii features and the spectra of SN 2019neq start to resemble a SN Ic BL spectrum at maximum luminosity.In addition, at this phase Si ii was identified by Könyves-Tóth et al. (2020a) with a spectral fit obtained with syn++.After +19 days, the continuum gradually cools down, and up to +80 days the spectral features do not evolve significantly, except for their intensity and velocity (see Sec. 4.3).At late phases, when the forbidden emission lines start to be seen in the spectrum, the Mg ii feature is expected to be contaminated or replaced by the semiforbidden Mg i]  4571 feature; however, this region is not covered by the latest spectrum (+231 days).Finally, we used the positions of the narrow H, [O iii] host-galaxy emission lines measured in the +68, +80 and +231 days spectra to derive a host galaxy redshift of  = 0.10592 ± 0.00005 (where the uncertainty is given by the dispersion of the measurements).
We also estimated the metallicity at the site of SN 2019neq by means of the narrow emission lines from the host galaxy which show up at late/nebular phases.We measured the flux emitted within these lines after optimizing the extraction from the two-dimensional spectrum for the host-galaxy and calibrating it against the host photometry at the SN location (see also Lunnan et al. 2014).To derive the metallicity, we used the tool pymcz (Bianco et al. 2016).In short, pymcz computes the metallicity from the flux measurements of those narrow emission lines via a number of metallicity diagnostics.To associate the errorbars to the metallicty measurements, pymcz randomly samples a Gaussian distribution whose mean and standard deviation are given by the flux measurements and their uncertainties, respectively.In the case of SN 2019neq, we measured the flux emitted within the 6717 emission lines.The nebular spectrum of SN 2019neq at 231 days after maximum does not show the [N ii]  6548, 6583 narrow emission lines.The only diagnostic which could be suitably used is then the m08_o3o2 (Maiolino et al. 2008), giving a metallicity 12 + log 10 (O/H) ≃ 8.3, which corresponds to  ≃ 0.4 Z ⊙ assuming a solar metallicity 12 + log 10 (O/H) = 8.69 (Asplund et al. 2009).
Finally, we estimated the star formation rate (SFR) of the host galaxy of SN 2019neq based on the measurements of the flux emitted by the reddening corrected narrow H using equation 2 of Kennicutt (1998).The derived SFR is ≃ 2.6 M ⊙ yr −1 , similar to the SFRs measured by Chen et al. (2017) for a sample of galaxies hosting SLSNe I.In addition, we estimated the specific SFR, i.e. the SFR expressed per unit of stellar mass content of the host galaxy, sSFR ≡ SFR/ * .To estimate  * , we used equation 8 of Taylor et al. (2011) after having measured the  host = 20.4 ± 0.1 mag and  host = 19.7 ± 0.1 mag for the host galaxy.We adopted Galactic extinction corrections from Schlafly & Finkbeiner (2011) and we obtained K-corrections scaling a starburst-galaxy template from Calzetti et al. (1994) to the  host and  host magnitudes, thus obtaining log 10  * / ⊙ ≃ 9.1 and sSFR ≃ 2.3 Gyr −1 .These values are in a good agreement with the average properties of a sample of 31 SLSNe I host galaxies studied by Lunnan et al. (2014).We stress that the SFR and sSFR values deduced can be highly affected by the error on the K-corrections.

Photometric and spectroscopic comparisons
In order to highlight its similarities/differences with other SLSNe I, we selected a sample including both very slow-and fast-evolving SLSNe I to compare their K-corrected -filter absolute LCs and spectra with those of SN 2019neq.The slow-SLSNe I subsample consists of SN 2020wnt (Gutiérrez et al. 2022), SN 2015bn (Nicholl et al. 2016a) and PTF12dam (Nicholl et al. 2013;Chen et al. 2015;Vreeswĳk et al. 2017), while the fast-SLSNe I one includes SN 2011ke and PTF11rks (Inserra et al. 2013;Quimby et al. 2018).In addition, we imported the best-sampled and non-heavily oscillating -filter ZTF LCs among those published in Chen et al. (2023a) (see Tab. 1).We compared the LCs of the SLSNe I of this sample with that of SN 2019neq.In addition, we investigated if SN 2019neq has a peculiar behaviour compared to the expected correlation between the SLSNe I maximum luminosity and the evolutionary timescales (see Inserra et al. 2018).4000 5000 6000 7000 8000 9000 Rest wavelength λ The LC comparison was done both in absolute and normalized units (see Fig. 4, left and right panel, respectively).When the Kcorrections were not available from the spectra, we approximate it as −2.5 log 10 (1 + ) (Hogg et al. 2002), which was shown to be a reasonable approximation for SLSNe I (Chen et al. 2023a).The LCs of the ZTF sub-sample are coloured according to their evolutionary velocity: we used the normalized magnitude at a phase of -15 rest-frame days,  norm,−15 , as a proxy for it.Hence, faster-evolving SLSNe I have greater values of  norm,−15 and vice versa.To mea-sure  norm,−15 , for each SN we fitted a polynomial to the -filter LC around the maximum luminosity.The degree of the polynomial and the fit domain were varied in order to minimize the root-mean squared of the fit.This quantity was also used as errorbar on  norm,−15 , although this choice may overestimate it.The values of  norm,−15 and  max of the selected comparison sample are listed in Tab. 1 too.We then plot the  max values versus  norm,−15 (Fig. 5).Our sample confirms the above-mentioned correlation between peak luminosities and the evolutionary timescales, with some exceptions nonetheless: SN 2019neq is actually one of the outliers as it evolves faster than the other SLSNe I with comparable peak luminosities; SN2018bym and the peculiar SN 2020wnt are slower than the other SLSNe I having comparable peak luminosities.On the contrary, SN 2015bn, PTF12dam, SN 2010gx and SN 2011ke better agree with the general trend followed by the ZTF subsample (PTF11rks was excluded due to the lack of pre-maximum data).However, much larger samples are required to draw general conclusions on this.Furthermore, we compared three spectra of SN 2019neq (at phases −6, +10, +231 days) with the available spectra for the comparison sample (including some of those belonging to the ZTF sample, see Fig. 6) at similar phases.In particular, we included pre-maximum spectra of SN 2019nhs and SN 2019unb and post maximum spectrum of SN 2019cdt.We also added a spectrum of the type-Ic BL SN 2003jd at about its maximum luminosity given its resemblance to SLSNe I spectra after maximum (see Sec. 1) .The absorptions on the blue part of the early spectrum of SN 2019neq are different from those of SN 2015bn, PTF11rks and SN 2020wnt, and hence due to other transitions (see also Könyves-Tóth & Vinkó 2021), while they possibly share the broad features from Si ii and O ii on the red side.In particular, the spectrum of SN 2019neq at −6 days looks more similar to the almost coeval one of SN 2019unb, whose location in the ( norm,−15 ,  max ) space is also peculiar (see Tab. 1).At later phases, the spectrum of SN 2019neq (+10 days after maximum) is similar to that of PTF11rks at the same phase and shows several broad features which nearly reproduce the spectral behaviour of SNe Ic BL at maximum luminosity.To see this, we compared the spectrum of SN 2019neq at 10 days after maximum with that of the SN 2003jd at maximum luminosity, which was given as best-match template by gelato (Harutyunyan et al. 2008).Also the post-maximum spectrum of SN 2020wnt, PTF12dam and SN 2019cdt are similar to the other spectra shown at comparable phases, but the prominent feature seen in these spectra at about 4700 Å cannot be easily explained by the Mg ii 4571 only.
Finally, the spectrum of SN 2019neq observed 231 days after maximum was compared with the spectra of PTF12dam and SN 2020wnt.At these epochs, the spectra look different.In particular, the nebular spectra of PTF12dam and SN 2020wnt still show some broad features between 5000 and 6000 Å and more prominent [Ca ii]  7291, 7323 and O i  7774 features.The nebular spectrum of SN 2019neq looks instead pretty much featureless, except for the [O i]  6300,6364; however, its lower signal-to-noise ratio makes a more thorough comparison difficult.

SED, colours, photospheric temperature and radius evolutions
We discuss the SED evolution by analyzing the time variation of properties deduced from the multicolor photometry presented in Sec. 2.
In detail, we present the rest-frame colours (Fig. 7, 8), photospheric temperature and radius evolution.Photospheric temperatures were estimated by fitting a blackbody curve to the SED and the photospheric radius was inferred via the Stefan-Boltzmann law.To exclude any possible UV/NIR deviations from the blackbody, we repeated the calculation excluding the 2, 2, 1, , ,  s -fluxes from the SED (Figs. 9, 10).In addition, for each epoch  =  + 24 (hence the phase from explosion, assuming an explosion phase of −24 days, see Sec. 2.2) we estimated the photospheric temperature and radius from spectra: to do this, we fitted blackbody curve to spectra and computed for each epoch  the radius as  phot () ×  (see Dessart et al. 2015, and discussion in Sec. 4.3).However, we found no major temperature differences between the three estimates and conclude that, at least at photospheric post-maximum phases (see later), the SED is reasonably well described by a blackbody over UV/optical/NIR wavelengths.
We compare the color evolution of SN 2019neq with other SLSNe I with different evolutionary timescales.After the initial nearly-constant phase, the colours of SN 2019neq rapidly evolve towards red similar to the fastest-evolving SLSNe I of this sample (see e.g.PTF11rks, SN 2010gx and SN 2020ank, Inserra et al. 2013;Pastorello et al. 2010;Kumar et al. 2020).In contrast, colours of slower evolving SLSNe I redden at a slower pace compared to the fastest ones; the behaviour of SN 2017gci is odd, but its atypical color evolution is likely ascribable to the ∼ 0.6-mag-wide bumps in its LC (Fiore et al. 2021).However, in all cases SLSNe I display blue colours around the pre-maximum/early-post-maximum phases (Pastorello et al. 2010;Quimby et al. 2011;Chomiuk et al. 2011;Leloudas et al. 2012;Inserra et al. 2013).In particular, at this phase the rest-frame  −  colour remains almost constant at ∼ −0.25 mag for every SLSNe I considered in Fig. 7; which likely means that the SED does not significantly evolve.This is true also for SN 2019neq: to verify whether this behaviour is seen in other color indexes, we extend the time coverage of the observed photometry by performing synthetic photometry on the spectra (see dotted lines in Fig. 8), this is particularly useful in  band (although the flux calibration of the spectra is less precise and synthetic-photometry measurements are in part extrapolated, see Fig. 8, right panel).As can be seen in Fig. 8, the initial flattening appears present in  −  and  −  too, although very early photometric measurements can be done only in  and  filters.Subsequently, we see that after ∼ 30 − 40 days from the maximum luminosity the  −  colour curve flattens again.
The nearly constant initial behaviour of the color could reflect the early photospheric-temperature evolution; however, due to the lack of UV data at pre-maximum phases, the initial plateau phase cannot be fully trusted.Soon after maximum, the photospheric temperature starts to decline, reaching a new constant value after ∼ 30 days (similar to other SLSNe I, see also Nicholl et al. 2017).Interestingly, at comparable phases, the photospheric radius reaches a maximum (Fig. 10) and then starts to recede: a similar trend of the photospheric-radius evolution is described by a simple model (Liu et al. 2018) assuming homologous expansion and constant opacity.However, more detailed calculations of the photospheric radius with a proper time-varying opacity are needed to constrain the ejectadensity profile of SN 2019neq, and this will be the subject of future investigations.

Photospheric velocity
To estimate the photospheric velocity of SN 2019neq, the minima of the O i  7774 P-Cygni profiles were measured with a gaussian fit (see Fig. 11) after having normalized and continuum-subtracted the spectra.These measurements were performed on the spectra at phases −6 to 80 days with respect to the maximum luminosity.By measuring the Doppler shift compared to the rest-frame wavelength of the corresponding maximum, we obtained the expansion velocity.This was   used as a proxy for the photospheric velocity (Dessart et al. 2015) (see Fig. 12).Errorbars were estimated by changing the continuum level multiple times before performing the fit.The comparison between the radii deduced from photometry and spectroscopy confirms that the O i is a good tracer of the photosphere for SN 2019neq for  between ∼ −5 and 18 days, where we inferred an average photospheric velocity of ∼ 12500 km s −1 .Overall, the Doppler shift measured with respect to the rest-frame wavelength of the emissions corresponds to a photospheric velocity (O I) ≈ 12500 − 15000 km s −1 .We also compared the photospheric velocities of SN 2019neq with those of other SLSNe I whose velocity measurements are available in literature: iPTF13ehe (Yan et al. 2015), iPTF15esb (Yan et al. 2017), SN 2015bn (Nicholl et al. 2016a), SN 2017gci (Fiore et al. 2021) and SN 2018hti (Lin et al. 2020;Fiore et al. 2022).In addition, given the unusually high velocities of SN 2019neq, we included in the comparison sample also two SNe Ic BL, namely SN 2003jd (Valenti et al. 2008) and SN 2007ru (Sahu et al. 2009) (see Fig. 12).The high velocities of SN 2019neq do not seem to be correlated with a steeper average velocity gradient compared to slower SLSNe I (like SN 2015bn or SN 2018hti), as suggested by Inserra et al. (2018),(Könyves-Tóth & Vinkó 2021), whereas this remains true for iPTF13ehe and iPTF15esb.The same appears to be also valid for the SNe Ic BL SN 2007ru and SN 2003jd, the latter reaching even higher velocities at maximum epochs.We notice that the velocities of SN 2019neq deduced by us are in tension with those deduced by  2021) used a different method to infer the  phot values: they cross-correlated the observed spectra of SN 2019neq with a template syn++ spectrum computed with  phot = 10000 km s −1 and found a velocity difference Δ X .The Δ X values were then used to obtain the  phot measurements after having applied a correction (see Section 4.1 and equations 10, 11 in Könyves-Tóth & Vinkó 2021).

Modelling SN 2019neq
We used the nebular spectrum of SN 2019neq at 231 days to infer an estimate of the mass of the oxygen zone of the ejecta.The absolute integrated luminosity emitted within the O i  7774 feature can be used to infer the O-zone mass (see equations 7,8 of Jerkstrand et al. 2017).In addition, we used the sumo models for the nebular spectra6 of SNe Ic (Jerkstrand et al. 2017).The broad constraint obtained on the ejecta mass will be used in Sec.4.4.2 as prior to favour one of the possible best-fit bolometric LCs obtained with the mosfit (Guillochon et al. 2017b,a) tool.For this same SN, Könyves-Tóth et al. ( 2020b) estimated an ejecta mass of ∼ 23 M ⊙ assuming an ejecta opacity  = 0.1 cm 2 g −1 .In the following, we will more carefully investigate this comparing the photometric best-fit value with the analysis of nebular spectroscopy.

Nebular-spectrum modelling
The sumo nebular spectra are single-zone models computed for three different compositions (pure-O, C-burning ashes and OMg).For each composition, spectra calculations assume a phase of 400 days after the explosion, the homologous expansion of the ejecta at a constant velocity  = 8000 km s −1 and  = 100 random clumps for different ejecta masses7  sumo ejecta = 3, 10, 30, energy deposition  dep = 2.5 × with which sumo models are computed as  sumo ejecta to distinguish them from other ejecta-mass estimates.Before measuring the O i  7774 absolute luminosity, we estimated the residual contribution of the host-galaxy emission using the hostgalaxy spectrum adjacent to that of SN 2019neq.Then this was scaled and subtracted from the spectrum of SN 2019neq until the narrow galaxy emission was removed from the SN spectrum.Then, we measured an integrated luminosity of the O i  7774 line of  7774 ≲ 2 × 10 40 erg s −1 .Due to the relatively low signal-to-noise ratio between 7000-8000 Å, we assume this value as an upper limit.Similar to the sumo models, we assumed a maximum velocity of the fluid elements  = 8000 km s −1 .Aside this parameter, the only tunable quantities are the clumping factor  and the electron fraction   .We computed the O-zone mass for the  values with which sumo models are computed, i.e.  = 0.1, 0.01, 0.001, a phase from the explosion of 231 + 24 days = 255 days (see Sec. 2.2) and we consider a range of   = 0.05 − 0.1 (this range encompasses the typical values used in sumo models, Jerkstrand et al. 2017).Results are shown in Tab. 2. The case  = 0.1 leads to very high O-zone masses ( () > 100 M ⊙ ): based on the rapid LC-evolution timescales, in the following we will not consider this case anymore.The  = 0.01, 0.001 cases allow for lower O-zone masses permitted by lowermass ejecta.In addition, our mass estimates are strongly affected by our assumption on the electron fraction   , the clumping factor  and the maximum velocity of the clumps .
We then consider sumo models computed for  = 0.01, 0.001,  sumo ejecta = 3, 10, 30 for pure-Oxygen, OMg and C-burning abundances, and different values of the energy deposition  dep = 2 × 10 41 − 2 × 10 42 erg s −1 (see Figs. B1, B2, B3, B4, B5, B6).[O i] seems to be better described by models with  sumo ejecta =10 and  = 0.01 for higher energy deposition, thus possibly favouring higher values of the electron fraction   .Other models underpredict the [O i]  6300,6364 feature, except the C-burning models with  sumo ejecta = 30 and  dep > 1 × 10 42 erg s −1 .Furthermore, OMg-composition models seem less suitable to explain the O i features seen in the nebular spectrum of SN 2019neq, in particular for  dep = 2 × 10 41 erg s −1 (see Figs. B5,B6).However, the temporal discrepancy between the observed nebular spectrum of SN 2019neq and the sumo solutions makes the sumo spectra by a factor (400/255) 3 ∼ 3.9 less dense than the observed, if ejecta mass and expansion velocity are the same: this makes their spectral comparisons less relevant to discriminate among the different ejecta configurations.In addition, we warn the reader that a more careful estimate of the physical parameters of SN 2019neq would need independent constraints to break some degeneracies (e. g. between   and  or between  in and  ejecta ); this is however beyond the scope of the present work.

Light-curve modelling
We modelled the LCs of SN 2019neq to estimate the physical parameters of the explosion assuming the magnetar and the ejecta-CSM interaction.We used the Modular Open Source Fitter for Transients (mosfit) tool (Guillochon et al. 2017b,a) which offers a set of modules to fit the observed multicolour LCs for different kind of transients.For SN 2019neq, we ran mosfit using the slsn and the csm modules, computing synthetic LCs powered by magnetar spindown (Kasen & Bildsten 2010;Inserra et al. 2013) and CSM (Chatzopoulos et al. 2012) interaction, respectively.We fixed  = 0.1 cm 2 g −1 for both models which is suitable for SNe Ic (e.g.Nagy 2018).
The csm fit was set up adopting a broken power-law ejecta density profile as in Chatzopoulos et al. (2012) with fixed exponents  = 12 and  = 2, suitable for a steady wind-like CSM.A shell-like csm ( = 0) fit is not well constrained by the data, hence we excluded it.The  = 0 csm best-fit parameters are a CSM mass of  CSM ≃ 0.7 M ⊙ , a progenitor radius  0 ≃ 1.4 × 10 14 cm, a CSM density  ≃ 1.8 × 10 −10 g cm −3 , an ejecta mass  ejecta ≃ 21.9 M ⊙ , an average ejecta velocity  ejecta ≃ 4100 km s −1 and a temperature floor  min ≃ 3700 K.This corresponds to a kinetic energy of ∼ 3.7 × 10 51 erg.The ejecta mass for the csm fit broadly agrees with the sumo models computed  sumo ejecta = 30.The slsn fit predicts a lower ejecta mass (∼ 12 M ⊙ ), thus in better agreement with sumo spectra with  sumo ejecta = 10−30.This is especially true for sumo models computed with  sumo ejecta = 10 and  = 0.01.To verify whether higher ejecta-mass models return reasonable results, we performed another slsn fit, this time fixing an ejecta mass  ejecta = 25 M ⊙ and fitting the opacity: not surprisingly, in this last case the best-fit opacity is lower ( = 0.05 cm 2 g −1 ), but however reasonable in the case of H/He-poor events (Nagy 2018).Both the  = 0.05 cm 2 g −1 and the  = 0.1 cm 2 g −1 slsn fits provide ejecta-mass estimates that are broadly consistent with what can be inferred via the nebular-spectrum interpretation presented in the previous Section.The other best-fit parameters for the two slsn fits are a magnetic field  ≃ 6 × 10 14 G, an initial spin period  init ∼ 1.1 − 1.5 ms, a temperature floor  min ≃ 6000 K and an average ejecta velocity  ejecta ≃ 10000 km s −1 .This correspond to a kinetic energy of 1.3 − 2.7 × 10 52 erg.The best-fit parameters for the csm and slsn are listed in Tab.3,4 and the corner plots showing their posterior distributions are shown in Fig. B7, B8, respectively.The slsn and csm LCs are shown in Fig. 13, 14 and are able to capture the pre-maximum/maximum luminosity epochs.The early post-maximum photometric detections in 2 and 2 bands cannot be properly described by both models; moreover, red and NIR-bands photometric points at a phase later than ∼ +60 days are better described by the csm fit.However, the UV/NIR photometric coverage is too scarce to use the previous argument as discriminating factors between the two models.We notice also that the  min and the  ejecta values predicted by the slsn model are in better agreement with the measured data (Figs.9, 12) compared to those predicted in the csm case.
Moreover, the spectra of SN 2019neq do not display the prominent narrow/multicomponent emission features seen e. g. in the luminous type- IIn SN 2006gy and SN 2006tf (Smith et al. 2007, 2008;Kiewe et al. 2012).Based on these reasons, we disfavour the CSM-interaction scenario for SN 2019neq and suggest that a millisecond magnetar endowed with a magnetic field  ≃ 6 × 10 14 G is the engine driving the observed luminosity of SN 2019neq.CSM interaction might have contributed during pre-/post-maximum epochs and resulted in the possible complexities of the LCs of SN 2019neq.However, we cannot rule out that more peculiar (e. g.  disk-like) CSM configurations might have reprocessed the UV/X-ray photons before reaching the photosphere, thus quenching the narrow/multicomponent line formation (Andrews & Smith 2018).

CONCLUSIONS
We analyzed the spectrophotometric observations of SN 2019neq, which is among the fastest-evolving SLSNe I ever observed.Its LCs rise towards maximum in about 24 rest-frame days and decline at a ∼ 2.6 times slower rate.Similar to other SLSNe I, SN 2019neq displays approximately constant blue colours during the pre-maximum/maximum epochs.At these phases, the spectra display W-shaped absorptions attributed to O ii, and at a phase 16-19 days remind the spectra of SNe Ic/Ic BL at maximum luminosity.From the narrow lines attributed to the host galaxy, we inferred a metallicity  ∼ 0.4  ⊙ , a SFR ∼ 2.6 M ⊙ yr −1 and a sSFR of ∼ 2.3 Gyr −1 at the SN site.A late nebular spectrum of SN 2019neq was scaled in absolute specific luminosity and compared with different models computed with the sumo single-zone, radiative-transfer code designated for nebular spectra of SNe Ic; we also used a set of analytical relations  (Jerkstrand et al. 2017) to estimate the ejecta O mass.In addition, multicolor LCs of SN 2019neq were modelled with the mosfit tool using the csm and slsn modules to test the ejecta-CSM interaction and the magnetar scenarios, respectively.Based on the best-fit results and using the inferred nebular properties as a broad prior, we conclude that the spindown radiation of a millisecond, young magnetar with a magnetic field of ∼ 6 × 10 14 G likely powered the luminosity of SN 2019neq at around the maximum luminosity.Given the degeneracy between the ejecta mass  ejecta and the ejecta opacity , it is hard to retrieve reasonable estimates, hence we provide for them the following ranges:  ejecta ≃ 10 − 30 M ⊙ and  ≃ 0.05 − 0.1 cm 2 g −1 .
In addition, we cannot rule out that the ejecta of SN 2019neq have interacted with a previously-ejected CSM, but either this requires a peculiar CSM topology, or that its contribution is subdominant in powering the SN around maximum luminosity.
Due to their rapid photometric evolution, SLSNe I like SN 2019neq are challenging to follow up, in particular at nebular epochs.Their careful study warrants a strong observational effort which future wide-field surveys like the Legacy Survey of Space and Time at the Vera Rubin Observatory (Villar et al. 2018) will allow for.In addition, a larger dataset is useful to improve the statistical significance of SLSNe sample studies.These are crucial to understand the diversity of SLSNe and to unravel the physical reasons for their luminous nature.Sec.4.2) and propagated in our analysis as an additional uncertainty due to the non-standard instrumental photometric systems.

A2 K-correction
We computed K-corrections to account for the effect of the cosmological redshift on the magnitude measured in the observer frame band-pass filters.Similar to S-corrections, we computed the Kcorrections performing synthetic magnitudes measurements on the spectra of SN 2019neq (see Sec. 3) via pysynphot.In detail, for each band-pass filter, we derived a synthetic magnitude both for the rest-frame spectrum (for which we computed a synthetic magnitude  ,rest ) and for the observed one (for which we computed a synthetic magnitude  ,obs ).For each epoch and filter, the K-correction was computed as  ,rest −  ,obs .The resulting K-corrections are listed in Tab.C10.Similar to the Δ corr , the K-corrections for the 2, 2, 1, , ,  s -filter magnitudes were estimated as before but using blackbody fits to the SED.

APPENDIX B: FIGURES APPENDIX C: TABLES
This paper has been typeset from a T E X/L A T E X file prepared by the author.
MNRAS 000, 1-16 ( 2022) The superluminous SN 2019neq 17 Figure1.S-corrected lightcurves for SN 2019neq in 2, 2, 1, , , , ,  ,  , , ,  , ,  s bands plotted in black, brown, cyan, pink, green, dark blue, red, purple, blue, magenta, orange, dark silver, silver, and yellow, respectively.Photometry obtained with different instruments are plotted with different symbols, as labelled in the upper right-hand corner.The -filter LC was fitted with a high-order polynomial (red solid line) to obtain the epoch of maximum luminosity (see text).Arrows correspond to 2.5 detection limits.Magnitudes are in the AB system.The upper inset on the right-hand side shows the period in which the late bump is visible (additional data fromHosseinzadeh et al. (2022) are represented by black dots); the lower-left and -right insets show the pre-maximum phases in  and  filters, respectively

Figure 2 .
Figure 2. Pseudo-bolometric LC of SN 2019neq.For comparison also the 56 Co-decay slope is also shown (black-dashed line).

Figure 3 .
Figure3.The spectral evolution of SN 2019neq reported to the rest-frame.The spectra (light blue solid lines) were also smoothed with a Savitzky-Golay filter(Savitzky & Golay 1964, black solid lines) to make the features more distinguishable.Spectra were offset and scaled.Identified spectral features are marked with black dashed lines at their corresponding rest-frame wavelength.Numbers on the right-hand side of each spectrum indicate their phase.

Figure 4 .Figure 5 .
Figure 4. Comparisons of the absolute -filter LC of SN 2019neq (black dots) with those of two fast-evolving SLSNe I, SN 2011ke (cyan triangles) and PTF11rks (purple hexagons) (data from Inserra et al. 2013) and two slow-evolving SLSNe I, SN 2015bn (green diamonds) (Nicholl et al. 2016a) and PTF12dam (orange diamonds) (Nicholl et al. 2013; Chen et al. 2015; Vreeswĳk et al. 2017).Right panel: LCs comparison normalized to the same maximum luminosity.Left panel: LCs comparison without rescaling.The absolute LCs were computed assuming the cosmological parameters used in this work (see Introduction).

Figure 6 .
Figure6.Comparison of three spectra of SN 2019neq (blue) with those of other SLSNe at broadly coeval phases with respect to maximum.The phases are reported in square brackets.Line identifications are marked with blue shaded areas and for each of them the corresponding ion is labelled on the top.For the last spectrum of SN 2019neq we plot both the original (light-blue) and its smoothed version (with a Savitzky-Golay filter, blue).

Figure 7 .
Figure 7. Rest-frame  −  evolution of SN 2019neq (blue filled dots) compared to that of a sample of fast-evolving (blue-edged markers), intermediate (magenta-edged markers) and slow-evolving SLSNe I (green-edged markers).The solid-blue line joins the rest-frame  −  points computed via synthetic photometry on the spectra of SN 2019neq.

Figure 8 .Figure 9 .
Figure 8. Left panels: rest-frame color evolution of SN 2019neq (red dots are photometric points and green dots are obtained from synthetic photometry on the observed spectra).Right panel: four representative spectra of SN 2019neq (see Fig. 3) within ∼ 80 days from the maximum luminosity.The dashed curves represent the throughput of the , , ,  and -passband filters.

Figure 11 .
Figure 11.Normalized and continuum subtracted spectra of SN 2019neq, with the photospheric velocity being measured by fitting the minimum of the O i 7774 absorption feature with a gaussian (black solid lines).
Photospheric velocity v phot [km s

Figure A1 .
Figure A1.S-corrections for SN 2019neq for the different instrumental configurations used for the photometric follow-up of SN 2019neq.

Figure B1 .
Figure B1.Comparison of the observed nebular spectrum of SN 2019neq (red line) with different sumo spectra computed for the C-burning model with a filling factor  = 0.001 for different ejecta masses  sumo ejecta = 10, 30 and different energy depositions  dep .Different colours correspond to different  sumo ejecta , as labelled in the top-right corner.Solid lines correspond to an energy deposition  dep = 2.5 × 10 41 erg s −1 , dotted lines to  dep = 5 × 10 41 erg s −1 , dashed-dotted lines to  dep = 1 × 10 42 erg s −1 , dashed lines to  dep = 2 × 10 42 erg s −1 .

Figure B8 .
Figure B8.Same as Fig. B7, but for the csm fit.

Table 1 .
The SLSNe I sample used in this work (including SN 2019neq) with the values of  norm,−15 and  max and the reference used for the data (errors in parenthesis).Absolute magnitudes in AB system, corrected for foreground extinction and K-corrected, and referred to the cosmological model assumed in the present work.

Table 3 .
mosfit csm physical parameters (the opacity, marked with ( * ) is a fixed parameter; the other ones are best-fit values).

Table 4 .
mosfit slsn physical parameters (the one with ( * ) is a fixed parameter, the other ones are best-fit values).

Table C7 .
S-corrections for IO:O filters.

Table C10 .
K-corrections expressed in magnitudes.