power-law index of the CSM | $\mathcal{U}\left(0,2\right)$ | … | … | … | … | … | $0.6\pm 0.3$ | $1.8\pm 0.1$ | $0.5\pm 0.2$ | …
density profile $s$ |  |  |  |  |  |  |  |  |  |
slope of the outer SN ejecta | $\mathcal{U}\left(8,12\right)$ | … | … | … | … | … | $8.8^{+0.6}_{-0.4}$ | $9.7^{+0.8}_{-0.9}$ | $10.0^{+0.6}_{-0.7}$ | …
density profile $n$ |  |  |  |  |  |  |  |  |  |
slope of the inner SN ejecta | fixed | … | … | … | … | … | 0 | 0 | 0 | …
density profile $\delta$ |  |  |  |  |  |  |  |  |  |
progenitor radius $R_{0}$ (AU) | $\log\mathcal{U}\left(0.1,1000\right)$ | … | … | … | … | … | $617^{+160}_{-149}$ | $59^{+67}_{-35}$ | $10^{+6}_{-4}$ | …
Fit quality
log Bayesian evidence ($\log\leavevmode\nobreak\ Z$) |  | 516 | 640 | 639 | 541 | 497 | 644 | 642 | 638 | 287
Number of free parameters |  | 11 | 12 | 10 | 7 | 10 | 11 | 14 | 13 | 7
Derived properties
$\gamma$-ray escape time $t_{0}$ (day) |  | $320^{+30}_{-20}$ | $9200^{+76200}_{-8400}$ | $630\pm 80$ | $570\pm 20$ | $290\pm 20$ | … | $23700^{+64800}_{-17400}$ | $300\pm 30$ | $690^{+60}_{-50}$
nickel mass $M_{\rm Ni}$ $\left(M_{\odot}\right)$ |  | … | $32^{+22}_{-18}$ | $32^{+11}_{-9}$ | $42\pm 3$ | … | … | $31^{+4}_{-3}$ | $20^{+6}_{-5}$ | $63^{+31}_{-19}$
kinetic energy $E_{\rm kin}$ $\left(10^{51}\rm erg\right)$ |  | $21\pm 2$ | $17^{+9}_{-7}$ | $45^{+22}_{-19}$ | $14\pm 1$ | $23\pm 2$ | $14^{+3}_{-4}$ | $12\pm 1$ | $12\pm 2$ | $14^{+1}_{-2}$
rotational energy $E_{\rm rot}$ $\left(10^{51}\rm erg\right)$ |  | $33^{+5}_{-4}$ | $0.09^{+0.05}_{-0.04}$ | $0.11^{+0.03}_{-0.02}$ | … | … | … | … | … | …

262626The model ‘56Ni (red)’ only fitted the data in the $r$ and redder bands.
We used uniform ($\mathcal{U}$) and log uniform ($\log\mathcal{U}$) priors.
The uncertainties of the marginalised posteriors are quoted at $1\sigma$
confidence. The explosion date is measured with respect to the date of the
first detection. All marginalised posteriors are reported in linear units. The
Bayesian evidence is reported in log units. The kinetic energy of the ejecta
was computed via $E_{\rm kin}=1/2\,M_{\rm ej}\,v_{\rm scale}^{2}$ and the
rotational energy of the magnetar via $E_{\rm rot}=2\times
10^{52}\,\left(M_{\rm NS}/1.4\,M_{\odot}\right)^{3/2}\,\left(P_{0}/1\,\rm
ms\right)^{-2}\leavevmode\nobreak\ \rm erg$.

The CSM + 56Ni model might be viable. However, the inferred nickel fraction
and CSM parameters are very sensitive to the a priori unknown effective
opacities $\kappa$ and $\kappa_{\gamma}$. The two parameters are different for
each radiation component and come with their own uncertainties. Furthermore,
the contribution of CSM interaction and 56Ni to the observed light curve
changes with time resulting in a time-variable $\kappa$ and $\kappa_{\gamma}$.
Capturing these complexities is not possible with that particular CSM + 56Ni
model. Furthermore, the Chatzopoulos et al. (2012) CSM model used in MOSFiT is
debated (see the discussion in Sorokina et al. 2016; Moriya et al. 2018b),
limiting the interpretability of the results. More sophisticated CSM + 56Ni
models are needed to disentangle the two sources of energy (e.g. Chevalier &
Fransson 2017; Tolstov et al. 2017a; Suzuki & Maeda 2021; Takei et al. 2022).

The pure magnetar model can be rejected on physical and statistical grounds.
The aforementioned increased over-prediction of observed flux at late time
makes this model unviable. Furthermore, the magnetar model requires extreme
conditions ($M_{\rm NS}\sim 2.1\leavevmode\nobreak\ M_{\odot}$,
$P_{0}=1\leavevmode\nobreak\ \rm ms$ and $M_{\rm ej}\sim
81\leavevmode\nobreak\ M_{\odot}$) to squeeze out as much energy as possible
from the neutron star. The lower limit on the progenitor mass ($M_{\rm
progenitor}>M_{\rm ej}+M_{\rm NS}$) exceeds $83\leavevmode\nobreak\
M_{\odot}$. Explosion models predict that such a massive star leaves behind a
black hole but not a neutron star (Heger & Woosley 2002). From an empirical
point of view, the H-poor SLSNe, which are thought to be powered by a
magnetar, have ejecta masses of only $\sim 5\leavevmode\nobreak\ M_{\odot}$.
The most massive ejecta reach a few $10\leavevmode\nobreak\ M_{\odot}$ but
never exceed $50\leavevmode\nobreak\ M_{\odot}$ (Nicholl et al. 2017;
Blanchard et al. 2020; Tinyanont et al. 2023; Chen et al. 2023b; West et al.
2023), in stark contrast to the $\sim 81\leavevmode\nobreak\ M_{\odot}$
required to power SN 2018ibb. Previously, Eftekhari et al. (2021) reported
that SN 2018ibb is powered by a magnetar (using the same magnetar model) but
using a very sparse data set. We caution against this practice. Even with our
comprehensive data set, only the data after $t_{\rm max}$+600 days broke the
degeneracy between the central-engine models and the nickel-powered models.
Furthermore, the best-fit magnetar properties presented here and in Eftekhari
et al. (2021) are significantly different, demonstrating that datasets with a
large wavelength coverage and a wide time span are required to determine the
powering mechanism of SLSNe. Our results echo the conclusions from Moriya et
al. (2017), who performed a parameter study of magnetar and nickel models,
that the two models could produce indistinguishable light curves if the data
set covers only a limited time interval. Very late-time observations, such as
those presented here, are critical to break the degeneracy in the light curve
modelling (Moriya et al. 2017).

The fallback model has implausible parameters. The total energy input from the
fallback accretion required to power the light curve is given by
$2.5\,L_{1}\times t^{-2/3}_{\rm tr}$ (Moriya et al. 2018a). For the best-fit
values (Table 26), we infer a total energy input of
0.9–$4.4\leavevmode\nobreak\ M_{\odot}\,c^{2}$. Assuming a realistic
conversion efficiency from the fallback accretion disk to the large-scale
outflow of $10^{-3}$ (i.e. the SN luminosity; Dexter & Kasen 2013), the mass
of the accretion disk would have to be 900–$4400\leavevmode\nobreak\
M_{\odot}$, which is unphysical. In the case of a jet forming during the
accretion process, the conversion efficiency could be $\sim 10\%$ (McKinney
2005; Kumar et al. 2008; Gilkis et al. 2016) and the mass of the accretion
disk would only have to be 9–$44\leavevmode\nobreak\ M_{\odot}$. Irrespective
of the accretion efficiency, the fallback model always over-predicts the
brightness by $\sim 2$ mag at late times. To mitigate this issue, the
accretion rate would also have to be fine-tuned to match observations. Owing
to these issues, we do not deem the fallback model to be viable for SN
2018ibb.

To strengthen our conclusions on the model selection, we also fit the
observations using the software package Redback (Sarin et al. 2023), which
implements these different models. We fit the multi-band data in magnitude
space with a Gaussian likelihood function and the exact same priors, and
utilise the nestle272727http://kylebarbary.com/nestle/ sampler implemented in
bilby (Ashton et al. 2019; Romero-Shaw et al. 2020). The fit parameters of the
pure magnetar and nickel models as well as the magnetar + nickel model are
consistent with the results from MOSFiT (Table 61). The pure nickel model is
the only viable model. Furthermore, the necessity for a large ejecta is solely
determined by the long rise. It does not depend on the availability of the
data in the blue bands. We verified that by fitting the data in the $r$ band
and in redder filters with MOSFiT and Redback (Tables 26 and 61).

In summary, the multi-band light curve is consistent with being powered by
radioactive 56Ni. The undulations, visible in several bands, are not captured
by this nickel model. This could point to an additional source of energy, for
example, interaction with circumstellar material. The required 56Ni mass of
$\sim 35\leavevmode\nobreak\ M_{\odot}$ can only be produced during a PISN
explosion. Central engine models can be excluded with high confidence. The
necessity for a large ejecta and nickel mass is solely determined by the long
rise, the high peak luminosity and the slow decline. It does not depend on the
availability of the data in particular filters.

Due to this extraordinary result, we want to briefly comment on whether the
derived ejecta mass ($M_{\rm ej}\sim 83\leavevmode\nobreak\ M_{\odot}$),
nickel fraction ($\sim 50\%$), velocity ($\sim 4000\leavevmode\nobreak\ \rm
km\,s^{-1}$) and kinetic energy ($14\times 10^{51}\leavevmode\nobreak\ \rm
erg$) are sensible for PISN models (values taken from Table 26). The PISN
models by Kozyreva et al. (2017), Gilmer et al. (2017) and Heger & Woosley
(2002) predict for a progenitor with a nickel yield of $\sim
35\leavevmode\nobreak\ M_{\odot}$ an ejecta mass of $\sim
130\leavevmode\nobreak\ M_{\odot}$, a nickel fraction of $\sim 30\%$, an
ejecta velocity of 9000–11000 km s-1 and a kinetic energy of 80–$90\times
10^{51}\leavevmode\nobreak\ \rm erg$. The model-predicted ejecta velocities,
kinetic energy and nickel fractions appear to be in tension with our results,
but that is not critical for the following reasons. i) The nickel model used
in the light-curve fitting is an analytical one-dimensional model, but not a
PISN model. ii) The velocity of the nickel model is a scaling velocity, which
is not the ejecta velocity. In Section 4.3.2, we showed that the ejecta
velocity is $8500\leavevmode\nobreak\ \rm km\,s^{-1}$. The fastest portions of
the ejecta move at $\sim 12,500\leavevmode\nobreak\ \rm km\,s^{-1}$. These
velocities yield a more realistic kinetic energy of 60–$130\times
10^{51}\leavevmode\nobreak\ \rm erg$, consistent with PISN model predictions.
iii) Different light-curve fitting codes give slightly different results. iv)
PISN models have large uncertainties in predicting the energy release and the
amount of ${}^{56}\mathrm{Ni}$. They depend on the nuclear network and
assumptions of a given stellar evolution code, such as convection, details of
the solution of the stellar structure equations and the dynamical phase of the
pair-instability explosion, etc. For example, ${}^{56}\mathrm{Ni}$ masses
calculated by Takahashi (2018) and Umeda & Nomoto (2002) differ significantly
from Heger & Woosley (2002). In Umeda & Nomoto (2002), a
$270\leavevmode\nobreak\ M_{\odot}$ star would produce $10\leavevmode\nobreak\
M_{\odot}$ of ${}^{56}\mathrm{Ni}$ and a similar star in Takahashi (2018)
would produce between a few $M_{\odot}$ masses and $>30\leavevmode\nobreak\
M_{\odot}$ of ${}^{56}\mathrm{Ni}$. In contrast to that, a
$250\leavevmode\nobreak\ M_{\odot}$ star in Heger & Woosley (2002) yields
$\sim 55\leavevmode\nobreak\ M_{\odot}$ of ${}^{56}\mathrm{Ni}$. As shown by
Farmer et al. (2019, 2020) and Kawashimo et al. (2023), the choice of the
${}^{12}\mathrm{C}(\alpha,\gamma)^{16}\mathrm{O}$ reaction rate leads to the
scatter in the final ${}^{56}\mathrm{Ni}$ mass between $10M_{\odot}$ and
$70\leavevmode\nobreak\ M_{\odot}$ for the He-core stellar model of
$130\leavevmode\nobreak\ M_{\odot}$, and a factor of 2 in the amount of the
released explosion energy.

####  5.2.3 Matching the light curve with PISN templates

Motivated by the light-curve fits, we compare the bolometric light curve to
the PISN templates from Kasen et al. (2011), Gilmer et al. (2017) and Kozyreva
et al. (2017). The grid of models from Kasen et al. (2011) are the metal-free
helium models from Heger & Woosley (2002). The Gilmer et al. (2017) and
Kozyreva et al. (2017) models assume a metallicity of 7% solar. Very massive
stars in low-metallicity environments ($Z\sim 0.07\leavevmode\nobreak\
Z_{\odot}$) lose their hydrogen envelopes during the early evolution, assuming
up-to-date wind mass-loss rates. The details about the used mass-loss rates
are described in Ekström et al. (2012) and Yusof et al. (2013). Therefore,
these stars are hydrogen-free by the time of the pair-instability episode, and
the helium-core models from Heger & Woosley (2002) are a good representation
of these explosions. This is in agreement with the models from Gilmer et al.
(2017). Their suite of models, which were computed self-consistently, are
initially hydrogen-rich models; however, owing to mass loss, the highest-mass
models become hydrogen-free by the time of the explosion.

Figure 25: Comparison of SN 2018ibb with the PISN models P250Ni25 and P250Ni34
from Kozyreva et al. (2017) (left panel) and the He100, He125, and He130 from
Heger & Woosley (2002) (right panel). Templates with nickel masses of
34–$44\leavevmode\nobreak\ M_{\odot}$ are required to describe the entire
bolometric light curve from $t_{\rm max}$$-93$ to $t_{\rm max}$+706 days.
Models with $M({\rm Ni})=25\leavevmode\nobreak\ M_{\odot}$ systematically
underestimate the observed bolometric light curve, but they could still be
viable if CSM interaction contributes significantly throughout the evolution.

Among the suitable models, we chose the P250Ni25 and P250Ni34 templates from
Gilmer et al. (2017), and the He100, He120, He125 and He130 templates from
Kasen et al. (2011) (where the number stands for the helium core mass in
$M_{\odot}$). The models have nickel yields between 5.8 and
$40\leavevmode\nobreak\ M_{\odot}$. The vital properties of these models are
presented in Table 9. The P250Ni25 model starts with an initial mass of
$250\leavevmode\nobreak\ M_{\odot}$. At the time of the explosion, a helium
core of $127\leavevmode\nobreak\ M_{\odot}$ has formed, which is similar to
the helium models He125 and He130. We note that the P250Ni25 model not only
loses its hydrogen envelope but also most of its helium layer (total mass of
$2.6\leavevmode\nobreak\ M_{\odot}$ before the loss of the He layer) and ends
up as a bare carbon-oxygen core with a tiny helium fraction by the time of the
pair-instability explosion. In contrast to that, the He125 and He130 models
evolve without mass losses and retain $2.4\leavevmode\nobreak\ M_{\odot}$ and
$2.8\leavevmode\nobreak\ M_{\odot}$ of helium, respectively.

Table 9: Summary of PISN model parameters

Name | $M(\rm ZAMS)$ | $M(\rm He)$ | $M(\rm Ni)$ | Metallicity | $v_{\rm ejecta}$ | $E_{\rm kin}$
---|---|---|---|---|---|---
| $\left(M_{\odot}\right)$ | $\left(M_{\odot}\right)$ | $\left(M_{\odot}\right)$ | $\left(Z/Z_{\odot}\right)$ | $\left(\rm km\,s^{-1}\right)$ | $\left(10^{51}\,\rm erg\right)$
He100 | 205 | 100 | 6 | 0.01 | 8400 | 42
He120 | 242 | 120 | 26 | 0.01 | 10000 | 71
He125 | 251 | 125 | 34 | 0.01 | 10300 | 79
He130 | 260 | 130 | 44 | 0.01 | 10600 | 87
P250Ni25 | 250 | 127 | 25 | 0.07 | 7500 | 82
P250Ni34 | 250 | 127 | 34 | 0.07 | 8850 | 82
282828The mass in column (3) lists the mass of the He core before the
progenitor explodes. The values of the He100–He130 models were taken from
https://2sn.org/DATA/HW01/bulk_yields.txt and are based on Heger & Woosley
(2002). The values of the P250 models are from Kozyreva et al. (2017) and
Gilmer et al. (2017). The ZAMS masses of the He100–He130 models are computed
with $M({\rm He})=13/24\times\left[M({\rm ZAMS})-20\,M_{\odot}\right]$ (Heger
& Woosley 2002).

To build the bolometric light curves of the PISN models, we use the
hydrodynamics radiative transfer code STELLA (Blinnikov et al. 2006). The
slight difference of the P250Ni34 light curve between our calculation and that
in Gilmer et al. (2017) and Kozyreva et al. (2017) is caused by the different
versions of STELLA used in the two studies. A relevant discussion can be found
in Kozyreva et al. (2020). The re-calculated light curves of the helium models
are consistent with those calculated with the spectral synthesis code SEDONA
(Kasen et al. 2011). In Figure 25, we compare the bolometric light curve of SN
2018ibb to those computed for a series of PISN models. The PISN templates with
nickel yields between 34 and $44\leavevmode\nobreak\ M_{\odot}$ (He125, He130,
P250Ni34) provide excellent matches to the rise, the peak, and the fading
parts of the bolometric light curve of SN 2018ibb. While the He125 and
P250Ni34 models describe the rise and peak well, they systematically
underestimate the late-time flux. Such a deviation at the late epochs may not
necessarily refute these two models since the observed bolometric flux of SN
2018ibb may include a time-varying contribution from CSM interaction. As we
show in Section 5.2.4, this contribution is not negligible and could boost the
luminosity by a few 0.1 dex.

Figure 26: Late-time spectra of SN 2018ibb at 287 and 637 days after its
maximum. Overlaid are the computed PISN spectra from Jerkstrand et al. (2016)
scaled to these epochs. The shaded region indicates the uncertainty of the
explosion time. The He130 model provides an adequate description of the
emission redwards of 6000 Å at $t_{\rm max}$+286.7 days, but a worse match for
the second epoch. The observed spectra show a considerable excess at shorter
wavelengths that is not expected from the model spectra. We argue that the
blue excess is due to the interaction of the SN ejecta with circumstellar
material, which is not included in existing PISN models. The He100 model
matches the observation of neither epoch. The vertical bars at the top of each
panel indicate the location of telluric features.  Figure 27: Nebular spectrum
of SN 2018ibb (black) at $t_{\rm max}$+637.3 days and its decomposition into a
CSM interaction (blue) and PISN component (red). This decomposition reveals
that the shape of the spectrum in the blue is similar to the pseudo-continuum
seen in interaction-powered SNe (Type Ia-CSM, Ibn, Icn and IIn). The emission
lines in the blue arise either from CSM interaction or from material in the
CSM shell excited by the SN light. The dotted vertical lines indicate the
location of strong galaxy emission lines.

The templates with $M({\rm Ni})\sim 25\leavevmode\nobreak\ M_{\odot}$ (He120
and P250Ni25) also provide reasonable matches to the data, even though they
exhibit a faster rise and produce peak luminosities that are $\approx 0.3$ dex
fainter. In the case that CSM interaction contributes at all times, the
apparent tension might be alleviated. The PISN model He100, which produces the
smallest amount of Ni in our set, generates a light curve that is incompatible
with observations. The peak bolometric luminosity of such a model is 0.8 dex
lower compared to the observed value, implying that a different energy source
must account for $>84\%$ of the observed peak luminosity.

####  5.2.4 Late-time spectra of SN 2018ibb compared to PISN models

Jerkstrand et al. (2016) computed spectra of the He100 [$M({\rm
Ni})=5.8\leavevmode\nobreak\ M_{\odot}$] and He130 [$M({\rm
Ni})=44\leavevmode\nobreak\ M_{\odot}$] PISN models at 400 and 700 days after
the explosion.292929The model spectra of the P250 templates will be presented
in a forthcoming paper by Kozyreva et al. (in prep.). To compare these spectra
with the observations, we need to constrain the poorly measured explosion date
of SN 2018ibb. The bolometric light curve of the He100 and He130 models peaked
at $\approx 130$ rest-frame days. Assuming that SN 2018ibb’s bolometric light
curve peaked up to 20 days before the peak in the $g$—$r$ band, we can scale
the computed spectra to the epochs of the observed spectra via
$\exp\left(\Delta t/t_{1/2,\rm Co}\right)$, where $\Delta t$ is the phase
difference between the observed and computed spectra and $t_{1/2,\rm Co}$ is
the mean lifetime of 56Co.

Figure 26 shows the observed spectra of SN 2018ibb at $t_{\rm max}$+286.7 (top
row) and $t_{\rm max}$+637.3 days (bottom row) in black. The upper left and
the bottom left panels compare the earlier and the later spectra with the
phase-adjusted He100 model (red) at $t_{\rm max}$+400 days and $t_{\rm
max}$+700 days, respectively. The right column presents the same comparison to
the He130 model spectra. The phase-adjusted spectra have a shaded band to
indicate the impact of the uncertain peak time of the bolometric light curve
on the model flux. We selected the observed spectra at these specific epochs
to minimise the phase correction and cover a wide wavelength range.

The He100 model fails to match the spectra of SN 2018ibb at both $t_{\rm
max}$+286.7 and $t_{\rm max}$+637.3 days. The predicted emission lines are
significantly weaker compared to the data, and the relative strength of the
features does not match the shape of the observed spectra. In addition, the
model spectra also exhibit lines that are significantly narrower compared to
the observed spectra since the He100 model yields a lower ejecta velocity
(Table 9). The He130 model provides a better match. At $t_{\rm max}$+286.7
days, the model spectrum describes the observed spectrum redwards of 6000 Å
well, in terms of the absolute and relative strength of the features as well
as the line widths. The computed spectrum also matches the observed NIR
spectrum, albeit the strongest predicted feature at 1.2 $\mu$m (Fe i and Si i)
is redshifted to a region that is strongly affected by atmospheric absorption
(indicated by the black-shaded region in the upper half of the figure). The
match at $t_{\rm max}$+637.3 days appears to be less plausible compared to
that at $t_{\rm max}$+286.7 days. While the model reproduces the [O ii]+[Ca
ii] at 7300 Å, the observed spectrum shows an elevated continuum level and
stronger [O i] $\lambda\lambda$ 6300, 6364 in emission. The observed spectrum
also shows prominent O i $\lambda$ 7773 in emission that is not generated by
the model. However, in Section 5.1 we showed that O i $\lambda$ 7773 is
produced by the CSM interaction.

Bluewards of 6000 Å, the discrepancy between the observed and computed spectra
is considerable in both epochs. A similar excess in the blue part of the
spectrum was observed in other slow-evolving SLSNe (e.g. Jerkstrand et al.
2017), and it was used as a critical piece of evidence against the PISN
interpretation (e.g. Dessart et al. 2013; Nicholl et al. 2013). Mazzali et al.
(2019) performed detailed modelling of a nebular phase spectrum of the
candidate PISN 2007bi and confirmed these conclusion from Dessart et al.
(2013), Nicholl et al. (2013) and Jerkstrand et al. (2017). In addition, these
authors put forward the idea that if the core of the explosion was mixed, Ca,
Mg and O would become efficient coolants and would produce conspicuous
emission lines at 4570 Å (Mg), 6300,6363 Å (O), and 7291, 7324 Å (Ca).
However, PISN calculations do not support the notion of mixing in the inner
ejecta (Mazzali et al. 2019; Chen et al. 2020).

This blue excess is not a critical piece of evidence against the PISN
interpretation for SN 2018ibb. In Sections 4.3.3, 4.3.4 and 5.1, we showed
that SN 2018ibb is not solely powered by 56Ni. SN 2018ibb’s progenitor had an
eruptive mass-loss episode shortly before the explosion. The interaction of
the SN ejecta with CSM contributes to the observed light curve via discrete
emission lines, and it could even produce a blue pseudo-continuum similar to
that seen in interaction-powered SNe (Silverman et al. 2013; Hosseinzadeh et
al. 2017; Perley et al. 2022)303030The pseudo-continuum in Type Ibn SNe is the
product of the blending of thousands of iron emission lines (e.g. Dessart et
al. 2022).. This raises the questions of whether the blue excess in SN 2018ibb
is similar to that seen in interaction-powered SNe and how large the
contribution of CSM interaction is to the bolometric light curve.

In Figure 27, we further inspect the spectrum of SN 2018ibb at $t_{\rm
max}$+637 days against the phase-adjusted spectrum of the He130 model. We
attempt to decompose the spectrum of SN 2018ibb into two elements, namely a
PISN and an ejecta-CSM interaction component. The CSM component is represented
by a spectrum of the Type Icn SN 2021csp obtained at $\sim 52.7$ days after
the explosion from Perley et al. (2022).313131Using spectra of Type Ia-CSM or
Ibn SNe would give similar results. We decided to use an Icn template because
it shows no emission lines from H and He, which are present in Type Ia-CSM and
Ibn SNe. Its flux scale is scaled so that the sum of the PISN and CSM
components (green) matches the shape of SN 2018ibb’s pseudo-continuum. This
approach is similar to that in Ben-Ami et al. (2014), where these authors used
a spectrum of a Type IIn SN to deduce that the ejecta of the Type Ic SN 2010mb
interacted with a large amount of H-free circumstellar material. Indeed, this
toy model captures the general shape of SN 2018ibb, suggesting that a
considerable fraction of the flux bluewards of 6000 Å is produced by the CSM
interaction. Most of the emission lines in the blue and O i $\lambda$ 7773
feature were not observed in the spectrum of SN 2021csp. However, we have
shown that some of the observed lines in SN 2018ibb, for example, O i
$\lambda$ 7773,9262, [O ii] $\lambda\lambda$ 7320,7330 and [O iii]
$\lambda\lambda$ 4959,5007, are generated by the CSM interaction. Others, such
as [O i] $\lambda\lambda$ 6300,6364 and Mg i] $\lambda$ 4571, are likely
formed in the unshocked SN ejecta or the contact discontinuity (cool-dense
shell) between the SN ejecta and the CSM (Section 5.1).

Assuming SN 2018ibb’s progenitor is similar to the He130 star model, we can
roughly estimate the fractions of the observed bolometric flux that have been
produced by the nickel decay and the CSM interaction. The bolometric
luminosity calculated at $t_{\rm max}$+286.7 days covers the wavelength range
from 3020 to 14,250 Å (rest-frame). The phase-adjusted model spectrum from
Jerkstrand et al. (2016) accounts for 70% of the bolometric flux. In other
words, the nickel-powered light curve would be 0.1 dex fainter than the
observed bolometric light curve. At $t_{\rm max}$+637.3 days, the observed
bolometric luminosity covers the range from 3930 Å to 8500 Å (rest-frame). The
phase-adjusted PISN spectrum accounts for only 21% of the observed bolometric
flux. This means that the Ni-powered light curve would be 0.7 dex fainter. To
illustrate that, we show in Figure 28 the observed bolometric light curve
(solid blue lines) and the fraction of the observed bolometric light curve
that can be attributed to the He130 model (dashed red lines).

Figure 28: Observed late-time bolometric light curve (solid blue lines) and
the fraction of light that could be attributed to 56Ni after accounting for
CSM interaction (dashed, red). The dotted green curves show the 56Ni light
curve after adding the missing IR flux (up to $5\leavevmode\nobreak\ \mu$m).
The IR correction pushes the light curves back to the regime of PISN models
that produce 25–$44\leavevmode\nobreak\ M_{\odot}$ of 56Ni. Even in the case
of a substantial contribution from CSM interaction, a total amount of
25–$44\leavevmode\nobreak\ M_{\odot}$ of 56Ni appears to be essential to power
SN 2018ibb.

However, there is a critical detail that we need to take into account before
drawing a conclusion. The bolometric light curves of the P250Ni25—34 and
He100–He130 models extend to 50,000 Å. Therefore, our observed bolometric
light curve could miss a substantial fraction of the true bolometric flux. The
Jerkstrand et al. (2016) model spectra cover the wavelength range from the far
UV to 25,000 Å, and figure 13 in Jerkstrand et al. (2016) shows the fraction
of light emitted between 25,000 Å and 50,000 Å, allowing us to estimate the
missing IR fractions. At $t_{\rm max}$+286.7 days, the missing IR fraction is
0.14 dex. This fraction increases to $\sim 0.39$ dex at $t_{\rm max}$+637.3
days, due to the shorter wavelength coverage of the observed bolometric light
curve and an increased mid-IR contribution from the PISN model. The dotted
green line in Figure 28 shows the estimated Ni-powered light curve of SN
2018ibb after correcting for the missing IR flux. The IR correction
fortuitously compensates for most of the observed bolometric flux lost to CSM
interaction, corroborating that even with significant CSM interaction a total
mass of 25–44 M⊙ of 56Ni is still needed to power the light curve and spectra.

A progressively increasing contribution from CSM interaction to the bolometric
flux is not a contrived scenario. If the shock is radiative, as is expected
for a high metallicity and dense CSM, then the luminosity from the shock is
$\sim\dot{M}\Delta v_{\rm rel}^{3}/2$, where $\Delta v_{\rm rel}$ is the
relative velocity of the ejecta and the CSM. If the density gradient, $n$, of
the ejecta is steep, the shock velocity is only slowly decreasing, $v_{\rm
s}\propto t^{-1/(n-s)}$ for a CSM with $r^{-s}$ density profile, and the shock
luminosity will only be a slowly decreasing function of time (Chevalier &
Fransson 2017). Because the radioactive input decreases exponentially, it is
expected that the shock contribution will increase relative to the
radioactively powered input.

####  5.2.5 [Co ii] $\lambda$ 1.025 $\mu$m

The NIR spectra of SN 2018ibb after $t_{\rm max}$+300 days reveal an emission
line at 1.025 $\mu$m that we interpret as [Co ii] (a triplet of individual
lines at 1.019, 1.025 and 1.028 $\mu$m, which result from the 9-1, 10-2, and
11-3 transitions as sorted from higher to lower energies, respectively; Figure
29). The line luminosities are $(2.9\pm 0.8)\times 10^{40}\leavevmode\nobreak\
\rm erg\,s^{-1}$ and $(5.4\pm 1.4)\times 10^{40}\leavevmode\nobreak\ \rm
erg\,s^{-1}$ at $t_{\rm max}$+352.6 and $t_{\rm max}$+377.5 days,
respectively. Assuming optically thin LTE, we can convert the line luminosity
to a (temperature-dependent) Co ii mass. The line luminosity of the [Co ii]
$\lambda$ $1.025\leavevmode\nobreak\ \mu$m multiplet can be written as the sum
of the individual transitions

$L({\rm{Co\textsc{ii}}})=N_{9}\,A_{9-1}\,E_{9-1}+N_{10}\,A_{10-2}\,E_{10-2}+N_{11}\,A_{11-3}\,E_{11-3},$

where $N_{u}$ is the total number of ions in the upper state $u$, $A_{u-l}$
the transition rate for spontaneous emission from the upper state $u$ to the
lower state $l$, and $E_{u-l}$ the energy level of the transition. We use
Quinet (1998) for the values of the Einstein coefficients and the energy
levels. The partition function can be taken as 20, with an error less than a
factor of 2 for reasonable temperatures. Then, using equation 42 in Jerkstrand
(2017), the cobalt mass is given by

$M\left({\rm{Co\textsc{ii}}}\right)\gtrsim 0.5\leavevmode\nobreak\
M_{\odot}\,\times\,\left[\frac{L\left(\rm{Co\textsc{ii}}\right)}{5\times
10^{40}\leavevmode\nobreak\ \rm
erg/s}\right]\,\times\,\frac{\exp{\left(15410/T\right)}}{\exp{\left(15410/5000\right)}},$

where $T$ is the temperature of the ejecta, in units of $K$. The temperature
factor (the ratio of exponentials) varies from 0.2 at $T=10,000$ K to 100 at
$T=2000$ K.

To calculate the initial nickel mass, we need to account for the amount that
has decayed over time. The initial nickel mass is a factor of
$\exp\left(t/t_{1/2,\rm Co}\right)\gtrsim\exp\left(450/111\right)\simeq 60$
larger, where $t$ is the time since explosion and $t_{1/2,\rm Co}$ is the mean
lifetime of 56Co. Averaging over the line luminosities of the two epochs, the
inferred 56Ni mass is $\gtrsim 30\leavevmode\nobreak\ M_{\odot}$ if $T\approx
5000$ K. This estimate is consistent with the inferred 56Ni mass from the
light-curve modelling (Sections 5.2.1, 5.2.2, 5.2.3). However, lower values of
the nickel mass would be expected if the temperature is higher
($6\leavevmode\nobreak\ M_{\odot}$ at $T=10,000$ K). For temperatures below
$\sim 3500$ K, the nickel mass becomes unphysically large,
$>100\leavevmode\nobreak\ M_{\odot}$. We note that [Co ii] $\lambda$ 1.025
$\mu$m can be blended with S ii $\lambda$ 1.032 $\mu$m323232This feature
consists of six lines between 1.0287 and 1.0370 $\mu$m., indicated by the
hatched region in Figure 29.

This Ni-mass estimate assumes that the transitions are optically thin. The
Sobolev optical depth of the 9-1 transition line in LTE is (Jerkstrand et al.
2017, ignoring stimulated emission):

$\displaystyle\tau_{9,1}$
$\displaystyle=A_{9,1}\,\lambda_{9,1}^{3}\,\frac{1}{8\pi}\frac{g_{9}}{g_{1}}\,n_{1}t$
$\displaystyle\approx
0.08\times\left(\frac{M\left({\rm{Co\textsc{ii}}}\right)[450d]}{1\leavevmode\nobreak\
M_{\odot}}\right)\frac{x_{1}}{f},$

where $\lambda_{9,1}$ is the wavelength of the emitted photon, $g_{n}$ is the
multiplicity of the $n$th state, $n_{1}$ is the number density of atoms in the
ground state, $x_{1}$ is the fraction of Co II ions in the ground state, and
$f$ is the filling factor for the 56Ni zone. In LTE at 5000 K $x_{1}$ is
$\approx 0.5$, whereas at lower temperatures and/or in NLTE $x_{1}$ is
typically higher (towards unity). A typical CCSN has a characteristic filling
factor of $f\sim 0.1$ for any given zone, which means that 0.5 $M_{\odot}$ of
56Co are optically thin at $\sim$ 450 days. For SLSNe, filling factors for the
oxygen zones have been derived to be $f\approx 10^{-3}-10^{-2}$ (Jerkstrand et
al. 2017). If these filling factors also hold for the 56Ni zone of SN 2018ibb,
then the Co II lines would be optically thick at $\sim$ 450 days, and
determining a mass from them would be impossible at that time (Jerkstrand
2017). One may note that numerical simulations of PISNe show little clumping
or mixing of the inner material (Chen et al. 2020), and the low filling
factors derived for other SLSNe may be due to mixing from the central engine
(e.g. Suzuki & Maeda 2021) or compression by circumstellar interaction (van
Marle et al. 2010).

Figure 29: Zoom-in of the region from 9000 Å to 10,500 Å at $t_{\rm
max}$+377.5 days. Cobalt has its strongest feature at 1.025 $\mu$m. Its
tentative detection translates to 56Ni mass of $\gtrsim 30\leavevmode\nobreak\
M_{\odot}$, consistent with the light curve modelling. The second strongest
cobalt feature is at 9340 Å. Its location is indicated by a fiducial Gaussian
centred 9340 Å. Its integrated luminosity is expected to be between 50 and
100% of [Co ii] $\lambda$ 1.025 $\mu$m. The absence of [Co ii]$\lambda$ 9340
is not an argument against the discovery of [Co ii] $\lambda$ 1.025 $\mu$m
(see Section 5.2.5 for details). Lines from other elements that could blend
with the [Co ii] lines are marked. Regions of strong atmospheric absorption
are indicated by vertical bars at the top of the figure.

In the stripped-envelope-supernova models from Jerkstrand et al. (2015), [Co
ii] $\lambda$ 1.025 $\mu$m is the strongest predicted line from cobalt. The
second strongest Co feature is a blend of two lines at 9338 and 9344 Å. This
[Co ii] feature could be blended with the red wing of the O i $\lambda$ 9263
recombination line. In optically thin LTE, the expected line ratio between [Co
ii] $\lambda$ 9340 and [Co ii] $\lambda$ 1.025 $\mu$m is between $0.5-1$ for a
wide range of plausible temperatures. To examine whether [Co ii] $\lambda$
9340 could be present in the spectrum at $t_{\rm max}$+387 days, we show in
Figure 29 a Gaussian centred at 9340 Å that has either the same integrated
luminosity as [Co ii] $\lambda$ 1.025 $\mu$m or a luminosity that is 50%
smaller. Clearly, [Co ii] $\lambda$ 9340 is not present in our data at those
luminosities.

In the He130 model of Jerkstrand et al. (2016) at 400 days after the
explosion, neither of the [Co ii] lines are present in any significant
strength as they are absorbed by line blocking extending into the NIR. Under
conditions with less line blocking (as in the Jerkstrand et al. 2015 CCSN
models), the [Co ii] $\lambda$ 1.025 $\mu$m line can still be visible, also as
iron has few strong emission lines around this particular wavelength. The same
cannot be said about the 9340 Å region where iron is stronger. In a PISN
ejecta, the densities are about 100-times higher for a given epoch, and the
NIR region is still largely opaque at 400 days after the explosion. To explain
the observed [Co ii] $\lambda$ 1.025 $\mu$m line but the absence of the [Co
ii] $\lambda$ 9340 feature, we need to call upon absorption of the 9340 Å line
but not the 1.025-$\mu$m line. In other words, the He130 model reproduces the
spectral shape near 9340 Å but not at 1.025 $\mu$m if SN 2018ibb is a PISN.

While the association of [Co ii] $\lambda$ 1.025 $\mu$m could be the smoking
gun signature that SN 2018ibb is a PISN, we caution that the interpretation
hinges on the detection of a single line. An IR spectrum with NIRSpec
(Jakobsen et al. 2022) aboard the James Webb Space Telescope could resolve
such ambiguity. It is the only instrument that can provide an uncensored view
from 1 to $5\leavevmode\nobreak\ \mu$m. Such a spectrum could reveal, for
instance, Co, Fe and Ni lines at $>2.7\leavevmode\nobreak\ \mu$m as seen in
the Type Ia SN 2021aefx (Kwok et al. 2023).

###  5.3 Velocity evolution

In Section 4.3.2, we showed that the photospheric velocity is $\approx
8500\leavevmode\nobreak\ \leavevmode\nobreak\ \rm km\,s^{-1}$ and remained
constant between $t_{\rm max}$ and $t_{\rm max}$+100 days. Here, we contrast
our observations with predictions of PISN models.

From the observational point of view, the photospheric velocity is defined via
the minimum of the Fe ii $\lambda\lambda$ 4924,5018,5169 P-Cygni profiles,
which serve as a good proxy of the photospheric velocity (Dessart & Hillier
2005). In contrast to that, there is no specific location of the photosphere
from the radiative transfer point of view, because its location is wavelength
dependent. An average photosphere could be estimated, for instance, via the
Rosseland mean optical depth (Rybicki & Lightman 1986) and electron
scattering.

The photospheric velocity of the light curves simulated with STELLA is defined
as the velocity of the mass shell where the integrated optical depth in the
$B$ band is equal to 2/3. The integrated optical depth accounts for the total
opacity, which includes contributions from continuum opacity (photoionisation,
free-free absorption, and electron scattering processes assuming local
thermodynamical equilibrium in the plasma) and line interactions. The location
of the photosphere defined with this method coincides to a large extent with
the location of the electron-scattering photosphere (Kozyreva et al. 2022) and
permits a direct comparison with our observations.

The velocity evolution of the P250Ni25, P250Ni34, He100, and He130 PISN models
together with the observed Fe ii velocities are shown in Figure 30. The match
between the observed and predicted velocity tracks of the P250Ni25—34 and
He130 models is very good in both the predicted velocities and the velocity
evolution, in particular the He130 model. Small deviations exist in the
absolute values of the velocities (He130: $|v_{\rm obs}-v_{\rm
predicted}|\lesssim 1000\leavevmode\nobreak\ \leavevmode\nobreak\ \rm
km\,s^{-1}$ and P250Ni25—34: $|v_{\rm obs}-v_{\rm
predicted}|\sim+2000\leavevmode\nobreak\ \leavevmode\nobreak\ \rm km\,s^{-1}$)
but not the velocity evolution. These discrepancies are in part due to
differences in the definitions of the synthetic velocity and the physical
model used in STELLA. Furthermore, none of the parameters of the PISN models
were tuned to fit the observations. The velocity track of the He100 model is
inconsistent with the observations ($|v_{\rm obs}-v_{\rm
predicted}|\gtrsim+3000\leavevmode\nobreak\ \leavevmode\nobreak\ \rm
km\,s^{-1}$). As shown in Section 5.2.3, the predicted bolometric light curve
of this model also does not match the observations.

Figure 30:  Photospheric velocity evolution of SN 2018ibb (probed with Fe ii)
and the predicted evolution from different PISN models. The synthetic
velocities of the P250Ni25—34 and He130 models are very similar to the
observed velocities and velocity evolution. These are also the models that
match the observed bolometric light curve well, in particular the He130 model.

###  5.4 Comparison to other slow-evolving SLSNe

Among the $\gtrsim 200$ H-poor SLSNe known, only eight objects belong to such
a phenomenological subclass of slow-evolving SLSNe: SN 1999as (Hatano et al.
2001), SN 2007bi (Gal-Yam et al. 2009), PS1-11ap (McCrum et al. 2014),
PTF12dam (Nicholl et al. 2013), LSQ14an (Inserra et al. 2016), PS1-14bj
(Lunnan et al. 2016), SN 2015bn (Nicholl et al. 2016b) and SN
2018ibb.333333The ZTF-I SLSN sample contains a further possible slow-evolving
SLSN. SN 2018lzx has a rise ($t_{1/e,\rm rise}=60.5^{+8.2}_{-6.8}$ days) and a
decline ($t_{1/e,\rm decline}=108.8^{+10.0}_{-13.2}$ days) time scale
comparable to SN 2018ibb (Table 5). The peak absolute magnitude is 0.2 mag
brighter than that of SN 2018ibb (Chen et al. 2023a). Owing to its high
redshift of $z=0.44$, the light curve spans a short time interval before SN
2018lxz faded below the detection threshold, and the quality of the spectra is
significantly lower compared to that of SN 2018ibb. Therefore, we exclude this
SLSN from the comparison. In the following sections, we compare the
photometric and spectroscopic properties of SN 2018ibb to those of the
historical slow-evolving SLSNe to comprehensively examine its exceptional
properties. We omit SN 1999as from this analysis because its light curve and
spectra were never published.

We utilise the multi-band and bolometric light curves and host-subtracted
spectra of LSQ14an presented in Inserra et al. (2017) and Jerkstrand et al.
(2017), PS1-11ap from McCrum et al. (2014), PS1-14bj from Lunnan et al.
(2016), PTF12dam from Nicholl et al. (2013), Chen et al. (2015) and Quimby et
al. (2018), SN 2007bi from Gal-Yam et al. (2009), Young et al. (2010) and
Jerkstrand et al. (2017), and SN 2015bn from Nicholl et al. (2016b, a, 2018)
and Jerkstrand et al. (2017). Furthermore, we use the Fe ii velocity
measurements from Liu et al. (2017) and Lunnan et al. (2016). All light curves
and spectra were corrected for MW extinction. The spectrum of SN 2015bn in
Nicholl et al. (2018) is not corrected for any host contribution. In Appendix
F, we describe our approach to subtract the host contamination for SN 2015bn.

Table 10: Light-curve properties of slow-evolving SLSNe

| Redshift | $t_{1/e,\rm rise}$ | $t_{1/e,\rm decline}$ | $M_{g,\rm peak}$ | $(g-r)_{\rm peak}$
---|---|---|---|---|---
|  | (day) | (day) | (mag) | (mag)
SN 2018ibb | 0.166 | 68 | 102 | $-21.8$ | $-0.12$
LSQ14an1 | 0.163 | … | $\sim 100$ | $<-20.8$ | $-0.21$
PS1-11ap2 | 0.524 | $<25$ | 38 | $-21.8$ | …
PS1-14bj2 | 0.521 | $83$ | 130 | $-20.6$ | …
PTF12dam | 0.107 | $50$ | 56 | $<-21.7$ | $-0.20$
SN 2007bi3 | 0.128 | $<23$ | $<77$ | $-21.3$ | 0
SN 2015bn | 0.114 | $<31$ | 56 | $-22.0$ | $-0.17$
ZTF SLSNe |  | $29$ | 43 | $-21.5$ | $-0.12$
343434The peak absolute magnitudes and the colours at $t_{\rm peak}$ are
k-corrected and corrected for MW extinction. The uncertainty on the time
scales is of the order of a few days. The peak magnitudes and colours have a
statistical error of $\approx 0.1$ mag. The row ‘ZTF SLSNe’ summarises the
median values of the homogeneous ZTF SLSN sample (Chen et al. 2023a). Chen et
al. (2023a) do not report the median rise and decline $1/e$ time-scale. We
inferred that value from their median $1/10$ rise time-scale and the
relationship between the rise and decline time scales reported in that paper.
$1$$1$footnotetext: The light curve of LSQ14an covers only the post-peak
evolution.
$2$$2$footnotetext: The redshifts of PS1-11ap and PS1-14bj are so high that
the rest-frame $g-r$ colour cannot be inferred from spectra published in
Lunnan et al. (2014) and McCrum et al. (2014).
$3$$3$footnotetext: The light curve of SN 2007bi is only well-observed in $R$
band. We use the $R$-band data as an upper limit on the $g$-band time scale.

####  5.4.1 Light curves

Following the methodology of Chen et al. (2023a), we measure for each slow-
evolving SLSN the k-corrected peak absolute magnitude in the $g$ band, the
k-corrected rest-frame $g-r$ colour, and the $1/e$ rise and decline time-
scales of the $g$-band light curves.353535Owing to the high redshift of
PS1-11ap and PS1-14bj, we use their $i$-band light curves, which probe a rest-
frame wavelength interval similar to that of the $g$ band of SN 2018ibb. All
measurements are summarised in Table 10. We also report in that table the
measurements of SN 2018ibb and, for a broader comparison, the median values of
the homogeneous ZTF SLSN sample (Chen et al. 2023a).

Slow-evolving SLSNe, including SN 2018ibb, have peak absolute magnitudes
between $\sim-20.8$ and $-22$ mag in the $g$ band and k-corrected $g-r$
colours between $-0.2$ and 0 at peak (Table 10). Both their absolute peak
magnitudes and the peak $g-r$-colours are comparable to the median values of
the ZTF SLSN sample (median values being $M_{g,\rm peak}-21.5$ mag and
$g-r=-0.12$ mag; Table 10). The rising parts of the light curves of the
historical slow-evolving SLSNe are not well sampled, limiting the comparison
with SN 2018ibb and the ZTF SLSN-I sample. Only PS1-14bj has a rise time that
is at least as long as that of SN 2018ibb and even 30 days longer than that of
SN 2018ibb. The decline time scales of the historical slow-evolving SLSNe are
well measured. They vary between 38 and 130 days, placing those events above
the average of the ZTF-I sample (Table 10). Yet, only one historical slow-
evolving SLSN had a decline time-scale as extreme as SN 2018ibb. With a
decline time scale of 130 days, PS1-14bj evolves even slower than SN 2018ibb
but its peak luminosity in the rest-frame $g$-band was 1.2 mag fainter than
that of SN 2018ibb. This makes SN 2018ibb an unprecedented case even among the
most extreme SLSNe known. LSQ14an also has a decline time scale of 100 days,
but the observed light curve only covers the declining light curve, adding an
unknown systematic error to its time scale measurement.

Figure 31: SN 2018ibb in the context of the phenomenological sub-class of
slow-evolving SLSNe. Even among this rare sub-class of SLSNe, SN 2018ibb with
its exceptionally broad light curve and high peak luminosity is an extreme
object with unprecedented properties.

Figure 31 shows the $r$-band absolute magnitude light curves of all slow-
evolving SLSNe. The supernovae 2015bn and 2018ibb are the only SLSNe with
observations extending beyond 500 rest-frame days after maximum light. The
light curve of SN 2015bn faded much faster than that of SN 2018ibb. At about
400 days after peak, the decline slowed down and became very gradual. In
contrast, SN 2018ibb’s light curve faded linearly with a decline slope of
$\sim 1.1\leavevmode\nobreak\ \rm mag\,(100\leavevmode\nobreak\ days)^{-1}$
that steepened to $\sim 1.5\leavevmode\nobreak\ \rm
mag\,(100\leavevmode\nobreak\ days)^{-1}$ at 500 days after maximum. These
differences translate into differences in the powering mechanisms. Magnetars
lose their rotational energy efficiently through dipole radiation, which
scales as $\dot{E}_{\rm rot}\propto t^{-2}$. The energy deposition (and
henceforth the SN luminosity) evolves as a power law. Therefore, the light
curve is expected to flatten at later times (in the time vs. magnitude space).
Radioactive material has an exponentially declining energy deposition rate,
which results in a linear decline in the time vs. magnitude space. The loss of
$\gamma$-ray trapping accelerates the fading independently of the powering
mechanism, but it only modifies the light curve without altering its general
shape, i.e. 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). Therefore, the power-law-shaped decline of SN 2015bn could point
to magnetar powering, as concluded in Nicholl et al. (2018). In return, SN
2018ibb’s continued linear decline excludes powering by a magnetar.

Figure 32: The bolometric light curves of SN 2018ibb and historical slow-
evolving SLSNe in the context of PISN models with nickel masses between 5 and
$44\leavevmode\nobreak\ M_{\odot}$. SN 2018ibb is the only SLSN whose entire
light curve from $t_{\rm max}$$-93$ to $t_{\rm max}$+706 days is consistent
with PISN templates. The other SLSNe have either too fast declining light
curves, light curve shapes inconsistent with PISN models, or their light
curves are poorly sampled, hindering a comparison with PISN templates. The
grey-shaded region indicates the $1\sigma$ uncertainty.

Nicholl et al. (2017) fitted the multi-band data of the slow-evolving SLSNe
with the slsn magnetar model in MOSFiT. This model provides an adequate
description of the observations, even to the data of SN 2015bn at 1000 rest-
frame days after maximum (Nicholl et al. 2018). The best-fit parameters cover
the range from 5.3 to $14\leavevmode\nobreak\ M_{\odot}$ for the ejecta mass
$M_{\rm ej}$ (median being $6.3\leavevmode\nobreak\ M_{\odot}$), 0.1 to
$0.8\times 10^{14}$ G for the orthogonal component of the magnetic field
strength $B$ (median being $0.3\times 10^{14}$ G), and 2.3 to 3.9 ms for the
initial spin period $P_{0}$ (median being 2.8 ms). These values are typical
for SLSN light curves fitted with that particular magnetar model (Nicholl et
al. 2017; Chen et al. 2023b). SN 2018ibb has starkly different values (Section
5.2.2; Table 26). The best fit requires a magnetar with an initial spin period
of 1 ms and an ejecta mass of $86\leavevmode\nobreak\ M_{\odot}$, to squeeze
out as much energy as possible from the magnetar model. As alluded to in
Section 5.2.2, such massive stars do not have neutron star remnants.
Furthermore, the magnetar model overpredicts the late-time flux significantly
due to its power-law-shaped energy deposition.

In Figure 32, we compare the bolometric light curves of the slow-evolving
SLSNe to the suite of PISN models used for SN 2018ibb in Section 5.2.3. The
bolometric light curves of all historical slow-evolving SLSNe are either
inconsistent with the PISN models or the comparison is inconclusive: PTF12dam
evolves too fast, the light curves of PS1-11ap and SN 2015bn have a different
shape to PISN templates, PS1-14bj shows a flattening at late times, and the
bolometric light curves of SN 2007bi and LSQ14an have no pre-max data. The
lack of an estimate of the rising bolometric light curve for the latter two
objects precludes concluding whether the two SLSNe could be PISNe or not.
Dedicated studies on PS1-11ap, PS1-14bj, PTF12dam and SN 2015bn revealed that
the magnetar model provides an adequate description of the light curves
(Nicholl et al. 2013; McCrum et al. 2014; Lunnan et al. 2016; Nicholl et al.
2018; Vurm & Metzger 2021).

In conclusion, SN 2018ibb is the _only_ SLSN among the hundreds of SLSNe known
whose entire light curve is consistent with PISN models. This result is even
more revelatory considering that the bolometric light curve covers an
exceptionally wide time interval from $t_{\rm max}$$-93$ to $t_{\rm max}$+706
days.

####  5.4.2 Spectra

In this section, we compare the spectroscopic properties of SN 2018ibb to
those of other slow-evolving SLSNe. First, we compare the photospheric
velocities measured with the Fe ii $\lambda$ 5169 region. The top panel of
Figure 33 displays the photospheric velocities, measured from the Fe ii
$\lambda$5169 region, of SN 2018ibb and other slow-evolving SLSNe363636LSQ14an
is omitted from this comparison. Conspicuous galaxy lines from [O iii]
$\lambda\lambda$ 4959,5007 contaminate the Fe ii $\lambda$5169 region. (in
colour) and of the ZTF-I sample (kernel density estimate). The slow-evolving
SLSNe have velocities between 8000 and $12,000\leavevmode\nobreak\ \rm
km\,s^{-1}$ at peak, lower than the median value of the ZTF-I sample
($14,800\leavevmode\nobreak\ \rm km\,s^{-1}$). PTF12dam and SN 2015bn have the
fastest expanding ejecta ($\sim 12,000\leavevmode\nobreak\ \rm km\,s^{-1}$ at
peak), but their ejecta rapidly decelerate to $\sim 6000\leavevmode\nobreak\
\rm km\,s^{-1}$ in $\sim 60$ rest-frame days. Their velocities and velocity
evolution is similar to those of other SLSNe (Liu et al. 2017). In stark
contrast to that, SN 2018ibb has a velocity of merely
$8500\leavevmode\nobreak\ \rm km\,s^{-1}$, comparable to those of PS1-14bj and
SN 2007bi. Furthermore, the velocity of SN 2018ibb remains constant for 100
rest-frame days, which has not seen for any other SLSN before. Though the
velocities of PS1-14bj and SN 2007bi are very similar to SN 2018ibb, the
spectroscopic sequences of these two events are limited, precluding comparing
their velocity evolution to that of SN 2018ibb.

Next, we explore the spectroscopic properties of slow-evolving SLSNe during
their photospheric and nebular phases. Panel A in Figure 34 shows the
photospheric spectra at the time of maximum light. PTF12dam and SN 2015bn
sustained a hot photosphere with a temperature of $\gtrsim 12,000$ K (Nicholl
et al. 2016b; Vreeswijk et al. 2014). One of the strong features in their
spectra is a comb of O ii absorption lines, a characteristic feature of SLSNe,
which are only seen in photospheres with $T>15,000$ K (Quimby et al. 2018) and
probably also require non-thermal excitation (Mazzali et al. 2019). The
spectra of PS1-11ap, PS1-14bj and SN 2018ibb are cooler (blackbody
temperatures of 10,000 to 12,000 K). Their spectra do not show O ii absorption
lines but instead absorption lines from Ca, Fe, Mg, O and Si (see Figure 7 for
the locations). Common to PTF12dam, SN 2015bn and SN 2018ibb is the presence
of [Ca ii] $\lambda\lambda$ 7291, 7323 in emission. It is one of the strongest
features seen in nebular phase spectra of SNe (Filippenko 1997; Gal-Yam 2017)
but is only seen during the photospheric phase in slow-evolving SLSNe (Gal-Yam
et al. 2009; Inserra et al. 2017; Nicholl et al. 2019). It is also seen in SN
2007bi and LSQ14an, but these SLSNe lack spectra at peak.

Figure 33: Fe ii ejecta velocities of slow-evolving SLSNe (in colour) and
general SLSNe samples (grey) at the time of maximum (top panel) and as a
function of time (bottom panel). SN 2018ibb has a markedly low velocity at the
time of maximum and a flat velocity evolution, which is in stark contrast to
the bulk of the SLSN population. Its velocity at peak is similar to the slow-
evolving SLSNe PS1-14bj and SN 2007bi. However, both comparison objects lack
spectra at earlier and later times.

Around the time of maximum light (Panel B), all objects have similar spectra.
Due to the differences in the ejecta velocities, features appear sharper in
LSQ14an, PS1-14bj and SN 2018ibb that in PTF12dam and SN 2015bn. Some clear
differences are well visible though. LSQ14an, PS1-14bj and SN 2018ibb reveal
[O iii] in emission. As we concluded in Section 5.1, the 7300 Å feature in SN
2018ibb is not dominated by [Ca ii], but [O ii] $\lambda\lambda$ 7320,7330. In
the other objects, the centre of the 7300 Å feature is consistent with [Ca
ii]. Moreover, the 7300 Å feature is well developed in SN 2007bi, LSQ14an and
SN 2018ibb but still very weak in PTF12dam and SN 2015bn. The line profiles
also differ. In SN 2018ibb the line profile is flat-topped but triangular and
skewed to the blue for the other objects.

During the early nebular phase ($t/t_{\rm decl.}\sim 2$; Panel C), all objects
show a blue pseudo-continuum with superimposed forbidden and allowed emission
lines from calcium, magnesium and oxygen (for the line identifications see
Figure 7). SN 2018ibb and PS1-14bj are spectroscopically indistinguishable,
though their overlap in wavelength is limited and PS1-14bj is significantly
fainter than SN 2018ibb (Figure 31). The other objects reveal an increasing
level of dissimilarities (LSQ14an $\rightarrow$ SN 2015bn $\rightarrow$ SN
2007bi). LSQ14an has a similar blue pseudo-continuum but its emission lines
are not well developed. This is best seen in Ca ii $\lambda\lambda$ 3933,
3968, [O iii] $\lambda$ 4363 and [Ca ii] $\lambda\lambda$ 7291,7324 + [O ii]
$\lambda\lambda$ 7320,7330. The SNe 2007bi and 2015bn have redder pseudo-
continua and significantly weaker [Ca ii]+[O ii]. Moreover, SN 2007bi has only
a few features bluewards of 5000 Å.

Figure 34: Comparison of the spectra of SN 2018ibb to those of other slow-
evolving SLSNe between $t_{\rm max}$+30 days and $t_{\rm max}$+1000 days
(darker colour: 5 Å binning; light shade: unbinned spectra). Photospheric
phase (Panels A–B): Around $t_{\rm max}$ (Panel A), the spectra of PTF12dam
and SN 2015bn are characterised by a hot continuum with superimposed O ii
absorption lines as seen in many SLSNe at a similar epoch. SN 2007bi, LSQ14an
and SN 2018ibb have cooler photospheres, and their spectra exhibit absorption
lines from Ca, Fe, Mg, O and Si (see Figure 7 for their locations) but not O
ii. At around $t_{\rm max}$+60 days (Panel B), all spectra appear similar,
though differences exist. LSQ14an, PS1-14bj, and SN 2018ibb are the only SLSNe
showing [O iii] in emission. Furthermore, SN 2007bi, LSQ14an and SN 2018ibb
exhibit strong [Ca ii] + [O ii] in emission. This feature is also present in
PTF12dam and SN 2015bn but is less pronounced. Nebular phase (Panels C–E):
Differences start to emerge during the early nebular phase and become stronger
with time. SN 2018ibb, LSQ14an and PS1-14bj continue to show conspicuous [O
iii] in emission, in contrast to PTF12dam, SN 2007bi and SN 2015bn that have
very strong [O i] and O i in emission. SNe 2015bn and 2018ibb are the only
SLSNe with spectra at $\sim$$t_{\rm max}$+1000 days (Panel E). SN 2018ibb
continues to show intermediate-width [O iii], whereas the spectrum of SN
2015bn exhibits [O i]. The elevated noise in the SN 2018ibb spectrum at
$t_{\rm max}$+989.2 days at $\lambda>6000\leavevmode\nobreak\ \rm\AA$ is due
to residuals of the skyline subtraction. The dashed vertical lines indicate
the expected locations of emission lines commonly seen from H ii regions.

These differences develop further with time. During the late nebular phase
($t/t_{\rm decl.}\sim 5$; Panel D), the pseudo-continuum of all objects fades.
SN 2018ibb and LSQ14an are characterised by a weaker [O i] $\lambda\lambda$
6300,6364 than SNe 2007bi and 2015bn. The ratio between [Ca ii]+[O ii] and [O
i] is 2–3:1. Intriguingly, the emission lines of SN 2018ibb evolved much
slower than for LSQ14an. Now, LSQ14an exhibits more conspicuous emission lines
than SN 2018ibb, best seen in [O iii] and [Ca ii]+[O ii]. The [O ii] feature
of SN 2018ibb has a Lorentzian profile, whereas the profile of LSQ14an is
double-peaked. In contrast to LSQ14an and SN 2018ibb, PTF12dam, SN 2007bi and
SN 2015bn have exceptionally strong [O i]. It is, in fact, their strongest
feature. Moreover, the [O i] is markedly narrower than for SN 2018ibb and
LSQ14an: 6000–9000 km s-1 vs. 16,000 km s-1. The [Ca ii]+[O ii] to [O i] ratio
is 1:2–3 and inverted compared to LSQ14an and SN 2018ibb.

Panel E shows spectra of SN 2018ibb and SN 2015bn at 1000–1100 rest-frame days
after maximum ($t/t_{\rm decline}=9$–13). These are the only two SLSNe with
such extensive spectroscopic observations. Despite the low signal-to-noise
ratio, their spectra exhibit well-defined SN features. SN 2018ibb continues to
show intermediate-width [O iii] with a similar width as in the spectral epochs
before, whereas SN 2015bn exhibits [O i] like in the previous epochs.

Figure 35 presents NIR spectra of LSQ14an, SN 2015bn and SN 2018ibb at
3–4-times their respective decline time scales. All spectra reveal only very
few features beyond 1 $\mu$m, which is expected for models of PISNe
(Jerkstrand et al. 2016), SLSNe (Jerkstrand et al. 2017), and regular
stripped-envelope supernovae (Jerkstrand et al. 2015). Some of the brightest
expected features are redshifted to regions of strong atmospheric absorption
at the average redshift of SLSNe. A feature that has been commonly seen among
all known SLSNe is O i $\lambda$ 1.13$\mu$m. SN 2018ibb reveals an emission
feature at 1.025 $\mu$m, which we identified as [Co ii] (Section 5.2.5). [Co
ii] is not present in any of the other spectra. The data quality of the
spectra of LSQ14an and SN 2015bn is higher compared to that of SN 2018ibb,
suggesting that if a substantial amount of 56Ni was also formed in these
supernovae, the [Co ii] line should have been visible. Instead, SN 2015bn
reveals Mg i $\lambda$ 1.50 $\mu$m that is not visible in SN 2018ibb but
possibly in LSQ14an (Jerkstrand et al. 2017).

In conclusion, SN 2018ibb is spectroscopically similar to other SLSNe,
including slow-evolving SLSNe. During the photospheric phase, SN 2018ibb
stands out by its low ejecta velocity and flat velocity evolution. The early
nebular phase does not differ from other SLSNe. Very late-time observations
($t/t_{\rm decl.}>5$) show clear differences between SN 2018ibb and other
SLSNe, for example, the weak and broad [O i] that stays optically thick
throughout the entire evolution. Late-time NIR spectroscopy revealed the
tentative detection of [Co ii] in SN 2018ibb. This feature is unprecedented
for a SLSN and could be the smoking gun that SN 2018ibb is powered by the
decay of 56Ni. In Sections 4.3.3, 5.1, 4.3.4 and 5.2.4, we argued that the
blue pseudo-continuum in SN 2018ibb is produced by the interaction of the SN
ejecta with CSM. The prevalence of this feature in the other slow-evolving
SLSNe raises the question of whether CSM interaction is also present in these
objects. If this is the case, it is necessary to treat nebular spectra of
SLSNe as the sum of at least two powering mechanisms, for example, magnetar +
CSM or 56Ni + CSM, necessitating more complex SLSN models than the ones that
currently exist. This also means that distinguishing between different
powering mechanisms is more difficult and requires comprehensive data sets.

Figure 35: Late-time NIR spectra of SN 2018ibb and the slow-evolving SLSNe
LSQ14an and SN 2015bn. The strongest features are labelled. SN 2018ibb is the
only SLSN that shows cobalt in emission, which has its strongest optical-NIR
feature at 1.025 $\mu$m. Its luminosity translates to a nickel mass of
$\gtrsim 30\leavevmode\nobreak\ M_{\odot}$, consistent with the light curve
modelling. All spectra are scaled so that O i $\lambda$ 7773 has the same
amplitude in all objects. Regions of strong atmospheric absorption are
cropped.

###  5.5 Is SN 2018ibb a pair-instability supernova?

Models of H-poor PISNe make very clear predictions for PISNe in the regime of
SLSNe ($M_{\rm peak}\lesssim-20$ mag) for their light curves (Kasen et al.
2011; Dessart et al. 2013; Gilmer et al. 2017; Kozyreva et al. 2017), ejecta
velocities, spectra (Dessart et al. 2013; Jerkstrand et al. 2016) and the
environments in which their progenitors are formed (Langer et al. 2007). In
Sections 4.6 and 5.2, we tested the most critical predictions of the PISN
models on the light curves, spectra and host galaxy. SN 2018ibb passes most
tests of PISN models with a nickel yield of 25–44 $M_{\odot}$. However, SN
2018ibb did not comply with two predictions, although it could pass these
tests with the interpretations that we propose. Table 11 summarises all tests.

The tentative detection of [Co ii] $\lambda$ 1.025$\mu$m is unprecedented for
a PISN candidate. However, existing PISN models do not predict significant
emission from [Co ii] because of line blocking. In Section 5.2.5, we proposed
that line blocking could be less severe than predicted by existing models.

The shape and the relative line intensities of the nebular spectra of SN
2018ibb are compatible with those predicted by PISN models. Our observations
reveal a significant excess at wavelengths shorter than 5000 Å, which should
not be present if 25–44 $M_{\odot}$ of iron group elements were formed. As we
concluded in Section 5.2.4, we propose that CSM interaction may account for
some, if not all, of the excess. PISN models consider mass loss (Kasen et al.
2011; Gilmer et al. 2017; Kozyreva et al. 2017; Dessart et al. 2013) in the
evolution of the progenitor star. However, their light curves and spectra are
computed assuming that any interaction between the SN ejecta and the
circumstellar material is negligible. Kasen et al. (2011) pointed out that the
CSM interaction could actually have a non-negligible contribution.
Furthermore, the CSM might not only be produced by stellar winds but also by
eruptions similar to that seen in Eta Carinae in 1843. That this is indeed a
non-negligible effect is corroborated by recent findings in Chen et al.
(2023b). These authors studied the light curves of 77 events from the
homogenous ZTF SLSN-I sample, and concluded that CSM is common around H-poor
SLSNe (in at least 25–44% of the events) and that it contributes to the
observed emission, albeit finding spectroscopic evidence in the spectra is
difficult. Owing to the lack of predictions of PISN models on CSM interaction,
we cannot firmly conclude that SN 2018ibb is a PISN.

Our observations demonstrate that interactions between the SN ejecta and the
ambient CSM play a non-negligible effect in the observed photometric and
spectroscopic properties (Sections 4.3.3, 4.3.4, 5.1, 5.2.4). PISN models of
H-poor progenitors with CSM are urgently needed. In the coming years, the
Rubin Observatory, and the James Webb, Euclid and Roman Space Telescopes will
systematically explore the high-redshift Universe. Since PISNe require metal-
poor stars and the early Universe was less chemically enriched than today,
PISNe are thought to be more abundant at higher redshifts. Several teams have
proposed search strategies to find PISNe with these new observing facilities
(e.g. Wang et al. 2017; Regős et al. 2020; Moriya et al. 2022b, a). However,
their search strategies are based on PISN models that, for instance, do not
include CSM interaction. Considering that these new observing facilities
either just started or will commence their science operations in the next
years, it is critical to expand the suite of existing PISN models in order to
find high-$z$ PISNe in real-time.

###  5.6 Could SN 2018ibb be a pulsational pair-instability supernova?

Table 11: Summary of the PISN tests applied on SN 2018ibb

Test | Condition | Observation | Section | Pass | Reference
---|---|---|---|---|---
Light curve
Rise time | 120 – 150 days |  $>93$ days | 4.2 | ? | 1, 2
Decline rate |  1 mag (100 day)-1 |  1.1 mag (100 day)-1 | 4.2 | ✓ | 1, 2
Peak absolute magnitude $M_{\rm bol}$ |  $-20$ – $-22.5$ mag |  $<-21.8$ mag | 4.2 | ✓ | 1, 2
Nickel mass |  10 – 40 $M_{\odot}$ |  25 – 40 $M_{\odot}$ |  5.2.1, 5.2.2, 5.2.3 | ✓ | 1, 2
PISN template | He100 – He130, | He120 – He130, | 5.2.3 | ✓ | 1, 2, 3
| P200 – P250 | P250, P250Ni34 |  |  |
Spectra
Photospheric velocity |  7000–11,000 km s-1 |  8500 km s-1 | 4.3.2 | ✓ | 3, 4
Photospheric velocity evolution | flat | flat | 5.3 | ✓ | here
Nebular spectra | He100, He130 | He130, but blue excess | 5.2.4 | ✘ | 5
[Co ii] $\lambda$ 1.025$\mu$m | not predicted | detected | 5.2.5 | ✘ | 5
Contribution from CSM interaction to the light curve and spectra
CSM interaction | not explored | observed |  5.1, 4.3.4, 5.2.4 | ? |
Host galaxy
Metallicity | $<Z_{\odot}/3$ | $0.25^{+0.07}_{-0.06}\leavevmode\nobreak\ Z_{\odot}$ | 4.6 | ✓ | 6

1) Kasen et al. (2011); 2) Kozyreva et al. (2017); 3) Gilmer et al. (2017); 4)
https://2sn.org/DATA/HW01/bulk_yields.txt based on Heger & Woosley (2002) 5)
Jerkstrand et al. (2016); 6) Langer et al. (2007)

The massive eruptions in a pulsational pair-instability supernova (PPISN) with
a large kinetic energy can, under the right conditions, be an ideal case for a
luminous interacting SN, as demonstrated in several studies (e.g. Woosley et
al. 2007; Yoshida et al. 2016; Woosley 2017; Leung et al. 2019; Marchant et
al. 2019; Renzo et al. 2020). While the PPI mechanism is difficult to avoid
for a He core in the mass range of 40–65 $M_{\odot}$ (Woosley et al. 2007),
the number of pulses and the interval between these, as well as the mass
ejected and their kinetic energies, are more uncertain and differ between
various studies. For a bright event to take place, the relative velocities
between the shells of the different ejections, as well as their relative
masses are important. The brightest event would result from the collision
between a very fast, massive shell and a shell of low or zero velocity. The
first shell must also be dense enough for the shocks to be radiative and
massive enough for the second shell to be completely decelerated. Finally, the
collision has to take place close enough to the star, on the order of
$10^{15}-10^{16}$cm, so that it will radiate the energy on a timescale of
approximately a year. This means that the interval between the pulses should
not be more than approximately a year. However, a collision at a very small
radius, and short time interval, will result in a very optically thick shell
where most of the released energy will go into adiabatic expansion. In
summary, there are a number of conditions which have to be fulfilled for a
bright SN to result. This has been illustrated in detail by the different
radiation-hydrodynamical models, for example, Woosley (2017).

Below, we discuss the most extreme models in order to judge whether a pure
PPISN could explain the large total radiated energy we find for SN 2018ibb.
For a pure He core, Woosley (2017) finds an upper limit to the kinetic energy
of $\sim 2\times 10^{51}$ erg, with the highest energy from the highest He
core mass, if no additional power source (e.g. magnetar or black hole
accretion) is involved. The most extreme model with a $62\leavevmode\nobreak\
M_{\odot}$ He core resulted in a 36 $M_{\odot}$ ejecta with a total kinetic
energy of $2.1\times 10^{51}$ erg. This is distributed over several pulses,
with most of the energy being dissipated in the first pulse. Without any
previous strong mass loss this will, however, not be converted into radiation
over a timescale of approximately a year. This is also confirmed by the light
curve models in Woosley (2017).

Brighter light curves could be obtained for models with a remaining hydrogen
envelope. The most extreme, T130D in Woosley (2017), had three pulses,
ejecting the 70 $M_{\odot}$ hydrogen-rich envelope with a kinetic energy of
$1.5\times 10^{51}$ erg. About 3300 years later a second pulse ejected a 7.7
$M_{\odot}$ shell with He, C and O and energy $1.1\times 10^{51}$ erg, and
after another 8 months a 13.5 $M_{\odot}$ shell and energy $1.5\times 10^{51}$
erg. The last two shells were ejected close enough in time to collide and
create a luminous SN with a total radiated energy of $4.5\times 10^{50}$ erg.

Similar calculations have been done by Marchant et al. (2019) and Leung et al.
(2019), using the MESA code (while Woosley 2017 used the Kepler code).
Qualitatively, these models agree, especially in the higher energies and mass
ejected, as well as the number of pulses with increasing He core mass. In
particular, Leung et al. (2020) find a maximum kinetic energy of $2.8\times
10^{51}$ erg for the highest PPI He core mass, similar to the corresponding
model by Woosley et al. (2007). However, as discussed by Leung et al. (2019),
there are also substantial quantitative differences between the models,
including ejected masses and time interval between the pulses. Some of the
differences can be traced back to the treatment of shocks, and convection in
both the hydrostatic and hydrodynamic phases.

We note that a large kinetic energy in PPISN models has also been invoked to
explain the light curves of other luminous SNe. For the FBOT AT2018cow, Leung
et al. (2020) invoked a kinetic energy of $5\times 10^{51}$ erg from a 42
$M_{\odot}$ He-core interacting with an ejected shell with mass of 0.5
$M_{\odot}$. An obvious solution to supply the extra energy is a hybrid model
with a combination of a PPISN and the energy from a magnetar or accretion.
This has been discussed by Woosley (2017) and for other energetic SLSNe
including PTF12dam (Tolstov et al. 2017a), Gaia16apd (Tolstov et al. 2017b)
and iPTF16eh (Lunnan et al. 2018a). However, it remains unclear how a magnetar
can be formed from the core collapse of the very massive He core in a PPISN.

In summary, a pure PPISN, close to the upper He core mass limit, may
potentially explain the observed radiated energy of $>3\times
10^{51}\leavevmode\nobreak\ \rm erg$ (Section 4.2.2). The conversion of
kinetic energy to radiative energy, however, requires rather special
conditions in terms of pulse intervals, ejecta mass and velocities. The
uncertainties in the models are, unfortunately, large, and it is difficult to
draw any firm conclusions. Additional energy sources can not be excluded, such
as a magnetar or a black hole. A contribution from a magnetar would result in
the flattening of the late-time light curve which is in stark contrast to our
observations. In the case a black hole was formed during the gravitational
collapse of the progenitor star, the accretion rate would need to be well-
tuned to be consistent with the exponentially declining light curve, making
the PPISN scenario less likely.

###  5.7 PISN rate constraint

Assuming that SN 2018ibb is a genuine PISN and that the PISN and SLSN
luminosity functions are roughly the same, we can measure the PISN rate as a
function of SLSNe. As of 25 July 2023, the spectroscopically complete ZTF
Bright Transient Survey (BTS; Fremling et al. 2020; Perley et al. 2020)
includes 24 genuine SLSNe-I with an absolute peak magnitude of $<-21$ mag and
an apparent peak magnitude of $<19$ mag.373737The BTS has a spectroscopic
completeness of 95% down to $m=18.5$ mag
(https://sites.astro.caltech.edu/ztf/bts/explorer.php). Although, the
completeness is lower at $m=19$, we can go down to that fainter magnitude
limit to get a minimum number of spectroscopically confirmed SLSNe and place
an upper limit on the PISN-to-SLSN rate. None of these objects are
photometrically similar to SN 2018ibb. Using Binomial statistics, we place an
upper limit of $14\%$ ($2\sigma$ confidence) on the PISN-to-SLSN rate. The
latest version of the BTS catalogue reports a volumetric SLSN-I rate of $5\pm
2\leavevmode\nobreak\ \rm Gpc^{-3}\,yr^{-1}$ between $z=0.1$ and $z=0.2$ (D.
A. Perley; priv. comm.). This translates the PISN-to-SLSN rate to a volumetric
rate of $<0.7\leavevmode\nobreak\ \rm yr^{-1}\,Gpc^{-3}$ in the same redshift
interval ($2\sigma$ confidence).

With a peak magnitude of $r\sim 17.7$ mag (Table 5), SN 2018ibb passes the
brightness cut of the BTS survey. However, it failed the BTS quality cuts as
it contaminated all ZTF reference images. Assuming that SN 2018ibb would have
passed the BTS quality cuts, we can place a lower limit on the volumetric
rate. Following the methodology in Perley et al. (2020), the volumetric rate
is

$R=\frac{N}{T\times 4/3\,\pi\,D^{3}_{L}\times f_{\rm sky}\times f_{\rm
retention}}\leavevmode\nobreak\ \\\ $

where $N$ is the number of events, $T$ is the BTS survey duration of 4.9
years, $D_{L}$ is the maximum distance at which a transient with a given
magnitude can be detected, $f_{\rm sky}$ is the average nightly BTS sky
coverage of 35%, and $f_{\rm retention}$ is the retention factor ($=1$, as SN
2018ibb is assumed to be detected). For $N$ and $D_{L}$, we apply the
following considerations. Firstly, a Poisson distribution with a mean of 1
($=\rm one$ PISN detection in the entire survey) has a lower $2\sigma$
confidence interval of 0.025 objects (Gehrels 1986). Secondly, at the BTS
completeness limit of 18.5 mag, SN 2018ibb would have been detected up to a
distance of 0.8 Gpc which we assume for $D_{L}$. This yields a lower limit of
$>9\times 10^{-3}\leavevmode\nobreak\ \rm yr^{-1}\,Gpc^{-3}$ at $2\sigma$
confidence (i.e. $0.2\%$ of the SLSN rate).

Expected PISN rates vary between $10^{-5}$ and $10^{-2}$ of the core-collapse
supernova rate (e.g. Langer et al. 2007; Briel et al. 2022). They are very
sensitive to the assumed shape of the IMF, the metallicity of the progenitor
star (Langer et al. 2007), stellar evolution models (e.g. nuclear reaction
rates, rotation, mixing, mass-loss history; Takahashi 2018; du Buisson et al.
2020; Farmer et al. 2020), the cosmic star-formation history (Briel et al.
2022), redshift (Langer et al. 2007; Briel et al. 2022), and whether the
progenitor system evolves as an isolated object or in a binary system (Briel
et al. 2022; Tanikawa et al. 2023). This precludes a detailed comparison of
the observed and expected rates. Our rate measurement of
$0.009$–$0.7\leavevmode\nobreak\ \rm yr^{-1}\,Gpc^{-3}$ is significantly
smaller than all but one prediction. Briel et al. (2022) computed the PISN
rate using various prescriptions for the cosmic star-formation history. For
one of the prescriptions, these authors deduced a rate of $\sim
0.8\leavevmode\nobreak\ \rm yr^{-1}\,Gpc^{-3}$, comparable with our upper
limit of $<0.7\leavevmode\nobreak\ \rm yr^{-1}\,Gpc^{-3}$. A more in-depth
comparison is beyond the scope of this paper.

##  6 Conclusion

In this paper, we have presented observations of the slow-evolving H-poor SLSN
2018ibb covering an exceptionally long time interval from $-93$ to +989 rest-
frame days after maximum. SN 2018ibb shares many similarities with H-poor
SLSNe, but its properties are extreme even for SLSNe. It is one of the slowest
evolving SLSNe known. The slow evolution is apparent through the long rise of
$>93$ rest-frame days from 10% peak flux to peak, the slow decline of merely
1.1 mag (100 days)-1, and the low photospheric velocity of 8500 km s-1 that
remains constant between the time of maximum and the following 100 rest-frame
days. At peak, SN 2018ibb reached an absolute magnitude of $M_{r}=-21.7$ mag,
comparable to the bulk of the SLSN population. The bolometric light curve had
a peak luminosity of $>2\times 10^{44}\leavevmode\nobreak\ \rm erg\,s^{-1}$.
During its lifetime, SN 2018ibb radiated $>3\times 10^{51}$ erg. The peak
luminosity and total radiated energy are strict lower limits.

We compared SN 2018ibb with PISN and SLSN models. SN 2018ibb complies with
most tests of PISN models with peak luminosity $<-20$ mag, and possibly all
tests with the interpretations that we propose, making SN 2018ibb the best
PISN candidate, to date. Specifically, SN 2018ibb passes the following tests:

  1. i)

a rise time of $>93$ days (expected: 120–150 days)

  2. ii)

a decline time scale of 1.1 mag (100 day)-1 (expected: 1.1 mag (100 day)-1)

  3. iii)

the modelling of the multi-band light curves with physical SLSN models and the
Katz et al. (2013) method point to the production of
25–$44\leavevmode\nobreak\ M_{\odot}$ 56Ni (expected:
10–$44\leavevmode\nobreak\ M_{\odot}$)

  4. iv)

the bolometric light curve is consistent with PISN templates that produce 25
and $44\leavevmode\nobreak\ M_{\odot}$ of 56Ni

  5. v)

a low ejecta velocity of $8500\leavevmode\nobreak\ \rm km\,s^{-1}$ (expected:
7000–$11,000\leavevmode\nobreak\ \rm km\,s^{-1}$)

  6. vi)

a flat evolution of the velocity after maximum light during the photospheric
phase (expected: flat evolution)

  7. vii)

a low metallicity of $\sim 0.25$ solar (expected: $<1/3$ solar)

  8. viii)

none of the $>200$ SLSNe has properties similar to SN 2018ibb (expected: PISNe
are rare).

Such a huge amount of nickel of 25–$44\leavevmode\nobreak\ M_{\odot}$ can only
be produced in a pair-instability-supernova explosion of a star with a He-core
mass of 120–130 $M_{\odot}$ at the time of the explosion (ZAMS mass of
approximately 240–$260\leavevmode\nobreak\ M_{\odot}$). However, SN 2018ibb
does not comply with the following tests:

  1. i)

the tentative detection of [Co ii] $\lambda$ 1.025 $\mu$m in emission,
implying $M(^{56}{\rm Ni})\gtrsim 30\leavevmode\nobreak\ M_{\odot}$ (expected:
no [Co ii] in emission)

  2. ii)

the nebular spectra are similar to the He130 [$M(^{56}{\rm
Ni})=44\leavevmode\nobreak\ M_{\odot}$] PISN model but show a substantial
excess bluewards of 5000 Å due to CSM interaction.

The tentative detection of [Co ii] is unprecedented for a PISN candidate and
any SLSN. It could be the smoking-gun evidence of SN 2018ibb being a PISN,
though the line identification hinges on the detection of a single line. PISN
models predict no significant [Co ii] $\lambda$ 1.025 $\mu$m in emission
because of line blocking extending by iron to the NIR. We propose that the
line blocking might be over-estimated in existing models.

While the late-time spectra are similar to PISN models, they also exhibit a
blue excess that should not be present due to the massive line-blanketing of
25–44 $M_{\odot}$ iron-group elements. A similar blue excess was also observed
in previous PISN candidates. Its presence was used as a critical piece of
evidence against the PISN interpretation. We argue that this is not the case
for SN 2018ibb. Three lines of evidence reveal that SN 2018ibb is not solely
powered by radioactivity and that CSM interaction is also at play: i) the
detection of a slow-moving CSM shell around the progenitor star; ii) the
presence of similarly slow O i, [O ii], [O iii] emission lines; and iii) a
blue pseudo-continuum similar to that of interaction-powered SNe. This
suggests that some, if not all, of the blue excess is produced by CSM
interaction. We stress that even after accounting for a substantial
contribution of CSM interaction to the bolometric flux,
25–$44\leavevmode\nobreak\ M_{\odot}$ are still required to power the entire
bolometric light curve.

PISN models consider mass-loss episodes (winds and to some level eruptions) to
evolve their progenitors to the point of explosion. However, the SN light
curves and spectra are computed in sterile environments, assuming that any
interaction between the SN ejecta and the circumstellar material is
negligible. Our observations demonstrate that CSM interaction is an important
non-negligible effect that needs to be systematically explored in PISN models.
The lack of such PISN models is the reason why we cannot conclusively argue
for SN 2018ibb being a PISN.

Our data set disfavours central engine models (magnetar powering and fallback
accretion onto a black hole), the magnetar+56Ni model and pure CSM models. The
continued linear decline out to $t_{\rm max}$+706 days and the absence of any
light curve flattening, expected for magnetar models, are in conflict with
existing analytical prescriptions of magnetar models. Furthermore, the
inferred values of the physical parameters of the magnetar and magnetar+56Ni
models are in conflict with existing stellar evolution models. A model with a
simple-power-law-shaped fallback accretion rate, the default assumption in
fallback models, would also result in a flattening of the light curve in
contradiction with our observations. Analytical CSM models did not provide an
adequate description either.

The extensive, high-quality dataset of SN 2018ibb is predestined to perform
definitive tests with SLSN and PISN models, and to explore rare explosion
mechanisms, such as axion-instability supernovae (AISNe; Sakstein et al.
2022). Simulations by Mori et al. (2023) suggest that AISNe evolve faster and
are bluer than PISNe for a given He-core mass. AISNe might also be more
abundant than PISNe. Therefore, revealing the powering mechanism of SN 2018ibb
will have immediate consequences not only for SN science but also for stellar
evolution theory. The final confirmation of a PISN would also have
ramifications for the interpretation of the observed drop in the black hole
mass function and, therefore, gravitational wave astronomy.

In the coming years, the Rubin Observatory, and the James Webb, Euclid and
Roman Space Telescopes will be used to search for SLSNe, PISNe, and the
explosions of Population III stars in the high-redshift Universe. To make this
leap forward, the community requires a significantly improved understanding of
the powering mechanisms and the progenitors of SLSNe. This can be accomplished
with i) comprehensive data sets of low-$z$ SLSNe similar to the one presented
here and ii) more complex theoretical models with clear predictions for light
curves and spectra. The James Webb Space Telescope could be transformative for
studying low-redshift SLSNe. Its IR spectrograph NIRspec has the sensitivity
to provide an uncensored view from 1 to 5 $\mu$m. Such an IR spectrum of a SN
2018ibb-like event could reveal strong emission lines from cobalt, nickel and
iron between 2 and 5 $\mu$m during the nebular phase, which would be the
smoking-gun evidence for powering by 56Ni.

###### Acknowledgements.

We thank the referee for a careful reading of the manuscript and for helpful
comments that improved this paper. We thank Stéphane Blondin (Laboratoire
d’Astrophysique de Marseille, France), Luc Dessart (Sorbonne Université,
France), Sebastian Gomez (Space Telescope Science Institute, USA), Ryosuke
Hirai (Monash University, Australia), Boaz Katz (Weizmann Institute of
Science, Israel), Keiichi Maeda (Kyoto University, Japan), Ilya Mandel (Monash
University, Australia), and Kanji Mori (Fukuoka University, Japan) for
fruitful discussions. U.C. Berkeley undergraduate students Nachiket Girish,
Andrew Hoffman, Evelyn Liu, Shaunak Modak, Jackson Sipple, Samantha Stegman,
Kevin Tang, and Keto Zhang helped obtain data with the Lick/Nickel telescope.
Z. Chen acknowledges support from the China Scholarship Council. A. V.
Filippenko’s supernova group at U.C. Berkeley received financial support from
the Christopher R. Redlich Fund, Gary & Cynthia Bengier, Clark & Sharon
Winslow, Sanford Robertson, Frank and Kathleen Wood (T. G. Brink is a Wood
Specialist in Astronomy), Alan Eustace (W. Zheng is a Eustace Specialist in
Astronomy), and numerous other donors. C. Fransson acknowledges support from
the Swedish Research Council and the Swedish National Space Board. J. P. U.
Fynbo acknowledges support from the Carlsberg Foundation. The Cosmic Dawn
Center (DAWN) is funded by the Danish National Research Foundation under grant
No. 140. M. Gromadzki is supported by the EU Horizon 2020 research and
innovation programme under grant agreement No. 101004719. A. Jerkstrand
acknowledges support from the European Research Council (ERC) under the
European Union’s Horizon 2020 Research and Innovation Programme (ERC Starting
Grant No. [803189]). H. Kuncarayakti was funded by the Academy of Finland
projects 324504 and 328898. G. Leloudas and M. Pursiainen are supported by a
research grant (19054) from VILLUM FONDEN. R. Lunnan is supported by the
European Research Council (ERC) under the European Union’s Horizon Europe
research and innovation programme (grant agreement No. 10104229 - TransPIre).
T. E. Müller-Bravo and L. Galbany acknowledge financial support from the
Spanish Ministerio de Ciencia e Innovación (MCIN), the Agencia Estatal de
Investigación (AEI) 10.13039/501100011033, the European Social Fund (ESF)
‘Investing in your future’, and the European Union Next Generation EU/PRTR
funds under the PID2020-115253GA-I00 HOSTFLOWS project, the 2019 Ramón y Cajal
program RYC2019-027683-I, the 2021 Juan de la Cierva program FJC2021-047124-I,
and from Centro Superior de Investigaciones Científicas (CSIC) under the PIE
project 20215AT016, and the program Unidad de Excelencia María de Maeztu
CEX2020-001058-M. M. Nicholl is supported by the European Research Council
(ERC) under the European Union’s Horizon 2020 research and innovation
programme (grant agreement No. 948381) and by a Fellowship from the Alan
Turing Institute. D. Polishook is grateful for the Wise Observatory staff. A.
Rossi acknowledges support from Premiale LBT 2013. M. Rigault has received
funding from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation programme (grant agreement No. 759194 -
USNAC). N. Sarin is supported by a Nordita Fellowship. Nordita is funded in
part by NordForsk. S. Schulze acknowledges support from the G.R.E.A.T.
research environment, funded by Vetenskapsrådet, the Swedish Research Council,
project number 2016-06012. L. J. Shingles acknowledges support by the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation program (ERC Advanced Grant KILONOVA No. 885281). L. Tartaglia
acknowledges support from MIUR (PRIN 2017 grant 20179ZF5KS). Y. Yang
acknowledges support from a Benoziyo Prize Postdoctoral Fellowship and the
Bengier-Winslow-Robertson Fellowship. This work was funded by ANID, Millennium
Science Initiative, ICN12_009.
Based in part on observations at the European Southern Observatory, Program
IDs 199.D-0143, 0105.D-0380, 0106.D-0524, 1103.D-0328, 2102.D-5026, and
2104.D-5006 (PIs C. Inserra, S. Schulze, and S. J. Smartt); Gemini-South,
Program ID 2021B-Q-901 (PI A. Gal-Yam); Hubble Space Telescope, Program ID
GO-16657 (PI C. Fremling); Keck, Program IDs C323, U023, U025 (PIs S. R.
Kulkarni, A. V. Filippenko); Large Binocular Telescope, Program ID
DDT_2019B_13 (PI E. Palazzi); Las Cumbres Observatory, Program IDs
FTPEPO2017AB-001, KEY2017AB-001, SUPA2019A-001, SUPA2019A-002, SUPA2019B-007,
and NOAO2020B-012 (PIs P. J. Brown, K. De); Liverpool Telescope, Program IDs
JL18B06, JL18B07, JL19A24, JL19B11, and JL20B15 (PI D. A. Perley); Nordic
Optical Telescope, Program IDs 57-502, 58-802, and 61-606, (PIs G. Leloudas,
J. Sollerman); P200 (PI L. Yan); and XMM-Newton, Program ID 08221501 (PI R.
Margutti).
We thank the staffs of the many observatories at which we conducted
observations.
This work has made use of data from the European Space Agency (ESA) mission
Gaia383838https://www.cosmos.esa.int/gaia, processed by the Gaia Data
Processing and Analysis
Consortium393939https://www.cosmos.esa.int/web/gaia/dpac/consortium (DPAC).
Funding for the DPAC has been provided by national institutions, in particular
the institutions participating in the Gaia Multilateral Agreement.
Part of the funding for GROND (both hardware as well as personnel) was
generously granted from the Leibniz-Prize to Prof. G. Hasinger (DFG grant HA
1850/28-1).
This work is based in part on observations made with the Large Binocular
Telescope (LBT). The LBT is an international collaboration among institutions
in Italy, the United States, and Germany. LBT Corporation partners are
Istituto Nazionale di Astrofisica, Italy; The University of Arizona on behalf
of the Arizona university system; LBT Beteiligungsgesellschaft, Germany,
representing the Max Planck Society, the Astrophysical Institute Potsdam, and
Heidelberg University; The Ohio State University; and The Research Corporation
on behalf of The University of Notre Dame, University of Minnesota, and
University of Virginia.
Some of the observations with the Las Cumbres Observatory data have been
obtained via OPTICON proposals and as part of the Global Supernova Project.
The OPTICON project has received funding from the European Union’s Horizon
2020 research and innovation programme under grant agreement No 730890.
This work made use of data supplied by the UK Swift Science Data Centre at the
University of Leicester.
CRTS is supported by the U.S. National Science Foundation (NSF) under grants
AST-0909182, AST-1313422, and AST-1413600. The Catalina Sky Survey (CSS) is a
NASA-funded project supported by the Near Earth Object Observation Program
(NEOO) under the Planetary Defense Coordination Office (PDCO).
This publication makes use of data products from the Two Micron All Sky
Survey, which is a joint project of the University of Massachusetts and the
Infrared Processing and Analysis Center/California Institute of Technology,
funded by NASA and the U.S. NSF.
Based in part on observations obtained with the Samuel Oschin Telescope
48-inch and the 60-inch Telescope at the Palomar Observatory as part of the
Zwicky Transient Facility project. ZTF is supported by the U.S. NSF under
grant AST-1440341 and a collaboration including Caltech, IPAC, the Weizmann
Institute of Science, the Oskar Klein Center at Stockholm University, the
University of Maryland, the University of Washington, Deutsches Elektronen-
Synchrotron and Humboldt University, Los Alamos National Laboratories, the
TANGO Consortium of Taiwan, the University of Wisconsin at Milwaukee, and
Lawrence Berkeley National Laboratories. Operations are conducted by COO,
IPAC, and UW. The SED Machine is based upon work supported by the U.S. NSF
under grant 1106171.
Partially based on observations made with the Nordic Optical Telescope, owned
in collaboration by the University of Turku and Aarhus University, and
operated jointly by Aarhus University, the University of Turku and the
University of Oslo, representing Denmark, Finland and Norway, the University
of Iceland and Stockholm University at the Observatorio del Roque de los
Muchachos, La Palma, Spain, of the Instituto de Astrofisica de Canarias. This
work makes use of observations from the Las Cumbres Observatory network. The
Las Cumbres Observatory team is supported by NSF grants AST-1911225 and
AST-1911151.
Some of the data presented herein were obtained at the W. M. Keck Observatory,
which is operated as a scientific partnership among the California Institute
of Technology, the University of California, and NASA; the observatory was
made possible by the generous financial support of the W. M. Keck Foundation.
KAIT, and its ongoing operation were made possible by donations from Sun
Microsystems, Inc., the Hewlett-Packard Company, AutoScope Corporation, the
Lick Observatory, the U.S. NSF, the University of California, the Sylvia & Jim
Katzman Foundation, and the TABASGO Foundation. A major upgrade of the Kast
spectrograph on the Shane 3 m telescope at Lick Observatory was made possible
through generous gifts from William and Marina Kast as well as the Heising-
Simons Foundation. Research at Lick Observatory is partially supported by a
generous gift from Google.

## References

  * Ageorges et al. (2010) Ageorges, N., Seifert, W., Jütte, M., et al. 2010, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7735, Ground-based and Airborne Instrumentation for Astronomy III, ed. I. S. McLean, S. K. Ramsay, & H. Takami, 77351L
  * Ambikasaran et al. (2015) Ambikasaran, S., Foreman-Mackey, D., Greengard, L., Hogg, D. W., & O’Neil, M. 2015, IEEE Transactions on Pattern Analysis and Machine Intelligence, 38, 252
  * Anderson (2018) Anderson, J. et al. 2018, EUROPEAN SOUTHERN OBSERVATORY, Doc. No. VLT-MAN-ESO-13100-1543
  * Andrews & Martini (2013) Andrews, B. H. & Martini, P. 2013, ApJ, 765, 140
  * Angus et al. (2019) Angus, C. R., Smith, M., Sullivan, M., et al. 2019, MNRAS, 487, 2215
  * Aoki et al. (2014) Aoki, W., Tominaga, N., Beers, T. C., Honda, S., & Lee, Y. S. 2014, Science, 345, 912
  * Appenzeller et al. (1998) Appenzeller, I., Fricke, K., Fürtig, W., et al. 1998, The Messenger, 94, 1
  * Arabsalmani et al. (2018) Arabsalmani, M., Møller, P., Perley, D. A., et al. 2018, MNRAS, 473, 3312
  * Ashton et al. (2019) Ashton, G., Hübner, M., Lasky, P. D., et al. 2019, ApJS, 241, 27
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481
  * Barkat et al. (1967) Barkat, Z., Rakavy, G., & Sack, N. 1967, Phys. Rev. Lett., 18, 379
  * Bauer et al. (2008) Bauer, F. E., Dwarkadas, V. V., Brandt, W. N., et al. 2008, ApJ, 688, 1210
  * Becker (2015) Becker, A. 2015, HOTPANTS: High Order Transform of PSF ANd Template Subtraction, Astrophysics Source Code Library, record ascl:1504.004
  * Belczynski et al. (2020) Belczynski, K., Klencki, J., Fields, C. E., et al. 2020, A&A, 636, A104
  * Bellm et al. (2019a) Bellm, E. C., Kulkarni, S. R., Barlow, T., et al. 2019a, PASP, 131, 068003
  * Bellm et al. (2019b) Bellm, E. C., Kulkarni, S. R., Graham, M. J., et al. 2019b, PASP, 131, 018002
  * Ben-Ami et al. (2014) Ben-Ami, S., Gal-Yam, A., Mazzali, P. A., et al. 2014, ApJ, 785, 37
  * Bertin (2012) Bertin, E. 2012, in Astronomical Society of the Pacific Conference Series, Vol. 461, Astronomical Data Analysis Software and Systems XXI, ed. P. Ballester, D. Egret, & N. P. F. Lorente, 263
  * Bertin & Arnouts (1996) Bertin, E. & Arnouts, S. 1996, A&AS, 117, 393
  * Blagorodnova et al. (2018) Blagorodnova, N., Neill, J. D., Walters, R., et al. 2018, PASP, 130, 035003
  * Blanchard et al. (2020) Blanchard, P. K., Berger, E., Nicholl, M., & Villar, V. A. 2020, ApJ, 897, 114
  * Blanchard et al. (2019) Blanchard, P. K., Nicholl, M., Berger, E., et al. 2019, ApJ, 872, 90
  * Blanton & Roweis (2007) Blanton, M. R. & Roweis, S. 2007, AJ, 133, 734
  * Blinnikov et al. (2006) Blinnikov, S. I., Röpke, F. K., Sorokina, E. I., et al. 2006, A&A, 453, 229
  * Blinnikov & Sorokina (2010) Blinnikov, S. I. & Sorokina, E. I. 2010, arXiv:1009.4353
  * Bradley et al. (2020) Bradley, L., Sipőcz, B., Robitaille, T., et al. 2020, astropy/photutils: 1.0.0, Zenodo
  * Branch et al. (2005) Branch, D., Baron, E., Hall, N., Melakayil, M., & Parrent, J. 2005, PASP, 117, 545
  * Brands et al. (2022) Brands, S. A., de Koter, A., Bestenlehner, J. M., et al. 2022, A&A, 663, A36
  * Breeveld et al. (2011) Breeveld, A. A., Landsman, W., Holland, S. T., et al. 2011, in American Institute of Physics Conference Series, Vol. 1358, American Institute of Physics Conference Series, ed. J. E. McEnery, J. L. Racusin, & N. Gehrels, 373–376
  * Briel et al. (2022) Briel, M. M., Eldridge, J. J., Stanway, E. R., Stevance, H. F., & Chrimes, A. A. 2022, MNRAS, 514, 1315
  * Bromm & Larson (2004) Bromm, V. & Larson, R. B. 2004, ARA&A, 42, 79
  * Burrows et al. (2005) Burrows, D. N., Hill, J. E., Nousek, J. A., et al. 2005, Space Sci. Rev., 120, 165
  * Buzzoni et al. (1984) Buzzoni, B., Delabre, B., Dekker, H., et al. 1984, The Messenger, 38, 9
  * Byler et al. (2017) Byler, N., Dalcanton, J. J., Conroy, C., & Johnson, B. D. 2017, ApJ, 840, 44
  * Calzetti et al. (2000) Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
  * Chambers et al. (2016) Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv:1612.05560
  * Chandra et al. (2012) Chandra, P., Chevalier, R. A., Chugai, N., et al. 2012, ApJ, 755, 110
  * Chandra et al. (2015) Chandra, P., Chevalier, R. A., Chugai, N., Fransson, C., & Soderberg, A. M. 2015, ApJ, 810, 32
  * Chatzopoulos & Wheeler (2012) Chatzopoulos, E. & Wheeler, J. C. 2012, ApJ, 760, 154
  * Chatzopoulos et al. (2012) Chatzopoulos, E., Wheeler, J. C., & Vinko, J. 2012, ApJ, 746, 121
  * Chen et al. (2020) Chen, K.-J., Woosley, S. E., & Whalen, D. J. 2020, ApJ, 897, 152
  * Chen et al. (2018) Chen, T. W., Inserra, C., Fraser, M., et al. 2018, ApJ, 867, L31
  * Chen et al. (2013) Chen, T.-W., Smartt, S. J., Bresolin, F., et al. 2013, ApJ, 763, L28
  * Chen et al. (2015) Chen, T. W., Smartt, S. J., Jerkstrand, A., et al. 2015, MNRAS, 452, 1567
  * Chen et al. (2017) Chen, T.-W., Smartt, S. J., Yates, R. M., et al. 2017, MNRAS, 470, 3566
  * Chen et al. (2023a) Chen, Z. H., Yan, L., Kangas, T., et al. 2023a, ApJ, 943, 41
  * Chen et al. (2023b) Chen, Z. H., Yan, L., Kangas, T., et al. 2023b, ApJ, 943, 42
  * Chevalier & Blondin (1995) Chevalier, R. & Blondin, J. M. 1995, ApJ, 444, 312
  * Chevalier & Fransson (1992) Chevalier, R. A. & Fransson, C. 1992, ApJ, 395, 540
  * Chevalier & Fransson (1994) Chevalier, R. A. & Fransson, C. 1994, ApJ, 420, 268
  * Chevalier & Fransson (2017) Chevalier, R. A. & Fransson, C. 2017, in Handbook of Supernovae, ed. A. W. Alsabti & P. Murdin, 875
  * Chornock & Filippenko (2008) Chornock, R. & Filippenko, A. V. 2008, AJ, 136, 2227
  * Conroy et al. (2009) Conroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486
  * Coppejans et al. (2018) Coppejans, D. L., Margutti, R., Guidorzi, C., et al. 2018, ApJ, 856, 56
  * Coppejans et al. (2021) Coppejans, D. L., Matthews, D., Margutti, R., et al. 2021, Transient Name Server AstroNote, 73, 1
  * Corsi et al. (2014) Corsi, A., Ofek, E. O., Gal-Yam, A., et al. 2014, ApJ, 782, 42
  * Crowther (2007) Crowther, P. A. 2007, ARA&A, 45, 177
  * Dark Energy Survey Collaboration et al. (2016) Dark Energy Survey Collaboration, Abbott, T., Abdalla, F. B., et al. 2016, MNRAS, 460, 1270
  * De Cia et al. (2018) De Cia, A., Gal-Yam, A., Rubin, A., et al. 2018, ApJ, 860, 100
  * de Ugarte Postigo et al. (2012) de Ugarte Postigo, A., Fynbo, J. P. U., Thöne, C. C., et al. 2012, A&A, 548, A11
  * Dekany et al. (2020) Dekany, R., Smith, R. M., Riddle, R., et al. 2020, PASP, 132, 038001
  * Dessart & Audit (2018) Dessart, L. & Audit, E. 2018, A&A, 613, A5
  * Dessart & Hillier (2005) Dessart, L. & Hillier, D. J. 2005, A&A, 439, 671
  * Dessart et al. (2022) Dessart, L., Hillier, D. J., & Kuncarayakti, H. 2022, A&A, 658, A130
  * Dessart et al. (2013) Dessart, L., Waldman, R., Livne, E., Hillier, D. J., & Blondin, S. 2013, MNRAS, 428, 3227
  * Dexter & Kasen (2013) Dexter, J. & Kasen, D. 2013, ApJ, 772, 30
  * Dey et al. (2019) Dey, A., Schlegel, D. J., Lang, D., et al. 2019, AJ, 157, 168
  * Drake et al. (2011) Drake, A. J., Djorgovski, S. G., Mahabal, A., et al. 2011, ApJ, 735, 106
  * Drlica-Wagner et al. (2018) Drlica-Wagner, A., Sevilla-Noarbe, I., Rykoff, E. S., et al. 2018, ApJS, 235, 33
  * du Buisson et al. (2020) du Buisson, L., Marchant, P., Podsiadlowski, P., et al. 2020, MNRAS, 499, 5941
  * Eftekhari et al. (2019) Eftekhari, T., Berger, E., Margalit, B., et al. 2019, ApJ, 876, L10
  * Eftekhari et al. (2021) Eftekhari, T., Margalit, B., Omand, C. M. B., et al. 2021, ApJ, 912, 21
  * Ekström et al. (2012) Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146
  * Elbaz et al. (2007) Elbaz, D., Daddi, E., Le Borgne, D., et al. 2007, A&A, 468, 33
  * Ellison et al. (2004) Ellison, S. L., Ibata, R., Pettini, M., et al. 2004, A&A, 414, 79
  * Evans et al. (2009) Evans, P. A., Beardmore, A. P., Page, K. L., et al. 2009, MNRAS, 397, 1177
  * Evans et al. (2007) Evans, P. A., Beardmore, A. P., Page, K. L., et al. 2007, A&A, 469, 379
  * Farmer et al. (2020) Farmer, R., Renzo, M., de Mink, S. E., Fishbach, M., & Justham, S. 2020, ApJ, 902, L36
  * Farmer et al. (2019) Farmer, R., Renzo, M., de Mink, S. E., Marchant, P., & Justham, S. 2019, ApJ, 887, 53
  * Fesen et al. (1999) Fesen, R. A., Gerardy, C. L., Filippenko, A. V., et al. 1999, AJ, 117, 725
  * Filippenko (1982) Filippenko, A. V. 1982, PASP, 94, 715
  * Filippenko (1997) Filippenko, A. V. 1997, ARA&A, 35, 309
  * Filippenko et al. (2001) Filippenko, A. V., Li, W. D., Treffers, R. R., & Modjaz, M. 2001, in Astronomical Society of the Pacific Conference Series, Vol. 246, IAU Colloq. 183: Small Telescope Astronomy on Global Scales, ed. B. Paczynski, W.-P. Chen, & C. Lemme, 121
  * Finkbeiner et al. (2016) Finkbeiner, D. P., Schlafly, E. F., Schlegel, D. J., et al. 2016, ApJ, 822, 66
  * Fiore et al. (2021) Fiore, A., Chen, T. W., Jerkstrand, A., et al. 2021, MNRAS, 502, 2120
  * Foreman-Mackey et al. (2014) Foreman-Mackey, D., Sick, J., & Johnson, B. 2014, Python-Fsps: Python Bindings To Fsps (V0.1.1)
  * Fowler & Hoyle (1964) Fowler, W. A. & Hoyle, F. 1964, ApJS, 9, 201
  * Fransson (1984) Fransson, C. 1984, A&A, 132, 115
  * Fransson et al. (2002) Fransson, C., Chevalier, R. A., Filippenko, A. V., et al. 2002, ApJ, 572, 350
  * Fransson et al. (2014) Fransson, C., Ergon, M., Challis, P. J., et al. 2014, ApJ, 797, 118
  * Fremling et al. (2018a) Fremling, C., Dugas, A., & Sharma, Y. 2018a, Transient Name Server Classification Report, 2018-1877, 1
  * Fremling et al. (2018b) Fremling, C., Dugas, A., Sharma, Y., Yaron, O., & Gal-Yam, A. 2018b, Transient Name Server Classification Report, 2018-2180, 1
  * Fremling et al. (2020) Fremling, C., Miller, A. A., Sharma, Y., et al. 2020, ApJ, 895, 32
  * Fremling et al. (2016) Fremling, C., Sollerman, J., Taddia, F., et al. 2016, A&A, 593, A68
  * Fynbo et al. (2009) Fynbo, J. P. U., Jakobsson, P., Prochaska, J. X., et al. 2009, ApJS, 185, 526
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1
  * Gal-Yam (2012) Gal-Yam, A. 2012, Science, 337, 927
  * Gal-Yam (2017) Gal-Yam, A. 2017, in Handbook of Supernovae, ed. A. W. Alsabti & P. Murdin, 195
  * Gal-Yam (2019a) Gal-Yam, A. 2019a, ApJ, 882, 102
  * Gal-Yam (2019b) Gal-Yam, A. 2019b, ARA&A, 57, 305
  * Gal-Yam et al. (2022) Gal-Yam, A., Bruch, R., Schulze, S., et al. 2022, Nature, 601, 201
  * Gal-Yam et al. (2009) Gal-Yam, A., Mazzali, P., Ofek, E. O., et al. 2009, Nature, 462, 624
  * Gehrels (1986) Gehrels, N. 1986, ApJ, 303, 336
  * Gehrels et al. (2004) Gehrels, N., Chincarini, G., Giommi, P., et al. 2004, ApJ, 611, 1005
  * Gerosa & Fishbach (2021) Gerosa, D. & Fishbach, M. 2021, Nature Astronomy, 5, 749
  * Gilkis et al. (2016) Gilkis, A., Soker, N., & Papish, O. 2016, ApJ, 826, 178
  * Gilmer et al. (2017) Gilmer, M. S., Kozyreva, A., Hirschi, R., Fröhlich, C., & Yusof, N. 2017, ApJ, 846, 100
  * Gimeno et al. (2016) Gimeno, G., Roth, K., Chiboucas, K., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9908, Ground-based and Airborne Instrumentation for Astronomy VI, ed. C. J. Evans, L. Simard, & H. Takami, 99082S
  * Goldoni et al. (2006) Goldoni, P., Royer, F., François, P., et al. 2006, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 6269, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, ed. I. S. McLean & M. Iye, 62692K
  * Gomez et al. (2019) Gomez, S., Berger, E., Nicholl, M., et al. 2019, ApJ, 881, 87
  * Graham et al. (2019) Graham, M. J., Kulkarni, S. R., Bellm, E. C., et al. 2019, PASP, 131, 078001
  * Greiner et al. (2008) Greiner, J., Bornemann, W., Clemens, C., et al. 2008, PASP, 120, 405
  * Guillochon et al. (2018) Guillochon, J., Nicholl, M., Villar, V. A., et al. 2018, ApJS, 236, 6
  * Hatano et al. (2001) Hatano, K., Branch, D., Nomoto, K., et al. 2001, in American Astronomical Society Meeting Abstracts, Vol. 198, American Astronomical Society Meeting Abstracts #198, 39.02
  * Heger et al. (2003) Heger, A., Fryer, C. L., Woosley, S. E., Langer, N., & Hartmann, D. H. 2003, ApJ, 591, 288
  * Heger & Woosley (2002) Heger, A. & Woosley, S. E. 2002, ApJ, 567, 532
  * HI4PI Collaboration et al. (2016) HI4PI Collaboration, Ben Bekhti, N., Flöer, L., et al. 2016, A&A, 594, A116
  * Hillebrandt & Niemeyer (2000) Hillebrandt, W. & Niemeyer, J. C. 2000, ARA&A, 38, 191
  * Hjorth & Bloom (2012) Hjorth, J. & Bloom, J. S. 2012, in “Gamma-Ray Bursts”, ed. S. E. W. C. Kouveliotou, R. A. M. J. Wijers (Cambridge University Press), 169–190
  * Hodgkin et al. (2021) Hodgkin, S. T., Harrison, D. L., Breedt, E., et al. 2021, A&A, 652, A76
  * Höflich (1991) Höflich, P. 1991, A&A, 246, 481
  * Hook et al. (2004) Hook, I. M., Jørgensen, I., Allington-Smith, J. R., et al. 2004, PASP, 116, 425
  * Hosseinzadeh et al. (2017) Hosseinzadeh, G., Arcavi, I., Valenti, S., et al. 2017, ApJ, 836, 158
  * Hosseinzadeh et al. (2022) Hosseinzadeh, G., Berger, E., Metzger, B. D., et al. 2022, ApJ, 933, 14
  * Hsiao et al. (2019) Hsiao, E. Y., Phillips, M. M., Marion, G. H., et al. 2019, PASP, 131, 014002
  * Huber et al. (2015) Huber, M., Chambers, K. C., Flewelling, H., et al. 2015, The Astronomer’s Telegram, 7153, 1
  * Inserra (2019) Inserra, C. 2019, Nature Astronomy, 3, 697
  * Inserra et al. (2016) Inserra, C., Bulla, M., Sim, S. A., & Smartt, S. J. 2016, ApJ, 831, 79
  * Inserra et al. (2018) Inserra, C., Nichol, R. C., Scovacricchi, D., et al. 2018, A&A, 609, A83
  * Inserra et al. (2017) Inserra, C., Nicholl, M., Chen, T. W., et al. 2017, MNRAS, 468, 4642
  * Inserra et al. (2013) Inserra, C., Smartt, S. J., Jerkstrand, A., et al. 2013, ApJ, 770, 128
  * Jakobsen et al. (2022) Jakobsen, P., Ferruit, P., Alves de Oliveira, C., et al. 2022, A&A, 661, A80
  * Jansen et al. (2001) Jansen, F., Lumb, D., Altieri, B., et al. 2001, A&A, 365, L1
  * Jeffery (1999) Jeffery, D. J. 1999, arXiv:9907015
  * Jerkstrand (2017) Jerkstrand, A. 2017, in Handbook of Supernovae, ed. A. W. Alsabti & P. Murdin, 795
  * Jerkstrand et al. (2015) Jerkstrand, A., Ergon, M., Smartt, S. J., et al. 2015, A&A, 573, A12
  * Jerkstrand et al. (2016) Jerkstrand, A., Smartt, S. J., & Heger, A. 2016, MNRAS, 455, 3207
  * Jerkstrand et al. (2017) Jerkstrand, A., Smartt, S. J., Inserra, C., et al. 2017, ApJ, 835, 13
  * Johnson et al. (2021) Johnson, B. D., Leja, J., Conroy, C., & Speagle, J. S. 2021, ApJS, 254, 22
  * Junde (1999) Junde, H. 1999, Nuclear Data Sheets, 86, 315
  * Kangas et al. (2022) Kangas, T., Yan, L., Schulze, S., et al. 2022, MNRAS, 516, 1193
  * Kasen et al. (2011) Kasen, D., Woosley, S. E., & Heger, A. 2011, ApJ, 734, 102
  * Kasen (2004) Kasen, D. N. 2004, PhD thesis, University of California, Berkeley
  * Katz et al. (2013) Katz, B., Kushnir, D., & Dong, S. 2013, arXiv e-prints, arXiv:1301.6766
  * Kawashimo et al. (2023) Kawashimo, H., Sawada, R., Suwa, Y., et al. 2023, arXiv:2306.01682
  * Kennicutt (1998) Kennicutt, Robert C., J. 1998, ARA&A, 36, 189
  * Könyves-Tóth & Vinkó (2021) Könyves-Tóth, R. & Vinkó, J. 2021, ApJ, 909, 24
  * Kozyreva & Blinnikov (2015) Kozyreva, A. & Blinnikov, S. 2015, MNRAS, 454, 4357
  * Kozyreva et al. (2014) Kozyreva, A., Blinnikov, S., Langer, N., & Yoon, S. C. 2014, A&A, 565, A70
  * Kozyreva et al. (2017) Kozyreva, A., Gilmer, M., Hirschi, R., et al. 2017, MNRAS, 464, 2854
  * Kozyreva et al. (2022) Kozyreva, A., Janka, H.-T., Kresse, D., Taubenberger, S., & Baklanov, P. 2022, MNRAS, 514, 4173
  * Kozyreva et al. (2018) Kozyreva, A., Kromer, M., Noebauer, U. M., & Hirschi, R. 2018, MNRAS, 479, 3106
  * Kozyreva et al. (2020) Kozyreva, A., Shingles, L., Mironov, A., Baklanov, P., & Blinnikov, S. 2020, MNRAS, 499, 4312
  * Kramida et al. (2018) Kramida, A., Yu. Ralchenko, Reader, J., & and NIST ASD Team. 2018, NIST Atomic Spectra Database (ver. 5.5.6), [Online]. Available: https://physics.nist.gov/asd [2018, April 6]. National Institute of Standards and Technology, Gaithersburg, MD.
  * Krühler et al. (2008) Krühler, T., Küpcü Yoldas, A., Greiner, J., et al. 2008, ApJ, 685, 376
  * Krühler et al. (2015) Krühler, T., Malesani, D., Fynbo, J. P. U., et al. 2015, A&A, 581, A125
  * Krumholz et al. (2019) Krumholz, M. R., McKee, C. F., & Bland-Hawthorn, J. 2019, ARA&A, 57, 227
  * Kumar et al. (2008) Kumar, P., Narayan, R., & Johnson, J. L. 2008, MNRAS, 388, 1729
  * Kuncarayakti et al. (2022) Kuncarayakti, H., Maeda, K., Dessart, L., et al. 2022, ApJ, 941, L32
  * Kwok et al. (2023) Kwok, L. A., Jha, S. W., Temim, T., et al. 2023, ApJ, 944, L3
  * Lacy et al. (2020) Lacy, M., Baum, S. A., Chandler, C. J., et al. 2020, PASP, 132, 035001
  * Lang et al. (2010) Lang, D., Hogg, D. W., Mierle, K., Blanton, M., & Roweis, S. 2010, AJ, 139, 1782
  * Langer et al. (2007) Langer, N., Norman, C. A., de Koter, A., et al. 2007, A&A, 475, L19
  * Larson et al. (2003) Larson, S., Beshore, E., Hill, R., et al. 2003, in AAS/Division for Planetary Sciences Meeting Abstracts, Vol. 35, AAS/Division for Planetary Sciences Meeting Abstracts #35, 36.04
  * Law et al. (2019) Law, C. J., Omand, C. M. B., Kashiyama, K., & et al. 2019, ApJ, 886, 24
  * Leloudas et al. (2012) Leloudas, G., Chatzopoulos, E., Dilday, B., et al. 2012, A&A, 541, A129
  * Leloudas et al. (2017) Leloudas, G., Maund, J. R., Gal-Yam, A., et al. 2017, ApJ, 837, L14
  * Leloudas et al. (2015) Leloudas, G., Schulze, S., Krühler, T., et al. 2015, MNRAS, 449, 917
  * Leung et al. (2020) Leung, S.-C., Blinnikov, S., Nomoto, K., et al. 2020, ApJ, 903, 66
  * Leung et al. (2019) Leung, S.-C., Nomoto, K., & Blinnikov, S. 2019, ApJ, 887, 72
  * Levan et al. (2013) Levan, A. J., Read, A. M., Metzger, B. D., Wheatley, P. J., & Tanvir, N. R. 2013, ApJ, 771, 136
  * Li & McCray (1992) Li, H. & McCray, R. 1992, ApJ, 387, 309
  * Li et al. (2003) Li, W., Filippenko, A. V., Chornock, R., et al. 2003, PASP, 115, 453
  * Liu et al. (2017) Liu, Y.-Q., Modjaz, M., & Bianco, F. B. 2017, ApJ, 845, 85
  * Lunnan et al. (2018a) Lunnan, R., Chornock, R., Berger, E., et al. 2018a, ApJ, 852, 81
  * Lunnan et al. (2014) Lunnan, R., Chornock, R., Berger, E., et al. 2014, ApJ, 787, 138
  * Lunnan et al. (2016) Lunnan, R., Chornock, R., Berger, E., et al. 2016, ApJ, 831, 144
  * Lunnan et al. (2018b) Lunnan, R., Fransson, C., Vreeswijk, P. M., et al. 2018b, Nature Astronomy, 2, 887
  * Mackey et al. (2003) Mackey, J., Bromm, V., & Hernquist, L. 2003, ApJ, 586, 1
  * Madau & Dickinson (2014) Madau, P. & Dickinson, M. 2014, ARA&A, 52, 415
  * Maoz et al. (2014) Maoz, D., Mannucci, F., & Nelemans, G. 2014, ARA&A, 52, 107
  * Marchant et al. (2019) Marchant, P., Renzo, M., Farmer, R., et al. 2019, ApJ, 882, 36
  * Margalit et al. (2018) Margalit, B., Metzger, B. D., Berger, E., et al. 2018, MNRAS, 481, 2407
  * Margutti et al. (2018) Margutti, R., Chornock, R., Metzger, B. D., et al. 2018, ApJ, 864, 45
  * Marino et al. (2013) Marino, R. A., Rosales-Ortega, F. F., Sánchez, S. F., et al. 2013, A&A, 559, A114
  * Masci et al. (2019) Masci, F. J., Laher, R. R., Rusholme, B., et al. 2019, PASP, 131, 018003
  * Matthews et al. (2021) Matthews, D., Margutti, R., Coppejans, D., et al. 2021, Transient Name Server AstroNote, 139, 1
  * Mazzali et al. (2019) Mazzali, P. A., Moriya, T. J., Tanaka, M., & Woosley, S. E. 2019, MNRAS, 484, 3451
  * McCrum et al. (2014) McCrum, M., Smartt, S. J., Kotak, R., et al. 2014, MNRAS, 437, 656
  * McKee & Ostriker (2007) McKee, C. F. & Ostriker, E. C. 2007, ARA&A, 45, 565
  * McKinney (2005) McKinney, J. C. 2005, ApJ, 630, L5
  * Medford et al. (2020) Medford, M. S., Lu, J. R., & Schlafly, E. F. 2020, Research Notes of the American Astronomical Society, 4, 38
  * Milisavljevic et al. (2012) Milisavljevic, D., Fesen, R. A., Chevalier, R. A., et al. 2012, ApJ, 751, 25
  * Milisavljevic et al. (2009) Milisavljevic, D., Fesen, R. A., Kirshner, R. P., & Challis, P. 2009, ApJ, 692, 839
  * Milisavljevic et al. (2018) Milisavljevic, D., Patnaude, D. J., Chevalier, R. A., et al. 2018, ApJ, 864, L36
  * Miyazaki et al. (2002) Miyazaki, S., Komiyama, Y., Sekiguchi, M., et al. 2002, PASJ, 54, 833
  * Modigliani et al. (2010) Modigliani, A., Goldoni, P., Royer, F., et al. 2010, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7737, Observatory Operations: Strategies, Processes, and Systems III, ed. D. R. Silva, A. B. Peck, & B. T. Soifer, 773728
  * Modjaz et al. (2016) Modjaz, M., Liu, Y. Q., Bianco, F. B., & Graur, O. 2016, ApJ, 832, 108
  * Moorwood et al. (1998) Moorwood, A., Cuby, J. G., & Lidman, C. 1998, The Messenger, 91, 9
  * Mori et al. (2023) Mori, K., Moriya, T. J., Takiwaki, T., et al. 2023, ApJ, 943, 12
  * Moriya et al. (2017) Moriya, T. J., Chen, T.-W., & Langer, N. 2017, ApJ, 835, 177
  * Moriya et al. (2022a) Moriya, T. J., Inserra, C., Tanaka, M., et al. 2022a, A&A, 666, A157
  * Moriya et al. (2018a) Moriya, T. J., Nicholl, M., & Guillochon, J. 2018a, ApJ, 867, 113
  * Moriya et al. (2022b) Moriya, T. J., Quimby, R. M., & Robertson, B. E. 2022b, ApJ, 925, 211
  * Moriya et al. (2018b) Moriya, T. J., Sorokina, E. I., & Chevalier, R. A. 2018b, Space Sci. Rev., 214, 59
  * Moriya et al. (2018c) Moriya, T. J., Terreran, G., & Blinnikov, S. I. 2018c, MNRAS, 475, L11
  * Moriya et al. (2019) Moriya, T. J., Wong, K. C., Koyama, Y., et al. 2019, PASJ, 71, 59
  * Murase et al. (2016) Murase, K., Kashiyama, K., & Mészáros, P. 2016, MNRAS, 461, 1498
  * Murase et al. (2021) Murase, K., Omand, C. M. B., Coppejans, D. L., et al. 2021, MNRAS, 508, 44
  * Nadyozhin (1994) Nadyozhin, D. K. 1994, ApJS, 92, 527
  * Nicholl et al. (2019) Nicholl, M., Berger, E., Blanchard, P. K., Gomez, S., & Chornock, R. 2019, ApJ, 871, 102
  * Nicholl et al. (2016a) Nicholl, M., Berger, E., Margutti, R., et al. 2016a, ApJ, 828, L18
  * Nicholl et al. (2016b) Nicholl, M., Berger, E., Smartt, S. J., et al. 2016b, ApJ, 826, 39
  * Nicholl et al. (2018) Nicholl, M., Blanchard, P. K., Berger, E., et al. 2018, ApJ, 866, L24
  * Nicholl et al. (2017) Nicholl, M., Guillochon, J., & Berger, E. 2017, ApJ, 850, 55
  * Nicholl et al. (2013) Nicholl, M., Smartt, S. J., Jerkstrand, A., et al. 2013, Nature, 502, 346
  * Nicholl et al. (2015a) Nicholl, M., Smartt, S. J., Jerkstrand, A., et al. 2015a, MNRAS, 452, 3869
  * Nicholl et al. (2015b) Nicholl, M., Smartt, S. J., Jerkstrand, A., et al. 2015b, ApJ, 807, L18
  * Oke et al. (1995) Oke, J. B., Cohen, J. G., Carr, M., et al. 1995, PASP, 107, 375
  * Oke & Gunn (1982) Oke, J. B. & Gunn, J. E. 1982, PASP, 94, 586
  * Omand & Jerkstrand (2023) Omand, C. M. B. & Jerkstrand, A. 2023, A&A, 673, A107
  * Omand et al. (2018) Omand, C. M. B., Kashiyama, K., & Murase, K. 2018, MNRAS, 474, 573
  * Omand & Sarin (2023) Omand, C. M. B. & Sarin, N. 2023, arXiv:2308.12997
  * Osterbrock (1989) Osterbrock, D. E. 1989, Astrophysics of gaseous nebulae and active galactic nuclei
  * Osterbrock & Ferland (2006) Osterbrock, D. E. & Ferland, G. J. 2006, Astrophysics of gaseous nebulae and active galactic nuclei
  * Ouchi et al. (2004) Ouchi, M., Shimasaku, K., Okamura, S., et al. 2004, ApJ, 611, 660
  * Perley (2019) Perley, D. A. 2019, PASP, 131, 084503
  * Perley et al. (2020) Perley, D. A., Fremling, C., Sollerman, J., et al. 2020, ApJ, 904, 35
  * Perley et al. (2016a) Perley, D. A., Quimby, R. M., Yan, L., et al. 2016a, ApJ, 830, 13
  * Perley et al. (2022) Perley, D. A., Sollerman, J., Schulze, S., et al. 2022, ApJ, 927, 180
  * Perley et al. (2016b) Perley, D. A., Tanvir, N. R., Hjorth, J., et al. 2016b, ApJ, 817, 8
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13
  * Pogge (2019) Pogge, R. W. 2019, rwpogge/modsCCDRed: v2.0.1
  * Pogge et al. (2010) Pogge, R. W., Atwood, B., Brewer, D. F., et al. 2010, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7735, Ground-based and Airborne Instrumentation for Astronomy III, 77350A
  * Poidevin et al. (2023) Poidevin, F., Omand, C. M. B., Könyves-Tóth, R., et al. 2023, MNRAS, 521, 5418
  * Prochaska et al. (2020a) Prochaska, J., Hennawi, J., Westfall, K., et al. 2020a, The Journal of Open Source Software, 5, 2308
  * Prochaska et al. (2007) Prochaska, J. X., Chen, H. W., Bloom, J. S., et al. 2007, ApJS, 168, 231
  * Prochaska et al. (2020b) Prochaska, J. X., Hennawi, J., Cooke, R., et al. 2020b, pypeit/PypeIt: Release 1.0.0
  * Pursiainen et al. (2018) Pursiainen, M., Castro-Segura, N., Smith, M., et al. 2018, The Astronomer’s Telegram, 12300, 1
  * Pursiainen et al. (2023) Pursiainen, M., Leloudas, G., Cikota, A., et al. 2023, A&A, 674, A81
  * Pursiainen et al. (2022) Pursiainen, M., Leloudas, G., Paraskeva, E., et al. 2022, A&A, 666, A30
  * Quimby et al. (2018) Quimby, R. M., De Cia, A., Gal-Yam, A., et al. 2018, ApJ, 855, 2
  * Quimby et al. (2011) Quimby, R. M., Kulkarni, S. R., Kasliwal, M. M., et al. 2011, Nature, 474, 487
  * Quinet (1998) Quinet, P. 1998, A&AS, 129, 147
  * Rakavy et al. (1967) Rakavy, G., Shaviv, G., & Zinamon, Z. 1967, ApJ, 150, 131
  * Regős et al. (2020) Regős, E., Vinkó, J., & Ziegler, B. L. 2020, ApJ, 894, 94
  * Renzo et al. (2020) Renzo, M., Farmer, R. J., Justham, S., et al. 2020, MNRAS, 493, 4333
  * Rigault et al. (2019) Rigault, M., Neill, J. D., Blagorodnova, N., et al. 2019, A&A, 627, A115
  * Rodrigo et al. (2012) Rodrigo, C., Solano, E., & Bayo, A. 2012, SVO Filter Profile Service Version 1.0, IVOA Working Draft 15 October 2012
  * Romero-Shaw et al. (2020) Romero-Shaw, I. M., Talbot, C., Biscoveanu, S., et al. 2020, MNRAS, 499, 3295
  * Roming et al. (2005) Roming, P. W. A., Kennedy, T. E., Mason, K. O., et al. 2005, Space Sci. Rev., 120, 95
  * Rybicki & Lightman (1986) Rybicki, G. B. & Lightman, A. P. 1986, Radiative Processes in Astrophysics
  * Sakstein et al. (2022) Sakstein, J., Croon, D., & McDermott, S. D. 2022, Phys. Rev. D, 105, 095038
  * Salvadori et al. (2019) Salvadori, S., Bonifacio, P., Caffau, E., et al. 2019, MNRAS, 487, 4261
  * Sarin et al. (2023) Sarin, N., Hübner, M., Omand, C. M. B., Setzer, C. N., & et al. 2023, arXiv:2308.12806
  * Schlafly & Finkbeiner (2011) Schlafly, E. F. & Finkbeiner, D. P. 2011, ApJ, 737, 103
  * Schulze et al. (2018) Schulze, S., Krühler, T., Leloudas, G., et al. 2018, MNRAS, 473, 1258
  * Schulze et al. (2021) Schulze, S., Yaron, O., Sollerman, J., et al. 2021, ApJS, 255, 29
  * Seifert et al. (2010) Seifert, W., Ageorges, N., Lehmitz, M., et al. 2010, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7735, Ground-based and Airborne Instrumentation for Astronomy III, ed. I. S. McLean, S. K. Ramsay, & H. Takami, 77357W
  * Selsing et al. (2019) Selsing, J., Malesani, D., Goldoni, P., et al. 2019, A&A, 623, A92
  * Serkowski et al. (1975) Serkowski, K., Mathewson, D. S., & Ford, V. L. 1975, ApJ, 196, 261
  * Shahbandeh et al. (2022) Shahbandeh, M., Hsiao, E. Y., Ashall, C., et al. 2022, ApJ, 925, 175
  * Sharon & Kushnir (2020) Sharon, A. & Kushnir, D. 2020, MNRAS, 496, 4517
  * Shingles et al. (2021) Shingles, L., Smith, K. W., Young, D. R., et al. 2021, Transient Name Server AstroNote, 7, 1
  * Silverman et al. (2012) Silverman, J. M., Foley, R. J., Filippenko, A. V., et al. 2012, MNRAS, 425, 1789
  * Silverman et al. (2013) Silverman, J. M., Nugent, P. E., Gal-Yam, A., et al. 2013, ApJS, 207, 3
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163
  * Smartt et al. (2015) Smartt, S. J., Valenti, S., Fraser, M., et al. 2015, A&A, 579, A40
  * Smith et al. (2020) Smith, K. W., Smartt, S. J., Young, D. R., et al. 2020, PASP, 132, 085002
  * Smith et al. (2016) Smith, M., Sullivan, M., D’Andrea, C. B., et al. 2016, ApJ, 818, L8
  * Smith (2006) Smith, N. 2006, ApJ, 644, 1151
  * Sobolev (1957) Sobolev, V. V. 1957, Sov. Ast., 1, 678
  * Sorokina et al. (2016) Sorokina, E., Blinnikov, S., Nomoto, K., Quimby, R., & Tolstov, A. 2016, ApJ, 829, 17
  * Speagle (2020) Speagle, J. S. 2020, MNRAS, 493, 3132
  * Stahl et al. (2019) Stahl, B. E., Zheng, W., de Jaeger, T., et al. 2019, MNRAS, 490, 3882
  * Stalder et al. (2017) Stalder, B., Tonry, J., Smartt, S. J., et al. 2017, ApJ, 850, 149
  * Steele (2004) Steele, I. A. 2004, Astronomische Nachrichten, 325, 519
  * Stetson (1987) Stetson, P. B. 1987, PASP, 99, 191
  * Stritzinger et al. (2002) Stritzinger, M., Hamuy, M., Suntzeff, N. B., et al. 2002, AJ, 124, 2100
  * Strüder et al. (2001) Strüder, L., Briel, U., Dennerl, K., et al. 2001, A&A, 365, L18
  * STScI Development Team (2013) STScI Development Team. 2013, pysynphot: Synthetic photometry software package, Astrophysics Source Code Library, record ascl:1303.023
  * Suzuki & Maeda (2021) Suzuki, A. & Maeda, K. 2021, ApJ, 908, 217
  * Swartz et al. (1995) Swartz, D. A., Sutherland, P. G., & Harkness, R. P. 1995, ApJ, 446, 766
  * Takahashi (2018) Takahashi, K. 2018, ApJ, 863, 153
  * Takei et al. (2022) Takei, Y., Tsuna, D., Kuriyama, N., Ko, T., & Shigeyama, T. 2022, ApJ, 929, 177
  * Tanikawa et al. (2023) Tanikawa, A., Moriya, T. J., Tominaga, N., & Yoshida, N. 2023, MNRAS, 519, L32
  * Terreran et al. (2017) Terreran, G., Pumo, M. L., Chen, T. W., et al. 2017, Nature Astronomy, 1, 713
  * The Dark Energy Survey Collaboration (2005) The Dark Energy Survey Collaboration. 2005, arXiv:0510346
  * The LIGO Scientific Collaboration et al. (2020) The LIGO Scientific Collaboration, the Virgo Collaboration, Abbott, R., et al. 2020, arXiv:2010.14533
  * The LIGO Scientific Collaboration et al. (2021) The LIGO Scientific Collaboration, the Virgo Collaboration, the KAGRA Collaboration, et al. 2021, arXiv:2111.03634
  * Theodosiou & Federman (1999) Theodosiou, C. E. & Federman, S. R. 1999, ApJ, 527, 470
  * Tinyanont et al. (2023) Tinyanont, S., Woosley, S. E., Taggart, K., et al. 2023, ApJ, 951, 34
  * Tody (1993) Tody, D. 1993, in Astronomical Society of the Pacific Conference Series, Vol. 52, Astronomical Data Analysis Software and Systems II, ed. R. J. Hanisch, R. J. V. Brissenden, & J. Barnes, 173
  * Tolstov et al. (2017a) Tolstov, A., Nomoto, K., Blinnikov, S., et al. 2017a, ApJ, 835, 266
  * Tolstov et al. (2017b) Tolstov, A., Zhiglo, A., Nomoto, K., et al. 2017b, ApJ, 845, L2
  * Tonry et al. (2018) Tonry, J., Denneau, L., Heinze, A., et al. 2018, Transient Name Server Discovery Report, 2018-1722, 1
  * Tonry (2011) Tonry, J. L. 2011, PASP, 123, 58
  * Turner et al. (2001) Turner, M. J. L., Abbey, A., Arnaud, M., et al. 2001, A&A, 365, L27
  * Umeda & Nomoto (2002) Umeda, H. & Nomoto, K. 2002, ApJ, 565, 385
  * van Dokkum (2001) van Dokkum, P. G. 2001, PASP, 113, 1420
  * van Marle et al. (2010) van Marle, A. J., Smith, N., Owocki, S. P., & van Veelen, B. 2010, MNRAS, 407, 2305
  * Vergani et al. (2015) Vergani, S. D., Salvaterra, R., Japelj, J., et al. 2015, A&A, 581, A102
  * Vernet et al. (2011) Vernet, J., Dekker, H., D’Odorico, S., et al. 2011, A&A, 536, A105
  * Vreeswijk et al. (2014) Vreeswijk, P. M., Savaglio, S., Gal-Yam, A., et al. 2014, ApJ, 797, 24
  * Vurm & Metzger (2021) Vurm, I. & Metzger, B. D. 2021, ApJ, 917, 77
  * Wang et al. (2017) Wang, L., Baade, D., Baron, E., et al. 2017, arXiv:1710.07005
  * Wang et al. (1997) Wang, L., Wheeler, J. C., & Höflich, P. 1997, ApJL, 476, L27
  * Wang et al. (2015) Wang, S. Q., Wang, L. J., Dai, Z. G., & Wu, X. F. 2015, ApJ, 799, 107
  * West et al. (2023) West, S. L., Lunnan, R., Omand, C. M. B., & et al. 2023, A&A, 670, A7
  * Woosley (2017) Woosley, S. E. 2017, ApJ, 836, 244
  * Woosley et al. (2007) Woosley, S. E., Blinnikov, S., & Heger, A. 2007, Nature, 450, 390
  * Wygoda et al. (2019) Wygoda, N., Elbaz, Y., & Katz, B. 2019, MNRAS, 484, 3941
  * Xing et al. (2023) Xing, Q.-F., Zhao, G., Liu, Z.-W., et al. 2023, Nature, 618, 712
  * Yan et al. (2017) Yan, L., Lunnan, R., Perley, D. A., et al. 2017, ApJ, 848, 6
  * Yan et al. (2015) Yan, L., Quimby, R., Ofek, E., et al. 2015, ApJ, 814, 108
  * Yaron & Gal-Yam (2012) Yaron, O. & Gal-Yam, A. 2012, PASP, 124, 668
  * Yoshida et al. (2016) Yoshida, T., Umeda, H., Maeda, K., & Ishii, T. 2016, MNRAS, 457, 351
  * Young et al. (2010) Young, D. R., Smartt, S. J., Valenti, S., et al. 2010, A&A, 512, A70
  * Yusof et al. (2013) Yusof, N., Hirschi, R., Meynet, G., et al. 2013, MNRAS, 433, 1114

##  Appendix A Photometric observations

###  A.1 All-sky surveys

##### ATLAS

— Between 2018 and 2019, the ATLAS survey utilised two 0.5 m telescope systems
on Haleakala and Mauna Loa, Hawaii (USA). Each of the telescopes has a field
of view of 29 square degrees. The two telescopes work in tandem to survey the
entire visible sky from $-45^{\circ}<\delta<90^{\circ}$ with a 2-day cadence.
ATLAS observes in two wide filters, called ‘cyan’ or ‘c’, which roughly covers
the SDSS/Pan-STARRS $g$ and $r$ filters (4200–6500 Å), and ‘orange’ or ‘o’,
which covers the SDSS/PanSTARRS $r$ and $i$ (5600–8200 Å) to a depth of $\sim
19$ mag ($5\sigma$, averaged over both telescope sites and weather conditions;
Tonry 2011). All data immediately go through an automatic data-processing
pipeline, described in Stalder et al. (2017).

After the field of SN 2018ibb appeared from behind the sun, the ATLAS survey
obtained the first image on 7 August 2018 ($t_{\rm max}$$-$107 days). On 10
September, SN 2018ibb got brighter than the detection threshold of the ATLAS
pipeline. To recover emission from SN 2018ibb between 7 August and 10
September and boost the quality of all data, we ran the ATLAS Forced
Photometry service404040https://fallingstar-data.com/forcedphot (Shingles et
al. 2021) on all images starting from summer 2018. We clipped and binned the
data using the python script
plot_atlas_fp.py414141https://gist.github.com/thespacedoctor/
gi86777fa5a9567b7939e8d84fd8cf6a76. We used a bin size of 2 days for the first
observing season and 5 days for the second observing season. Measurements with
a significance of $<3\sigma$ were removed from the final light curve. To