A Precursor Plateau and Pre-Maximum [O II] Emission in the Superluminous SN2019szu: A Pulsational Pair-Instability Candidate

We present a detailed study on SN2019szu, a Type I superluminous supernova at $z=0.213$, that displayed unique photometric and spectroscopic properties. Pan-STARRS and ZTF forced photometry shows a pre-explosion plateau lasting $\sim$ 40 days. Unlike other SLSNe that show decreasing photospheric temperatures with time, the optical colours show an apparent temperature increase from $\sim$15000 K to $\sim$20000 K over the first 70 days, likely caused by an additional pseudo-continuum in the spectrum. Remarkably, the spectrum displays a forbidden emission line even during the rising phase of the light curve, inconsistent with an apparently compact photosphere. We show that this early feature is [O II] $\lambda\lambda$7320,7330. We also see evidence for [O III] $\lambda\lambda$4959, 5007, and [O III] $\lambda$4363 further strengthening this line identification. Comparing with models for nebular emission, we find that the oxygen line fluxes and ratios can be reproduced with $\sim$0.25 M$_{\odot}$ of oxygen rich material with a density of $\sim10^{-15} \rm{g cm}^{-3}$. The low density suggests a circumstellar origin, but the early onset of the emission lines requires that this material was ejected within the final months before the terminal explosion, consistent with the timing of the precursor plateau. Interaction with denser material closer to the explosion likely produced the pseudo-continuum bluewards of $\sim$5500 \AA. We suggest that this event is one of the best candidates to date for a pulsational pair-instability ejection, with early pulses providing the low density material needed for the forbidden emission line, and collisions between the final shells of ejected material producing the pre-explosion plateau.


INTRODUCTION
Superluminous supernovae (SLSNe) are a class of supernovae (SNe), initially categorised as events with absolute magnitudes exceeding  < −21 mag (Quimby et al. 2011;Gal-Yam 2012).In addition to their high luminosities, these events radiate ∼ 10 51 erg when integrated over their broad light curves (Gal-Yam 2012; Lunnan et al. 2018b).The luminous nature of these events means we can observe them out to redshifts  > 4 (Cooke et al. 2012), and so even though the volumetric rate of these events is ∼ 1 in a few thousand SNe (Quimby et al. 2013;Frohmaier et al. 2021), they make up roughly 1% of the SNe discovered today (Fremling et al. 2020).This is aided by the wide-field surveys available which can probe the entire night sky instead of targeting only nearby massive galaxies.SLSNe were ★ E-mail: aa@star.sr.bham.ac.uk (AA) missed by earlier surveys due to their preference for metal-poor dwarf galaxies as hosts (Chen et al. 2013;Lunnan et al. 2014;Leloudas et al. 2015;Perley et al. 2016;Schulze et al. 2018).As more of these events have been discovered, the strict magnitude cut off for classification has since been replaced by spectral classification around peak luminosity (Gal-Yam 2019a; Quimby et al. 2018), driven by events with SLSN-like spectra but intermediate luminosities (De Cia et al. 2018;Angus et al. 2019;Gomez et al. 2022).
SLSNe can further be classified into Type I and Type II SLSNe, analogous to their less luminous counterparts.Type II SLSNe often resemble lower luminosity SNe IIn, with narrow hydrogen lines and a small subset displaying broad hydrogen lines (Kangas et al. 2022).Interaction with circumstellar material (CSM) is thought to the main power source for Type II SLSNe (Ofek et al. 2014;Inserra et al. 2018).Type I SLSNe (often simply called SLSNe) lack hydrogen in their spectra.These spectra are characterised by a steep blue con-tinuum indicative of high temperatures, and often show prominent O II absorption lines at early times (Quimby et al. 2011), eventually evolving to be similar to SNe Ic when at comparable temperatures (Pastorello et al. 2010).However, a small fraction of these events, show evidence of H at late times (≳ 30 days; Pursiainen et al. 2022;Yan et al. 2018) but this is not necessarily related directly to the power source at maximum light, and instead may be a product of interaction with its environment at scales ≳ 10 16 cm (Yan et al. 2015).A handful of events have also shown evidence for helium in their photospheric spectra (Quimby et al. 2018;Yan et al. 2020).
One of the biggest questions remaining about SLSNe pertains to their powering mechanism.Typical hydrogen-poor SNe are powered by the decay of radioactive nickel ( 56 Ni), but this explanation does not seem to work for most SLSNe for a number of reasons.Powering the peak luminosities of SLSNe with nickel decay would require ∼ 5 − 20 M ⊙ of 56 Ni, which is too high compared to the inferred ejected mass from the light curves (Inserra et al. 2013;Blanchard et al. 2018).This amount of nickel could likely only be produced in Pair-Instability SNe (PISN) of stars with initial masses  ≳ 140 M ⊙ (Heger & Woosley 2002;Gal-Yam et al. 2009).One of the best PISN candidates to date SN2018ibb is thought to be powered by 25−44 M ⊙ of freshly synthesised 56 Ni produced by a star with a helium core mass of 120 − 130 M ⊙ (Schulze et al. 2023).Although stars in this mass range required for PISNe have been observed, mass-loss on the main sequence makes PISN formation challenging except at very low metallicities (Yusof et al. 2013).However, some models suggest a magnetic field at the surface of the star could quench the mass loss for stars at solar metallicity, allowing enough mass to remain for the PI mechanism (Georgy et al. 2017).More problematic for SLSNe, any model producing this amount of nickel would result in a spectrum dominated by iron-group elements, which is inconsistent with the blue spectra of SLSNe (Dessart et al. 2012;Nicholl et al. 2013;Jerkstrand et al. 2016).However, these models often assume interaction between the SN ejecta and CSM is negligible which would only likely be the case if observed early enough that explosive nucleosynthesis is not affecting the layers of ejecta (Kasen et al. 2011).
Some theories suggest an internal power source such as the spin down energy of a magnetar or an accreting black hole could power SLSNe.However the accreted mass required in the latter scenario (≫ 100 M ⊙ ) often exceeds the mass of any reasonable star (Moriya 2018).In the magnetar scenario, the remnant is a fast rotating neutron star with a very strong magnetic field  ∼ 10 13 − 10 14 G (Kasen & Bildsten 2010).Nearly 10% of newly born neutron stars have B-fields in the range 10 13 − 10 15 G lasting over 1000 years after their birth, and so it is plausible that these could exist to produce SLSNe (Kouveliotou et al. 1998;Woods & Thompson 2006).This mechanism could explain the long duration of the light curves as the magnetar releases its rotational energy at the dipole spin-down rate, which remains high for days to weeks in this range of magnetic field strengths (Ostriker & Gunn 1971;Kasen & Bildsten 2010).
Another possible explanation to power SLSNe is interaction with CSM.This theory proposes a core collapse supernova that has large amounts of CSM created through stellar winds and ejections throughout the life of the progenitor star.The SN ejecta is able to catch up to this material because it has much higher velocities, and is rapidly decelerated if the CSM is massive enough.This creates a shock that deposits energy in the ejecta and CSM, the cooling of which can produce a bright and long-lived light curve (Chevalier & Irwin 2011;Smith et al. 2007).The mass of CSM required to efficiently power a bright light curve must be comparable to the ejecta mass (Chevalier & Irwin 2011;Ginzburg & Balberg 2012), ranging from a few solar masses up to a few tens (Nicholl et al. 2014).Getting this much mass close to the star just before explosion is difficult to explain using stellar winds, even in Wolf Rayet stars with high mass loss rates (Mauron & Josselin 2011;Sander et al. 2022).An alternative way to produce massive CSM is through discrete outbursts.Stars with masses in the range 70 − 140 M ⊙ are thought to undergo pulsational pair instability PPI eruptions (Woosley 2017).These stars are not massive enough to experience terminal pair instability, instead the star violently expels up to tens of solar masses worth of material towards the end of its life due to this mechanism (Woosley 2017).This has been suggested as a way to get sufficiently massive CSM to power SLSNe (Woosley et al. 2007).
The CSM model has been questioned as the main power source for hydrogen-poor SLSNe.We would naively expect to see narrow lines in the spectra from slow moving material if interaction was at play, but this is not seen in all SLSNe (Nicholl et al. 2015a).Instead, there is evidence that CSM interaction may play a role in powering some SLSNe including late time interaction producing H emission (Yan et al. 2018;Pursiainen et al. 2022), post peak bumps in the light curves (Hosseinzadeh et al. 2022), blue pseudo-continua, and early appearances of forbidden emission lines such as [O II] and [O III] (Schulze et al. 2023).Light curve and spectral modelling suggest that both central engines and CSM interaction may help to power this class of events (Chen et al. 2017(Chen et al. , 2023b)).Spectra taken at different phases of the SN evolution allow us to probe different regions of the ejecta.At early times the ejected material from the explosion is still optically thick and obscures the view of the inner layers.As the ejecta expands it becomes less dense, leading to more states leaving local thermal equilibrium (LTE) and lower populations of excited states, reducing the number of optically thick lines and bound-free continua.This spectral transition from photospheric to nebular is also driven by decreasing temperatures which results in fewer lines that are capable of significant cooling, and fewer excited states with enough population to provide opacity.This transition occurs typically on the timescale of hundreds of days and results in a spectrum dominated by low-lying forbidden transitions (Jerkstrand 2017).This contradicts observations in which some SLSNe have shown these forbidden emission lines early on during their photospheric phase.This includes SN2018ibb which displayed signs of a possible [Ca II] 7291,7323 at -1.4 days before peak, becoming prominent by 30 days later (Schulze et al. 2023).Other SLSNe have also shown signs of early forbidden emission lines including the earliest spectra of SN2007bi (Gal-Yam et al. 2009), and LSQ14an (Inserra et al. 2017), in which this line ∼ 50 days after peak was also attributed to [Ca II].This suggests some lower density regions exist in the ejecta or their surroundings.In principle these lines could have appeared even earlier if earlier spectra were obtained, challenging our understanding of the structure of this massive ejecta.
In this paper we present and analyse SN2019szu, a slowly evolving SLSN that showed forbidden emission lines remarkably soon after the time of explosion, at least 16 days before maximum light.We identify these lines as singly-and doubly-ionized oxygen, arising in a low density, hydrogen-poor CSM, and use this to place important constraints on the progenitor of this event.The structure of this paper is as follows.Section 2 outlines the data collected for this object.Section 3 covers the analysis of the host galaxy, and the photometric and spectroscopic data collected for the target.In Section 4, we discuss spectral models to fit this event as well as mosfit models of the light curve.We then discuss the implications of these results and how they fit into our understanding of SN2019szu in Section 5. Lastly, in Section 6 we present our conclusions.
Light curve of SN2019szu.All magnitudes are in given in the AB system and are not corrected for Milky Way extinction.The grey vertical line indicates time of explosion based on light curve fits from mosfit.Phase given in rest frame days with respect to maximum light in the -band.Includes data from ATLAS, Pan-STARRS, ZTF, Swift, LCO, and NTT.MJD is in the observer frame and 3 upper limits are indicated via inverted triangles.The ATLAS  and  bands are plotted without upper limits for clarity and the o band data points are binned to a 2 day cadence.Vertical markers on the top axis correspond to the spectra in Figure 6 and indicate when they were obtained.

Discovery and Classification
The Asteroid Terrestrial-impact Last Alert System (ATLAS) project (Tonry et al. 2018) discovered SN2019szu with the designation AT-LAS19ynd on 2019-10-21.The ATLAS-HKO (Haleakala) unit de-tected the supernova in the  band at 19.4 mag following a shallow non-detection 2 days prior at a limiting magnitude of  = 17.7 mag (Tonry et al. 2019).The transient was identified on multiple images and as it was coincident with a faint host galaxy (see Section 3.1), the ATLAS Transient server reported it as a SN candidate (Smith et al. 2020).An earlier detection was made by the Zwicky Transient Fa-  tion et al. 2020).The absolute magnitude coupled with the very blue shape of the spectral continuum, the dwarf nature of the host galaxy and its strong emission lines similar to other SLSN hosts (Leloudas et al. 2015) cemented the SLSN designation.However, since the initial spectrum did not cover H, later spectra were needed to confirm its lack of hydrogen and type I designation.SN2019szu was also included as part of a large population study by Chen et al. (2023a,b).This sample consisted of 78 H-poor SLSNe detected by ZTF over the span of 3 years.

Photometry
Observations of this target were obtained from a number of telescopes.Follow-up observations in gri were obtained with Las Cumbres Observatory (LCO) using a number of their 1m telescopes across multiple observatories in the network.After 350 days, deeper images were obtained with the ESO New Technology Telescope (NTT) in gri with the ESO Faint Object Spectrograph and Camera (EFOSC2) as part of the Extended Public ESO Spectroscopic Survey of Transient Objects (ePESSTO; Smartt et al. 2015).These images were reduced using the ePESSTO pipeline for the NTT images (Smartt et al. 2015), and the BANZAI2 pipeline for the LCO images.Photometry on these images was performed with the use of photometry sans frustration, a python wrapper for point-spread function (PSF) photometry using Astropy and Photutils (Nicholl et al. 2023).Zeropoints were calculated by cross matching sources in the field with the Pan-STARRS catalog (Flewelling et al. 2020).The photometry in  was template subtracted using archival PS1 images as templates.This was especially important in the late time photometry where we believe the host plays a significant contribution to the flux detected.
The Neil Gehrels Swift Observatory (hereafter Swift; Gehrels et al. 2004) began observations of the field of SN2019szu on the 21st November 2019 using the UV Optical Telescope (UVOT; Roming et al. 2005).SN2019szu was detected in all 6 optical/UV UVOT filters.Summed images were created by combining individual exposures taken during observations.Source counts were extracted from these summed images using a source region of 5" radius.Background counts were extracted using a circular region of radius 20" located in a source-free region.The count rates were obtained from the image lists using the Swift tool uvotsource.The count rates were then converted to AB magnitudes using the UVOT photometric zero points (Poole et al. 2008;Breeveld et al. 2011).
Data were also collected from the ZTF forced photometry server (Masci et al. 2019) and the ATLAS forced photometry server (Tonry et al. 2018;Smith et al. 2020;Shingles et al. 2021), both of which are performed on difference images.ZTF , , and ATLAS  band data points were binned together in a daily cadence whereas the  band data were binned together every 2 days to reduce noise.
Following the discovery of SN2019szu, we examined forced photometry at the location of the SN in Pan-STARRS1 and Pan-STARRS2 images obtained in survey mode (Chambers et al. 2016) from MJD 57362 (2015-12-06) onwards.Typically, 4 × 45 second exposures are obtained in survey mode in one of ,  or  filters on any given night, and photometric calibration and difference imaging is performed via the image processing pipeline (IPP; Magnier et al. 2020).The  filter is a broadband composite  +  + , with measured AB magnitudes roughly equivalent to those in the  band.The individual measurements for each nightly quad were stacked in order to improve signal to noise, and to obtain deeper upper limits in case of a non-detection.
All absolute magnitudes are calculated using the distance modulus and a simple K-correction of 2.5 log(1 + ).

Polarimetry
An epoch of polarimetry was also obtained on 17-01-2020 (31 days after peak in rest frame) using the Alhambra Faint Object Spectrograph and Camera (ALFOSC) instrument at the Nordic Optical Telescope (NOT) in the V band.The reduction and analysis is described in Pursiainen et al. (2023).The signal-to-noise ratio (S/N) is quite low compared to values traditionally needed for linear polarimetry.We find S/N ≳ 100 only for a small aperture size of ≤ 9 pixels, and falls quickly to 35 for aperture sizes above 20 pixels (larger apertures are necessary to account for any difference in point spread function between the ordinary and extraordinary beams).Although we measure an overall polarisation of  = 3.1 ± 1.3% for SN2019szu, compared to a an interstellar polarisation  ISP = 0.70±0.21%measured from a bright nearby star, the low S/N of the observation precludes a strong claim of polarized emission from SN2019szu.This is described in detail by (Pursiainen et al. 2023).

Spectroscopy
An initial spectrum of SN2019szu was obtained using the Spectrograph for the Rapid Acquisition of Transients (SPRAT) instrument on the Liverpool Telescope (LT).Spectroscopic follow-up observations of this target were then undertaken by ePESSTO using the NTT with EFOSC2 (Smartt et al. 2015).Most of these spectra were obtained with Gr#13 except for one epoch with Gr#16 on 2020-01-17 to extend the wavelength coverage redwards.The latter was averaged in the overlapping region with the Gr#13 spectrum obtained one day prior.Another spectrum was obtained on 2020-08-21 from MMT using the Binospec spectrograph covering a similar wavelength range to Gr#13 (Fabricant et al. 2019).A full breakdown of observations is given in Table 1.
All data were reduced using dedicated instrument-specific pipelines that apply de-biasing, flat-fielding, trace extraction, wavelength calibration and flux calibration using standard stars observed with the same setup.Spectra were then also flux corrected using both  and -band photometry.

Host Galaxy
The host galaxy of SN2019szu is a faint dwarf detected in the Pan-STARRS catalogue (PSO J002.5548-19.6923).There is no catalogued redshift or distance information and so estimates for the redshift were derived from the [OIII] 4959,5007 narrow host lines observed in the SN spectra.This gave a redshift of  = 0.213±0.0003which is also in agreement with H and H measurements from latter spectra.
The host is detected only in the Pan-STARRS  and  bands with a Kron magnitude  = 21.75 ± 0.08 mag (Flewelling et al. 2020).It is not catalogued in the NASA Extragalactic Database (Helou et al. 1991).In Schulze et al. (2018), the median SLSN host galaxy had   = −17.10± 1.45 mag.The host for SN2019szu has an absolute magnitude   = −18.38 ± 0.08 mag (at this redshift, the -band is similar to rest-frame ).This indicates the SN2019szu host galazy is well within the normal range of host luminosities (a proxy for masses) found in Schulze et al. (2018).
Hydrogen line ratios were measured based on the narrow emission lines observed in the late-time spectra of SN2019szu in order to estimate any reddening due to the host galaxy.This gave a ratio of H/H = 3.53 ± 0.13 for the spectrum at +211 days, and H/H = 4.52 ± 0.45 for the spectrum at +262 days.The H/H value could only be calculated for the +211 day spectrum and resulted in a value of H/H = 0.73 ± 0.17.Both of these values are above the expected ratios of H/H = 2.86, and H/H = 0.47 (Osterbrock & Ferland 2006).While the H/H ratio supports negligible extinction, the H/H ratio indicates significant reddening from the host, which is unexpected for a galaxy of this size.We can quantify the relation between Balmer decrement and colour excess  ( − ) described by Domínguez et al. (2013) giving  ( −) = 0.29 ± 0.05 and optical extinction of   = 1.2±0.3.This is much larger than typically inferred for SLSN host galaxies, which is generally <0.5 mag and averages ∼0.1 mag (Schulze et al. 2018).It is therefore likely that our measured Balmer decrement is unreliable, due to contamination from the SN spectrum.A low host extinction is also supported by light curve models (Section 4.2) and the lack of NaI D 5890, 5896 absorption which is used as a indicator of dust extinction (Poznanski et al. 2012).We therefore neglect host extinction in our analysis, as applying more host extinction didn't affect our spectral measurements significantly.The Milky Way extinction in the direction of SN2019szu is  ( − ) = 0.018 (Schlafly & Finkbeiner 2011), and this correction was applied to all spectra.
The metallicity of the galaxy was calculated from the  23 method outlined in Kobulnicky et al. (1999), which uses the fluxes of the [O II] 3727, [O III] 4959,5007, and H lines.As the metallicity appeared to be in the region between the metal rich and the metal poor branches of this relation, it was calculated for both branches.The metal rich branch yields a value of 12 + log(O/H) = 8.36 ± 0.07, and the metal poor branch a value of 12 + log(O/H) = 8.36 ± 0.05.Both of these values are in agreement with one another, and also consistent with work by Schulze et al. (2018) who found SLSNe host galaxies tended to have values of 12 + log(O/H) < 8.4.Limits on the metallicity can also be placed by measuring the ratio of [N II] 6583 to H. [N II] is not observable in SN2019szu and so using a 2 and 3 limit on this detection yields 12 + log(O/H) < 8.06 and 12 + log(O/H) < 8.14 respectively (Marino et al. 2013).This reinforces the low metallicity nature of the host.
The host for SN2019szu shows strong H and [O III] 5007 emission lines.In the +211 day Binospec sepctrum, these lines have equivalent widths of 148 Å and 127 Å respectively which are lower limits due to contamination from the SN continuum.This is similar to the sample of SLSN host galaxies studied by Leloudas et al. (2015), where ∼50% of the sample occurred in extreme emission line galaxies.
Using the star formation rate (SFR) diagnostics in Kennicutt (1998) provides two different measures of the SFR.The first uses the strength of H and gives a value of SFR = 0.4 M ⊙ yr −1 .The second uses [O II] 3727 and gives SFR = 0.3 − 0.4 M ⊙ yr −1 .Comparing to the sample of SFRs found in Leloudas et al. (2015) which ranged from 0.01 − 6.04 M ⊙ yr −1 , this is a typical star formation rate for SLSNe hosts.

Light Curve
The multi-band light curve of SN2019szu is shown in Figure 1 with a range spanning over 600 days in the observer frame.The event reaches a peak magnitude of  ,peak = −21.59± 0.06, which is close to the volume corrected median peak magnitude  peak = −21.31± 0.73 mag described by Lunnan et al. (2018b), and  ,peak = −21.14±0.75mag calculated by De Cia et al. (2018) for SLSNe, where the error represents the 1 spread.A more recent study by Chen et al. (2023a) found  ,peak = −21.54+1.12  −0.61 mag (not corrected for Malmquist bias).The main rising light curve is preceded by a plateau lasting 40 days in the rest frame.This commences at MJD 58700, and hovers around  = 21.6 mag before beginning to rise after MJD 58750.This is equivalent to an absolute magnitude of M  ∼ −18.7 mag and a luminosity ∼ 10 43 erg s −1 .This plateau was also observed in the -band.Figure 2 shows historical photometry in the -band showing deep upper limits down to 22.5 mag indicating this plateau is a emergent feature.Other SLSNe have shown signs of early excesses such as SN2006oz, in which a precursor plateau lasting 10 days was observed before the full monotonic rise.This was thought to be due to a recombination wave in the surrounding CSM consistent with the transition from O III to O II (Leloudas et al. 2012).Similarly SN2018bsz showed a slowly rising plateau lasting ∼ 30 days (Anderson et al. 2018).In LSQ14bdq the precursor peak was suggested to be caused by the cooling of extended stellar material (Nicholl et al. 2015b).In both of these cases the precursor events may have occurred after explosion, unlike in SN2019szu where the long plateau clearly precedes the explosion date inferred from the rising light curve.This feature in SN2019szu is also unusual due to its very flat nature over a long timescale which is not consistent with cooling material and may require an additional source of energy injection.
This light curve rise is captured well by ZTF and ATLAS.We estimated the date of maximum light to be MJD 58826 in the  band by fitting a low-order polynomial, giving a rise time of around 80 days.We take this to be the time of peak throughout.We caution that our fits to the peak are somewhat limited by SN2019szu entering solar conjunction around 50 days after discovery.Observations resumed when it became visible again around 90 days later.
The light curve appears to peak in the bluer bands first and has a rather flat shape or possible plateau at peak, which lasts longer in redder bands.In the  band this flattening lasts approximately 80 days.This is similar to other events such as SN2020wnt which also showed this plateau behaviour (Gutiérrez et al. 2022;Tinyanont et al. 2023).However, SN2020wnt also showed indications of an initially faster decline in bluer bands, which is not apparent in this event (Figure 1).The UV bands all show rising light curves until the solar conjunction, indicating a peak later in the evolution for these bands.
To parameterise the light curve peak, the exponential rise and decline timescales were determined giving an e-folding rise and decline of  −rise ∼ 48 days, and  −decline ∼ 100 days respectively.These timescales were determined by fitting low order polynomials to the g band light curve.In general SLSNe tend to have  rise ∼  decline /2 with slower evolving events also having slower rise times (Nicholl et al. 2015a).The rise and decline timescales of SN2019szu are consistent with this expectation.Some events however show a more skewed relation such as SN2017egm, which had a fast rise time  −rise ∼ 20 days, and slow decline with an estimated e-folding decline time of  −decline ∼ 60 days (Bose et al. 2018).
A small bump in the light curve can be seen around MJD 59100, corresponding to ∼200 days post peak.This is not unusual for SLSNe and a large fraction of SLSNe-I show these undulations (Hosseinzadeh et al. 2022;Chen et al. 2023b).This undulation is observed in the gr bands, which have sufficient coverage to observe variations at this phase.As discussed in Nicholl et al. (2016b); Inserra et al. (2017); Li et al. (2020), these undulations can be the result of collisions of the ejecta with shells or clumps of material.This interaction with circumstellar material can produce bluer colours ( − ) during the interaction due to heating of the ejecta (Chen et al. 2023b).An alternative theory is variation in the power output from a central engine such as the energy output from a magnetar (Metzger et al. 2018;Chen et al. 2023b).Although this would produce variability on timescales shorter than the observed bumps, this can be smoothed at early times if the variation is shorter than the photon diffusion time through the ejecta.This also implies that the undulations are more likely visible at later times as the ejecta becomes more transparent -this is supported by the fact that 73% of undulations are found post peak (Chen et al. 2023b).A central engine can also produce variations in the ejecta opacity via increased ionisation and hence more electron scattering even with constant energy input.It does this by creating an ionisation front that propagates outwards and breaks out from the front of the ejecta leading to a rebrightening (Metzger et al. 2014;Omand & Jerkstrand 2023).As the ejecta cools it can recombine, leading to a change in opacity again which could result in an observed undulation.Chen et al. (2023b) prefer this latter explanation for SN2019szu as this mechanism allows for a higher UV flux due to the decrease in bound-bound transitions for ionised metals.
Late time observations of SN2019szu (Figure 2 show deep upper limits in both  and , at levels below the precursor plateau observed.This supports the idea that this plateau is related to the SN event.

Colour
The colour evolution of SN2019szu is shown in Figure 3.The colour was calculated using Superbol, a python package that interpolates light curves in order to perform spectral energy distribution (SED) fits (Nicholl 2018).The  −  colour was calculated using our band data, and interpolated -band points from Superbol using a polynomial fit.The pre-peak colour shows a dramatic evolution to the blue, from an initial value of  −  = 0.05 ± 0.13 mag, dropping down to  −  = −0.51± 0.15 mag just before maximum light.This is not behaviour exhibited by other SLSNe, which tend to show a Colour evolution of SN2019szu in  −  compared a subset of SLSNe including SN2007bi, PTF12dam, and SN2015bn which all showed early emission of the 7300 Å line in their spectra (Gal-Yam et al. 2009;Nicholl et al. 2013Nicholl et al. , 2016b)).PS1-14bj is also included due to its extremely slow evolution (Lunnan et al. 2016).Other well observed SLSNe with published Superbol fits and more typical colour evolutions have also been plotted including SN2010gx, SN2011ke, LSQ12dlf, and SN2013dg (Pastorello et al. 2010;Inserra et al. 2013;Nicholl et al. 2014).All photometry has been K-corrected.
general colour increase, becoming redder over time shown by the sample of events in Figure 3.The colour of SN2019szu, exhibits only a very gradual change in colour after peak, hovering around  −  ∼ −0.3 mag.This is similar to other SLSNe which show a steady blue colour around peak before evolving dramatically towards the red, consistent with fast cooling after peak.Slow SLSNe such as SN2007bi, PTF12dam, and SN2015bn show a much more gradual colour evolution to the red (Gal-Yam et al. 2009;Nicholl et al. 2013Nicholl et al. , 2016b)), similar to the post peak evolution of SN2019szu.PS1-14bj is another event showing near constant colour but is overall a much redder event (Lunnan et al. 2016).
After the break in data, the colour appears to redden, consistent with slow cooling as seen in other SLSNe.At around 200 days post peak the colour once again begins to decrease.This turning point corresponds to the bump in the light curve mentioned in Section 3.2.This observation could be explained by an ionisation front breaking out of the ejecta (Metzger et al. 2014;Omand & Jerkstrand 2023), or interaction with circumstellar material (Inserra et al. 2017).Both of these mechanisms would heat the ejecta and therefore create a bluer colour.This is not something seen in other bumpy light curves, for example SN2015bn does not have a dramatic colour change around its bump at +50 days (Nicholl et al. 2016b).We also caution that the late time  −  colour could be affected by the strong emission lines from the host found in the -band wavelength range.

Bolometric Luminosity
The bolometric luminosity of SN2019szu was calculated using Superbol (Nicholl 2018), as shown in Figure 4. To do this the light curve in each band was interpolated to epochs with -band data.A constant colour relation was assumed for bands with fewer data points, or by fitting a low-order polynomial to capture the general shape of the light curve.The flux was corrected for time-dilation assuming a redshift of  = 0.213, and extinction corrected assum-ing a value of  ( − ) = 0.018 (Schlafly & Finkbeiner 2011).As discussed in Section 3.1, we assume negligible extinction from the host galaxy.The data was also K-corrected to shift the fluxes and effective filter wavelengths to their rest-frame values.The resulting spectral energy distribution (SED) was fit with a modified blackbody function suppressed below a cut-off wavelength (Nicholl et al. 2017b;Yan et al. 2018): Where   is the wavelength dependent flux,  is the wavelength,  is a nominal index for which we used a value of 3, and a cutoff wavelength of  0 =3000 Å was used.These values were chosen based on fitting Eq. 1 to each SED and averaging the best fits.Multicolour information was not available for the pre-explosion plateau and so this data was not included in the SED fitting.Fitting this equation using data from all bands produced the blackbody (BB) temperature (T) and radius (R).Superbol calculates the bolometric luminosity ( bol ) by integrating numerically under the observed SED points, and extrapolating the missing flux outside the wavelength range using the best-fitting absorbed BB model (Figure 4).The initial peak is fit using data from the uvw1, uvm2, uvw2, and Ugcwroi filters; the  and  bands were ignored due to their very sparse data points.Both bands were also very noisy and so did not provide much extra information compared to the much cleaner  band.After the first break due to solar conjunction, only gri data was used in order not to extrapolate too far in time in the other bands, instead opting to extrapolate further in wavelength.As we will discuss in section Section 3.3.2, the SED shape appears flatter than a blackbody in the redder bands.We therefore performed additional fits excluding the i band.These experiments showed no significant difference to the best-fit  bol , but did marginally affect  and  (Figure 4).We also caution that a late times (>200 days) the SED does not resemble a blackbody as evidenced by the nebular spectra in Figure 6, and so the bolometric luminosity should be treated with a degree of caution.
,  and  bol were calculated for each epoch, however, it is important to note that at later times the blackbody fits are more contaminated by nebular lines and as the event transitions from photopheric to nebular, the blackbody fit becomes less reliable.Beyond 300 days, we have detections in only one band, which is insufficient to measure a temperature and radius.However, assuming that the colour does not change dramatically, we are able to estimate the bolometric luminosity at the time of this last detection.
In Figure 5, we can see our event compared to some other well observed SLSNe with published Superbol data.Even compared to these other events, the slow nature of this event is apparent with only PS1-14bj having a comparable gradual decline.It is much harder to compare the rise to peak of these events as they are not all as well sampled, but we note that SN2019szu has a much faster rise to peak than PS1-14bj.
Integrating the bolometric light curve gives us an energy of  = 2.6×10 51 erg.This is a lower limit on the energy radiated by this event due to our finite sampling in time and wavelength.It also does not consider the energy released outside of time span covered by our photometry.The energy is consistent with other SLSNe which typically radiate ∼ 10 51 erg over their lifetimes (Gal-Yam 2012; Lunnan et al. 2018b).Purple stars represent the luminosity derived by integrating under the +211 and +262 day spectra.Middle: Blackbody temperature of SN2019szu calculated with black stars using all bands and blue points excluding .Bottom: Blackbody radius of SN2019szu calculated with black stars using all bands and blue points excluding .The radius is fit (purple line) to the initial rise up to -20 days relative to peak and indicates and expanding photosphere at ∼1200 km s −1 .After solar conjunction, points are unfilled to represent sparse photometry and large uncertainties.

Blackbody Temperature and Radius
Figure 5 shows the blackbody temperature and radius compared to other SLSNe, calculated using the method discussed in 3.2.2.Although this method allowed us to extrapolate the SED outside of the observed bands and measure the total luminosity, the continuum shape visible in our spectra clearly deviates from a blackbody even at early times as seen in Figure 6.This will be discussed in Section 3.3.2,but for our purposes here this means that forcing a blackbody fit on this data may result in temperature and radius measurements that are not physically meaningful.As the deviations become apparent at restframe wavelengths longer than 5500 Å, these SED fits were repeated with the i band removed.Removing this band did not produce any changes in temperature or radius larger than the 1 uncertainties on these quantities at early times, however some late time points between 200-300 days show significant deviations.The evolution excluding the i band shows significant fluctuations in short timescales for both temperature and radius, as well as much larger uncertainties.For this reason the fits including the  band were chosen for further analysis.After solar conjunction, points are unfilled to represent sparse photometry and large uncertainties.(Nicholl et al. 2016b(Nicholl et al. ,a, 2013;;Lunnan et al. 2016;Lin et al. 2023;Bose et al. 2018;Nicholl et al. 2017aNicholl et al. , 2014;;Gutiérrez et al. 2022) By looking at Figure 5, it is apparent that the temperature evolution is not consistent with other SLSNe, as it appears to increase over the first ∼100 days.This is in stark contrast with most other events, which tend to show a decreasing temperature (Chen et al. 2023a).However, it is consistent with the colour evolution found in Section 3.2.1 which indicated the SN became bluer with time around peak.At later times the temperature of SN2019szu has larger errors which makes it harder to constrain but it appears as though the SLSN stays much hotter than the other events, and also remains roughly constant rather than increasing or decreasing drastically.SN2020wnt and PS1-14bj both show increases in temperature post explosion (Gutiérrez et al. 2022;Lunnan et al. 2016).In SN2020wnt this increase is only before the peak of the light curve and lasts ∼50 days post explosion before decreasing in temperature.PS1-14bj instead shows a steady increase over the entire time frame, however both events stay cooler than SN2019szu.Lunnan et al. (2016) suggest that late-time heating, due to X-ray to UV breakout from a central engine, could explain this increase.We also note that the large discrepancy in the temperature evolution of SN2019szu compared to other events might suggest that indeed a blackbody is not an accurate representation of its SED. in a population paper of 78 SLSNe-I observed from March 17, 2018 to October 31, 2020.That paper highlighted the anomalous nature of SN2019szu (designated ZTF19acfwynw).They presented the same increasing temperature profile that we found in Figure 4. Chen et al. (2023a) provide possible explanations of CSM interaction providing an additional heating source, or ejecta being ionised by a central engine such as a magnetar.Alternatively the apparently rising tem-perature may be due to mismatch between the true SED and the assumption of a thermal spectrum.These possibilities will be discussed in Section 5.
Looking at Figure 4, we can see the radius appears to peak at  = (2.13 ± 0.22) × 10 15 cm at around 20 days before maximum light.Afterwards, the radius appears to either decrease very slowly or remain relatively constant up to the break in the data.Fitting this Phase (day) +211 initial rise up to -20 days with a linear function appears to show the photosphere expanding at  ∼ 1700 km s −1 .After the break the radius appears to have decreased significantly, where it continues a slow, gradual decline down to  = (4.31± 1.26) × 10 14 cm by 300 days.Overall the blackbody radius of SN2019szu appears quite compact compared to other SLSNe in our comparison sample, by a factor of few.Slower moving ejecta for SN2019szu could be one possible explanation for this difference.However, it may also be that the deviations of the spectrum from a true blackbody, as indicated by the apparently increasing temperature (and discussed further in the next section), lead to an underestimate of the true radius.This is supported by SN2018ibb which had an ejecta velocity of 8500 km s −1 but displayed a steady photophere radius of ∼ 5 × 10 15 cm over the course of 100 days (Schulze et al. 2023), comparable to the radius of other events.

Spectra
Figure 6 shows the evolution of the spectra of SN2019szu, starting at −30 days pre-peak, up to 262 days post peak.The first spectrum shows a featureless blue continuum and although quite noisy, the narrow [O III] doublet at 4959, and 5007 from the host galaxy is visible.Later spectra were obtained using NTT and MMT, a full breakdown of this is given in Table 1.The transition between a photospheric spectrum and a nebular spectrum occurs in the gap between the spectra at +31 and +211 days.Unfortunately this transition was not observed due to constraints from the Sun for ground-based telescopes.

Photospheric Spectra
The 'early' spectra here are defined as those taken before the break and showing photospheric absorption lines.In this case they run from −30 days to 30 days post peak in the rest frame.Narrow host-galaxy lines are seen with the [O III] emission lines at 4959, 5007, and H at 6563.The characteristic steep blue continuum and O II absorption lines associated with SLSNe are apparent in the 3000-5000 Å range, as well as typical Fe II and Fe III lines around 5000 Å blended with the O II absorption lines.Figure 7 shows the most prominent lines that can be seen in the spectra.An approximate SN velocity can be determined by measuring the blueshift of the O II absorption lines (Gal-Yam 2019b) yielding an ejecta velocity of  ej ∼ 4500 km s −1 .Chen et al. (2023b) found in their sample of events that the median O II derived velocity was 9700 km s −1 around time of maximum light.SN2019szu has a considerably slower velocty but still within the range measured by their sample.These spectra also include a broad emission line at 7300 Å, an unusual feature at such an early phase for any SN, though a similar feature has been seen as early as ∼50 days post-peak in a handful of slowly evolving SLSNe (Gal-Yam et al. 2009;Nicholl et al. 2016b;Inserra et al. 2017).There is only one other event SN2018ibb, in which this feature has been observed around peak (Schulze et al. 2023).We will discuss this in detail in Section 3.3.3.
In Figure 6 one can clearly see a lack of evolution in SN2019szu between −16 and +30 days.Other SLSNe have shown similar periods with minimal spectral evolution around maximum light; for example, SN2015bn maintained a constant spectrum dominated by O II and Fe III between at least −27 and +7 days (Nicholl et al. 2016b), but its spectrum then cooled and evolved quickly between 7 and 30 days.SN2019szu shows little evolution until at least 30 days, and if it did undergo any period of rapid cooling this was unfortunately unobserved while the object was in solar conjunction.This could be due to all of these spectra being obtained during the plateau at peak.The lack of line velocity evolution over this time also may favours the magnetar central engine for this event (Mazzali et al. 2016).The only feature that appears to change in this time period is the small bump that increases around 4340 Å.This wavelength is consistent with [O III] 4363, indicated by the first dashed grey line on Figure 6.This line has not been identified previously in such an early spectrum of a SLSN, and we explore this possible identification in Section 3.3.3.Any variation in [O III] 4959,5007 over this time frame is hard to determine due to its blend with Fe II lines.There is also an unidentified emission line that varies in flux around 3850 Å.If this was caused by Ca II H & K, it would require a blueshift of 7500 km s −1 , and if it was caused by [O II] 3727 it would require a redshift of 10500 km s −1 , both of which seem unlikely due to the lack of other features with similar blue/redshifts.

Continuum Shape
The shape of the continuum in the early time spectra shows a relatively flat region between 5500-7000 Å combined with a steep blue shape below ∼5500 Å.We attempted to fit the continuum with a single blackbody (Figure 8) and with the sum of two blackbodies of independent temperature.Our best fits indicate a singular source at ∼13000 K which captures most of the shape of the flat region.Adding a second component for a double blackbody did not improve this fit and in fact preferred assigning both components the same temperature.Moreover, the continuum shape between 5500-7000 Å does not resemble other SLSNe at a similar epoch, most of which can be well approximated by a simple blackbody.
A possible explanation for this could be an additional pseudocontinuum.In this case interaction with CSM produces a forest of narrow Fe II lines blended together, bluewards of ∼5500 (Inserra et al. 2016).This has been observed in various types of interacting supernovae including Type Ia-CSM and Type Ibn such as SN2014av and SN2006jc (Pastorello et al. 2016).Figure 8 shows a spectrum for SN2021csp, a SN Icn that has been argued to have exploded within H and He-poor CSM.This event shows a strong pseudo-continuum bluewards of ∼ 5000 Å attributed to a forest of Fe lines produced due to interaction with the surrounding medium (Perley et al. 2022;Fraser et al. 2021).Figure 8 Jerkstrand et al. (2017).Bottom: The photospheric spectrum of SN2015bn near maximum light, and a composite model where the 104_400 Cburn model has been added to this spectrum.The SN2015bn spectrum is used to represent a characteristic SLSN spectrum and the model spectrum is used to represent a separate emitting region producing the forbidden emission lines.This is compared with a photospheric spectrum of SN2019szu.
summing a blackbody at 20000 K to an arbitrarily scaled spectrum of SN2021csp.A hotter blackbody was needed as the best fit at 13000 K under predicted in the bluer wavelengths.This approximately recreates the unique continuum shape of SN2019szu, with a flat red region combined with the steep continuum in the blue.In reality, the component originating from the SN ejecta is more complicated than the simple blackbody used here (e.g. the SN spectrum contains O II absorption lines).However we use this composite model purely to show we can achieve a similar continuum shape to SN2019szu with an additional interaction component.
The pseudo-continuum in SN2019szu is apparent in the spectrum obtained at -16 days relative to peak.This suggests the SN was already interacting with CSM by this phase.

The 7300 Å Line
SN2019szu shows a prominent broad emission line at ∼7300 Å throughout its photospheric phase.While this line has been observed in some SLSNe in the late photospheric phase, it is already apparent in the earliest spectrum of SN2019szu that covers this wavelength range, meaning it is present at least 16 days before maximum light (Figure 9).The line itself does not evolve much over the course of our observations, with a similar full width at half maximum (FWHM) of ∼ 7000 km s −1 during the photospheric phase.This velocity drops to ∼ 5000 km s −1 FWHM in the late-time spectra which could be caused by a velocity gradient or change in optical depth.The total line flux stays relatively consistent during the early spectra at around 2 × 10 −14 erg cm −2 s −1 and decreases to 3 × 10 −15 erg cm −2 s −1 in the late-time spectra.
Other slow-evolving events such as SN2007bi, PTF12dam, SN2015bn, and LSQ14an have shown early emission of forbidden and semi forbidden lines ranging from 50-70 days post peak (Gal-Yam et al. 2009;Young et al. 2010;Nicholl et al. 2013Nicholl et al. , 2016b;;Inserra et al. 2017).One of the slowest events SN2018ibb displayed forbidden emission lines even earlier, apparent in its earliest spectrum at −1.4 days relative to peak (Schulze et al. 2023).Forbidden lines are formed when the radiative de-excitation dominates rather than collisional de-excitation.The conditions needed to form these forbidden  emission lines are generally not seen until the SN reaches its nebular phase, when material is much more diffuse, temperatures are lower, and the energy deposition from the power source is lower (Jerkstrand 2017).
In previous SLSNe, the line at 7300 Å was usually identified as [Ca II] 7291,7323.SN2018ibb is an interesting case as it is one the earliest examples of a SLSN showing nebular emission lines at 1.4 days before peak (Schulze et al. 2023).This was the earliest spectrum obtained of this object but still showed evidence for an emission line at 7300 Å that strengthened over time.The line profile shifted from a top-hat shape to bell shaped and also shifted by a few Å redwards over the first 100 days.Schulze et al. (2023) attribute this to the profile originally displaying [Ca II] and slowly becoming dominated by [O II] 7320,7330.In the case of LSQ14an, the line was strong in the earliest spectrum obtained 55 days after peak.Jerkstrand et al. (2017) attribute this line to emission from [O II] 7320,7330 as opposed to [Ca II] as identified in other SLSNe spectra, due to the strength of this line and the lack of [O I] 6300 emission throughout the nebular phase, suggesting oxygen was primarily in higher ionisation states throughout the ejecta.Inserra et al. (2017) suggest this line is a combination of both [O II] and [Ca II] in order for the line to match the widths of other [O III] lines present.The line also appears to be slightly asymmetric with a small narrow peak on the blue side of the profile which is similar to the shape of the line in SN2019szu (Figure 9).This similar asymmetry could indicate the presence of [Ca II] in the SN2019szu spectra as the blue peak is close to the rest-frame wavelength of [Ca II] 7291.However the low spectral resolution and signal-to-noise ratio makes it difficult using the line Other effects could also play a role in the asymmetry seen for the 7300 Å line in SN2019szu.Continuous scattering can be caused by free electrons or dust resulting in a profile with a blueshifted peak with a longer red tail (Jerkstrand 2017).Considering the first epoch in Figure 9, the base velocity of the line is ∼ 6000 km s −1 with a peak at ∼ −1500 km s −1 .This would put us in the regime where the electron scattering optical depth   = 2 − 3 (Jerkstrand 2017).The evolution of the line then would be consistent with a declining   .In the nebular phase   ≲ 1, leading to a diminished impact on the perceived blueshift of the line (Jerkstrand 2017).But seminebular lines can still emerge even with electron scattering depths of a few, and so the true velocity of the line may be slower.In some Type Ia SNe a similar effect of a blueshifted profile is thought to be caused by an aspherical explosion.Dong et al. (2018) show that in SN1991bg-like events, the central peak can be shifted off-center by ∼1000 km s −1 .Alternatively the profile could be explained by the SN blocking the view of the fastest receeding material, thus resulting in a slight blueshift.
Although the 7300 Å line appears to be double peaked in the +30 day and +262 day spectra, we believe this is an artefact from the smoothing.This is apparent by looking at the line in the +211 day spectrum obtained using Binospec which has a better wavelength resolution and does not indicate a double peaked nature.However, double peaked features have been observed in SLSNe such as SN2018bsz which displayed a H profile with a strong redshifted peak, and a weaker blueshifted peak.An asymmetric, disk-like CSM structure was used to explain how this profile could form (Pursiainen et al. 2022).
SN2019szu is the first SLSN that shows the 7300 Å line in a signif-icantly pre-maximum spectrum, present alongside the characteristic O II absorption lines.At the early times seen in this event, the expected high radiation temperatures would imply that Ca II would be ionised and so lines from this species would not be observed.This is supported by the singly ionised oxygen lines evident in Figure 7; O II has a higher ionisation potential than Ca II.This suggests that [O II] may be a more likely explanation for an emission line at 7300 Å appearing so early in the SN evolution.We can test whether [O II] is a valid identification for this line by comparing to spectral models to predict other lines we expect to appear in these spectra.Nebular models from Jerkstrand et al. ( 2017 2014) used this method for SN2010mb, a SN with evidence for interaction with large amounts of hydrogen-free CSM.With this in mind we identify the emission line at 7300 Å as [O II] 7320,7330 with a net blueshift indicating a velocity ∼ 1500 km s −1 .We will explore the physical implications of this line in Section 4.1.

Late Spectra
At late times the spectrum has transitioned into its nebular phase, defined by the broad emission lines and lack of thermal continuum.The spectra of SN2019szu at these phases are much more consistent with other SLSNe, as can be seen in Figure 11.Prominent lines (labelled in Figure 7 One notable difference is the lack of [O I] 6300 emission in SN2019szu.This line is visible in most of the spectra of the other SLSNe at a similar phase and so reflects a high and persistent degree of oxygen ionisation in SN2019szu.In the sample of SLSN nebular spectra analysed by Nicholl et al. (2019a), three out of 12 events showed evidence for weak [O I] 6300 and with [O II] dominating the 7300 Å region.This could be due to runaway ionisation which is believed to occur in the magnetar central engine scenario for low ejecta masses and a high power pulsar wind nebula.The result of this is a sharp switch of the spectra from O I dominated to O II/O III dominated.In the O II/O III dominated space this leads to a suppression of [O I] 6300 emission as seen in SN2019szu (Jerkstrand et al. 2017;Omand & Jerkstrand 2023).Alternatively, this line could emerge at later phases even for events with massive ejecta.SN2015bn and LSQ14an are modelled to have ≳10 M ⊙ of ejecta and show little [O I] emission at early on in their nebular phase (+243 and +111 days respectively).But both events show much stronger emission after 300-400 days (Jerkstrand et al. 2017;Nicholl et al. 2016a).We do not have spectra for SN2019szu at such a late phase. The

MODELLING AND INTERPRETATION
We have demonstrated that SN2019szu is a bright and slowly evolving SLSN, which displays surprising and persistent [O II] and [O III] emission lines even in early spectra obtained 16 days before maximum light.Our goal for this section is to determine the physical parameters of the oxygen emitting region and the ejecta overall, in order to understand how forbidden emission can arise at a time when the density in the ejecta is typically expected to be too high.
The typical electron density of expanding ejecta can be written as: where  is the mean atomic weight,  ej and  ej are the ejecta mass and velocity respectively,   is the electron fraction,  is the time from explosion, and f is the filling factor (Jerkstrand 2017).The critical density for the [O II] line is  crit ≈ 10 7 cm −3 (Appenzeller & Oestreicher 1988).Above this density, collisional de-excitation dominates rather than radiative de-excitation, therefore suppressing the formation of forbidden lines.Assuming the ejecta is mostly oxygen, with a velocity of ∼ 10000 km s −1 , at 80 days after explosion this results in   ∼ 10 10 ( ej /M ⊙ ) cm −3 .For typical SLSN ejecta masses of around 10 M ⊙ (Blanchard et al. 2020), this is several orders of magnitude greater than the [O II] critical density.This suggests that [O II] emission is more likely to arise from lower density material, presumably expelled before the first detection of SN2019szu.The early appearance of these lines requires material close to the SN site, therefore excluding material other than that originating from the SN or its progenitor.This argument still holds even if the 7300 Å line has been misidentified, since the critical density for [Ca II] is similar to that for [O II].As discussed in Section 3.3.2, the unusual shape in the spectral continuum may be attributed to interaction, providing further evidence for a (presumably H-poor) CSM close to the explosion site.

Models for the Oxygen Emission Lines
Jerkstrand et al. ( 2017) investigated the emission from long-duration SLSNe during their nebular phase using models from their SUMO code.These models explore oxygen-rich compositions and single zones, as this allowed exploration of parameters agnostic to the powering mechanism (Jerkstrand et al. 2017).Different compositions included a pure oxygen zone (pureO), and carbon burning ashes (Cburn).Here we choose the latter series to explore fully, as it is more physically motivated.Element abundances for this model were taken from the ONeMg zone in the Woosley et al. (2007) models assuming a star  ZAMS = 25 M ⊙ collapsing into a supernova.The model assumes 100 randomly distributed spherical clumps with vacuum in between, at 400 days post explosion.The ejecta is also assumed to have been travelling at a constant velocity of 8000 km s −1 .The energy deposited is assumed to be from high energy sources such as gamma rays.Each model had varying ejecta mass ( ej ), filling factor (  ), and energy deposition ( dep ).
A small subset of models were chosen from the grid of parameters to explore based on the line strengths of [O I] 6300, 6364 compared to [O II] 7320, 7330, where any models with visible [O I] emission compared to [O II] were rejected.This is motivated by Figure 6, where these particular lines are isolated from other lines in the late-time spectra of SN2019szu, and the lack of [O I] emission is clear.In the models, a lack of [O I] can occur due to near-complete oxygen ionisation.A full breakdown of models passing our selection is given in Table 2.The 104_400 model also creates [O III] lines at 4363 Å, 4959 Å, and 5007 Å which when scaled and added to an early spectrum of SN2015bn in Figure 10 recreates the line profile seen in SN2019szu around 5000 Å.The 4363 Å line also lines up with the emission feature seen in Figure 6 that appears to increase in strength over time.This provides solid evidence that these additional features can be explained by oxygen rich material.All of these oxygen lines also have a similar blueshift derived from the peak of the [O II] 7320,7330 and [O III] 4363 lines, which indicates a velocity of  = 1500 km s −1 .
We now seek to determine the heating rate and the density of oxygen needed to explain the [O II] line luminosity in SN2019szu.In the sample of four models, the luminosity of the 7300 Å line is roughly proportional to the energy deposited.From this we can estimate the energy deposition that gives the correct line luminosities for the SN2019szu spectra.The luminosity of the [O II] line in SN2019szu is ∼ 1 × 10 42 erg s −1 at early times and ∼ 4 × 10 41 erg s −1 in our nebular spectra, which would require  dep, early ∼ 4 × 10 42 erg s −1 and  dep, late ∼ 1 × 10 42 erg s −1 assuming the [O II] luminosity to be directly proportional to deposition.The required deposition at early times is also outside the range explored by Jerkstrand et al. (2017) but only by a factor of 2. Extrapolating to higher energies is not trivial, as changing the deposition does not only affect the line luminosities; it will also influence the temperature and ionisation of the ejecta, and hence the line ratios.In order to keep the ionisation state constant for different energies and masses, the following relation must hold: where  is the recombination rate, and   is the electron number density. is the ionisation rate per particle, defined as: Here  ion,OII is the fraction of deposited energy used in ionising O II to O III,  dep is the energy deposition per unit volume, and  is the ionisation potential.We can approximate this as  ion,OII ∼  ion  OII , where  OII is the fraction of O II in the gas.We assume   to be proportional to the oxygen number density   , which is then proportional to the density of the material  for our sample of models.This is motivated by the grid of models which show a ratio of   /  ∼ 1 − 1.3 for the subset used (Jerkstrand et al. 2017).If  is not too strong a function of temperature, we can assume it will be constant for all of these models with the same composition.This leads to the relation which can then be used to calculate the density of the oxygenemitting material needed in order to maintain a constant ionisation fraction for a given energy deposition.The fraction of energy going into ionisations,  ion , decreases with higher energy depositions, but the SUMO models show that this function varies very slowly and we may assume it to be constant in the regime we are working in.
To use the relation in Equation 5, we first identified models that produced lines consistent with SN2019szu.This was done by comparing the line ratios of [O II] 7320, 7330, [O III] 4363, and [O III] 4959, 5007 in the model series to our observed spectra.This is shown in Figure 12 for the spectrum at +262 days.This method was only applicable to the final two nebular spectra in Table 1, as the lines were not isolated enough in earlier spectra.Both spectra indicated  dep ≃ 1 × 10 42 erg s −1 provided roughly the correct ionisation balance for a 3 M ⊙ model with a filling factor  = 0.1.This model has  model = 6.7 × 10 −16 g cm −3 .We verified that applying a significant host galaxy extinction (Section 3.1) did not change the line ratios significantly and so we disregarded this effect.
We can thus scale the deposition to match the observed spectra, and determine the corresponding density to match the inferred ionisation state.Using  dep, early and  dep, late , this resulted in  early ∼ 1 × 10 −15 g cm −3 and  late ∼ 7 × 10 −16 g cm −3 .However, as the ionisation state cannot be constrained directly in the early spectra, the early measurements assume a constant ionisation state over the course of the SN lifetime, and so  early should be taken as a rough guide only.This density will be used further in Section 4.3 to constrain the CSM parameters.We also note that our measurement of the late time [O II] 7320, 7330 line may be contaminated by nebular emission from the SN ejecta.Based on Equation 5, a small change in  dep would result in an even smaller change in the density as it scales with √︁  dep , and so this calculation is not very sensitive to contamination from the ejecta.We can use the early time energy estimate as a check as we would not expect to see contamination during the photospheric stage.We can see that the density estimates at early and late times are within a factor of a few from each other.
Another check is the line flux remaining relatively constant over the course of its evolution and so we believe the line flux at times is dominated by CSM interaction.

Magnetar Model and Ejecta Mass Estimate
Modelling the light curve allows us to estimate the mass and velocity of the SN ejecta.Semi-analytic magnetar-powered models have been extensively tested in the literature and shown to able to reproduce the light curves of most observed SLSNe (e.g.Inserra et al. 2013;Nicholl et al. 2017b;Chen et al. 2023a).Although we have found evidence for CSM interaction in this event through the material needed to produce the 7300 Å line, the low density of this material suggests that it is likely not sufficient to power the full luminosity at the bright peak of the light curve, so multiple power sources may be at play.Nicholl et al. (2015a) modelled 24 SLSNe light curves using formulae detailed by Chatzopoulos et al. (2012), and showed that densities of 10 −12 g cm −3 were needed to match the rise-decline relation seen in SLSN light curves.Although these analytic CSM models are widely used (Chatzopoulos et al. 2012), the complicated geometry of the interaction and additional flexibility from separate ejecta and CSM density profiles can make the results difficult to interpret, and lead to mass estimates that are quite discrepant with hydrodynamical models (e.g.Sorokina et al. 2016;Nicholl et al. 2020;Suzuki et al. 2021).We therefore restrict our mass estimates to the magnetar model only, but note that the additional presence of CSM introduces an additional systematic uncertainty.We quantify the extent of the CSM contribution in Section 5.
To constrain the magnetar parameters that would be needed to power SN2019szu, and to estimate the ejecta mass and velocity, we employed the Modular Open Source Fitter for Transients (mosfit) (Guillochon et al. 2018;Nicholl et al. 2017b).mosfit is a fully Bayesian code that fits physical models to the multi-band light curves of transients.The plateau was excluded from this fitting as it is believed to be a result of pre-explosion activity and the magnetar model would not be able to fit such a feature.For the magnetar model the assumption is the energy input at time  is given by: Where  mag is the rotational energy of the magnetar,  mag is the timescale it spins down on. is the spin period of the magnetar and  ⊥ is the perpendicular component of the magnetic field to the spin axis.
The default mosfit priors were edited to account for the specifics of this event.In initial fits the value of A V,host railed against the upper end of the default prior, with A V,host approaching 0.5 mags, inconsistent with a dwarf host galaxy.The parameter fit in the models for this is the hydrogen column density ( H ), which is related to A V,host by  V =  H /1.8 × 10 21 .We therefore fixed  H to 10 16 cm −2 to ensure that an unrealistic extinction did not bias our fits.
The prior for scaling velocity ( ej ) was set to be a Gaussian distribution with mean = 5000 km s −1 and standard deviation = 3000 km s −1 , based on the blueshift of the O II lines in the photospheric spectra as described in Section 3.3.1.This approximation can be made for SLSNe around peak as the O II lines form in a region close to the photosphere.The shallow nature of the photosphere means it can be described using a single velocity (Gal-Yam 2019b), which can be used as a proxy for the ejecta velocity.The broad Gaussian also allows for a more typical SLSN velocity ∼ 10, 000 km s −1 as the lines may not have formed in the same location as the photosphere, and hence could have differing velocities.The prior on the minimum temperature ( min ) at late times was adjusted to cover a broader range from 10 3 to 10 5 K with a flat distribution in log space, the physical motivation being the unusual behaviour of the colour temperature of SN2019szu (Section 3.2.1).The prior for the magnetic field () was expanded to include lower values from 10 12 G to 10 15 G to allow for the slow evolution of SN2019szu.The optical opacity  was fixed to 0.15 g cm −2 , the median value found in the population study by Nicholl et al. (2017b).
Model fits are shown in Figure 13, and the derived posteriors are given in Figure 14.Overall the model fits cannot fully capture the the bumps and wiggles in the light curves which would require more nuanced models.It also under predicts the luminosity in some of the bluer bands such as , , and , whilst over predicting the  band luminosity at peak, and the  band at late times.This perhaps is a reflection of the unusual shape of the SED for this event.The magnetic field  = 0.37 +0.02 −0.02 × 10 14 G is at the low end for the population studied by Nicholl et al. (2017b) ( = 0.8 +1.1 −0.6 × 10 14 G), which is perhaps unsurprising given the slow decline together with Equation 8.A lower value of  results in energy deposition over a longer time frame, leading to a longer lived SLSN.We can also see that the spin period found for SN2019szu ( = 1.97 +0.16  −0.14 ms) is within the 1 range of the median value found Nicholl et al. (2017b) −1.2 ms).The spin period sets the luminosity scale, as  mag () ∝  mag / mag ∼ P −4 .As SN2019szu is a comparable luminosity to other SLSNe this is unsurprising.
Our estimated ejected mass ( ej = 30.2+4.5 −4.5 ) is also on the upper end of the general population, with a median and 1 distribution of  ej, median = 4.8 +8.1 −2.6 (Nicholl et al. 2017b).A larger population study by Blanchard et al. (2020) found this distribution of masses spanned 3.6 − 40 M ⊙ for the general population and so SN2019szu also lies on the upper end of this.The mass estimate is derived from the light curve width and diffusion timescale, and is not as dependent on the type of model chosen provided the power source is internal to the ejecta (e.g. 56Ni decay or magnetar engine).An interaction model is more complex, as the CSM can provide additional diffusive mass.However, since we have calculated the density of the CSM to be low relative to typical models for SLSNe, we would expect that the diffusion time is still dominated by the ejecta.For this reason, we expect our estimate of the total mass should be the right order of magnitude, even if interaction with the CSM is also contributing to the light curve.In SN2019szu, the large mass creates a longer diffusion time to contribute to the slow evolution.mosfit also provides an estimated explosion time relative to the first data point of  ∼ −37.96 +2.12  −2.33 days.Combining this with our time of maximum light, we estimate a total rise time of ≈82 days.

Constraints on Pre-explosion Mass Loss
From previous analysis in Section 4, we know that SN2019szu has a region of low density, oxygen rich material needed to produce the forbidden oxygen lines at early times.The density of this material was calculated to be  ∼ 10 −15 g cm −3 in Section 4.1.This density is too low to be from the SN itself and so this region must be comprised of pre-expelled circumstellar material.Figure 13.Magnetar model fits to SN2019szu using mosfit.Pre-explosion data is plotted but was excluded from the fits.Upper limits are indicated via inverted triangles and band offsets for display are: uvw2+5.5;uvm2+5; uvw1+4.5;U+4; B+2.5; g+2; V+1; c+0; w-1; r-2; o-3; i-4.
Using the velocities of the CSM and ejecta, we can constrain when this material was ejected before explosion, such that the ejecta catches the CSM prior to the first NTT spectrum.The width of the [O II] 7320, 7330 line did not decrease much over time, changing from 7000 km s −1 to 5000 km s −1 over the course of ∼300 days.Throughout this period, it maintains a constant blueshift of ≈ 1500 km s −1 .We assume that this blueshift corresponds to the CSM velocity and explore the consequences of this, though we acknowledge that it is also possible to achieve a blueshifted profile with electron scattering (Jerkstrand 2017).
The mosfit model provides a measurement for the scaling velocity  ej = 4700 km s −1 , which is consistent with the SN velocity estimated from the O II absorption lines measured in Section 3.3.1.Another velocity estimate comes from the increase in radius (Figure 4).Assuming a  =  relation gives  ej = 1700 km s −1 over the first 30 days which is similar to the velocity estimates from the [O II] emission lines.This could suggest that at early times there is a photosphere within the CSM, however it is important to remember the limitations of the SED fitting used to produce these radii.
The presence of the [O II] emission line at ≈66 days from the estimated explosion date indicates that the SN ejecta has already caught up with the CSM at this phase, giving an upper limit on the CSM radius of 2.6 × 10 15 cm.This places an upper limit on the ejection time which must be more recent than 200 days before the peak, or ∼120 days before explosion, for a CSM velocity of  CSM ∼ 1500 km s −1 and SN velocity  SN ∼ 4500 km s −1 .Mass ejection on a timescale of only months before the explosion is also supported by the precursor plateau in the light curve.
Using the estimated velocity and expansion time to determine the CSM radius at different times, we can also estimate the mass of emitting CSM.Using the approximate radius ≈ 2.6 × 10 15 cm at the light curve peak results in  early ∼ 0.1 M ⊙ , whereas in the late spectra at ≈250 days post peak we find  late ∼ 0.25 M ⊙ based on the densities found in Section 4.1 and assuming spherical, uniform CSM.However, contamination from nebular lines at late times could result in a slight increase in  late .Both masses being similar suggest that the the SN had interacted with the majority of this CSM by the time of our spectral observations.However if this slight increase in apparent mass is real, it could point to a thicker shell, meaning the SN is able to interact with more material as it expands.It is also important to note that the ionisation state derived from the late-time spectra may not strictly hold at earlier times.
Using our CSM mass estimate of around 0.25 M ⊙ and the velocity of this material we can calculate the total kinetic energy (KE).We obtain a value of KE ∼ 10 49 erg adopting the velocity from the [O II] and [O III] emission lines.This is a small fraction of the total energy obtained by integrating the bolometric light curve of E∼ 3 × 10 51 erg and so only a small fraction of the total energy would be released via these line emissions.

The progenitor of SN2019szu
We have found that the early emission lines in SN2019szu require 0.25 M ⊙ of CSM ejected less than 200 days before maximum light.Explaining the mass ejected from the star in such a short amount of time requires creative explanations, as stellar wind mass loss rates, even for Wolf Rayet stars, are too low to explain the entirety of this CSM (Mauron & Josselin 2011;Sander et al. 2022).Typical evolution of stars thought to undergo PPI and describes temperature and luminosity models for these events.Less massive cores lead to smaller PPI shell masses at times closer to core collapse.In the models with He core masses between 48-52 M ⊙ (CO cores of ≳ 40 M ⊙ ), the shells of material ejected in each pulse collide with one another and produce sustained luminosities above ∼ 10 43 erg s −1 from 40 days before explosion.This could provide an explanation for the pre-explosion plateau observed in SN2019szu and would be the first direct observation of PPI ejections if true.This subset of models also describe typical timescales for the onset of the first PPI ejections before core collapse of > 3.6 × 10 6 s.This is also consistent with the timescales obtained from the 7300 Å line, which constrained the time of mass ejection to within the last 120 days before time of maximum light.
However, this subset of models produce ejected shell masses in the region of 6-8 M ⊙ .This is significantly more massive than the ≈ 0.25 M ⊙ of CSM responsible for the oxygen emission lines in SN2019szu, which would instead be consistent with ≈ 30 M ⊙ cores.However, at lower masses the PPI eruptions occur only in the final days before explosion, inconsistent with the light curve plateau.Alternatively, the discrepancy in CSM mass could be explained by the loss of the He layers before the onset of PPI, resulting in less massive shells of CO material.In fact, this would likely favour the scenario where the outermost shells of the PPI ejections produce the forbidden oxygen emission lines.We suggest that a stripped ∼ 40 M ⊙ CO core is more consistent with the observations, and the low mass, low density material producing forbidden lines is the outermost material from the earliest pulse.This assumption is supported by Renzo et al. (2020), who show that the density of the PPI ejecta decreases as you move radially outwards for a 50 M ⊙ He core (∼ 40 M ⊙ CO core).Their work also suggests that these pulses may have velocities of a few thousand km s −1 , in alignment with our measured velocity from the 7300 Å line.Later successive pulses interact with each other to produce the plateau (Woosley 2017), and the ejecta interacting with these shells could be the origin of the pseudo-continuum.We also note that many factors could affect our calculated CSM mass, such as the clumpiness of the ejecta or how much of the material has been excited.We assumed spherical geometry and that the material is shock excited, but it could also be radiatively excited.Nonetheless, it is striking that our estimates for the energy and timescale of CSM ejection are consistent with existing PPI models.
It is also important to note that in all non-rotating PPI models the star collapses into a black hole, and none of these models reach the luminosity of SLSNe unless the initial star was rotating rapidly enough to create a magnetar (Woosley 2017).Indeed, Nicholl et al. (2015a) estimated that to power the observed long-lived SLSNe via CSM interaction, densities  CSM ≳ 10 −12 gcm −3 and masses comparable to the ejected  ej are probably needed.This is much higher than our inferred density and CSM mass for SN2019szu, and so it is not clear that this event can be solely powered by CSM interaction.PPI SN models reach peak luminosities of 3 × 10 43 erg s −1 from collisions with PPI shells (Sukhbold & Woosley 2016;Woosley 2017).Comparing this to the bolometric light curve of SN2019szu, this extreme would result in only a maximum of a 10% contribution at peak.If sustained at late times, this would have a more pronounced effect but falls within the error bars of the bolometric luminosity measured.
A more luminous interaction, ∼ 10 44 erg s −1 , may be possible if the final SN ejecta interacts with a massive PPI shell (in contrast to the fainter shell-shell collisions).The nebular emission lines from CO material suggest this is unlikely to occur in SN2019szu.In luminous PPI models (Helium cores ≥ 48 M ⊙ , or CO cores ≳ 40 M ⊙ ), the bulk of CSM mass is in the first, He-rich pulse, which may or may not be present in SN2019szu depending on whether the He layer was already stripped.The low-mass pulses of CO material come later.The low-density, oxygen-rich CSM producing the emission lines suggests that the CO shells have not been entirely overrun by the SN ejecta, so the SN ejecta cannot have reached any massive He-rich shell even if it was present.Interaction of a massive ejecta with CO shells of ≲ 1 M ⊙ will not contribute a large amount of luminosity, as the fraction of energy thermalised is roughly  ej /( ej +  CSM ) ≪ 1.We do however note that there are lots of uncertainties associated with the final stages of massive star evolution, and so other unknown scenarios involving eruptive mass loss from stripped stars could provide an alternative explanation for SN2019szu, where interaction could contribute a larger fraction of the luminosity.In any case, the oxygen-rich CSM, luminous early plateau, and consistency of the peak luminosity and our estimated ejecta mass with the engine-powered models of Woosley (2017), make PPI a compelling explanation for the massloss in SN2019szu.

Implications for other SLSNe and PPI candidates
The definitive presence of nearby CSM could explain other shared aspects of the SLSN population.The combination of a pre-explosion plateau and the very early appearance of nebular line emission in SN2019szu shows that nebular line emission during the photospheric phase of SLSNe can be an indicator of mass ejection shortly before explosion.The 7300 Å lines in the spectra of other SLSNe such as SN2007bi, PTF12dam, LSQ14an and SN2015bn may therefore reveal recent mass ejection, potentially driven by the PPI mechanism, in those events too.We note that these are all long-lived SLSNe, likely indicating large ejecta (and therefore progenitor core) masses.
Many SLSNe also show signs of a pre-maximum bump in their light curves (Leloudas et al. 2012;Nicholl et al. 2015b;Nicholl & Smartt 2016;Smith et al. 2016), which could potentially be explained by multiple shells of CSM produced by PPI eruptions (Woosley et al. 2007).However, if this was the case we might expect to see spectral evidence of this interaction.Early spectra of other SLSNe such as LSQ14bdq with pre-maximum bumps do not show evidence of broad emission lines (Nicholl et al. 2015b).Therefore, other explanations such as post-shock cooling of extended stellar material or a recombination wave in the ejecta may be more plausible explanations for some events (Nicholl et al. 2015b;Leloudas et al. 2012).However, no spectrum has been obtained during the pre-maximum bump of a SLSN, so it is also possible that spectroscopic signatures of CSM during the bump could be erased by the time of maximum light, e.g. if the ejecta has overrun the CSM.
Post-maximum bumps in the light curve could also indicate interaction with PPI shells.The long term light curve of SN2017egm shows multiple late time bumps as well as varying levels of decline (Figure 5) (Lin et al. 2023;Zhu et al. 2023).Lin et al. ( 2023) reproduce the evolution of this event using four distinct shells of CSM produced by PPI ejections from a star with a 48-51 M ⊙ He core.This is similar to the proposed He-core mass of SN2019szu and so we might expect to see similar features between the two events.SN2017egm also displays very little [O I] 6300 compared to [O II] 7320,7330 which is explained by the high temperatures and therefore ionised neutral oxygen for the event.However, the light curve shapes differ with SN2017egm having a sharper peak, and a faster overall decline.The authors attribute this sharp peak to forward shock between the SN ejecta and the nearest CSM shell.SN2017egm also displays a shorter rise time of ∼30 days from explosion to peak (Bose et al. 2018).Another key difference is the environments in which these events occurred.SN2017egm is unique in that it originated from a large galaxy with a high metallicty (Nicholl et al. 2017a).SN2019szu originated in small dwarf galaxy with much lower metallicity, an environment where a rapidly rotating magnetar may be more likely to form.The differences in these events could be explained by a difference in internal engine powering, as an engine is required to match the luminosity of SN2019szu (Woosley 2017), and cannot be explained by CSM interaction alone.
We can see clear observational evidence for circumstellar material in other SLSNe.In iPTF16eh, a Mg II resonance doublet was observed to change from blueshifted emission to redshifted emission over time (Lunnan et al. 2018a).This is explained by reflection of light from a detached shell of CSM surrounding the SN.This material had a velocity of 3300 km s −1 and was thought to have been ejected 32 years prior to the supernova explosion (Lunnan et al. 2018a), a much longer timescale than derived for SN2019szu.In the case of SN2019szu, we do not see this changing Doppler shift, as the ejecta has already collided with the CSM.However, the mechanism proposed to produce the CSM in iPTF16eh is also pair instability ejections, but in that case from a more massive progenitor with a He core mass of ∼ 51 − 53 M ⊙ , which experiences the PPI earlier before explosion.
Some SLSNe show evidence for interaction only at late times with the appearance of broad H emission.In SN2018bsz, this feature is multi-component and appears at ∼30 days, accompanied by other hydrogen lines.Pursiainen et al. (2022) explain this using highly aspherical CSM with several emitting regions.In iPTF15esb, iPTF16bad and iPTF13ehe, this feature emerged at +73, +97 and +251 days respectively, implying the progenitors lost their hydrogen envelope several decades before the SN explosion, leading to a neutral hydrogen shell (Yan et al. 2015(Yan et al. , 2018)).Yan et al. (2018) suggest this eruptive mass loss could be common in SLSN progenitors.iPTF15esb in particular had a triple-peaked light curve, which could be explained by shells of CSM with a total mass ∼ 0.01 M ⊙ .Collisions between these shells or between shells and ejecta could provide the excess luminosity to power the light curve undulations in iPTF15esb.Undulations during the declining phase in other SLSNe have also been attributed to interaction (Nicholl et al. 2016c;Inserra et al. 2017;Hosseinzadeh et al. 2022), though central engine flaring or ionisation fronts are an alternative explanation (Metzger et al. 2014;Hosseinzadeh et al. 2022).
Other energetic SNe have shown signs of interaction with oxygen rich material, including the unusual Type Ic SN2010mb (Ben-Ami et al. 2014).This event had a slowly declining light curve thought to be the result of interaction with ∼3 M ⊙ of CSM.This resulted in spectral features such as a blue quasi-continuum and a strong [O I]5577 emission line at late times.However, this emission line had a narrow core and and required densities of ∼ 10 −14 g cm −3 .This is more dense than the CSM surrounding SN2019szu but could be partially explained by the slower velocity of the CSM in SN2010mb at 800 km s −1 .Ben-Ami et al. ( 2014) also suggest PPI ejections could be the source of this material.Other events that fall into this emerging population of Ic-CSM include SN2022xxf and SN2021ocs, both of which show signs of interaction with H/He poor CSM (Kuncarayakti et al. 2023(Kuncarayakti et al. , 2022)).
Looking at the inferred properties of the progenitor can also provide contextual clues for pair-instability candidates.For example, the type I SN2016iet was estimated to have had a CO core mass of ∼ 55 − 120 M ⊙ prior to explosion (Gomez et al. 2019).This event was best modelled with interaction with ∼35 M ⊙ of CSM, ejected within the last decade before explosion.This leads to a mass loss rate of ∼7 M ⊙ yr −1 , much higher than the inferred rate for SN2019szu.These high masses coupled with the low metallicity host galaxy is within the regime of PPI or pair-instability supernovae and is consistent with PPI models by Woosley (2017).
In summary, the short timescale between explosion and observation of the 7300 Å line in SN2019szu supports the theory that the CSM producing this line was ejected by its progenitor very shortly before explosion.This is supported by the precedent set by other SLSNe which have shown evidence for eruptions close to explosion, albeit on a longer range of timescales.The relatively tight constraints on the CSM mass and timing of ejection make SN2019szu one of the strongest candidates for a PPI SN to date, and suggests that some of these events can form the engines required to reach superluminous magnitudes.

CONCLUSIONS
This paper presents extensive optical follow up of the SLSN SN2019szu.This includes spectra taken over nearly 300 days in rest frame, and photometry over 800 days.This event is one of the slowest evolving SLSNe to date with a rise time of ∼ 80 days from explosion to peak and an exponential decline time of ∼100 days.This event also displayed a pre-explosion plateau at an absolute magnitude   ∼ −18.7 mag lasting 40 days.
SN2019szu displayed a remarkably early forbidden emission line at 16 days before maximum light during its photospheric phase, the earliest we have ever seen such emission lines.Using models of nebular SN spectra from Jerkstrand et al. (2017), we were able to not only determine that this line at ∼ 7300 Å originated from [O II] 7320, 7330 but also deduce parameters of the material of origin.We found that SN2019szu had at least ∼ 0.25 M ⊙ of H-poor and O-rich material with a density of ∼ 10 −15 g cm −3 .
The spectra of this event also showed a steep continuum in the blue, combined with a relatively flat continuum redwards of ∼ 5500 Å.This unique spectral shape was not well fit by a simple blackbody.Instead we showed that this shape could be recreated by combining a hot blackbody spectrum with that of an interacting SN.Combined with the bumpy light curve, sustained blue colours, and the O-rich CSM needed to produce the [O II] line at early times, this provides strong evidence that SN2019szu was interacting with nearby CSM.In order for the interaction to occur by 16 days before maximum light, it must have been ejected less than 120 days before explosion (assuming a CSM velocity of 1500 km s −1 based on the blueshift of the emission lines), suggesting that mass ejection is also responsible for the light curve plateau.
We conclude that producing ∼ 0.25 M ⊙ of hydrogen-poor CSM close to the time of explosion is not feasible using known mechanisms, such as stellar winds or eruptions from luminous blue variables.Instead we suggest pulsational pair-instability (PPI) ejections are a promising possibility.The PPI mechanism also can explain the lack of H and He in this CSM.PPI models from a stripped ∼ 40 M ⊙ CO core are consistent with our estimated CSM energetics and ejection timescale, the duration and luminosity of the preexplosion plateau, and the estimated ejecta mass from the SN light curve.
The detailed study of SN2019szu introduces a new observational approach that can be used to find signatures of PPI interactions.Early observations of nebular emission lines alongside the characteristic O II absorption lines could be used to probe the structure of SLSNe.Obtaining these observations as soon as possible after explosion could help provide stricter constraints on when and how much mass is ejected during these PPI ejections.Observing preexplosion activity will also provide more information on the pro-genitors with spectroscopic observations during this time helping to unravel the composition and velocity of this material.As growing numbers of SLSNe show evidence for interaction with CSM, adopting this approach will also help answer questions about the explosion mechanisms involved.This will be especially useful in future survey telescopes such as the Vera Rubin Observatory, which will be able to detect precursor activity in time for more detailed follow-up.This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 4 .
Figure 4. Measurements from fitting blackbodies to SEDs constructed from Superbol.Top: Bolometric luminosity of SN2019szu constructed using uvw1, uvm2, uvw2, and Ugcwroi up until the initial break (gcroi thereafter).Purple stars represent the luminosity derived by integrating under the +211 and +262 day spectra.Middle: Blackbody temperature of SN2019szu calculated with black stars using all bands and blue points excluding .Bottom: Blackbody radius of SN2019szu calculated with black stars using all bands and blue points excluding .The radius is fit (purple line) to the initial rise up to -20 days relative to peak and indicates and expanding photosphere at ∼1200 km s −1 .After solar conjunction, points are unfilled to represent sparse

Figure 7 .
Figure 7. Spectra of SN2019szu with the main features labelled.Phase given in rest frame days with respect to maximum light in the -band.Top: Photospheric spectra with predominantly SLSNe absorption lines and a visible emission line attributed to [O II] at 7320.The narrow components of the [O III] lines and the narrow H are attributed to the host-galaxy.Dashed vertical lines indicate the O II absorption lines (Quimby et al. 2018), if blueshifted with a velocity of 4500 km s −1 .Bottom: Nebular spectra with SLSN emission lines.Narrow lines are attributed to the host-galaxy.Y axis cut for clarity.

Figure 8 .
Figure8.Comparisons between the spectrum of SN2019szu at -16 days as well a smoothed spectrum of interacting Type Icn SN021csp at +53 days, with a strong pseudo-continuum(Perley et al. 2022).In orange is the best fit blackbody to SN2019szu with a temperature component  ∼ 13000 K, shifted vertically to align with the flat red continuum.In blue is a scaled blackbody with a temperature component  ∼ 20000 K.A composite spectrum is also shown created by summing the 20000 K blackbody and SN2021csp spectrum (representing the interaction component) which can recreate both the flat red continuum and the steep blue continuum shape.

Figure 10 .
Figure10.Comparisons of the SN2019szu with model spectra.Top: Photospeheric and nebular spectra of SN2019szu compared with a scaled version of the 104_400 Cburn model fromJerkstrand et al. (2017).Bottom: The photospheric spectrum of SN2015bn near maximum light, and a composite model where the 104_400 Cburn model has been added to this spectrum.The SN2015bn spectrum is used to represent a characteristic SLSN spectrum and the model spectrum is used to represent a separate emitting region producing the forbidden emission lines.

Figure 11 .
Figure 11.Comparison spectra of SN2019szu compared to a selection of other late time, nebular spectra of SLSNe showing the 7300 Å feature.Vertical dashed lines indicate positions of [O III] 4363, [O III] 4959,5007, [O I] 6300,6364, and [O II] 7320,7330.Phases given with respect to time of maximum light.Numbers in brackets represent the phase divided by exponential decline timescale for a more direct comparison of phase.Spectra are plotted by relative strength of [O I] and the emission line 7300 Å.
profile alone to determine to what extent [Ca II] or [O II] may be contributing.SN2007bi also showed a strong emission line at 7300 Å in its earliest spectrum at 50 days post-peak (Gal-Yam et al. 2009), though in this case it exhibited a more typical nebular phase dominated by [O I] 6300.
) predict emission from [O III] 4959, 5007, and [O III] 4363 to also appear if strong [O II] is present.Looking at Figure 7, one can see evidence for these lines with the same blueshift as the [O II] emission line.This is even more apparent in Figure 10: adding the nebular model from Jerkstrand et al. (2017) showing the [O II] and [O III] lines to a photospheric spectrum of SN2015bn greatly improves the match with SN2019szu in the areas with the [O II] and [O III] lines.It also recreates the profile seen around 5000 Å due to the blend of [O III] with Fe II/III.This method of combining photospheric and nebular spectra assumes the nebular component is transparent enough to not impede our view of the SN, and has been used before in the literature.In particular Ben-Ami et al. ( ) include [O II] 3727, [O II] 7320,7330, Ca II 3934,3969, [O III] 4363, [O III] 4959,5007, Mg I] 4571, and a blend of iron lines around 5200 Å including [Fe II] 5250.
[O II] and [O III] emission lines, previously blended with broad absorption lines of O II and Fe III, are now much more easily isolated, confirming our earlier identification of [O II] 7320,7330.

Figure 12 .
Figure 12.Line ratios of [O II] and [O III] lines in model Cburn spectrum with M=3 M ⊙ and f=0.1(Jerkstrand et al. 2017).Dashed horizontal lines indicate line ratios found in the +262 day spectrum for an unknown energy deposition.Our data matches well with all the line ratios at  dep = 1 × 10 42 erg s −1 .
The early light curve of SN2019szu in  and  from Pan-STARRS forced photometry, and  and  from ZTF forced photometry.All magnitudes are in given in the AB system and 3 upper limits are indicated via inverted triangles.Phase given in rest frame days with respect to maximum light in the -band.
Nicholl et al. (2019b)0-19 under the name ZTF19acfwynw(Bellm et al. 2019), with the data visible in the Lasair broker 1(Smith et al. 2019).Gaia also detected this transient on 2019-11-02 with an internal name Gaia19fcb(Wyrzykowski 2016).It was later classified as a SLSN-I byNicholl et al. (2019b)as part of the C-SNAILS survey at the Liverpool Telescope (LT).It was initially given this classification using the host galaxy redshift of  = 0.213 (based on using the [O III] doublet emission at 4959 Å and 5007 Å from the SLSN spectrum), and therefore an absolute magnitude  = −21 mag indicating a very luminous event.This redshift corresponds to a distance of  = 1060 Mpc assuming a Planck cosmology (Planck Collabora-

Table 1 .
Spectroscopic observations of SN2019szu.Phase is given in rest frame frame days with respect to the time of maximum light in the -band.* Identical arms.
Spectral evolution of SN2019szu, phase is given in rest frame days with respect to maximum light in the -band.Dashed vertical lines indicate blue shifted positions of ionised oxygen lines seen in the spectra.The first spectrum at −30 is presented as both a smoothed and unsmoothed version for clarity.Edges of of this spectrum are clipped to remove noisy edges.Spectrum at 30 days is a composite of spectra taken on 30 and 31 days.The spectrum at 211 days has been telluric corrected and also smoothed using a Savitzky-Golay filter to reduce noise.
Chen et al. (2023a)also analysed the ZTF light curve of SN2019szu shows a composite spectrum created by

Table 2 .
SUMO models used.All models assume ejecta 400 days post explosion with a velocity of 8000 km s −1 , and 100 randomly distributed spherical clumps.M is the mass ejected in each model, f is the filling factor, and E dep is the energy injected into the system.

Table A1 :
Photometric observations of SN2019szu.Upper limits are indicated by a 1 in the upper limit column.