11institutetext:  The Oskar Klein Centre, Department of Physics, Stockholm
University, Albanova University Center, 106 91 Stockholm,
Sweden22institutetext: The Oskar Klein Centre, Department of Astronomy,
Stockholm University, Albanova University Center, 106 91 Stockholm,
Sweden33institutetext: Heidelberger Institut für Theoretische Studien,
Schloss-Wolfsbrunnenweg 35, 69118 Heidelberg, Germany44institutetext:
Technische Universität München, TUM School of Natural Sciences, Physik-
Department, James-Franck-Straße 1, 85748 Garching, Germany55institutetext:
Max-Planck-Institut für Astrophysik, Karl-Schwarzschild Straße 1, 85748
Garching, Germany66institutetext: Max-Planck-Institut für Extraterrestrische
Physik, Giessenbachstraße 1, 85748, Garching, Germany77institutetext:
Department of Particle Physics and Astrophysics, Weizmann Institute of
Science, 234 Herzl St, 76100 Rehovot, Israel88institutetext: The Caltech
Optical Observatories, California Institute of Technology, Pasadena, CA 91125,
USA99institutetext: Finnish Centre for Astronomy with ESO (FINCA), 20014
University of Turku, Finland1010institutetext: Tuorla Observatory, Department
of Physics and Astronomy, University of Turku, 20014,
Finland1111institutetext: DTU Space, National Space Institute, Technical
University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby,
Denmark1212institutetext: Department of Physics, University of Oxford, Denys
Wilkinson Building, Keble Road, Oxford OX1 3RH, UK1313institutetext:
Astrophysics Research Centre, School of Mathematics and Physics, Queen’s
University Belfast, BT7 1NN, UK1414institutetext: Department of Astronomy,
University of California, Berkeley, CA 94720-3411, USA1515institutetext:
Birmingham Institute for Gravitational Wave Astronomy and School of Physics
and Astronomy, University of Birmingham, Birmingham B15 2TT,
UK1616institutetext: Nordita, Stockholm University and KTH Royal Institute of
Technology Hannes Alfvéns väg 12, 106 91 Stockholm, Sweden1717institutetext:
Division of Physics, Mathematics and Astronomy, California Institute of
Technology, Pasadena, CA 91125, USA1818institutetext: INAF - Osservatorio di
Astrofisica e Scienza dello Spazio, via Piero Gobetti 93/3, 40129 Bologna,
Italy1919institutetext: Physics Department and Tsinghua Center for
Astrophysics (THCA), Tsinghua University, Beijing, 100084,
China2020institutetext: Gemini Observatory / NSF’s NOIRLab, Casilla 603, La
Serena, Chile2121institutetext: MIT-Kavli Institute for Astrophysics and Space
Research, 77 Massachusetts Ave., Cambridge, MA 02139, USA2222institutetext:
Department of Astrophysics, University of Vienna, Türkenschanzstraße 17, 1180
Vienna, Austria2323institutetext: Cosmic Dawn Center (DAWN), Denmark; Niels
Bohr Institute, Copenhagen University, Jagtvej 128, 2200, Copenhagen N,
Denmark2424institutetext: Department of Astronomy, Cornell University, Ithaca,
NY 14853, USA2525institutetext: Cardiff Hub for Astrophysics Research and
Technology, School of Physics & Astronomy, Cardiff University, Queens
Buildings, The Parade, Cardiff, CF24 3AA, UK2626institutetext: Tuorla
Observatory, Department of Physics and Astronomy, 20014 University of Turku,
Finland2727institutetext: Finnish Centre for Astronomy with ESO (FINCA), 20014
University of Turku, Finland2828institutetext: Astrophysics Research
Institute, Liverpool John Moores University, Liverpool Science Park, 146
Brownlow Hill, Liverpool L3 5RF, UK2929institutetext: Large Binocular
Telescope Observatory, University of Arizona, 933 N. Cherry Ave., Tucson, AZ
85721, USA3030institutetext: George Mason University, Department of Physics &
Astronomy, MS 3F3, 4400 University Dr., Fairfax, VA 22030,
USA3131institutetext: GSI Helmholtzzentrum für Schwerionenforschung,
Planckstraße 1, 64291 Darmstadt, Germany3232institutetext: Department of
Astronomy and Astrophysics, University of California, Santa Cruz, CA
950643333institutetext: INAF - Osservatorio Astronomico di Padova, vicolo
dell’Osservatorio 5, 35122 Padova, Italy3434institutetext: INAF - Osservatorio
Astronomico d’Abruzzo, Via M. Maggini snc, 64100 Teramo,
Italy3535institutetext: European Southern Observatory, Alonso de Córdova 3107,
Casilla 19, Santiago, Chile3636institutetext: Millennium Institute of
Astrophysics MAS, Nuncio Monsenor Sotero Sanz 100, Off. 104, Providencia,
Santiago, Chile3737institutetext: INAF - Istituto di Astrofisica Spaziale e
Fisica cosmica Milano (IASF), via Alfonso Corti 12, 20133 Milano,
Italy3838institutetext: Institute of Space Sciences (ICE, CSIC), Campus UAB,
Carrer de Can Magrans, s/n, 08193 Barcelona, Spain3939institutetext: Institut
d’Estudis Espacials de Catalunya (IEEC), 08034 Barcelona,
Spain4040institutetext: Astronomical Observatory, University of Warsaw, Al.
Ujazdowskie 4, 00-478 Warszawa, Poland4141institutetext: IPAC, California
Institute of Technology, 1200 E. California Blvd, Pasadena, CA 91125,
USA4242institutetext: Center for Astrophysics, Harvard & Smithsonian, 60
Garden Street, Cambridge, MA 02138-1516, USA4343institutetext: Las Cumbres
Observatory, 6740 Cortona Dr. Suite 102, Goleta, CA, 93117,
USA4444institutetext: Department of Physics, University of California, Santa
Barbara, CA 93106-9530, USA4545institutetext: Department of Physics,
University of California, 1 Shields Avenue, Davis, CA 95616-5270,
USA4646institutetext: Université de Lyon, Université Claude-Bernard Lyon 1,
CNRS/IN2P3, IP2I Lyon, 69622 Villeurbanne, France4747institutetext: Lunar and
Planetary Laboratory, University of Arizona, Tucson, AZ 85721,
USA4848institutetext: The NSF AI Institute for Artificial Intelligence and
Fundamental Interactions, USA

#  1100 days in the life of the supernova 2018ibb —
The best pair-instability supernova candidate, to date

Steve Schulze  11<EMAIL_ADDRESS>Claes Fransson 22 Alexandra
Kozyreva  33 Ting-Wan Chen  445566 Ofer Yaron 77 Anders Jerkstrand  22 Avishay
Gal-Yam 77 Jesper Sollerman  22 Lin Yan 88 Tuomas Kangas  991010 Giorgos
Leloudas  1111 Conor M. B. Omand  22 Stephen J. Smartt  12121313 Yi Yang (杨轶)
1414 Matt Nicholl  15151313 Nikhil Sarin  161611 Yuhan Yao  1717 Thomas G.
Brink  1414 Amir Sharon 77 Andrea Rossi 1818 Ping Chen  77 Zhihao Chen  1919
Aleksandar Cikota  2020 Kishalay De NASA Einstein Fellow2121 Andrew J. Drake
1717 Alexei V. Filippenko 1414 Christoffer Fremling 88 Laurane Fréour 2222
Johan P. U. Fynbo  2323 Anna Y. Q. Ho  2424 Cosimo Inserra 2525 Ido Irani  77
Hanindyo Kuncarayakti  26262727 Ragnhild Lunnan  22 Paolo Mazzali 282855 Eran
O. Ofek 77 Eliana Palazzi  1818 Daniel A. Perley  2828 Miika Pursiainen  1111
Barry Rothberg 29293030 Luke J. Shingles  3131 Ken Smith 1313 Kirsty Taggart
3232 Leonardo Tartaglia  33333434 WeiKang Zheng 1414 Joseph P. Anderson
35353636 Letizia Cassara 3737 Eric Christensen 4747 S. George Djorgovski 1717
Lluís Galbany  38383939 Anamaria Gkini 22 Matthew J. Graham 1717 Mariusz
Gromadzki  4040 Steven L. Groom  4141 Daichi Hiramatsu  42424848 D. Andrew
Howell 43434444 Mansi M. Kasliwal  1717 Curtis McCully 4343 Tomás E. Müller-
Bravo  38383939 Simona Paiano 3737 Emmanouela Paraskeva 4545 Priscila J. Pessi
22 David Polishook  77 Arne Rau 66 Mickael Rigault  4646 Ben Rusholme  4141

(Received 10 May 2023; accepted 12 October 2023)

Stars with zero-age main sequence masses between 140 and
$260\leavevmode\nobreak\ M_{\odot}$ are thought to explode as pair-instability
supernovae (PISNe). During their thermonuclear runaway, PISNe can produce up
to several tens of solar masses of radioactive nickel, resulting in luminous
transients similar to some superluminous supernovae (SLSNe). Yet, no
unambiguous PISN has been discovered so far. SN 2018ibb is a hydrogen-poor
SLSN at $z=0.166$ that evolves extremely slowly compared to the hundreds of
known SLSNe. Between mid 2018 and early 2022, we monitored its photometric and
spectroscopic evolution from the UV to the near-infrared (NIR) with 2–10 m
class telescopes. SN 2018ibb radiated $>3\times 10^{51}\leavevmode\nobreak\
\rm erg$ during its evolution, and its bolometric light curve reached
$>2\times 10^{44}\leavevmode\nobreak\ \rm erg\,s^{-1}$ at its peak. The long-
lasting rise of $>93$ rest-frame days implies a long diffusion time, which
requires a very high total ejected mass. The PISN mechanism naturally provides
both the energy source (56Ni) and the long diffusion time. Theoretical models
of PISNe make clear predictions as to their photometric and spectroscopic
properties. SN 2018ibb complies with most tests on the light curves, nebular
spectra and host galaxy, and potentially all tests with the interpretation we
propose. Both the light curve and the spectra require 25–44 $M_{\odot}$ of
freshly nucleosynthesised 56Ni, pointing to the explosion of a metal-poor star
with a helium core mass of 120–130 $M_{\odot}$ at the time of death. This
interpretation is also supported by the tentative detection of [Co ii]
$\lambda$ 1.025$\mu$m, which has never been observed in any other PISN
candidate or SLSN before. We observe a significant excess in the blue part of
the optical spectrum during the nebular phase, which is in tension with
predictions of existing PISN models. However, we have compelling observational
evidence for an eruptive mass-loss episode of the progenitor of SN 2018ibb
shortly before the explosion, and our dataset reveals that the interaction of
the SN ejecta with this oxygen-rich circumstellar material contributed to the
observed emission. That may explain this specific discrepancy with PISN
models. Powering by a central engine, such as a magnetar or a black hole, can
be excluded with high confidence. This makes SN 2018ibb by far the best
candidate for being a PISN, to date.

###### Key Words.:

supernovae: individual: SN 2018ibb, ATLAS18unu, Gaia19cvo, PS19crg,
ZTF18acenqto

##  1 Introduction

Observations of stellar nurseries (e.g. Krumholz et al. 2019), as well as
massive stars (e.g. Crowther 2007) and their fates (e.g. Filippenko 1997; Gal-
Yam 2017) have led to stellar evolution models of an ever-increasing
complexity (e.g. McKee & Ostriker 2007). These models also predict the
existence of stars with $\gtrsim 100\leavevmode\nobreak\ M_{\odot}$ (e.g.
Heger & Woosley 2002; Heger et al. 2003), which may have no analogues in the
local Universe (Mackey et al. 2003; Bromm & Larson 2004; Langer et al. 2007,
but see Brands et al. 2022), and exotic types of stellar explosions (Fowler &
Hoyle 1964; Rakavy et al. 1967; Woosley et al. 2007; Sakstein et al. 2022).

One of those predicted, yet not securely discovered object classes, is pair-
instability supernovae (PISNe). This SN class is produced by the thermonuclear
runaway of metal-poor stars with zero-age main sequence (ZAMS) masses between
140 and $260\leavevmode\nobreak\ M_{\odot}$ (Fowler & Hoyle 1964; Barkat et
al. 1967; Rakavy et al. 1967). When such a massive star dies, its helium core
would have grown to 65–$130\leavevmode\nobreak\ M_{\odot}$ (Heger & Woosley
2002). The combination of a relatively low matter density and high temperature
leads to the production of $e^{-}e^{+}$ pairs, reducing the radiation pressure
that supports the star against the gravitational collapse. As a result,
implosive oxygen and silicon burning produce enough energy to revert the
collapse and obliterate the entire star, leaving no remnant behind.

During the past 15 years, PISNe have been a focus of fundamental physics and
SN science. Stars with helium-topped cores slightly less massive than $\sim
65\leavevmode\nobreak\ M_{\odot}$ presumably leave black holes behind, and
stars whose helium-topped cores exceed $\sim 130\leavevmode\nobreak\
M_{\odot}$ are thought to collapse directly into black holes. In this
paradigm, there should be a dearth of black holes with masses between $\sim
50$ and $\sim 120\leavevmode\nobreak\ M_{\odot}$ (Farmer et al. 2019; Renzo et
al. 2020). Observations by the Laser Interferometer Gravitational-Wave
Observatory (LIGO) and the Virgo interferometer found tentative evidence for
the existence of such a drop in the black-hole mass function (The LIGO
Scientific Collaboration et al. 2020). A more recent study by the LIGO-Virgo
collaboration using the larger Gravitational-Wave Transient Catalog 3 shows
that the evidence of a mass gap at $\sim 50\leavevmode\nobreak\ M_{\odot}$ is
inconclusive (The LIGO Scientific Collaboration et al. 2021). However, this
could be due to the inclusion of binary black holes that formed through
dynamical channels involving repeated mergers rather than evidence for the
lack of a mass gap (e.g. Belczynski et al. 2020; Gerosa & Fishbach 2021, and
references therein). A further consequence of the PISN model is its peculiar
nucleosynthetic pattern. In the case of an extremely metal-poor star being
formed from the gas of PISN ejecta, its chemical composition should show a
strong variance between elements with an odd and even nuclear charge (Heger &
Woosley 2002; Umeda & Nomoto 2002). Xing et al. (2023) recently reported the
discovery of a very metal-poor star with such a chemical signature (for
additional candidates see also Aoki et al. 2014; Salvadori et al. 2019),
lending support to the first stars having been very massive and exploding as
PISNe.

Finding PISNe is one of the main challenges in the SN field. PISN models
predict that up to $\sim 57\leavevmode\nobreak\ M_{\odot}$ of radioactive 56Ni
are produced during the thermonuclear runaway (Heger & Woosley 2002). Such
high Ni-yield PISNe are thought to produce long-lived (rise times $>80$ days),
luminous ($M_{\rm peak}<-21$ mag) transients (Kasen et al. 2011; Kozyreva et
al. 2017) in the regime of superluminous supernovae (SLSNe; Quimby et al.
2011; Gal-Yam 2012, 2019b). Although the powering mechanism of SLSNe is
debated (Gal-Yam et al. 2009; Blinnikov & Sorokina 2010; Inserra et al. 2013),
numerous studies of both H-poor and H-rich SLSNe have revealed that nickel is
not the primary source of energy (e.g. Chatzopoulos & Wheeler 2012; Chen et
al. 2013; Inserra et al. 2013; Nicholl et al. 2017; Inserra et al. 2018;
Moriya et al. 2018b; Gal-Yam 2019b; Inserra 2019; Kangas et al. 2022; Chen et
al. 2023b). Yet, a few SLSNe had markedly broad and luminous light curves
similar to predictions of PISN models, for example, SN 1999as, SN 2007bi,
PTF12dam, PS1-14bj, and SN 2015bn (Hatano et al. 2001; Gal-Yam et al. 2009;
Nicholl et al. 2013; Chen et al. 2015; Lunnan et al. 2016; Nicholl et al.
2016b; Kozyreva et al. 2017). However, the published candidates either had
incomplete datasets, not long enough rise times, too high ejecta velocities,
too blue spectra, or exploded in galaxies with a metallicity that is too high
to conclusively argue for the discovery of a PISN (e.g. Nicholl et al. 2013;
Jerkstrand et al. 2017).

Not all PISNe are expected to be superluminous. The stars at the low-mass end
of the PISN range ($M({\rm ZAMS})\gtrsim 140\leavevmode\nobreak\ M_{\odot}$
equivalent to a He-core mass of $\gtrsim 65\leavevmode\nobreak\ M_{\odot}$ at
the end of their evolution) are thought to produce a small-to-no-mass of 56Ni
(Kasen et al. 2011; Kozyreva et al. 2014; Gilmer et al. 2017). The bright
supernova OGLE-2014-SN-073 is a candidate for a low-mass PISN, as its light
curve would require only $\sim 1\leavevmode\nobreak\ M_{\odot}$ of 56Ni
(Terreran et al. 2017; Kozyreva et al. 2018, but see also Dessart & Audit 2018
and Moriya et al. 2018c for alternative interpretations). However, also in
that case, there is still some tension between observations and theoretical
models due to the unknown explosion date and the lack of sufficiently early
data to search for the luminous shock-breakout predicted by PISN models
(Kozyreva et al. 2018). The Type I SN 2016iet is another candidate for a
lower-mass PISN. Gomez et al. (2019) concluded that the light curve could be
powered by a few $M_{\odot}$ of 56Ni. However, the SN was discovered only
around maximum light. Furthermore, its nebular spectra are very different to
the PISN model spectra from Mazzali et al. (2019) as pointed out by Gomez et
al. (2019). While the inferred nickel mass of a few $M_{\odot}$ and ejecta
mass of $\sim 64\leavevmode\nobreak\ M_{\odot}$ could point to a PISN
explosion, other powering mechanisms are also possible (Gomez et al. 2019).

In June 2018, the Zwicky Transient Facility (ZTF; Bellm et al. 2019b; Graham
et al. 2019) started to survey the northern sky every 2–3 nights in two
filters and detects thousands of supernovae every year (Fremling et al. 2020;
Perley et al. 2020). Until autumn 2021, we carried out a systematic survey for
SLSNe in ZTF (Chen et al. 2023a, b). SN 2018ibb, the slowest evolving SLSN in
our sample, has several properties that match predictions of PISN models.
Between mid 2018 and early 2022, we built a comprehensive photometric and
spectroscopic dataset covering the evolution from $t_{\rm max}$$-$93 to
$t_{\rm max}$$+$1000 rest-frame days to scrutinise SLSN and PISN models. In
this paper, we present this dataset along with our conclusions on SN 2018ibb’s
source of energy and progenitor. The paper is structured as follows: we report
the SN discovery in Section 2 and describe the observations in Section 3. In
Section 4, we derive the properties of SN 2018ibb’s light curve, spectra and
host galaxy and in Section 5 we contrast SLSN and PISN models with our
dataset. Finally, in Section 6 we summarise our findings and present our
conclusions on the nature of SN 2018ibb and its connection to PISNe.

Throughout the paper, we provide all uncertainties at $1\sigma$ confidence.
The photometry is reported in the AB system. We assume $\Lambda$CDM cosmology
with $H_{0}=67.8\leavevmode\nobreak\ {\rm km\,s}^{-1}\,{\rm Mpc}^{-1}$,
$\Omega_{\rm M}=0.308$, and $\Omega_{\Lambda}=0.692$ (Planck Collaboration et
al. 2016). Phase information is reported in the rest-frame with respect to the
$g|r$-band maximum ($t_{\rm max}$) at MJD=58455.

##  2 Discovery

Figure 1: False-colour image of the field when SN 2018ibb was bright (left)
and after it had faded below the host level (right). The SN position, marked
by the crosshair, is located $\sim 1$ kpc from the centre of its star-forming
dwarf host galaxy ($M^{\rm host}_{r}\sim-15.4$ mag, $M_{\star}\sim
10^{7.6}\leavevmode\nobreak\ M_{\odot}$). For more information about the host,
see Section 4.6. The false-colour image was built with STIFF version 2.4.0
(Bertin 2012).

SN 2018ibb, located at $\alpha$ = 04:38:56.950, $\delta=-$20:39:44.10 (J2000),
was discovered by the Asteroid Terrestrial-impact Last Alert System (ATLAS;
Tonry 2011; Smith et al. 2020) survey as ATLAS18unu on 10 September 2018 with
an apparent magnitude of $o=18.89$ mag (wavelength range 5600–8200 Å; Tonry et
al. 2018). Later detections were reported by the public northern sky survey of
the Zwicky Transient Facility (Bellm et al. 2019a) on 16 November 2018
(internal name: ZTF18acenqto), the Pan-STARRS Survey for Transients (Huber et
al. 2015) on 8 January 2019 (internal name: PS19crg) and the Gaia Photometric
Science Alerts transient survey (Hodgkin et al. 2021) on 4 July 2019 (internal
name: Gaia19cvo). A false-colour image of the field when SN 2018ibb was bright
and after it had faded is shown in Figure 1. Fremling et al. (2018a) initially
classified SN 2018ibb as a Type Ia SN on 5 December 2018 but retracted this
classification on 6 December 2018 and set a new classification to ‘supernova’
on 6 December 2018 (Fremling et al. 2018b). Pursiainen et al. (2018) obtained
a spectrum with the 3.58 m New Technology Telescope at La Silla Observatory
(Chile) as a part of the Extended Public ESO Spectroscopic Survey of Transient
Objects (ePESSTO; Smartt et al. 2015) on 14 December 2018 and classified SN
2018ibb as a H-poor SLSN at $z=0.16$.

##  3 Observations and data reduction

###  3.1 Supernova photometry

Our imaging campaign had three tiers: i) all-sky surveys with sufficient depth
and cadence to monitor the evolution from $t_{\rm max}$$-$93 to $t_{\rm
max}$+306 days; ii) dedicated follow-up campaigns to expand the wavelength
coverage to the UV and near-IR and to extend the light curve coverage to
$t_{\rm max}$+1000 days; and iii) smaller targeted campaigns to mitigate data
gaps, expand the wavelength coverage to the near-IR, and ensure a good flux
calibration of the photometric and spectroscopic data. Owing to the large
number of facilities involved in this effort, we present the details of each
campaign and the data reduction in Appendix A.

The ground-based photometry was calibrated with field stars from PanSTARRS1
(Chambers et al. 2016, PS1), the Dark Energy Survey (DES; The Dark Energy
Survey Collaboration 2005), the Dark Energy Spectroscopic Instrument (DESI)
Legacy Imaging survey (LS; Dey et al. 2019), and the Two Micron All-Sky Survey
(2MASS; Skrutskie et al. 2006). We applied known colour equations between
PS1/DES and Bessell/GROND/SDSS/ZTF filters (Finkbeiner et al. 2016; Drlica-
Wagner et al. 2018; Greiner et al. 2008; Medford et al. 2020) and
Lupton111https://www.sdss.org/dr12/algorithms/sdssubvritransform, to account
for differences in the filter response function. We applied the offsets from
Blanton & Roweis (2007) to convert all measurements to the AB photometric
system. The Swift/UVOT data were calibrated with zeropoints from the Swift
pipeline and converted to the AB system following Breeveld et al. (2011). The
HST photometry was done with a custom-made aperture photometry tool, based on
the python package photutils (Bradley et al. 2020) version 1.5, using an
aperture with a diameter of 0$\aas@@fstack{\prime\prime}$5 and calibrated
against tabulated zeropoints in pysynphot version 2.0.0 (STScI Development
Team 2013).

SN spectra are characterised by strong absorption and emission features that
evolve with time. This can lead to time-dependent colour terms between similar
but not identical filters (e.g. Stritzinger et al. 2002) and add a non-
negligible systematic scatter to the light curves if these differences are not
corrected. To illustrate this issue, we compute the synthetic magnitude in
ZTF/$g$, GROND/$g$ and EFOSC2/$g$ at $t_{\rm max}$ and $t_{\rm max}$+210
days222The GROND, ZTF and EFOSC2 $g$-band filters have an effective wavelength
of $4504,\,4723,5104$ Å and width of $1373,1282,788$ Å, respectively
(retrieved from the Spanish Virtual; Rodrigo et al. 2012, and references
therein).. At $t_{\rm max}$, the colour term between the EFOSC2/GROND and ZTF
filters is $-0.01$ and $+0.04$ mag, respectively, but at $t_{\rm max}$+210
days the differences increased to $-0.13$ and $+0.12$ mag. Since the EFOSC2
and GROND data cover the late-time evolution, the differences in the filters
would be well visible in the final light curve if they remained uncorrected.

To calibrate the various datasets into the same photometric system, we defined
a set of reference filters consisting of the Swift filters, ZTF/$gr$,
GROND/$izJH$ and 2MASS/$K$. Then, we extracted synthetic photometry of all
ground-based filters used in our campaign from the Keck and VLT spectra
(Section 3.3), which were obtained in clear/photometric conditions, and
measured the expected colours with respect to our reference filter system as a
function of time. After applying this s-correction (Stritzinger et al. 2002),
we merged the different datasets to build a photometric sequence of SN 2018ibb
from $t_{\rm max}$$-93$ to $t_{\rm max}$+706 days. We omitted these
corrections for the $BVJHK$ data because most observations in these filters
were done with the same instrument. We also skipped applying an s-correction
on the $HST$ observations in F336W at $t_{\rm max}$+1645 days as the SN was
not detected anymore. Table 12 in Appendix A summarises the homogenised SN
photometry. The measurements are not corrected for Galactic extinction along
the line of sight [$E(B-V)=0.03$ mag; Schlafly & Finkbeiner 2011], but this
correction is applied to all derived properties and photometric data presented
in this paper.

The photometry is available on WISeREP333https://www.wiserep.org (Yaron & Gal-
Yam 2012). It is also available as a machine-readable table in the electronic
version of this paper.

###  3.2 Host galaxy photometry

We obtained additional photometry with the ESO VLT, the 3.58 m New Technology
Telescope and the Hubble Space Telescope approximately 1000 days after maximum
(Appendix A). The brightness of the host galaxy was measured with elliptical
apertures encircling the entire host galaxy and calibrated in the same way as
the SN photometry. The HST photometry was done with a custom-made aperture
photometry tool, based on the python package photutils, using an aperture
comparable in area to the ground-based images and calibrated against tabulated
zeropoints in pysynphot. In the $R$-band, we measure a brightness of $24.39\pm
0.05$ mag. The brightness in the other filters is reported in Table 1.

###  3.3 Spectroscopy

Table 1: Photometry of the host galaxy

Telescope | Instrument | Filter | Brightness
---|---|---|---
|  |  | (mag)
HST | WFC3 | $F336W$ | $>26.04$
NTT | EFOSC2 | $B$ | $24.94\pm 0.22$
VLT | FORS2 | $g\\_{\rm HIGH}$ | $24.95\pm 0.05$
VLT | FORS2 | $R\\_{\rm SPECIAL}$ | $24.39\pm 0.05$
VLT | FORS2 | $I\\_{\rm BESSELL}$ | $24.32\pm 0.10$
VLT | FORS2 | $z\\_{\rm SPECIAL}$ | $23.78\pm 0.14$
444All measurements are reported in the AB system and not corrected for
reddening. Non-detections are reported at $3\sigma$ confidence.

Table 2: Log of spectroscopic observations

MJD | Phase | Telescope/Instrument | Disperser | Slit | Wavelength | Spectral | Exposure
---|---|---|---|---|---|---|---
| (day) |  |  |  width (′′) | range (Å) | resolution | time (s)
58453.349 | -1.4 | Keck-I/LRIS | 400/3400 + 400/8500 | 1.0 | 3076 – 9350 | 600/1200 | 300/300
58461.248 | 5.4 | P60/SEDm | … | IFU | 4650 – 9200 | 100 | 2250
58464.228 | 7.9 | P60/SEDm | … | IFU | 3950 – 9200 | 100 | 2250
58465.246 | 8.8 | P200/DBSP | 600/316 | 1.5 | 3500 – 10,000 | 1000/1000 | 1200
58467.254 | 10.5 | NTT/EFOSC2 | Gr#13 | 1.0 | 3650 – 9250 | 350 | 3600
58480.217 | 21.6 | P60/SEDm | … | IFU | 5500 – 8850 | 100 | 1200
58483.292 | 24.3 | Lick/Kast | 600/4310 + 300/7500 | 2.0 | 3500 – 10,500 | 800 | 1960
58487.225 | 27.6 | Lick/Kast | 600/4310 + 300/7500 | 2.0 | 3500 – 10,500 | 800 | 2400
58490.934 | 30.8 | NOT/ALFOSC | Gr#4 | 1.3 | 3600 – 9600 | 280 | 1800
58491.226 | 31.1 | NTT/EFOSC2 | Gr#11 + Gr#16/OG530 | 1.0 | 3345 – 9995 | 460/460 | 1800/1800
58493.096 | 32.7 | VLT/X-shooter | … | 1.0/0.9/0.9 | 3000 – 24,800 | 5400/8900/5600 | 1800
58509.904 | 47.1 | NOT/ALFOSC | Gr#4 + WG345 | 1.3 | 3800 – 9450 | 280 | 600
58510.179 | 47.3 | Lick/Kast | 600/4310 + 300/7500 | 2.0 | 3500 – 10,500 | 800 | 3600
58515.100 | 51.5 | NTT/EFOSC2 | Gr#11 + Gr#16/OG530 | 1.0/1.0 | 3345 – 9995 | 460/460 | 1800/1800
58525.048 | 60.1 | VLT/X-shooter | … | 1.0/0.9/0.9 | 3000 – 24,800 | 5400/8900/5600 | 2400
58536.929 | 70.3 | NOT/ALFOSC | Gr#4 | 1.3 | 3900 – 9600 | 280 | 2400
58541.083 | 73.8 | NTT/EFOSC2 | Gr#11 + Gr#16/OG530 | 1.0/1.0 | 3345 – 9995 | 460/460 | 2200/2200
58550.037 | 81.5 | VLT/X-shooter | … | 1.0/0.9/0.9 | 3000 – 24,800 | 5400/8900/5600 | 3600
58559.153 | 89.3 | Lick/Kast | 600/4310 + 300/7500 | 2.0 | 3500 – 10,500 | 800 | 2400
58565.000 | 94.3 | VLT/X-shooter | … | 1.0/0.9/0.9 | 3000 – 24,800 | 5400/8900/5600 | 3600
58718.304 | 225.8 | NTT/EFOSC2 | Gr#13 | 1.0 | 3650 – 9250 | 350 | 5400
58724.585 | 231.2 | Keck-I/LRIS | 400/3400 + 400/8500 | 1.0 | 3076 – 9350 | 600/1200 | 300/300
58776.290 | 275.6 | NTT/EFOSC2 | Gr#13 | 1.5 | 3650 – 9250 | 230 | 5400
58776.907 | 276.1 |  LBT/MODS+LUCI\scriptsize a\scriptsize afootnotemark: $\scriptsize a$ | G400L/G670L/G200 | 1.2/1.2/1.0 | 3200 – 12,000 | 925/1150/1100 | 3049/3083/3800
58789.287 | 286.7 |  VLT/X-shooter\scriptsize b\scriptsize bfootnotemark: $\scriptsize b$ | … | 1.0/1.0/0.9 | 3000 – 20,700 | 5400/8900/5600 | 3600
58866.134 | 352.6 |  VLT/X-shooter\scriptsize b\scriptsize bfootnotemark: $\scriptsize b$ | … | 1.0/1.0/0.9 | 3000 – 20,700 | 5400/8900/5600 | 3600
58876.648 | 361.6 |  LBT/MODS+LUCI\scriptsize a\scriptsize afootnotemark: $\scriptsize a$ | G400L/G670L/G200/G200 | 1.2/1.2/1.0/1.0 | 3200 – 23,500 | 925/1150/1100/1100 | 7200/7200/3400/2800
58895.120 | 377.5 |  VLT/X-shooter\scriptsize b\scriptsize bfootnotemark: $\scriptsize b$ | … | 1.0/1.0/0.9 | 3000 – 20,700 | 5400/8900/5600 | 3600
59110.570 | 562.3 | Keck-I/LRIS | 600/4000 + 400/8500 | 1.0 | 3400 – 10,275 | 1000/1200 | 4935/4935
59113.780 | 565.0 | VLT/FORS2 | 300V | 1.0 | 3800 – 9600 | 440 | 6600
59198.044 | 637.3 | VLT/FORS2 | 300V | 1.0 | 3800 – 9600 | 440 | 14400
59607.089 | 988.1 | Gemini-S/GMOS | R150/GG455 | 1.0 | 5000 – 10,000 | 310 | 4800
59608.380 | 989.2 | VLT/FORS2 | 300V | 1.0 | 3800 – 9600 | 440 | 14400

555The modified Julian dates quote the beginning of each spectroscopic
observation. The phase is reported for the mid-exposure time in the rest-frame
with respect to the $g|r$-band maximum at MJD=58455. For the multi-arm
instruments Kast, LBT, LRIS and X-shooter, we report the exposure time and
spectral resolution of each arm. The wavelength ranges and the values of the
spectral resolutions are taken from instrument manuals and are reported in the
observer frame. The spectral resolutions refer to the unbinned data.
$a$$a$footnotetext: The optical spectra were obtained with MODS and the grisms
G400L and G670L. The NIR spectra were obtained with the G200 grating and the
$zJ$ ($t_{\rm max}$+276.1 days) and $zJ$+$HK$ ($t_{\rm max}$+361.6 days)
spectroscopy filters. The start time is the average of the start times of the
visual and NIR spectra.
$b$$b$footnotetext: The observation was done with a $K$-band blocking filter.

We collected a series of spectra spanning from the time of maximum to $t_{\rm
max}$+989.2 days. Similarly to the imaging campaign, we utilised a large
number of 2–10 m class telescopes. A brief summary of the observations is
provided in Table 2. The details of the observations and data reduction are
presented in Appendix B. All spectra were absolute-flux-calibrated with multi-
band photometry. Since the photometry was not obtained contemporaneously with
the spectroscopic observation, we linearly interpolated between adjacent
observations.

The spectra obtained after August 2021 have an increasing contribution from
the host galaxy. The host contamination was removed with the FORS2 spectrum
from January 2022 ($t_{\rm max}$+989.2 days). The slit did not cover the
entire host galaxy. We scaled the spectrum to the flux encircled by the slit.
We note that to determine whether SN 2018ibb contributed to the observed
spectrum from January 2022, we compared the observed spectrum to the fit of
the spectral energy distribution (SED) of the entire host galaxy. The
continuum level of the January 2022 spectrum is fully consistent with the best
fit to the host galaxy SED (Figure 37). The only remaining SN feature is broad
[O iii] $\lambda\lambda$ 4959,5007 in emission, produced by the interaction of
the SN ejecta with circumstellar material (Section 5.1). Owing to that, we
masked the region and estimated the host galaxy flux with linear
interpolation. To recover the host-subtracted spectrum of SN 2018ibb from the
January 2022 epoch, we utilised the best fit to the galaxy SED.

All data were also corrected for Milky-Way (MW) extinction. We note that a few
spectra were affected by adverse weather conditions. The absolute-flux
calibrated spectra without MW extinction correction are available on WISeREP.

###  3.4 Imaging polarimetry

Table 3: Log of polarimetric observations

MJD | Phase | Exposure | Mean airmass | Filter | $q$ | $u$ | $p$ | $\theta$
---|---|---|---|---|---|---|---|---
| (day) | time (s) |  |  | (%) | (%) | (%) | (∘)
58492.241 | 31.9 | $4\times 100$ | 1.60 |  $v$_HIGH | $0.14\pm 0.08$ | $-0.24\pm 0.08$ | $0.28\pm 0.08$ | $150.2\pm 7.9$
58512.121 | 49.0 | $4\times 100$ | 1.15 |  $v$_HIGH | $0.10\pm 0.09$ | $-0.31\pm 0.09$ | $0.33\pm 0.09$ | $140.0\pm 8.3$
58524.125 | 59.3 | $4\times 100$ | 1.35 |  $v$_HIGH | $0.11\pm 0.10$ | $-0.25\pm 0.10$ | $0.28\pm 0.10$ | $146.8\pm 10.2$
58565.007 | 94.4 | $4\times 250$ | 1.33 |  $v$_HIGH | $0.21\pm 0.07$ | $-0.09\pm 0.07$ | $0.23\pm 0.07$ | $168.5\pm 5.2$
58565.020 | 94.4 | $4\times 250$ | 1.43 |  $R$_SPECIAL | $0.45\pm 0.07$ | $-0.16\pm 0.07$ | $0.48\pm 0.07$ | $170.4\pm 4.0$

666The first four columns summarise the observations. The last four columns
report the debiased polarisation properties: Stokes parameters $q$ and $u$,
the debiased polarisation level $p=\sqrt{q^{2}+u^{2}}$, and the position angle
$\theta=1/2\,\arctan\left(u/q\right)$.

To measure the ejecta geometry, we acquired four epochs of imaging polarimetry
in the $v$_HIGH filter with VLT/FORS2 between $t_{\rm max}$+31.9 and $t_{\rm
max}$+94.4 days (Table 3). In addition, we got one epoch with the $R$_SPECIAL
filter at $t_{\rm max}$+94.4 days. Each polarisation measurement required four
exposures at four different retarder-plate angles: $0^{\circ}$,
$22\aas@@fstack{\circ}5$, $45^{\circ}$, and $67\aas@@fstack{\circ}5$. The beam
was split with a Wollaston prism into the ordinary (o) and the extraordinary
(e) ray. The o-ray and the e-ray were placed at the 7th and the 8th multi-
object spectroscopy (MOS) stripes, respectively.

We reduced the data in a standard manner using IRAF (Tody 1993) tasks. The
flux of the SN in the o-ray and e-ray were measured through aperture
photometry at all four retarder-plate angles using the DAOPHOT.PHOT package
(Stetson 1987). Stokes parameters and polarisation of the target were derived
based on the FORS2 manual (Anderson 2018), and the polarisation degrees were
corrected for polarisation bias, caused by the non-negativity nature of the
polarisation degree, following Wang et al. (1997). The extracted, debiased
polarisation properties are summarised in Table 3.

These values need to be corrected for polarisation induced by dichroic
extinction from non-spherical dust grains that aligned with the magnetic field
of the interstellar medium of the Milky Way (MW) and the host galaxy.
Following Serkowski et al. (1975), the polarisation level from the Milky Way
can be as high as $\lesssim 9\%\times E(B-V)$. With a Galactic extinction of
$E(B-V)=0.03$ mag towards SN 2018ibb, the MW polarisation level could be up to
0.26%. The determination of the interstellar polarisation from SN 2018ibb’s
host galaxy is not feasible. We note that the polarisation degree is only
$p\lesssim$0.3% in $v$_HIGH between $t_{\rm max}$+32 and $t_{\rm max}$+94 days
(see Table 3). Such a low level of polarisation is very unlikely to be caused
by a high intrinsic polarisation aligned and cancelled to a comparable level
of significant interstellar polarisation. Therefore, without correcting for
the polarisation from the host galaxy, the observations point to a high degree
of spherical symmetry of SN 2018ibb during the phase of our polarisation
measurement.

Table 4: Log of X-ray observations

MJD | Phase | Count rate | $F_{\rm X}$ | $L_{\rm X}$
---|---|---|---|---
| (day) |  ($10^{-3}\rm ct\,s^{-1}$) |  ($10^{-15}\rm erg\,s^{-1}\,cm^{-2}$) |  ($10^{42}\rm erg\,s^{-1}$)
Swift
58469.56 | $12.5^{+65.5}_{-4.1}$ | $<0.7$ | $<27.0$ | $<2.2$
58567.93 | $96.9^{+5.2}_{-6.9}$ | $<1.9$ | $<68.6$ | $<5.6$
58592.17 | $117.7^{+5.9}_{-6.0}$ | $<2.3$ | $<82.7$ | $<6.7$
58688.35 | $200.1\pm 24.2$ | $<0.6$ | $<22.9$ | $<1.9$
58741.72 | $245.9\pm 12.8$ | $<1.8$ | $<64.9$ | $<5.3$
XMM-Newton
58511.22 | $48.5\pm 0.3$ | $<9.1$ | $<17.1$ | $<1.4$
58561.70 | $91.8\pm 0.3$ | $<10.2$ | $<19.2$ | $<1.6$
58694.68 | $205.8\pm 0.3$ | $<19.9$ | $<37.4$ | $<3.0$
58723.37 | $230.4\pm 0.3$ | $<8.7$ | $<16.4$ | $<1.3$
777The phases report the mid-exposure time. For Swift data, the modified
Julian date reports the mid-exposure time. The phase error quotes the bin size
after dynamically rebinning the data. For XMM-Newton data, the modified Julian
date reports the beginning of the observation. The total phase error
corresponds to the on-source integration time. All limits are reported at
$3\sigma$ confidence. The measurements are corrected for MW absorption and
reported for the bandpass from 0.3 to 10 keV.

###  3.5 X-ray Observations

####  3.5.1 Swift/XRT

While monitoring SN 2018ibb with UVOT between $t_{\rm max}$+8.4 and $t_{\rm
max}$+224 days, Swift also observed the field with the X-ray telescope XRT
between 0.3 and 10 keV in photon-counting mode (Burrows et al. 2005). We
analysed these data with the online-tools of the UK Swift
team888https://www.swift.ac.uk/user_objects that use the software package
HEASoft version 6.26.1 and methods described in Evans et al. (2007, 2009).

SN 2018ibb evaded detection in all epochs. The median $3\sigma$ count-rate
limit of all observing blocks is 0.005 count s-1 (0.3–10 keV). Using the
dynamic rebinning option in the Swift online tools pushes the $3\sigma$ count-
rate limits to 0.002 count s-1 (median value). A list of the limits from the
stacking analysis is shown in Table 4. To convert the count-rate limits into a
flux, we used WebPIMMS999https://heasarc.gsfc.nasa.gov/cgi-
bin/Tools/w3pimms/w3pimms.pl and assumed a power-law spectrum with a photon
index101010The photon index is defined as the power-law index of the photon
flux density ($N(E)\propto E^{-\Gamma}$). of $\Gamma=2$ and a Galactic neutral
hydrogen column density of $1.97\times 10^{20}$ cm-2 (HI4PI Collaboration et
al. 2016). The average energy conversion factor for the unabsorbed flux is
$3.66\times
10^{-11}\,\rm\left(erg\,s^{-1}\,cm^{-2}\right)/\left(ct\,s^{-1}\right)$. The
median count-rate limit corresponds to an unabsorbed flux of $<7.4\times
10^{-14}\leavevmode\nobreak\ {\rm erg\,cm}^{-2}\,{\rm s}^{-1}$ between 0.3–10
keV and a luminosity of $<4.9\times 10^{42}\leavevmode\nobreak\ {\rm
erg\,s}^{-1}$. The flux and luminosity limits of the individual bins are shown
in Table 4.

####  3.5.2 XMM-Newton

The field of SN 2018ibb was also observed by XMM-Newton (Jansen et al. 2001,
Principal Investigator: R. Margutti, University of California, Berkeley, USA).
Four epochs were taken with the European Photon Imaging Camera (EPIC) with the
pn (Strüder et al. 2001) and MOS1—2 cameras (Turner et al. 2001) between 28
January 2019 and 28 August 2019 ($t_{\rm max}$+48.5 – $t_{\rm max}$+230.4). We
reduced the XMM-Newton/EPIC pn data using the XMM-Newton Science Analysis
System111111https://www.cosmos.esa.int/web/xmm-newton/sas (SAS) following
standard procedures. We extracted the source using a circular region with a
radius of $32^{\prime\prime}$, and the background from a source-free region on
the same CCD. The MOS data are shallower than the pn data, so we omit
reporting them in this paper.

All XMM-Newton observations led to non-detections with count rate limits
between 0.009 and 0.020 ct s-1 between 0.3 and 10 keV. Using the same spectral
model as for XRT and an energy conversion factor of $1.88\times
10^{-12}\,\rm\left(erg\,s^{-1}\,cm^{-2}\right)/\left(ct\,s^{-1}\right)$, these
limits translate to unabsorbed flux limits between 1.6 and $3.7\times
10^{-14}\leavevmode\nobreak\ {\rm erg\,cm}^{-2}\,{\rm s}^{-1}$. Table 4
summarises the measurements.

###  3.6 Radio observations

The field was observed by the VLA Sky Survey (Lacy et al. 2020) between 2 and
4 GHz on 27 October 2020 ($t_{\rm max}$+595 days). No source was detected. The
flux at the SN position is $-47\pm 223$ $\mu$Jy, translating to a $3\sigma$
flux limit of 622 $\mu$Jy and a luminosity of $2\times 10^{39}\,\rm
erg\,s^{-1}$. Eftekhari et al. (2021) presented sub-mm observations at 100 GHz
obtained with the Atacama Large Millimeter Array on 24 December 2019 ($t_{\rm
max}$+331 days). These authors also reported a non-detection with an r.m.s. of
19 $\mu$Jy, translating to a $3\sigma$ flux limit of 58 $\mu$Jy and a
luminosity of $5\times 10^{39}\,\rm erg\,s^{-1}$.

##  4 Results

###  4.1 Redshift

The X-shooter spectra between $t_{\rm max}$+32.7 and $t_{\rm max}$+94.3 days
show narrow absorption lines of Mg i $\lambda$ 2852 and Mg ii $\lambda\lambda$
2796,2803 from the host galaxy at a common redshift of $z=0.1660$ (Figure 2,
top panel). The low-resolution FORS2 spectrum obtained at $t_{\rm max}$+565.3
days, shown in the bottom panel of Figure 2, reveals narrow emission lines
from hydrogen and oxygen from the H ii regions in the host galaxy at the same
redshift as the absorption-line redshift. This redshift translates to a
luminosity distance of 822.6 Mpc using the cosmological parameters from Planck
Collaboration et al. (2016).

Figure 2: Galaxy absorption and emission lines at a common redshift of
$z=0.166$ in the supernova spectra at $t_{\rm max}$+32.7 days (top) and at
$t_{\rm max}$+565.3 days (bottom). The error spectrum of each epoch is shown
in grey. Figure 3: Multi-band light curve of SN 2018ibb from 1800 to 18,500 Å
(rest-frame) after correcting for the Galactic extinction. SN 2018ibb was
first detected by Gaia. The last non-detections before the first detection by
Gaia and ATLAS are shown by the downward-pointing triangles. With a rise time
of $>93$ rest-frame days, SN 2018ibb is one of the slowest evolving SLSN
known. The decline of 1.1 mag (100 days)-1 is similar to the decay time of
radioactive 56Co. After $t_{\rm max}$+575 days, the decline steepened to 1.5
mag (100 days)-1. The light curve shows undulations up to $t_{\rm max}$+100
days and a longer-lasting bump at $\sim 300$ rest-frame days. Vertical bars
represent the epochs of spectroscopy and imaging polarimetry. The absolute
magnitude is computed with $M=m-{\rm DM}(z)+2.5\,\log\left(1+z\right)$, where
DM is the distance modulus and $z$ the redshift.

###  4.2 Light curve

####  4.2.1 General properties

Figure 3 shows the evolution of SN 2018ibb from $t_{\rm max}$$-93$ to $t_{\rm
max}$+706 rest-frame days. The early-time evolution was captured by the ATLAS,
Gaia and ZTF surveys. Human scanners discovered SN 2018ibb shortly before
maximum light, which triggered our large monitoring campaign from UV to NIR
wavelengths. The $g$, $r$ and $o$ band light curves cover the evolution from
early to late times. We use these datasets to infer the time of maximum light
and the rise and decline time scales. Fitting the light curves with 3rd order
polynomials between ${\rm MJD}=58425$ and ${\rm MJD}=58485$ returns the time
of maximum light at ${\rm MJD}=58458\pm 2$, $58454\pm 2$ and $58452\pm 4$ in
$g$, $r$ and $o$, respectively. Throughout the paper, we adopt the weighted
mean MJD $58455\pm 2$ as the time of maximum light. At the time of peak, SN
2018ibb reached a brightness of $17.54\pm 0.02$, $17.65\pm 0.01$ and $17.92\pm
0.04$ in $g$, $r$ and $o$ band, respectively (all corrected for MW extinction;
Table 5).121212The host extinction is negligible (Section 4.6). Using the Keck
spectrum at $t_{\rm max}$$-$1.4 days, we infer k-corrected absolute magnitudes
of $-21.79\pm 0.02$, $-21.66\pm 0.01$ and $-21.43\pm 0.04$ mag in the
aforementioned bands (Table 5) and a k-corrected $g-r$ colour of $-0.12\pm
0.02$ mag at peak (corrected for MW extinction), a typical luminosity and
colour for a H-poor SLSN (Nicholl et al. 2015a; De Cia et al. 2018; Lunnan et
al. 2018a; Angus et al. 2019; Chen et al. 2023a).

Similar to Chen et al. (2023a), we measure the rise and decline time scales
from 10% and 50% peak flux to peak in all three bands. In the $g$ band, we
obtain $t_{1/2,\,\rm rise}=52\pm 1$ days, $t_{1/2,\,\rm decline}=88^{+1}_{-2}$
days, $t_{1/10,\,\rm rise}>79.3$ days, and $t_{1/10,\,\rm decline}=242\pm 1$
days, that is, 1 mag (100 days)-1 (all measured in the rest-frame). The light
curve parameters in the other bands are summarised in Table 5. Although the
Gaia light curve is poorly sampled, the data are of sufficient quality to
improve the lower limit on the rise timescale $t_{1/10,\,\rm rise}$. The Gaia
alert database reports the first detection on ${\rm MJD}=58346.11$ (16 August
2018), 11.5 and 13.3 rest-frame days before the first ZTF and ATLAS131313The
last ATLAS non-detection is from 17 August 2018, that is, 1.3 rest-frame days
after the first Gaia detection, reaching a limiting magnitude of $o\approx
19.9$ mag at 3 sigma confidence. detection, respectively. At the time of
discovery, SN 2018ibb had a brightness of 19.8 mag; around the time of maximum
light, the brightness reached 17.7 mag. This sets a lower limit of $>93.4$
days on $t_{1/10,\,\rm rise}$.

Table 5: Light curve properties

Property | $g$ | $r$ | $o$
---|---|---|---
Peak time (MJD) | $58458\pm 2$ | $58454\pm 2$ | $58452\pm 4$
$m_{\rm peak}$ (mag) | $17.54\pm 0.02$ | $17.65\pm 0.01$ | $17.72\pm 0.02$
$M_{\rm peak}$ (mag) | $-21.80\pm 0.02$ | $-21.66\pm 0.01$ | $-21.62\pm 0.02$
$t_{1/2,\rm rise}$ (day) | $52\pm 1$ | $60\pm 1$ | $64\pm 1$
$t_{1/2,\rm decline}$ (day) | $88^{+1}_{-2}$ | $93\pm 1$ | $95\pm 3$
$t_{1/e,\rm rise}$ (day) | $68.3\pm 0.4$ | $72.5\pm 0.5$ | $73.8\pm 0.5$
$t_{1/e,\rm decline}$ (day) | $102\pm 2$ | $117\pm 1$ | $>107$
$t_{1/10,\rm rise}$ (day) | $>79.3$ | $>82$ | $>80$
$t_{1/10,\rm decline}$ (day) | $242\pm 1$ | $248\pm 1$ | $235^{+1}_{-4}$
141414All magnitudes are corrected for Galactic extinction. The absolute
magnitudes include a k-correction inferred from the Keck spectrum at $t_{\rm
max}$$-1.4$ days. All time scales are reported in the rest-frame. The
uncertainties reflect the $1\sigma$ statistical errors.

Between July 2018 and the date of the first Gaia detection (16 August 2018),
the field was visible to observing facilities in the southern hemisphere. We
searched the data archives of the Australian Astronomical Observatory, the
European Southern Observatory, the Gemini Observatory, and the Las Cumbres
Observatory for serendipitous observations of this field but found no relevant
data. We conclude that SN 2018ibb’s progenitor exploded $>93$ rest-frame days
before the maximum light, but we have no firm constraint on the explosion
date.151515The Gaia alert database reports an observation on 5 July 2018 but
no measurement. This could either mean i) a non-detection (limiting magnitude
$G=20.7$ mag) and hence imposing an upper limit of $<129$ days on
$t_{1/10,\,\rm rise}$, ii) the observation was not performed, or iii) a
problem in the data processing (Hodgkin et al. 2021).

We obtained a final epoch of photometry with HST/WFC3 in $u$ band at $t_{\rm
max}$+1165 days as a part of an HST Snapshot programme to search for signs of
late-time CSM interaction in SNe (Fremling2021a). SN 2018ibb evaded detection.
We place an upper limit of $u=26.2$ ($3\sigma$ confidence), shown as a
downward pointing triangle in Figure 3.

SN 2018ibb’s light curve exhibits several peculiar properties. Figure 4
compares the absolute magnitude versus the rest-frame phase of SN 2018ibb to
the light curves of the 78 H-poor SLSNe from ZTF-I presented in Chen et al.
(2023a). The absolute magnitude of all objects is computed with $M=m-{\rm
DM}(z)+2.5\,\log\left(1+z\right)$, where DM is the distance modulus and $z$
the redshift. SN 2018ibb has the longest rise in the ZTF sample. The $g$-band
rise time $t_{1/10,\,\rm rise}$ exceeds the sample mean value (41.9 days) by a
factor of 2.1 times the sample standard deviation [$\sigma(t_{1/10,\,\rm
rise})=17.8$ days; Chen et al. 2023a]. This factor could increase to even
$4.9\sigma$ if the Gaia data are a good proxy of the rise time in ZTF/$g$. The
light curve fades by 1.1 mag (100 days)-1 for 500–600 days before the decline
steepens to 1.5 mag (100 days)-1. The decline time scale is slower than for
any of the other H-poor SLSNe from the ZTF-I sample. The rise is even slower
than any of the $>100$ H-poor SLSNe found by other surveys (as queried from
the Transient Name Server161616https://www.wis-tns.org and ADS Abstract
Service171717https://ui.adsabs.harvard.edu). Only the H-poor SLSN PS1-14bj
(Lunnan et al. 2016) had a longer rise. We discuss this in more detail in
Section 5.4.

A number of SNe have shown a pre-bump with observed peak luminosities between
$M_{g}\sim-18$ and $-23$ mag (Leloudas et al. 2012; Nicholl et al. 2015b;
Smith et al. 2016; Angus et al. 2019). Such a bump is not visible in the light
curve of SN 2018ibb (Figures 3 and 4). However, the progenitor of SN 2018ibb
exploded before the field came from behind the sun, precluding drawing a firm
conclusion on the absence or presence of a pre-bump.

Figure 4: $g$-band light curve of SN 2018ibb in the context of the homogeneous
ZTF-I SLSN sample. SN 2018ibb has a typical peak absolute magnitude. The rise
of $>93$ rest-frame days is significantly longer than of the average ZTF SLSN.
The long-lasting rise implies a long diffusion time, which requires a very
high total ejected mass. The high peak luminosity requires a very energetic
explosion. Both properties together hint to an explosion mechanism that might
be different from that of regular SLSNe.

####  4.2.2 Bolometric light curve

We compute the bolometric luminosity of SN 2018ibb over a wavelength range
from $\sim 1800$ to $\sim 14,300$ Å (rest-frame), which is defined by the
wavelength coverage of our photometric dataset. However, our dataset does not
have the same wavelength coverage throughout the entire duration of the
observations. In the following, we describe how the bolometric light curve is
constructed and discuss the bolometric corrections that we derived for time
intervals with incomplete spectral information.

Figure 5:  Bolometric light curve of SN 2018ibb from 1800 to 14,300 Å (rest-
frame). The dotted lines indicate time segments with partial wavelength
coverage. At peak SN 2018ibb reached a luminosity of $>2\times
10^{44}\leavevmode\nobreak\ \rm erg\,s^{-1}$. Integrating over the light curve
from $t_{\rm max}$$-93$ to $t_{\rm max}$+706 days yields a radiated energy of
$3>\times 10^{51}\leavevmode\nobreak\ \rm erg$. Both values are conservative
lower limits. The inset shows the evolution of the blackbody temperature and
radius of the photospheric phase where photometry has been carried out from
the $u$ to $H$ bands. The shaded regions indicate the $1\sigma$ statistical
uncertainties.

The bolometric light curve is built as follows: i) correcting all photometric
data for the MW extinction, ii) dividing the entire dataset into segments
defined by the observing seasons, iii) interpolating the light curve in each
band of each observing season with a Gaussian process with the python package
George version 0.4.0 (Ambikasaran et al. 2015)181818 We added a systematic
error of 5% to all optical and NIR filters and 10% to all UV filters in
quadrature to account for uncertainties in the flux calibration., iv)
constructing the spectral energy distributions for every time step, v)
calculating the bolometric flux by numerical integration of each SED, and vi)
multiplying the bolometric flux by $4\pi\,d^{2}_{L}$, where $d_{L}$ is the
luminosity distance, to obtain the bolometric luminosity.

Our dataset has the best spectral coverage between $t_{\rm max}$ and $t_{\rm
max}$+375 days: 1800 to 14,300 Å between $t_{\rm max}$ to $t_{\rm max}$+100
days and 3000 to 14,300 Å between $t_{\rm max}$+200 to $t_{\rm max}$+375 days.
The bolometric light curve for these time intervals are shown as solid lines
in Figure 5 and their $1\sigma$ confidence intervals as a shaded region. A
tabulated version can be found in Table 13. Based on the blackbody fits to the
data from $u$ to $H$ band, we estimate that $\lesssim 3\%$ of the observed
bolometric flux is emitted at longer wavelengths between $t_{\rm max}$ and
$t_{\rm max}$+100 days. Linearly extrapolating the observed SED from 1800 to
14,300 Å towards shorter wavelengths yields a missing UV contribution of $\ll
1\%$. Hence, we omit to correct the observed bolometric flux. At later phases,
the spectrum does not resemble a blackbody anymore (Figure 6), and we cannot
quantify the missing flux at longer and shorter wavelengths.

Figure 6: Spectroscopic sequence from 2500 Å to 10,000 Å and from the time of
maximum to $t_{\rm max}$+1000 days (rebinned to 5 Å and smoothed with a
Savitzky-Golay filter). Spectra up to $t_{\rm max}$+100 days (left panel) are
characterised by a blackbody continuum with superimposed absorption lines from
the SN ejecta, expanding with a velocity of $\sim 8500\leavevmode\nobreak\ \rm
km\,s^{-1}$. Between $t_{\rm max}$+100 and $t_{\rm max}$+225 days (while SN
2018ibb was behind the sun), the spectroscopic behaviour of SN 2018ibb evolved
drastically. The late-time spectra (right panel) are characterised in the blue
($<5000$ Å) by a pseudo-continuum and emission lines produced by the
interaction of the SN ejecta with circumstellar material and in the red
($>5000$ Å) by nebular emission lines from the 56Ni-heated SN core. The
regions with the fastest evolution are highlighted by the grey-shaded regions.
Figure 7 shows the identification of the most prominent features of the
photospheric and nebular phases. Regions affected by strong telluric
absorption were clipped. Their locations are indicated in Figure 7.

For the other epochs, we used these time intervals ($t_{\rm max}$ to $t_{\rm
max}$+100 days and $t_{\rm max}$+200 to $t_{\rm max}$+375 days) to estimate
bolometric corrections. The pre-max dataset consists of photometry in ZTF $g$
and $r$, and ATLAS $c$ and $o$ filters, and the Gaia white band. We only use
the ZTF data when computing the bolometric luminosity because the ATLAS and
Gaia filters are too broad for building SEDs. At the time of the first epoch
with coverage from $w2$ to $H$, $\sim 26\%$ of the bolometric flux was emitted
in $g+r$ band. We use this flux ratio as an estimate of the missing flux.
Since SN ejecta cool with time, such a universal correction will progressively
underestimate the bolometric flux towards earlier epochs. Between the first
and second observing seasons, we continued the follow-up with Swift/UVOT in
$ubv$ when SN 2018ibb was no longer visible from the ground. Similar to the
pre-max data, we chose time intervals with data from $w2$ to $H$ or $u$ to $H$
band to correct for the missing flux. At phases later than $t_{\rm max}$+500
days, photometric data are only available from $g$ to $z$ band. We omit to
apply any bolometric correction for this time interval because we have no good
estimate of the missing bolometric flux.

SN 2018ibb reached a peak luminosity of $L_{\rm bol,\,peak}\geq 2\times
10^{44}\leavevmode\nobreak\ {\rm erg\,s}^{-1}$. Integrating the light curve
from $t_{\rm max}$$-93$ to $t_{\rm max}$+706 days yields $\geq 3\times
10^{51}\leavevmode\nobreak\ {\rm erg}$ for the total radiated energy $E_{\rm
rad}$. We emphasise that both values are strict lower limits. Our multi-band
campaign only started when SN 2018ibb peaked in the $g$ and $r$ bands, which
was likely after the bolometric peak.

Between $t_{\rm max}$ and $t_{\rm max}$+100 days, the spectra of SN 2018ibb
are characterised by a cooling photosphere (Figure 6), and the spectral energy
distributions from the $u$ to $H$ band are adequately fitted with a Planck
function. The red and blue curves in the inset of Figure 5 show the evolution
of the blackbody temperature and radius (see also Table 13), respectively. The
photosphere has a temperature of 12,000 K at the time of maximum light and
cools by 3000 K in 100 rest-frame days. During the same time interval, the
location of the photosphere hardly changes from its mean value of $5\times
10^{15}$ cm. The values of the blackbody radius and temperature are comparable
to regular SLSNe (Chen et al. 2023a) and the slow-evolving SLSN 2015bn
(Nicholl et al. 2016b), which have observations in the UV. The blackbody
temperature of SN 2018ibb evolves slower than for regular SLSNe (Chen et al.
2023a), mirroring its slowly evolving light curve. We remark that including
data at shorter wavelengths would have led to lower temperatures ($\approx
0.1$ dex at $t_{\rm max}$) and larger radii ($\approx 0.12$ dex at $t_{\rm
max}$) due to absorption lines in the UV (Yan et al. 2017; Lunnan et al.
2018a; Angus et al. 2019). Owing to this, we omit these data to infer the
blackbody radius and temperature.

Figure 7: Line identification of the photospheric-phase spectrum (top) and
nebular-phase spectrum (bottom). Top: The photospheric phase spectrum was
fitted with the parameterised spectral synthesis code SYNOW (red curve). Most
of the spectral features can be attributed to O i, Mg ii, Si ii, Ca ii, and Fe
ii as seen in other SLSNe during their cool photospheric phase (Gal-Yam
2019b). In addition to the absorption lines in the SN ejecta, the photospheric
phase spectrum shows conspicuous [Ca ii] $\lambda\lambda$ 7291, 7324, a
feature that gets dominated by [O ii] $\lambda\lambda$7320,7330 at about
$t_{\rm max}$+30 days. Bottom: The spectrum of the nebular phase consists of a
blue pseudo-continuum and a series of allowed and forbidden emission lines
from singly and doubly ionised oxygen, calcium, magnesium and iron. Remarkable
is the presence of [O ii] and [O iii] in emission (as early as $t_{\rm
max}$+30 days), indicating ionising radiation from shock interactions (Section
5.1). SN absorption lines are indicated by dashed lines, and the locations
mark the absorption trough minima (blueshifted by $8500\leavevmode\nobreak\
\leavevmode\nobreak\ \rm km\,s^{-1}$ from their rest wavelengths). SN emission
lines are indicated by solid lines; their line centres are at the velocity
coordinate $v=0$. Regions of strong atmospheric absorption are grey-shaded.

###  4.3 Spectroscopy

####  4.3.1 Spectroscopic sequence

Figure 6 shows the spectral evolution between $\sim 2800$ Å to $\sim 10,000$ Å
from the time of maximum to $t_{\rm max}$+990 days (all rest-frame). The
spectra up to $t_{\rm max}$+100 days capture the photospheric phase. To
identify the elements and ions responsible for the most prominent features, we
model the spectrum at $t_{\rm max}$+32.7 days with the spectrum synthesis code
SYNOW (Branch et al. 2005). The SYNOW fit, shown in the top panel of Figure 7,
was obtained for a photospheric expansion velocity of 8000 km s-1 (Section
4.3.2) and for a blackbody temperature of 12,000 K (Section 4.2.2; a range in
the order of $\pm 500$ is applicable for both properties). The major ions that
are securely identified and match the spectrum well are those of: O i, Mg ii,
Si ii, Ca ii, and Fe ii (the Mg ii mainly improves the match of the feature
around 4400 Å, together with the Fe ii line), in agreement with Könyves-Tóth &
Vinkó (2021). Various additional iron group elements, such as Ti ii, clearly
help to lower the model flux on the blue side (3000–4000 Å). However, we do
not include those in the final SYNOW fit because the overall fit was not
convincingly improved. Intriguingly, the spectrum shows narrow absorption
lines from Fe II $\lambda\lambda$ 4924,5018,5169 (half-width at zero intensity
$\approx 1500\leavevmode\nobreak\ \rm km\,s^{-1}$), reminiscent of the H-poor
SLSNe 1999as and 2007bi. Modelling the photospheric-phase spectra of these two
SLSNe revealed that such narrow features are challenging for existing SLSN
models (Moriya et al. 2019). These narrow features could point to a velocity
cut in the density structure of the SN ejecta. This could be related to the
density structure of the progenitor or possibly point to the deceleration of
the outermost layer of the ejecta by the collision with dense circumstellar
material (Kasen 2004; Moriya et al. 2019). To differentiate between these
scenarios, spectral modelling is required. This is beyond the scope of this
paper. Absorption from O ii between 3700 Å and 4700 Å, as seen in many SLSN
spectra around peak (Quimby et al. 2018), is not present.

Owing to the limitations of the SYNOW approach, for instance, the simplifying
underlying assumptions such as spherical, homologous expansion, resonant
scattering line formation above a sharp blackbody spectrum-emitting
photosphere, we perform this modelling only for the identification and
verification of the major features. We avoid any fine-tuning of the different
ion parameters and assessing the elemental abundances or relative mass
fractions.

A complementary analysis with the National Institute of Standards and
Technology (NIST) Atomic Spectra Database (Kramida et al. 2018), following the
methodology described in Gal-Yam (2019a), which includes the same elements as
above for relative intensities $\geq 0.5$ in the range 2000–10,000 Å (and
$\geq 0.2$ in the range 3000–6000 Å for the Fe ii lines), reveals additional
possible identification of features that are not accounted for by the SYNOW
fit. For instance, lines of Mg ii and/or Si ii may contribute to the small dip
redwards of the $\sim 7773$ Å (rest-frame) O i triplet. Also, numerous Fe ii
lines may contribute to the valley around 3000–3200 Å as well as additional Mg
ii lines accounting for the dips around 4300 Å. Remarkably, in addition to
absorption lines from the SN ejecta, the first spectrum at $t_{\rm max}$$-1.4$
days shows conspicuous [Ca ii] $\lambda\lambda$ 7291,7323 in emission. This is
one of the strongest forbidden emission lines seen in nebular SN spectra
(Filippenko 1997; Gal-Yam 2017). The only SLSNe that show [Ca ii] during the
photospheric phase are slow-evolving SLSNe (e.g. SN 2007bi, LSQ14an and SN
2015bn; Gal-Yam et al. 2009; Nicholl et al. 2019; Inserra et al. 2017).

Figure 8: Spectroscopic sequence from $9500$ to $21,500$ Å. The spectral
sequence covers the evolution of the photospheric (top) and nebular (bottom)
phases. The NIR spectra at $>1$ $\mu$m show only a few features in contrast to
the optical spectra (Figure 6). The most prominent features are labelled. All
spectra were rebinned to 5 Å and smoothed with a Savitzky-Golay filter, except
the spectrum at $t_{\rm max}$+361.6 days that was rebinned to 10 Å. The grey
scale at the bottom of each panel displays the strength of telluric features
(white = transparent, black = opaque). In addition, regions of strong
atmospheric absorption are grey-shaded.

During the first seasonal observing gap, the photosphere recedes and we start
to see the core of the explosion. The nebular spectra (right panel in Figure
6) are dominated by emission lines with widths up to 10,000 km s-1 and a blue
pseudo-continuum, similar to that seen in SNe Ia-CSM, Ibn, and Icn and some
SNe IIn (e.g. Silverman et al. 2013; Hosseinzadeh et al. 2017; Gal-Yam et al.
2022; Perley et al. 2022). Following previous observations of slow-evolving
SLSNe (Nicholl et al. 2016a; Lunnan et al. 2016) and theoretical models by
Jerkstrand et al. (2016), we identify the most conspicuous emission lines as
allowed and forbidden transitions from neutral and ionised calcium, iron,
magnesium, and oxygen (Figure 7).

Common to both the photospheric and the nebular phase is that the evolution is
very slow with the exceptions of the regions at $\sim 4360$, $\sim 5000$ and
7300 Å (highlighted in grey in Figure 6). At about 30 days after maximum, the
region at $\sim 5000$ Å shows a rapidly growing emission feature. A weak
emission line at $\sim 4360$ Å also emerges and reveals a similar trend to the
$\sim 5000$ Å feature. Owing to this, we identify the two features as [O iii]
$\lambda$ 4363 and [O iii] $\lambda\lambda$ 4959,5007, respectively. Most
remarkably, the [O iii] $\lambda\lambda$ 4959,5007 emission lines are present
throughout the entire post-max evolution, even in the spectrum at $t_{\rm
max}$+989.2 days after all other SN features faded below the detection
threshold of the 4-hour VLT spectrum. This has never been observed in any SLSN
before. Simultaneous with the rise of [O iii], the centre of the 7300 Å
feature moves a few Å to longer wavelengths, the line profile changes from
roughly top-hat to bell-shaped, the width decreases and the peak flux
increases by a factor $\sim 2$ within $<60$ days (Figure 6, left panel). This
suggests that this line complex, commonly identified as [Ca ii]
$\lambda\lambda$ 7291,7324, gets dominated by [O ii] $\lambda\lambda$
7320,7330.

[O ii] and even more so [O iii] are not common features for SNe. [O iii] was
only observed in the slow-evolving H-poor SLSNe LSQ14an (Inserra et al. 2017)
and PS1-14bj (Lunnan et al. 2016) during the photospheric phase and in SN
2015bn in the nebular phase (Nicholl et al. 2016a; Jerkstrand et al. 2017).
Occasionally, it is also seen in regular H-poor and H-rich SNe predominantly
years after the explosion (e.g. SNe II 1979C and 1980K, Milisavljevic et al.
2009 and Fesen et al. 1999; SN IIb 1993J, Milisavljevic et al. 2012; SNe IIn
1995N, 1996cr, 2010jl, Fransson et al. 2002; Bauer et al. 2008; Milisavljevic
et al. 2012; Fransson et al. 2014; SN Ib 2012au, Milisavljevic et al. 2018),
and even more rarely during the photospheric phase of regular SNe (e.g. Type
Ic SN 2021ocs; Kuncarayakti et al. 2022). Possible mechanisms to produce [O
ii] and [O iii] are i) excitation by CSM interaction (Chevalier & Fransson
1994), ii) photoionisation by the interaction of the pulsar wind nebula with
the SN ejecta (Chevalier & Fransson 1992; Omand & Jerkstrand 2023), and iii)
radioactivity (for high ratios of deposited energy to O-density; Jerkstrand et
al. 2017). In Section 5.1, we show that [O ii] and [O iii] are produced by the
interaction of the SN ejecta with circumstellar material.

Our series of NIR spectra (shown in Figure 8) covers the photospheric phase
from $t_{\rm max}$+33 to $t_{\rm max}$+94 days, and the nebular phase from
$t_{\rm max}$+276 to $t_{\rm max}$+378 days. The NIR spectra show a limited
number of absorption and emission lines. The photospheric-phase spectra show
two features at 1.093 and 1.13 $\mu$m. Following Jerkstrand et al. (2015),
Hsiao et al. (2019) and Shahbandeh et al. (2022), we tentatively identify the
former as an absorption line of Mg ii $\lambda$ 1.092 $\rm\mu m$ blueshifted
by $\sim 8500\leavevmode\nobreak\ \rm km\,s^{-1}$, and the latter as the
recombination line O i $\lambda$ 1.13 $\mu$m. These features can also be
blended with emission lines from sulphur. The emission lines clearly stand out
in the nebular-phase spectra. Our NIR spectra at $t_{\rm max}$+378 days show a
prominent emission line at 1.025 $\rm\mu m$ that we tentatively identify as
[Co ii] $\lambda$ 1.025. This is the first time that a cobalt line has been
detected in a SLSN spectrum. In Section 5.2.5, we examine this detection in
more detail.

####  4.3.2 Ejecta velocity

The photospheric-phase spectra of SN 2018ibb show a large number of narrow
absorption lines, mirroring a low ejecta velocity and the slow light curve
evolution. The ejecta velocities are commonly measured from Fe ii $\lambda$
5169\. Owing to the high velocities of SLSNe (e.g. Liu et al. 2017; Chen et
al. 2023b), this line is usually blended with Fe ii $\lambda$ 4924 and Fe ii
$\lambda$ 5018, necessitating template matching techniques to extract the
velocities (Modjaz et al. 2016; Liu et al. 2017). However, the ejecta velocity
of SN 2018ibb is slow, and the Fe ii $\lambda$ 5169 region is not blended and
resolves into three absorption lines that we identify as Fe ii
$\lambda\lambda$ 4924, 5018 and 5169 (Figure 9). By measuring the minima of
the three absorption lines, we extract a photospheric velocity of $\approx
8500\leavevmode\nobreak\ \rm km\,s^{-1}$ that remains constant between $t_{\rm
max}$ and $t_{\rm max}$+100 days as demonstrated in Figure 9 (all measurements
are summarised in Table 6).

The maximum ejecta velocity is best determined from the blue edge of the
strong Mg ii $\lambda\lambda$ 2796,2803 and Ca ii $\lambda\lambda$ 3934, 3968
resonance lines. In Figure 10, we show the regions around the two features at
$t_{\rm max}$+32.7 days, centred on the blue doublet components. Because of
the complexity of line features, we omit to subtract any continuum. For
illustration purposes, we normalise the spectral regions so that the peak
intensity and maximum absorption approximately match both lines. The blue
components of the doublets exhibit complex profiles at low velocities because
of the superposition with the wings of the red doublet components. The highest
velocities are less affected by this. The Ca ii $\lambda$ 3934 line gives the
best estimate for the maximum ejecta velocity, $\sim
12,500\leavevmode\nobreak\ \rm km\,s^{-1}$. This is consistent with the extent
of the absorption component of Mg ii $\lambda$ 2796, which, however, is more
affected by other SN lines. The absorption minima of Mg ii $\lambda$ 2796 and
Ca ii $\lambda$ 3934 are at $\sim 8000\leavevmode\nobreak\ \rm km\,s^{-1}$,
but are affected by the doublet nature of the lines. Nonetheless, the
locations of the absorption minima are consistent with the photospheric
velocity determined from the absorption minima of Fe ii $\lambda\lambda$
4924,5018,5169.

Figure 9:  Zoom-in onto the Fe ii absorption lines from the SN ejecta at
selected epochs of the photospheric phase. The SN photosphere expands with a
velocity of merely $\approx 8500\leavevmode\nobreak\ \rm km\,s^{-1}$. There
are no signs of deceleration between $t_{\rm max}$ and $t_{\rm max}$+100 days.
Starting at about $t_{\rm max}$+30 days, emission from [O iii]
$\lambda\lambda$ 4959, 5007 (grey shaded region), produced by the interaction
of the SN ejecta with circumstellar material, contaminates the blue wing of Fe
ii $\lambda$ 5169\.

To put the measurements in the context of other SLSNe, we first compare the
velocity of SN 2018ibb at maximum light to those of SLSNe in the ZTF-I sample
(Chen et al. 2023b). The histogram in the top panel of Figure 11 shows a
kernel density estimate of the velocity distribution of the 27 SLSNe from the
ZTF-I sample with Fe ii velocities measured within $\pm 20$ rest-frame days
from maximum light. After bootstrapping the sample and propagating the
measurement uncertainties with a Monte Carlo simulation, the median velocity
of the ZTF-I sample is $14,800\leavevmode\nobreak\ \rm km\,s^{-1}$ and its
$1\sigma$ confidence region extends from $10,500$ to
$19,000\leavevmode\nobreak\ \rm km\,s^{-1}$. SN 2018ibb lies in the bottom 8%
of this sample, but its velocity is not unparalleled. SN 2019aamu had a lower
photospheric velocity at peak, but the measurement is poorly constrained
($5876^{+8110}_{-349}\leavevmode\nobreak\ \rm km\,s^{-1}$; Chen et al. 2023b).
In the bottom panel of Figure 11, we show the evolution of the Fe ii
velocities of SN 2018ibb together with those of H-poor SLSNe from Liu et al.
(2017) (in grey). Within 50 days after maximum, the ejecta usually decelerate
from $\sim 15,000\leavevmode\nobreak\ \rm km\,s^{-1}$ to $\lesssim
10,000\leavevmode\nobreak\ \rm km\,s^{-1}$, whereas SN 2018ibb shows no
evolution.

####  4.3.3 A CSM shell around the progenitor of SN 2018ibb

The X-shooter spectra between $t_{\rm max}$+32.7 days and $t_{\rm max}$+94.3
days show two Mg ii absorption line systems (Figure 12). The narrow component
is associated with the gas in the SLSN host galaxy (Section 4.6). The lines of
the broader component have a full width at half maximum (FWHM) of
$406\leavevmode\nobreak\ {\rm km\,s}^{-1}$ and are blueshifted by
$2918\leavevmode\nobreak\ {\rm km\,s}^{-1}$ (not varying between $t_{\rm max}$
and $t_{\rm max}$+90 days; upper panels in Figure 12). They are significantly
broader than expected for the interstellar medium in the dwarf host galaxy or
any intervening dwarf galaxy191919Based on the correlations between the
stellar mass of galaxies and the width of galaxy absorption and emission lines
(Krühler et al. 2015; Arabsalmani et al. 2018). but also significantly
narrower than the narrowest SN features ($\sim 1900\leavevmode\nobreak\ \rm
km\,s^{-1}$; measured from Fe ii). The equivalent widths are $2.00\pm 0.09$
and $1.27\pm 0.08$ Å for Mg ii $\lambda$ 2796 and Mg ii $\lambda$ 2803,
respectively. The observed line ratio is $1.57\pm 0.12$ in tension with the
predicted value of 2 for unsaturated lines. Assuming that the Mg ii lines are
unsaturated, we can convert their equivalent widths to a lower limit on the
column density of singly ionised magnesium in the CSM shell. The rest-frame
equivalent width ${\rm EW}_{\rm r}$ is related to the column density $N$, in
units of atoms per $\rm cm^{2}$, via $N=1.13\times 10^{20}\leavevmode\nobreak\
{\rm EW}_{\rm r}\,/\,\left(\lambda^{2}_{\rm r}\,f\right)$ where $\lambda_{\rm
r}$ is the rest-frame wavelength, in units of Å, and $f$ the oscillator
strength (Ellison et al. 2004). Using the oscillator strengths from Theodosiou
& Federman (1999) for Mg ii $\lambda$ 2796 and Mg ii $\lambda$ 2803, we derive
a lower limit of $N>5\times 10^{13}\,\rm cm^{-2}$.

Table 6: Fe ii absorption line velocities during the photospheric phase Phase | Velocity | Phase | Velocity
---|---|---|---
(day) | $\left(\rm km\,s^{-1}\right)$ | (day) | $\left(\rm km\,s^{-1}\right)$
-1.4 | $8489\pm 88$ | 51.5 | $8371\pm 126$
8.8 | $8610\pm 28$ | 60.1 | $8382\pm 211$
10.5 | $8349\pm 64$ | 70.3 | $8303\pm 198$
30.8 | $8453\pm 121$ | 73.8 | $8313\pm 198$
31.1 | $8426\pm 205$ | 81.5 | $8417\pm 200$
32.7 | $8637\pm 168$ | 94.3 | $8431\pm 201$
47.1 | $8433\pm 218$ |  |
Figure 10: Maximum ejecta velocity. The extent of the Ca ii $\lambda$ 3934
(black) absorption on the blue side can be traced to $\sim
12,500\leavevmode\nobreak\ \rm km\,s^{-1}$ at $t_{\rm max}$+32.7 days. Mg ii
$\lambda$ 2796 (dark grey) has a comparable maximum velocity, albeit this
region is affected by additional SN features. The location of the blue and red
doublet components of Mg ii $\lambda\lambda$ 2796, 2803 and Ca ii
$\lambda\lambda$ 3934, 3968 of the host galaxy ISM are indicated by brackets
in a darker shade at the top of the figure. We also mark the position of the
doublets of the CSM shell with brackets in a lighter shade. The CSM shell is
detected through an additional Mg ii absorption-line system blueshifted by
$2918\leavevmode\nobreak\ \rm km\,s^{-1}$. The CSM shell is not detected in Ca
ii.

The only other SLSN that showed such a blueshifted Mg ii component was the
H-poor SLSN iPTF16eh (Lunnan et al. 2018b). For that SLSN, the Mg ii doublet
was blueshifted by 3300 $\leavevmode\nobreak\ \rm km\,s^{-1}$. Lunnan et al.
(2018b) also detected Mg ii in emission between 100 and 300 days after maximum
light. Moreover, the line centre of the emission lines moved from -1600 to
+2900 $\leavevmode\nobreak\ \rm km\,s^{-1}$ during that time interval. These
authors attributed the blueshifted Mg ii absorption line system with a CSM
shell expelled decades before the explosion and the time and frequency
variable Mg ii emission lines with a light echo from that shell. How such a
light echo evolves depends mainly on its distance to the progenitor star. With
that in mind, we analyse the Keck and X-shooter spectra between $t_{\rm
max}$+230 and $t_{\rm max}$+378 days to constrain the properties of the CSM
shell. Rebinning the spectra reveals Mg ii in emission (Figures 6 and 7).
However, due to heavy rebinning, the information about the variability of the
line centre was lost. We can, therefore, not ascertain whether the Mg ii
emission is connected with illuminated magnesium in the CSM shell or produced
by the interaction of the SN ejecta with circumstellar material.

Motivated by the discovery of a CSM shell around SN 2018ibb, we next search
for corresponding Ca ii $\lambda\lambda$ 3934,3969 absorption in the X-shooter
spectrum from $t_{\rm max}$+32.7 (Figure 10). The search is aggravated by how
the two Ca ii doublets (CSM shell and SN ejecta) overlap in contrast to the Mg
ii doublets. Using the wavelength of Ca ii $\lambda$ 3934 as the velocity
reference, the blue doublet absorption of Ca ii should be at the same velocity
as the Mg ii doublet ($-2918\leavevmode\nobreak\ \rm km\,s^{-1}$). The red
component of the Ca ii doublet will, however, be displaced by 34.8 Å, or 2653
$\leavevmode\nobreak\ \rm km\,s^{-1}$, to $\sim-265\leavevmode\nobreak\ \rm
km\,s^{-1}$. The position of a possible Ca ii CSM component is marked with the
light grey bracket in Figure 10. We do indeed see a sharp drop in the Ca ii
profile at zero velocity, which could be the result of a red CSM absorption.
For the blue component, it is more difficult because we do not know the line
profile of the Ca ii absorption from the SN ejecta. Therefore, it is difficult
to assess the significance of this. However, we conclude that there is no
evidence for Ca ii absorption from the CSM shell.

####  4.3.4 Circumstellar interaction — bumps and undulations in the light
curve

The multi-band light curves show a series of bumps and wiggles throughout the
entire evolution of SN 2018ibb (Figures 3 and 5). Between $t_{\rm max}$ and
$t_{\rm max}$+100 days, the bumps are well visible from $u$ to $H$ band
(luminosity increases by a few 0.1 mag). The amplitudes of the bumps in SN
2018ibb are comparable to the bumps seen in light curves of the other SLSNe
(e.g. Nicholl et al. 2016b; Inserra et al. 2017; Fiore et al. 2021;
Hosseinzadeh et al. 2022; Chen et al. 2023b). Following the nomenclature in
Chen et al. (2023b), these bumps fall in the ‘weak’ category. The bumps in SN
2018ibb also introduce wiggles in the evolution of its blackbody radius and
temperature (Figure 5). These modulations are well within the measurement
uncertainties of the long-term trends of these parameters, hindering a more
in-depth analysis of these features.

The late-time photometric evolution of SN 2018ibb reveals an increase in
luminosity of 0.2 dex between $t_{\rm max}$+240 and $t_{\rm max}$+340 days
(Figures 3 and 5) that is well isolated allowing for a more in-depth analysis.
The bolometric light curve before and after this bump exhibits a decline rate
of 1.18 mag (100 days)-1. After subtracting the underlying fading light curve,
we conclude that the light curve bump lasted for $\sim 80$ days (measured
between zero intensity) and reached its highest luminosity at $t_{\rm
max}$+300 days. In total, $6.7\pm 0.8\times 10^{48}$ erg are radiated in
excess to the $8.1\times 10^{49}$ erg that SN 2018ibb would have emitted
without the bump during this time.

Figure 11: Fe ii ejecta velocities of SN 2018ibb and general SLSN samples
(grey) at the time of maximum (top panel) and as a function of time (bottom
panel). SN 2018ibb has a remarkably low velocity at the time of maximum and an
unprecedentedly flat velocity evolution, which is in stark contrast to known
SLSNe.  Figure 12: Normalised X-shooter spectra from $t_{\rm max}$+32.7 days
to $t_{\rm max}$+94.3 days (top panels) and their inverse-variance weighted
co-added spectrum (bottom panel). The individual and stacked spectra show
barely resolved, narrow absorption lines from the host ISM (marked by the
solid vertical lines). In addition, a blue-shifted (2918 $\leavevmode\nobreak\
\rm km\,s^{-1}$) absorption line system is visible (marked by the dashed
vertical lines). The FWHMs of the blue-shifted component are 406 km s-1,
significantly larger than the ISM lines but significantly smaller than the SN
lines. This blue-shifted absorption-line system is connected with a shell of
circumstellar material expelled by the progenitor star shortly before the
explosion. No significant evolution in the position or shape of the absorption
lines can be seen in the individual spectra (upper panels). The error spectrum
is shown in grey.

Figure 13 presents the spectroscopic evolution of SN 2018ibb during the bump
phase. Assuming that all spectral features fade on exponential timescales
similar to the multi-band and bolometric light curves, we use the spectra
obtained before (blue) and after (yellow) the light curve bump to interpolate
the spectrum at $t_{\rm max}$+286.7 days (black). Such an approach estimates
the spectroscopic behaviour of SN 2018ibb in the absence of the bump. The
bottom panel of Figure 13 shows the observed spectrum at $t_{\rm max}$+286.7
days in black and the estimated spectrum without the bump in red. The
difference spectrum (blue) reveals substantially enhanced line fluxes in [O
ii] and [O iii] but no change in [O i]. The light-curve bump might also have
increased the flux of the continuum level bluewards of 5000 Å. Its shape is
reminiscent of the blue pseudo-continuum seen in interaction-powered SNe
(Silverman et al. 2013; Hosseinzadeh et al. 2017; Gal-Yam et al. 2022; Perley
et al. 2022). Considering the similarity of the difference spectrum to the
spectrum before and after the bump raises the question of whether a larger
fraction of the emission bluewards of 5000 Å in all nebular spectra is due to
CSM interaction. We investigate that further in Section 5.2.4.

Figure 13: Impact of the light curve bump between $t_{\rm max}$+240 and
$t_{\rm max}$+340 days on the SN spectrum at $t_{\rm max}$+286.7 days. Top:
The spectra before, during and after the light curve bump. Bottom: The
observed spectrum at $t_{\rm max}$+286.7 days ($\approx 13$ days before the
peak of the bump) is shown in black. We estimate the ‘bump-free’ spectrum of
SN 2018ibb at $t_{\rm max}$+286.7 days (red) based on the spectra obtained
before and after the bump. The difference between the observed (black) and
interpolated (red) spectra at $t_{\rm max}$+286.7 days is shown in blue. It
reveals a series of emission lines that can be attributed to [O ii] and [O
iii]. An excess bluewards of 5000 Å is also visible, while no apparent
residual can be seen at the location of [O i].

###  4.4 Radio and X-ray emission

The interaction of the SN ejecta with circumstellar material and heating of
the SN ejecta by a central engine (e.g. magnetar or a black hole) can produce
thermal X-ray emission and non-thermal radio emission (Chevalier & Fransson
1992, 1994). SN 2018ibb was observed in the X-rays and radio between $t_{\rm
max}$+13 and $t_{\rm max}$+246 days (Sections 3.5 and 3.6). All observations
led to non-detections with detection limits between 1 and $6\times
10^{41}\,\rm erg\,s^{-1}$ in the X-rays and between $10^{39}$ and
$10^{40}\,\rm erg\,s^{-1}$ in the radio. To put those measurements in the
context of the UV-to-NIR bolometric light curve, we show the radio and X-ray
measurements together with the bolometric light curve in Figure 14. From that,
we conclude that $<2\%$ and $<10\%$ of the total emission are radiated in the
radio and X-rays, respectively.

Figure 14: Thermal and non-thermal emission of SN 2018ibb. Less than a few
percent of the total radiated energy is emitted in the radio and X-rays. The
luminosity limits lie in the ballpark of non-detections of other SLSNe, and
they are a factor of 50 larger than the luminosity of the four H-poor SLSNe
with either radio or X-ray detection. The limits of SN 2018ibb are larger than
the most luminous radio and X-ray SNe.

The non-detection limits are in the observed range of other SLSNe with X-ray
and radio observations (Levan et al. 2013; Coppejans et al. 2018; Margutti et
al. 2018; Law et al. 2019; Eftekhari et al. 2021; Murase et al. 2021). Only
four SLSNe were detected at X-ray or radio frequencies: PTF10hgi (radio;
Eftekhari et al. 2019; Law et al. 2019), PTF12dam (X-ray; Margutti et al.
2018; Eftekhari et al. 2021), SCP06P6 (X-ray; Levan et al. 2013), and SN
2020tcw (radio and X-ray; Coppejans et al. 2021; Matthews et al. 2021). Their
measurements202020PTF10hgi was detected in the radio $>7.5$ years after the SN
explosion. Owing to this, we omit to show PTF10hgi in that figure., shown in
Figure 14, are a factor of $>50$ smaller than the detection limits of SN
2018ibb.

To put the radio and X-ray properties of SN 2018ibb in the context of
interaction-powered SNe, we also show the light curves of the most luminous
X-ray and radio SNe in Figure 14. The Type IIn SNe 2006jd and 2010jl are the
most luminous X-ray SNe with absorption-corrected luminosities of $\sim
10^{42}\rm\,erg\,s^{-1}$ (Chandra et al. 2012, 2015). The radio-loudest SNe
(e.g. SN Ic-BL PTF11qcj Corsi et al. 2014) reached luminosities of $\sim
10^{38}\rm\,erg\,s^{-1}$, i.e. $\lesssim 10$ times fainter than the limits for
SN 2018ibb. Their observed luminosities before correcting for host absorption
can be significantly dimmer for hundreds of days (Chandra et al. 2015).

In conclusion, the non-detection of SN 2018ibb neither rules out CSM
interaction nor a central engine as the dominant powering mechanism.
Furthermore, the non-detection of SN 2018ibb also agrees with theoretical
models of magnetar- and interaction-powered SLSNe that predict no bright radio
and X-ray emission for years after the SN explosion (Murase et al. 2016;
Margalit et al. 2018; Omand et al. 2018).

###  4.5 Imaging polarimetry

Our polarimetric observations between $t_{\rm max}$+31.9 days and $t_{\rm
max}$+94.4 days revealed a polarisation signal of $0.27\pm
0.04\leavevmode\nobreak\ \%$ in $V$ (weighted average of all epochs) and
$0.48\pm 0.07\leavevmode\nobreak\ \%$ in the $R$ band (Table 3). Dust grains
in the Milky Way and the host galaxy could introduce a polarisation signal. As
detailed in Section 3.4, the polarisation level of the MW could be up to
0.26%. The level of polarisation from the SN host galaxy is unknown, meaning
that all reported measurements are upper limits.

Considering the observed low degree of polarisation and the consistent levels
of Stokes parameters measured from SN 2018ibb (Table 3), we conclude that the
continuum polarisation intrinsic to SN 2018ibb is $\lesssim 0.3\%$ in $V$ band
between $t_{\rm max}$+31.9 days and $t_{\rm max}$+94.4 days. To convert this
measurement into an asphericity of the ejecta, we assume an oblate ellipsoidal
ejecta with a Thomson scattering atmosphere and a number density distribution
of $N(r)\propto r^{-n}$, where $r$ is the ejecta radius and $n$ is the power-
law index. Adopting $p\lesssim 0.3\%$, we infer an axis ratio B/A (minor axis
vs. major axis) of $\gtrsim 0.9$ for an optical depth of $\tau=1$ and a power-
law index of $n=2$, and B/A of $\gtrsim 0.8$ for $\tau=5$ and $n=3$–5 (Höflich
1991). The degree of polarisation in the $R$ band is slightly higher
($p\approx 0.5\%$). Therefore, we cannot exclude that the continuum
polarisation is $p>0.3\%$. A polarisation degree $p\sim 0.5\%$ implies an axis
ratio B/A of $\sim 0.88$ for $\tau=1$ and $n=2$ (Höflich 1991).

Therefore, we suggest that SN 2018ibb’s photosphere exhibits a high degree of
spherical symmetry. Pursiainen et al. (2023) analysed the data of the 16
SLSNe-I with polarimetric observations, including SN 2018ibb. After correcting
the phases of all objects for the diverse photometric decline rates, the
properties of SN 2018ibb are well within the observed distribution. While some
of the events exhibit a non-zero level of polarisation at similar phases to SN
2018ibb (e.g. SN 2015bn and SN 2021fpl; Leloudas et al. 2017; Inserra et al.
2016; Poidevin et al. 2023), most SLSNe show a consistently low polarisation
degree at comparable normalised phases (see figure 6 in Pursiainen et al.
2023).

The presence of any component in the atmosphere of SN 2018ibb significantly
deviating from spherical symmetry is thus unlikely within the photospheric
phases covered by VLT polarimetry observations. Although Thomson scattering is
wavelength independent, broad emission lines (see spectra in Figure 6), which
are in general not polarised, may dominate the polarisation spectrum in the
$V$ band and produce the apparent low polarisation values. Furthermore, iron-
group elements in the ejecta (Figure 7) have a large number of bound-bound
transitions in the blue and UV part of the spectrum, which can also depolarise
the signal (e.g. Chornock & Filippenko 2008), accounting for the slightly
different polarisation levels measured in $V$ and $R$ bands.

###  4.6 Host galaxy

Table 7: Properties of the interstellar medium in the host galaxy

Transition | ${\rm EW}_{\rm r}$ | Flux
---|---|---
| (Å) | $\left(10^{-18}\,{\rm erg\,cm}^{-2}\,{\rm s}^{-1}\right)$
Absorption lines
Mn ii $\lambda$ 2594 | $0.18\pm 0.13$ | …
| ($<0.39$) |
Fe ii $\lambda$ 2600 | $0.07\pm 0.13$ | …
| ($<0.39$) |
Mn ii $\lambda$ 2606 | $-0.12\pm 0.15$ | …
| ($<0.45$) |
Mg ii $\lambda$ 2796 | $0.51\pm 0.04$ | …
Mg ii $\lambda$ 2804 | $0.46\pm 0.04$ | …
Mg i $\lambda$ 2852 | $0.14\pm 0.04$ | …
| ($<0.16$) |
Ca ii $\lambda$ 3934 | $0.03\pm 0.01$ | …
Ca ii $\lambda$ 3969 | $0.01\pm 0.01$ | …
| ($<0.03$) |
Emission lines
H$\beta$ | … | $3.68\pm 0.78$
[O iii] $\lambda$ 4363 | … | $0.25\pm 0.16$
| … | ($<0.48$)
[O iii] $\lambda$ 4959 | … | $1.96\pm 0.81$
| … | ($<2.43$)
[O iii] $\lambda$ 5007 | … | $12.96\pm 1.11$
H$\alpha$ | … | $10.43\pm 0.76$
[Nii] $\lambda$ 6584 | … | $0.12\pm 0.54$
| … | ($<1.62$)
212121We report rest-frame equivalent widths EWr for absorption lines and
fluxes for emission lines. For lines detected with a significance of
$<3\sigma$, we also report the $<3\sigma$ upper limits. The rest-frame
equivalent widths are measured by averaging the measurements from the
X-shooter spectra between $t_{\rm max}$+32.7 days and $t_{\rm max}$+94.3 days
of each line species. The emission lines are measured from the FORS2 spectrum
at $t_{\rm max}$+637 days. The emission-line fluxes are not corrected for
reddening.

SN 2018ibb’s host galaxy was detected in several optical broad-band filters
($m_{R}\sim 24.4$ mag; Table 1). A false colour image of the field is shown in
Figure 1. The SN explosion site, marked by the crosshair, is $\approx 1$ kpc
from the centre of its host galaxy, a common offset for SLSNe (Lunnan et al.
2014; Schulze et al. 2018, 2021). To infer the mass and star-formation rate of
the host, we model the observed spectral energy distribution (black data
points in Figure 15) with the software package Prospector version 1.1 (Johnson
et al. 2021).222222Prospector uses the Flexible Stellar Population Synthesis
(FSPS) code (Conroy et al. 2009) to generate the underlying physical model and
python-fsps (Foreman-Mackey et al. 2014) to interface with FSPS in python. The
FSPS code also accounts for the contribution from the diffuse gas based on the
Cloudy models from Byler et al. (2017). We use the dynamic nested sampling
package dynesty (Speagle 2020) to sample the posterior probability. We assume
a Chabrier IMF (Chabrier 2003) and approximate the star formation history
(SFH) by a linearly increasing SFH at early times followed by an exponential
decline at late times [functional form $t\times\exp\left(-t/t_{1/e}\right)$,
where $t$ is the age of the SFH episode and $t_{1/e}$ is the $e$-folding
timescale]. The model is attenuated with the Calzetti et al. (2000) model. The
priors of the model parameters are set identical to those used by Schulze et
al. (2021). The observed SED is adequately described by a galaxy model with a
stellar mass of $\log\,M_{\star}/M_{\odot}=7.60^{+0.19}_{-0.22}$ and star-
formation rate of $0.02^{+0.04}_{-0.01}\leavevmode\nobreak\ M_{\odot}\,{\rm
yr}^{-1}$ (grey curve in Figure 15).

Figure 15: Spectral energy distribution of the host galaxy from 1000 to 60,000
Å (black dots). The solid line displays the best-fitting model of the SED. The
red squares represent the model-predicted magnitudes. The fitting parameters
are shown in the upper-left corner. The abbreviation ‘n.o.f.’ stands for the
number of filters.

The mass and the star-formation rate of the host of SN 2018ibb agree with the
expected values of SLSNe-I host galaxies at $z<0.3$ (Leloudas et al. 2015;
Perley et al. 2016a; Chen et al. 2017; Schulze et al. 2018, 2021), although
both fall in the lower half of the distributions. The specific star-formation
rate (SFR normalised by the stellar mass of the host) is comparable to a
common star-forming galaxy of that stellar mass (grey band in Figure 16; Elbaz
et al. 2007) but in the lower half of the observed distribution of SLSN host
galaxies (Schulze et al. 2021). We caution that specific SFRs are notoriously
difficult to measure (e.g. see figure 3 in Schulze et al. 2021) as they rely
on well-sampled SEDs from the UV to the NIR.

The X-shooter spectra up until $T_{\rm max}+80$ days reveal narrow absorption
lines from Mg i and Mg ii from the interstellar medium in the host galaxy but
no absorption features from Ca ii, Fe ii, and Mn ii, which have prominent
features in the wavelength range accessible with X-shooter and are typically
seen in low-mass star-forming galaxies, e.g. Prochaska et al. (2007) and Fynbo
et al. (2009). The equivalent widths of the detected lines and the upper
limits of the strongest expected absorption features are reported in Table 7.
The measurements of Mg i $\lambda$ 2852 and Mg ii $\lambda\lambda$ 2796, 2804
are comparable to those of the SLSN host galaxies reported in Vreeswijk et al.
(2014). Following the methodology of de Ugarte Postigo et al. (2012), we infer
an absorption-line strength parameter of $\sim-3.5$ from Ca ii, Mg i and Mg
ii, putting the host of SN 2018ibb at the low-metallicity end of the
distribution (albeit the diagnostic is tailored to host galaxies of long-
duration gamma-ray bursts, which are also connected with the death of very
massive stars but which prefer galaxies with slightly higher metallicities and
slightly older stellar populations than SLSNe-I; Hjorth & Bloom 2012; Leloudas
et al. 2015; Vergani et al. 2015; Perley et al. 2016b; Schulze et al. 2018).

Figure 16: Star-formation rate and stellar mass of the host galaxy of SN
2018ibb in the context of SLSN-I host galaxies from the PTF survey (Schulze et
al. 2021). The host galaxy of SN 2018ibb lies in the expected parameter space
of SLSN host galaxies but in the lower half of the mass and SFR distributions
(kernel density estimates of the observed distributions are shown at the top
and to the right of the figure). Its specific star-formation rate (SFR / mass)
is comparable to the typical star-forming galaxies (grey band) but lower than
for an average SLSN host galaxy.

SN 2018ibb’s nebular spectra exhibit emission lines from hydrogen and oxygen
from H ii regions in the host galaxy. We measure their intensities by
integrating over their line profiles. To apply emission-line diagnostics for
measuring the oxygen abundance, we also need the flux of [N ii] $\lambda$
6584, which evaded detection. Using the H$\alpha$ line profile as a template
of the [N ii] $\lambda$ 6584 line profile, we measure the nominal flux and its
uncertainty. Table 7 summarises all measurements. Using the O3N2 metallicity
indicator with the calibration from Marino et al. (2013) yields a low oxygen
abundance of $12+\log\,\left({\rm O/H}\right)=8.06^{+0.07}_{-0.11}$ in
accordance with the low value from the absorption-line strength parameter. The
oxygen abundance is comparable to the mean of SLSN host galaxies at similar
redshifts (Leloudas et al. 2015; Perley et al. 2016a; Chen et al. 2017). The
non-detection of [N ii] $\lambda$ 6584 and [O iii] $\lambda$ 4363 adds a
systematic uncertainty to the inferred metallicity. However, the average
metallicity of a galaxy is correlated with its mass. Andrews & Martini (2013)
reported the mass-metallicity relation for star-forming galaxies between
$z=0.027$ and 0.25, a stellar mass of $\log\leavevmode\nobreak\
M_{\star}/M_{\odot}=7.4$–10.5, and metallicities measured with the $T_{e}$
method. Using this mass-metallicity relation yields an oxygen abundance of
$12+\log\,\rm O/H\approx 8.07^{+0.10}_{-0.12}$, identical to the value
inferred from the galaxy emission lines. Assuming a solar oxygen abundance of
8.67 (Asplund et al. 2009), the host galaxy metallicity is
$0.25^{+0.07}_{-0.06}$ solar.

The flux ratio between H$\alpha$ and H$\beta$ is $2.76\pm 0.62$, which is
consistent within $1\sigma$ with the theoretically expected value of 2.86 for
no extinction (assuming a temperature of $10^{4}$ K and an electron density of
$10^{2}\leavevmode\nobreak\ {\rm cm}^{-3}$ for Case B recombination;
Osterbrock 1989). We conclude that the host attenuation is negligible. The
H$\alpha$ flux translates to a star-formation rate of ${\rm SFR}=4.4\pm
0.3\times 10^{-3}\leavevmode\nobreak\ M_{\odot}\,{\rm yr}^{-1}$ using
Kennicutt (1998) and the relation from Madau & Dickinson (2014) to convert
from the Salpeter to the Chabrier IMF in the Kennicutt (1998) relation. This
value is lower than the SFR estimated from the host SED fitting but consistent
within $2\sigma$.

##  5 Discussion

###  5.1 SN ejecta emission vs. CSM interaction

In Section 4.3.3, we have shown that the progenitor of SN 2018ibb is embedded
in circumstellar material ejected shortly before the explosion. In this
section, we examine the line profiles and evolution of selected oxygen and
metal lines to infer the physical conditions of the SN ejecta and the CSM.

The line profiles are most clear in the nebular phase. Figure 17 shows the
continuum-subtracted Mg i] $\lambda$ 4571 line with the [O i] $\lambda\lambda$
6300, 6364 doublet at $t_{\rm max}$+286.7 days. Both lines extend to $\sim
10,000\leavevmode\nobreak\ \rm km\,s^{-1}$. Their maximum velocity hardly
changes up to the last well-observed epoch at $t_{\rm max}$+637.3 days. Its
similarity to the maximum velocity of the Ca ii $\lambda$ 3934 absorption line
(Figure 10) suggests Mg i] and [O i] are produced in the high-velocity ejecta.
The Mg i] line is well fitted with a parabolic line profile with similar
maximum velocity, shown by the dark-blue line in Figure 17. This indicates
emission from an optically thick shell with constant velocity (e.g. Fransson
1984).

A similar parabolic line profile is consistent with the red side of the [O i]
$\lambda$ 6364 doublet component. However, the blue side of the doublet,
dominated by the 6300 Å component, lacks most of the emission compared to the
Mg i] line. By $t_{\rm max}$+637.3 days (Figure 18), the blue doublet
component has grown and is now the stronger of the two lines. The evolution of
the [O i] line profile may be explained if both doublet components are
optically thick to at least $t_{\rm max}$+286.7 days. In that case, the blue
component will be scattered by the red component, which extends over most of
the blue component (the velocity difference between the two components is
$3016\leavevmode\nobreak\ \leavevmode\nobreak\ \rm km\,s^{-1}$). These photons
will either be thermalised or emerge on the red side of the 6364 Å doublet
component. Emission from the front side of the ejecta with velocities
$\lesssim-(v-3016\leavevmode\nobreak\ \rm km\,s^{-1}$) are only partially
scattered, and some of this emission may leak out, explaining the ‘bump’ at
$\sim-7500\leavevmode\nobreak\ \leavevmode\nobreak\ \rm km\,s^{-1}$. At
$t_{\rm max}$+565.0 days, the blue doublet component has grown, and the blue
wing is equally bright, or somewhat brighter, compared to the 6364 Å
component. This trend continues at $t_{\rm max}$+637.3 days. The expected 3:1
ratio is still not reached, indicating that the ejecta is not optically thin,
yet.

Figure 17: [O i] $\lambda\lambda$ 6300,6364 and Mg i] $\lambda$ 4571 lines at
$t_{\rm max}$+286.7 days. The Mg i] line is well fitted with a parabolic shape
(dark blue), expected from an optically thick expanding shell (e.g. swept up
CSM and unshocked SN ejecta), while the [O i] lines show a strong blue deficit
because the line-forming region is still optically thick. The [O i] doublet is
centred on the 6364 Å doublet component.

Using the Sobolev (1957) theory for the line formation, we can estimate the
optical depth $\tau$ for a given O i density $n({O\textsc{i}})$ in the ejecta
(e.g. Li & McCray 1992). Assuming LTE among the 3P ground state levels, the
optical depth of each line of the [O i] doublet is given by

$\tau=\frac{A\left({}^{1}{\rm D_{2}},{}^{3}{\rm
P}_{J}\right)\,\lambda\left({}^{1}{\rm D_{2}},{}^{3}{\rm
P}_{J}^{3}\right)^{3}\,g\left({}^{2}{\rm D}_{2}\right)}{8\pi\,g_{\rm tot}}\,\
n({\rm{O\textsc{i}}})\times t$

where $A({}^{1}{\rm D_{2}},{}^{3}{\rm P}_{J})$ with $J=2,1$ are the transition
probabilities for the 6300 Å and 6364 Å lines, respectively,
$\lambda\left({}^{1}{\rm D_{2}},{}^{3}{\rm P}_{J}^{3}\right)$ is the
wavelength of blue and red doublet component, respectively, $g_{\rm tot}=9$ is
the total statistical weight to the ground multiplet, and $t$ is the time
since the explosion, in units of day. Putting in the atomic constants, we get
an optical depth for the 6364 Å line of

$\tau=2.7\ \left(\frac{n({\rm{O\textsc{i}}})}{10^{10}\ {\rm
cm}^{-3}}\right)\left(\frac{t}{300\rm\ day}\right)$

and a depth that is a factor of 2.9 larger for the 6300 Å line. The typical O
i density needed to get an optically thick line is, therefore, $\ga 10^{9}$
cm-3. This can be compared to the mean oxygen density of the core. Assuming O
i is the dominant species of oxygen in the core, the number density is

$\displaystyle n({\rm O})$ $\displaystyle\approx 3\times
10^{7}\,f^{-1}\,\left(\frac{M({\rm
O})}{30\,M_{\odot}}\right)\,\left(\frac{v_{\rm ej}}{10^{4}\leavevmode\nobreak\
\rm km\,s^{-1}}\right)^{-3}\,\left(\frac{t}{300\ \rm day}\right)^{-3}\rm
cm^{-3}$ (1)

where $M(\rm O)$ is the mass of oxygen in the core, $f$ is the filling factor,
and $v_{\rm ej}$ the ejecta velocity. To get an optically thick 6364 Å line at
$t_{\rm max}$+286.7 days, in other words a density $\ga 10^{10}$ cm-3,
requires a very small oxygen filling factor, $\lesssim 10^{-3}$ (i.e. a highly
clumped medium), or an unphysically large oxygen mass of $\ga
10^{3}M_{\odot}$.

Figure 18: Late-time evolution of the [O i] $\lambda\lambda$ 6300,6364
doublet. Note the strong evolution on the blue side, while the red side of the
lines is evolving slower. This indicates a transition from optically thick to
optically thin [O i] lines, implying that the scattering in the absorption
part of the P-Cygni profile is decreasing. Regions of strong atmospheric
absorption are grey-shaded.

In a CSM and PPISN interaction scenario, the continued high optical depth at
late times could point to a highly compressed cool dense shell (CDS). A CDS
will form from the compression behind the shock, which results from the
interaction between the ejecta and the CSM (for a discussion see, e.g.
Chevalier & Fransson 2017). If the density of the CSM is large, the forward
shock will be radiative and dominate the emission, which will then have the
composition of the CSM. In the opposite case, the reverse shock dominates with
a composition typical of the outer ejecta. The latter case is more relevant
for lower mass loss rates. In both cases, the density enhancement behind the
cooling shock will be very large. Assuming an approximate pressure balance
behind the shock, the density enhancement will be of the order of $T_{\rm
shock}/T_{\rm ps}\approx 3/16\,\mu\,m_{\rm u}\,v_{\rm rel}^{2}/(k\,T_{\rm
ps})$, where $\mu$ is the mean molecular weight ($\sim 1.7$ for a fully
ionised oxygen gas), $m_{\rm u}$ the atomic mass unit, $v_{\rm rel}$ the
relative velocity between the CSM shell and the ejecta, $T_{\rm shock}$ the
temperature immediately behind the shock, and $T_{\rm ps}$ the post-shock
temperature in the CDS ($\sim 10^{4}$ K). With $v_{\rm rel}\approx
5000\leavevmode\nobreak\ \rm km\,s^{-1}$232323The photospheric velocity of the
ejecta is $<8500\leavevmode\nobreak\ \leavevmode\nobreak\ \rm km\,s^{-1}$ at
late times (Section 4.3.2), and the velocity of the CSM shell is
$2918\leavevmode\nobreak\ \leavevmode\nobreak\ \rm km\,s^{-1}$ (Section
4.3.3)., the compression is of the order of $10^{5}$. Both $v_{\rm rel}$ and
$T_{\rm ps}$ are uncertain and magnetic pressure could limit the compression.
The CDS is also most likely unstable (Chevalier & Blondin 1995), leading to
clumping of the shell and limiting of the compression. However, the estimate
shows that a very large density could result in the CDS, making the line
optically thick, equivalent to a low filling factor. In the PISN scenario,
strong clumping in the ejecta is needed. This is, however, not indicated from
simulations of PISN models without CSM by Chen et al. (2020).

We now turn to the origin of the higher ionisation [O ii] and [O iii] lines.
Figure 19 shows the line profiles of the [O i], [O ii] and [O iii] lines after
subtracting the continuum. Owing to the doublet nature of the lines, we centre
the line profiles on the blue component in the left panel and on the red
component in the right panel. These line widths can be compared to the
velocity of the CSM shell, the photospheric velocity, and the maximum velocity
of the ejecta (vertical lines in Figure 19). It is clear that the [O ii] and
[O iii] line widths are closer to the velocity of the CSM shell than to the
photospheric velocity of the SN ejecta; in contrast to the [O i] line, which
extends to the maximum velocity of the SN ejecta.

The differences in the origin of the forbidden oxygen lines are corroborated
by the O i-O iii line profiles (Figures 19, 20). The asymmetric [O iii]
$\lambda$ 4363 and [O iii] $\lambda\lambda$ 4959,5007 lines have little
emission in the red wings indicative of emission from a thin shell, where most
of the red emission is absorbed by the photosphere. Examples of this can be
seen in figure 4b in Fransson (1984). That scenario is also consistent with
the evolution of the [O ii] $\lambda\lambda$ 7320,7330 doublet. The 7300 Å
line, which may be a blend of the [O ii] $\lambda\lambda$ 7320,7330 lines and
[Ca ii] $\lambda\lambda$ 7291,7324, is shown in Figure 21, centred on the [O
ii] $\lambda$ 7320 line. Focusing on the [O ii] lines (left panel), the blue
wing has a nearly constant line profile between $t_{\rm max}$+231.2 and
$t_{\rm max}$+565.0 days. The red wing of the [O ii] $\lambda$ 7330 doublet
component gets considerably narrower during the same time interval. At $t_{\rm
max}$+637.3 days (right panel), the entire 7300 Å line profile changes quite
dramatically, becoming flat-topped and broader. This is a result of the [Ca
ii] $\lambda\lambda$ 7291,7324 lines becoming strong, while the [O ii] lines
get weaker.

Figure 19: [O i] $\lambda\lambda$ 6300,6364, [O ii] $\lambda\lambda$
7320,7330, and [O iii] $\lambda\lambda$ 4959,5007 line profiles at $t_{\rm
max}$+565.0 days. The velocity scale is centred on the blue doublet component
in the left panel and on the red component in the right panel. [O ii] and [O
iii] only reach out to approximately the velocity of the CSM shell ($v_{\rm
CSM}$), much less than the photospheric velocity ($v_{\rm FeII}$) and the
maximum velocity of the ejecta ($v_{\rm max}$). In contrast to that, the [O i]
profile extends to $\approx 12,500\leavevmode\nobreak\ \rm km\,s^{-1}$. This
points to [O ii] and [O iii] being produced close to the CSM shell whereas [O
i] is produced in the SN ejecta.

The increasing asymmetry of the [O ii] doublet may be qualitatively understood
by the CSM being occulted by the SN ejecta. Assuming that the optically thick
SN ejecta with the velocity $v_{\rm ej}$ slams into the CSM shell with a low
velocity (ideally $v\approx 0$) located at a distance $R_{\rm s}$ from the
progenitor star, the maximum velocity of the red wing $v_{\rm red}$ is (a pure
geometric effect)

$\begin{split}v_{\rm red}&=v_{s}\,\left(1-\left(v_{\rm ej}\,t/R_{\rm
s}\right)^{2}\right)\approx v_{\rm s}\,\left(1-\left(v_{\rm ej}\,t/v_{\rm
s}\left(t+\Delta t\right)\right)^{2}\right)\\\ &\approx v_{\rm
s}\left(1-\left(v_{\rm ej}\,t/v_{\rm s}\Delta t\right)^{2}\right)\end{split}$

where $\Delta t$ is the time between the shell ejection and the explosion and
$t$ the time since explosion. (We have assumed that $t\ll\Delta t$.) Because
the ejecta with the CDS, which may define the photosphere, expands with a much
higher velocity ($\sim 8,500\leavevmode\nobreak\ \rm km\,s^{-1}$ vs. $\sim
3000\leavevmode\nobreak\ \rm km\,s^{-1}$) a progressively increasing portion
of the dense CSM will be occulted by the photosphere and less of the
‘backside’ of the CSM will be seen. This would lead to the red side getting
narrower with time. At the same time, an increasing portion of the dense CSM
will be shocked, leading to a decreasing luminosity from the dense CSM,
including the [O ii-iii] emission.

The fact that the forbidden [O ii] and [O iii] lines are seen at about $t_{\rm
max}$+30 days adds additional constraints on the physical conditions where
they originate. The critical densities, above which collisional de-excitation
becomes important, are less than $\sim 2\times 10^{6}\,\rm cm^{-3}$ for [O ii]
and [O iii] (Osterbrock & Ferland 2006). This is much lower than the densities
expected in the ejecta (Equation 1). Therefore, the [O ii] and [O iii] lines
would be severely suppressed if they were coming from the ejecta. Not only do
[O ii] and [O iii] originate from the CSM, but also the recombination lines O
i $\lambda$ 7773 and O i $\lambda$ 9263\. The blue wings of their line
profiles are similar to [O ii] $\lambda$ 7320, extending to $\sim
5000\leavevmode\nobreak\ \rm km\,s^{-1}$ (Figure 22).

Figure 20: Comparison of the [O iii] $\lambda\lambda$ 4959,5007 lines with the
[O iii]$\lambda$ 4363 line. These lines are produced by the interaction of the
SN ejecta with circumstellar material. Due to the occulation of the CSM by an
optically thick SN ejecta, less of the ‘backside’ of the CSM is seen. The
velocities of the CSM shell ($v_{\rm CSM}$), the photospheric velocity
($v_{\rm FeII}$) and the maximum velocity of the ejecta ($v_{\rm max}$) are
indicated.

In summary, we propose a two-component scenario where the broad component,
seen in particular in the [O i] $\lambda\lambda$ 6300,6364 and Mg i] $\lambda$
4571 lines as well as the broad absorption in Mg ii $\lambda$ 2800 and Ca ii
$\lambda$ 3934 come from either the CDS or possibly the unshocked ejecta. The
low-velocity component seen in the [O iii] lines, as well as the [O ii] and O
i recombination lines come from the CSM shell at $\sim
3000\leavevmode\nobreak\ \leavevmode\nobreak\ \rm km\,s^{-1}$. The fact that
we see [O iii] $\lambda\lambda$ 4959,5007 emission even at $t_{\rm max}$+989.2
days means that the dense CSM must extend out to at least a few $\approx
10^{17}$ cm. The velocity width of the Mg ii absorption of the CSM shell of
$406\leavevmode\nobreak\ \rm km\,s^{-1}$ (Section 4.3.3; Figure 12) may
correspond to the velocity gradient over the CSM shell. That we do not see any
change in the width with time suggests that this gradient must be small enough
so that the velocity close to the shock is nearly constant. The origin of this
gradient is not clear, though. One explanation could be that this is the
result of a time-limited eruption, where a Hubble-like outflow is expected
after a few dynamical time scales. This has been observed, for instance, in
the Eta Carinae Homunculus nebula produced during the great eruption in 1843
(e.g. Smith 2006). The absence of H and He lines throughout the entire
evolution reveals (Figure 6) that the CSM shell must be processed gas from the
stripped progenitor. Any hydrogen and helium must have been lost before this
eruption and reside at much larger radii.

Among the $>200$ H-poor SLSNe known, SN 2018ibb is only the seventh object
with spectroscopic evidence of CSM interaction. In previous cases, CSM
interaction did not manifest itself via [O iii] in emission (a possible
candidate for CSM interaction with O-rich material is PS1-14bj; Lunnan et al.
2016). iPTF16eh revealed CSM interaction through a light echo produced in a
shell of H-poor and He-poor material (Lunnan et al. 2016). Late-time spectra
of iPTF10aagc, 13ehe, 15esb and 16bad (Yan et al. 2015, 2017) and SN 2018bsz
(Pursiainen et al. 2022) showed broad Balmer emission lines, suggesting that
their progenitors lost their hydrogen envelopes much closer to the time of the
terminal explosion than SN 2018ibb and iPTF16eh.

Figure 21: Evolution of the [O ii] $\lambda\lambda$ 7320,7330 + [Ca ii]
$\lambda\lambda$ 7291,7324 line complex. Left: Up to $t_{\rm max}$+565.0 days,
the line complex is dominated by [O ii]. The red wing narrows due to an
increasing occultation of the CSM shell by the optically thick expanding SN
photosphere. Right: Between $t_{\rm max}$+565.0 days and $t_{\rm max}$+637.2
days days, the line complex shifts to the blue, consistent with the [Ca ii]
line becoming more dominant. All profiles are centred on [O ii] $\lambda$
7320.

###  5.2 Constraints on the powering mechanism and progenitor

In the following, we contrast SLSN and PISN models with our photometric and
spectroscopic datasets and discuss the most likely powering mechanism and
progenitor of SN 2018ibb.

####  5.2.1 Modelling the bolometric light curve

We first analyse the bolometric light curve. Katz et al. (2013) proposed an
exact method for testing whether a light curve is powered by the decay of
radioactive material and, therefore, allows us to place an upper limit on any
56Ni produced during the explosion of SN 2018ibb’s progenitor. This method is
independent of details in the radiative transport, including the highly
uncertain opacity, the velocity distribution and the ejecta geometry. The
method is described in detail in Wygoda et al. (2019) and Sharon & Kushnir
(2020). In brief, the Katz integral is given by

Figure 22: Comparison of the [O ii] $\lambda\lambda$ 7320,7330 lines with the
O i $\lambda$ 7773 and O i $\lambda$ 9263 recombination lines. The similar
line profiles, extending to $\lesssim 5000\leavevmode\nobreak\ \rm km\,s^{-1}$
indicates an origin in a highly processed CSM shell.

$\displaystyle QT$ $\displaystyle=$ $\displaystyle LT+ET\leavevmode\nobreak\
\rm with$ $\displaystyle QT$ $\displaystyle=$
$\displaystyle\int^{t}_{0}dt^{\prime}\,t^{\prime}\,Q_{\rm
dep}\left(t^{\prime}\right),\leavevmode\nobreak\
LT=\int^{t}_{0}dt^{\prime}\,t^{\prime}\,L\left(t^{\prime}\right)$

and $ET$ is the integrated time-weighted luminosity that would be emitted if
no 56Ni were produced. Assuming that there is no additional source of energy,
$ET$ can be assumed to be negligible. The total energy deposition rate from
radioactive decay of 56Ni, $Q_{\rm dep}$, is given by (Jeffery 1999)

$\displaystyle Q_{\rm dep}(t)$ $\displaystyle\approx
Q_{\gamma}\left(1-e^{\left(-t_{0}/t\right)^{2}}\right)+Q_{e^{+}}(t)$

where $t_{0}$ is the $\gamma$-ray escape time. The deposition rates from
$\gamma$-ray photons and positrons are

$\displaystyle Q_{\gamma}$ $\displaystyle=\frac{M({\rm
Ni})}{M_{\odot}}\,\left(6.45\,e^{-t/t_{1/2,\rm Ni}}+1.38\,e^{-t/t_{1/2,\rm
Co}}\right)\times 10^{43}\leavevmode\nobreak\ \rm erg\,s^{-1}$ $\displaystyle
Q_{e^{+}}$ $\displaystyle=4.64\frac{M({\rm
Ni})}{M_{\odot}}\,\left(-e^{-t/t_{1/2,\rm Ni}}+e^{-t/t_{1/2,\rm
Co}}\right)\times 10^{41}\leavevmode\nobreak\ \rm erg\,s^{-1}$

where the mean lifetimes of 56Ni and 56Co are $t_{1/2,\rm Ni}=8.76$ days and
$t_{1/2,\rm Co}=111.4$ days, respectively (Junde 1999).

Since the explosion time is not well known, we vary the explosion time between
0 and 50 rest-frame days before the first detection and use the relation
$L/LT=Q/QT$ to determine the $\gamma$-ray escape time. We measure a range of
600 to 700 rest-frame days for $t_{0}$. After the $\gamma$-ray escape time is
determined, we infer the nickel mass by comparing the luminosity in the fitted
range to the deposited radioactive energy. The best fit for each point in the
$t_{\rm exp}$ grid is shown in Figure 23. Indeed, the declining light curve is
fully consistent with being powered by 24–$35\leavevmode\nobreak\ M_{\odot}$
of 56Ni. The upper bound could be even larger if the SN explosion happened
more than 50 rest-frame days before the detection by Gaia.

The rise time of 90–140 days, the range of $\gamma$-ray escape times and the
range of nickel masses are consistent with expectations from PISN models
(Kasen et al. 2011; Kozyreva & Blinnikov 2015; Kozyreva et al. 2017) and SN
Ia242424Type Ia supernovae are a class of supernovae that occurs in binary
systems in which one of the stars is a white dwarf (e.g. Hillebrandt &
Niemeyer 2000, and references therein). This explosion produces $\sim
0.6\leavevmode\nobreak\ M_{\odot}$ of 56Ni (e.g. Maoz et al. 2014, and
references therein). (another class of thermonuclear explosions; e.g. Wygoda
et al. 2019; Sharon & Kushnir 2020), after scaling their average nickel mass
to the nickel mass of SN 2018ibb. Both the excellent match with nickel
powering and coverage of the fading light curve for 706 days is unprecedented
for any of the $>200$ SLSNe known, suggesting that SN 2018ibb could indeed be
a PISN.

Figure 23:  The bolometric light curve of SN 2018ibb from 1800 to 14,300 Å
(rest-frame) and light-curve fits to the fading light curve using the Katz et
al. (2013) method. The entire fading light curve from up to 706 days after
peak is fully consistent by being powered 24–$35\leavevmode\nobreak\
M_{\odot}$ of 56Ni (dark red), suggesting that SN 2018ibb could be a pair-
instability supernova. At about 300 days after peak (i.e. $\gtrsim 400$ days
after the explosion), the $\gamma$-ray trapping decreases with time. The loss
of trapping is indicated by the difference between the light-red ($100\%$
trapping) and dark-red ($<100\%$ trapping) curves.

####  5.2.2 Modelling the broadband light curve

Next, we fit the multi-band light curve with the Modular Open-Source Fitter
for Transients MOSFiT software tool (Guillochon et al. 2018). In addition to
the MOSFiT nickel model that is based on the parameterisation by Nadyozhin
(1994), we also select the central-engine models slsn (describing powering by
a spin down of a magnetar; Nicholl et al. 2017) and fallback (describing
powering by a black hole accreting fallback material; Moriya et al. 2018a),
and the Chatzopoulos et al. (2012) model to characterise the powering by CSM
interaction. We also utilise the more complex models magni combining powering
by a magnetar and radioactive 56Ni (Blanchard et al. 2019) and csmni which
combines powering by CSM interaction and 56Ni. In all models, the photosphere
is assumed to have a blackbody spectral energy distribution at all times.
While this approximation is adequate during the photospheric phase, it is
inadequate at later times when the spectrum is dominated by emission lines and
an interaction-powered pseudo-continuum. The spectral energy distribution of
the model slsn is modified in the UV to account for absorption by the SN
ejecta. As our dataset covers a very long time span, the trapping of
$\gamma$-ray photons will eventually decrease, which accelerates the fading.
All chosen models include a component to account for the loss of trapping
(Nicholl et al. 2017).

The priors of the model parameters are chosen to cover a broad range of the
physically allowed parameter spaces. Their ranges and shapes are similar to
Nicholl et al. (2017), Kangas et al. (2022) and Chen et al. (2023a), and they
are summarised in Table 26. For the pure nickel model, we set the opacities
$\kappa$ and $\kappa_{\gamma}$ to 0.07 and $0.027\leavevmode\nobreak\ \rm
g\,cm^{-2}$, respectively (Swartz et al. 1995; Wang et al. 2015, and
references therein). For the models that include two sources of energy, in
which one of the energy sources is 56Ni, we perform fits with unconstrained
opacities as well as fits with the opacities set to
$\kappa=0.07\leavevmode\nobreak\ \rm g\,cm^{-2}$ and
$\kappa_{\gamma}=0.03\leavevmode\nobreak\ \rm g\,cm^{-2}$. The model
parameters are inferred in a Bayesian way using the nested sampler dynesty.
The fits of each model are shown in Figure 24. As the fit covers a wide time
interval, each panel in Figure 24 also contains a window zooming in onto the
region of maximum light. The marginalised posteriors of the model parameters
are summarised in Table 26.

Visually, all models capture the rise, peak and decline up to $t_{\rm
max}$+400 days. There are noticeable differences between the fits and the data
because of i) not all models can be correct, ii) the inherent assumptions of
each model, and iii) the assumption of a blackbody photosphere at all times.
Owing to this, none of the models can capture the bumps and undulations (see
inset in Figure 24). The significant deviation in the $z$ band at $>$$t_{\rm
max}$+200 days is due to the assumption of a blackbody photosphere. The late-
time spectra reveal a blue pseudo-continuum with super-imposed emission lines.
The luminescent [O ii] $\lambda\lambda$ 7320,7330 emission lines are
redshifted to the $z$ band and cause the apparent discrepancy between the data
and the models (Figure 6).

Besides these general caveats, differences in the fit qualities between the
models are visible. At epochs later than $t_{\rm max}$+500 days, the pure
central engine models fail to describe the data. The discrepancies grow with
time and reach $\sim 2$ mag per band at $t_{\rm max}$+706 days. In contrast to
that, all nickel and CSM models are able to describe the entire light curve
from $t_{\rm max}$$-93$ days to $t_{\rm max}$+706 days. The fundamental
difference in the late-time behaviour of the central-engine and nickel models
lies in how the models deposit energy into the ejecta (and hence the SN
luminosity). In the pure central engine models, the energy deposition evolves
as a power law in time. This leads to a decrease in the decline rate of the SN
brightness in the time vs. magnitude space. In contrast to that, radioactive
material has an exponentially declining energy deposition rate. This results
in a linear decline in the time vs. magnitude space, which fits our
observations well. The loss of $\gamma$-ray trapping accelerates the fading
independent of the powering mechanism, but it only modifies the light curve
without altering its general shape. In other words, the loss of gamma-ray
trapping cannot convert a power-law decline into an exponential decline (e.g.
Chen et al. 2015; Wang et al. 2015; Nicholl et al. 2018).

Figure 24:  Modelling of the light curves from the rest-frame UV to the NIR
with MOSFiT. All models provide an adequate description of the data up to
$t_{\rm max}$+400 days, though with differences in the fit quality. At later
times, the models diverge. The pure nickel model is the only model that
captures the evolution after $t_{\rm max}$+500 days and has physically
meaningful parameters. The central engine models (magnetar and fallback)
predict a flattening of the light curve due to a power-law-shaped heating rate
in contrast to powering by 56Ni that has an exponential energy deposition
rate. The magnetar+nickel model also captures the full evolution. However, the
inferred model parameters would be either physically implausible or require an
exotic star, which we deem not viable. We note that the Keck photometry at
$t_{\rm max}$+539 days and $t_{\rm max}$+562 days is not corrected for host
contamination owing to the lack of Keck reference images to perform image
subtraction. The expected host contribution is $\approx 10\%$.

The nickel models (56Ni and 56Ni+CSM—magnetar) provide an adequate fit to the
entire light curve. All models require $\approx 35\leavevmode\nobreak\
M_{\odot}$ (spread in the median values: 20–$63\leavevmode\nobreak\
M_{\odot}$; statistical error of the individual measurements: 7–67%; Table 26)
of freshly synthesised nickel, consistent with our conclusions on the
bolometric light curve (Section 5.2.1). The pure 56Ni has a worse fit
statistic because of the loss of $\gamma$-ray trapping. This leads to
underpredicting the expected brightness. This is not a critical issue. The
modelling assumes that the SED is still a blackbody, which is not the case
anymore (Section 5.2.4). Furthermore, fitting the bolometric light curve with
56Ni gave an excellent fit (Section 5.2.1. Such large nickel masses can only
be produced in a PISN explosion. PISN models predict no remnant after the
entire star is obliterated (Fowler & Hoyle 1964; Barkat et al. 1967; Rakavy et
al. 1967), eliminating the magnetar + 56Ni model. If we were to ignore stellar
evolution theory, the rotational energy of the magnetar, which defines how
much energy could be converted into radiation, would contribute $<1\%$ to the
total radiated energy ($\sim 1\times 10^{50}\leavevmode\nobreak\ \rm erg$),
whereas $>99\%$ of the radiated energy would come from the radioactive decay
of 56Ni and its daughter products (Table 26). Moreover, the inferred spin
period of $\sim 15$ ms is much larger than the median spin period of $\sim
2.6$ ms from the ZTF-I SLSN sample (Chen et al. 2023a, see also Nicholl et al.
2017 and Blanchard et al. 2020). Even the slowest spinning SLSN magnetars
never exceeded 6–7 ms (Nicholl et al. 2017; Blanchard et al. 2020; Chen et al.
2023b).252525All measurements are based on the fiducial assumption of dipole
spin-down radiation. For a discussion on the impact of non-dipole radiation
see Omand & Sarin (2023).

Table 8: Light-curve fits with MOSFiT: models, parameters, priors and
marginalised posteriors

Parameter | Prior | Magnetar | Magnetar + | Magnetar + |  56Ni | Fallback | CSM | CSM + | CSM + |  56Ni
---|---|---|---|---|---|---|---|---|---|---
|  |  |  56Ni |  56Ni |  |  |  |  56Ni |  56Ni | (red)
|  |  |  |  (fixed $\kappa$’s) |  (fixed $\kappa$’s) |  |  |  |  (fixed $\kappa$’s) |  (fixed $\kappa$’s)
Fitted properties
General
ejecta mass $M_{\rm ej}$ $\left(M_{\odot}\right)$ | $\log\mathcal{U}\left(1,300\right)$ | $81^{+12}_{-9}$ | $54^{+29}_{-23}$ | $141^{+33}_{-29}$ | $83\pm 4$ | $74\pm 5$ | $46^{+11}_{-13}$ | $34^{+3}_{-2}$ | $34\pm 6$ | $107^{+8}_{-9}$
explosion date $t_{\rm exp}$ (day) | $\mathcal{U}\left(-200,0\right)$ | $-19\pm 3$ | $-22\pm 2$ | $-22^{+2}_{-3}$ | $-70\pm 5$ | $-11\pm 2$ | $-27^{+4}_{-6}$ | $-17^{+2}_{-3}$ | $-14\pm 2$ | $-97^{+7}_{-8}$
‘$\gamma$-ray’ opacity $\kappa_{\gamma}$ $\left(\rm cm^{2}\,g^{-1}\right)$ | $\log\mathcal{U}\left(10^{-2},10^{4}\right)$ | $0.013\pm 0.002$ | $14^{+744}_{-14}$ | 0.03 | 0.03 | $0.010\pm 0.001$ | … | $155^{+1751}_{-143}$ | 0.03 | 0.03
optical opacity $\kappa$ $\left(\rm cm^{2}\,g^{-1}\right)$ | $\mathcal{U}\left(0.01,0.2\right)$ | $0.18\pm 0.02$ | $0.05\pm 0.02$ | 0.07 | 0.07 | $0.19\pm 0.01$ | … | $0.07\pm 0.01$ | 0.07 | 0.07
scaling velocity $v_{\rm scale}$ $\left({\rm km\,s}^{-1}\right)$ | $\mathcal{U}\left(1000,10000\right)$ | $5130^{+210}_{-200}$ | $5660^{+170}_{-160}$ | $5620^{+180}_{-170}$ | $4050\pm 130$ | $5560^{+240}_{-200}$ | $5510^{+200}_{-220}$ | $5950^{+230}_{-200}$ | $5850^{+190}_{-160}$ | $3600^{+260}_{-310}$
white noise parameter $\sigma$ | $\log\mathcal{U}\left(10^{-3},100\right)$ | $0.25\pm 0.01$ | $0.21\pm 0.01$ | $0.25\pm 0.01$ | $0.25\pm 0.01$ | $0.26\pm 0.01$ | $0.21\pm 0.01$ | $0.21\pm 0.01$ | $0.21\pm 0.01$ | $0.23\pm 0.01$
Magnetar model
magnetic field $B_{\perp}$ $\left(10^{14}\,{\rm G}\right)$ | $\log\mathcal{U}\left(0.01,20\right)$ | $0.72^{+0.05}_{-0.08}$ | $1.32^{+1.25}_{-1.24}$ | $1.05^{+0.47}_{-0.36}$ | … | … | … | … | … | …
neutron-star mass $M_{\rm NS}$ $\left(M_{\odot}\right)$ | $\mathcal{U}\left(1,2.2\right)$ | $2.1\pm 0.1$ | $1.4^{+0.5}_{-0.3}$ | $1.6^{+0.4}_{-0.3}$ | … | … | … | … | … | …
initial spin period $P_{0}$ (ms) | $\mathcal{U}\left(1,20\right)$ | $1.0^{+0.1}_{-0.0}$ | $15.5^{+3.0}_{-4.8}$ | $14.6^{+3.1}_{-3.5}$ | … | … | … | … | … | …
56Ni model
nickel fraction $f_{\rm Ni}$ | $\log\mathcal{U}\left(10^{-3},1\right)$ | … | $0.6\pm 0.2$ | $0.2\pm 0.1$ | $0.50\pm 0.02$ | … | … | $0.9\pm 0.1$ | $0.6\pm 0.1$ | $0.6^{+0.3}_{-0.2}$
Fallback model
luminosity $L_{1}$ $\left(10^{55}\rm erg\,s^{-1}\right)$ | $\log\mathcal{U}\left(10^{-4},10^{3}\right)$ | … | … | … | … | $4.5\pm 0.1$ | … | … | … | …
transition time $t_{\rm tr}$ (day) | $\log\mathcal{U}\left(10^{-4},10^{4}\right)$ | … | … | … | … | $0.003^{+0.014}_{-0.003}$ | … | … | … | …
CSM
CSM mass $M_{\rm CSM}$ $\left(M_{\odot}\right)$ | $\log\mathcal{U}\left(0.01,300\right)$ | … | … | … | … | … | $58^{+13}_{-11}$ | $1.2^{+0.8}_{-0.5}$ | $13^{+3}_{-2}$ | …
CSM density $\rho\left(10^{-14}\rm cm^{-3}\right)$ | $\log\mathcal{U}\left(10^{-14},10^{2}\right)$ | … | … | … | … | … | $19^{+9}_{-6}$ | $14^{+57}_{-10}$ | $36300^{+18700}_{-15900}$ | …