© The Authors 2025

1. Introduction

The final stages in the lives of massive stars are poorly known. In particular, key processes such as mass-loss mechanisms (e.g. through winds or eruptions) are difficult to characterise (e.g. Smith 2014). For example, Wolf-Rayet (WR) stars efficiently lose their external layers because of strong winds driven by radiation pressure (Abbott 1982), and when they explode, the resulting supernova (SN) lacks signatures of H, and in some cases they also lack He, such as in stripped-envelope SNe (SE SNe; Clocchiatti & Wheeler 1997). On the other hand, luminous blue variables (LBVs) are massive stars that undergo multiple eruptive mass-loss episodes (Humphreys & Davidson 1994). For example, Moriya et al. (2014) found mass-loss rates for SNe IIn above 10⁻³ M_⊙ yr⁻¹, while Dukiya et al. (2024) have found an astonishing mass-loss rate of 2 − 7 M_⊙ yr⁻¹ for the SN IIn ASASSN-14il. Typically, the mass loss is not steady with time. Each episode can shed a significant amount of mass from the progenitor star and is perhaps under the influence of a binary companion. Binarity is an important and perhaps even dominant factor in the production of SE SNe (Sana et al. 2012). In fact, massive stars are mostly found in binary systems, and the presence of a companion is likely to strip the donor star from its outer H layers after a common envelope (CE) phase (Podsiadlowski et al. 1992). This process can enhance the mass-loss rate even for lower-mass stars and generate a dense CSM around the donor (Chevalier 2012).

The net result of an enormous (0.01 − 0.1 M_⊙ yr⁻¹, Kiewe et al. 2012) mass loss is a massive circumstellar medium (CSM), which is revealed by narrow lines in the spectrum after the explosion (e.g. Type IIn SNe, Schlegel 1990, 1996; Filippenko 1997) and contributes to the luminosity through interaction with the SN ejecta. The luminosity of interacting SNe is mainly powered by the conversion of kinetic energy into radiation in the shock of the ejecta with the CSM. It is expected that the shock is stronger when a shell ejection occurs shortly before the explosion (less than a few years). In fact, the closer the CSM shell is to the progenitor, the denser it is, and thus the SN ejecta colliding into it give rise to a more intense shock wave. The density of the CSM and its spatial distribution affect the shape of the light curve (Khatami & Kasen 2024) and the strength of the emission lines, giving rise to asymmetries (Andrews et al. 2017) and peculiar line profiles such as double peaks (Andrews & Smith 2018). Not only that, the mass and density of the CSM can also affect the shape of the light curve (Khatami & Kasen 2024). When the CSM density is particularly high, the interaction can completely mask the internal power source to the extent that the underlying event can even be a non-terminal outburst (Vink 2015), which could be mistaken for a faint core-collapse (CC) SN or even a thermonuclear explosion (Silverman et al. 2013). Interestingly, Taddia et al. (2013) found that there does not appear to be a continuity of properties among SNe IIn but rather different subtypes, possibly hinting at the presence of different progenitors. Although the sample in Taddia et al. (2013) was quite small, this trend has been observed recently with bigger samples of SNe IIn (Hiramatsu et al. 2024; Ransome & Villar 2024).

Objects powered by strong interaction are rare. However, some extraordinary objects have been found in the past, and they present very well-sampled multiband light curves and spectra. A well-known SN that showed a huge interaction in the brightness and duration of the light curve is SN 2010jl (Smith et al. 2012b; Fransson et al. 2014; Ofek et al. 2014), a luminous Type IIn SN with a bright (∼ − 20 mag at peak) light curve that lasted more than 1000 days and spectra dominated by H Balmer lines with symmetric electron scattering-driven profiles. These characteristics are interpreted as signs of interaction with a massive (> 3 M_⊙) H-rich CSM possibly produced by an LBV progenitor (Fransson et al. 2014). The LBVs and hypergiant stars have been proposed as progenitors of several SNe IIn (Gal-Yam et al. 2007; Gal-Yam & Leonard 2009). Another example is SN 2013L (Andrews et al. 2017; Taddia et al. 2020), which also had an observed light curve lasting 1500 days and peaking at around −19 mag. The spectra are also slow-evolving and dominated by Balmer lines, but the profile of the Hα line is significantly asymmetric, with a blue shoulder that is interpreted as being due to the emergence of the shock wave not fully hidden by electron scattering (Taddia et al. 2020). We use these two SNe as reference throughout the paper.

Interaction can also occur for a limited time range, temporarily increasing the luminosity and changing the shape of the emission lines. This is the case, for example, of SN 1998S (Fassia et al. 2000, 2001), whose broad Balmer emission lines disappeared and reappeared again during the spectroscopic evolution due to the presence of at least two distinct CSM shells (Fassia et al. 2001).

Interaction is a powerful phenomenon that has been proposed to power many SNe, including the superluminous SN class (SLSNe; Smith et al. 2007; Gal-Yam 2012), and explain why they are so luminous and have such slow evolution. The mechanism powering SLSNe is still debated, but the slow rise seems to favour a central engine (particularly, a magnetar, Kasen & Bildsten 2010) in a massive ejecta with long diffusion time. However, we caution that Hiramatsu et al. (2024) have shown that there is no clear transition between SNe IIn and SLSNe and that the arbitrary cut should be removed.

In some cases, the CSM distribution may be strongly asymmetric. If it is also particularly close to the progenitor, the ejecta can quickly and completely engulf it. Narrow lines are then no longer visible in the spectra, and the broad lines only show P Cygni profiles typical of an expanding photosphere, but the interaction continues within the inner ejecta, providing energy to the light curve. In these cases, it is difficult to determine whether the luminosity is fully due to interaction or to a central engine such as a magnetar, but the extreme energy and duration of these events seem to point towards a combination of the two phenomena (Kangas et al. 2022; Pessi et al. 2023; Salmaso et al. 2023).

From a multimessenger point of view, interacting SNe IIn are particularly interesting because the shocked regions may provide a favourable environment for particle acceleration. These accelerated particles can decay into high-energy (HE) neutrinos, which in principle may be observed by neutrino detectors. There have been some tentative associations of interacting SNe with HE neutrino events. The first one was SN 2011fh, which was associated with a cascade event detected one day after the optical light curve peak. More recently, the Type Ibn SN 2023uqf has been found inside the errorbox of neutrino IC-231004A (Reusch et al. 2023; Stein et al. 2023) and more or less in time coincidence with the neutrino detection. However, given the cosmic SN rate, there is the possibility that this is a random association (Petropoulou et al. 2017), which is further supported by the fact that the neutrino errorbox is usually a couple of square degrees in size. Moreover, the current understanding is that even with state-of-the-art neutrino detectors, a SN would need to explode within a few Megaparsecs to produce a detectable flux (Valtonen-Mattila & O’Sullivan 2023). However, HE neutrino production depends on the strength of the interaction and the efficiency of the acceleration process. Therefore, it is important to provide empirical constraints on whether interacting SNe can produce HE neutrinos and to characterise their energetics and ejecta–CSM physical conditions.

In this paper, we focus on SNe that show clear narrow lines in their spectra and are thus classified as SNe IIn. We present a sample of four interacting SNe that display a high luminosity and slow spectro-photometric evolution. Analysis of the follow-up spectro-photometric data is used to derive global physical parameters such as kinetic energy, mass-loss rate, and mass of the CSM; to investigate their similarities and differences; and to explore the potential production of HE neutrinos. The paper is structured as follows: We give an overview of the objects and of their host galaxies in Sect. 2. In Sect. 3, we derive and analyse the bolometric light curves. In Sect. 4, we derive the CSM mass, radius, and mass-loss rate, and we compare the parameters with theoretical predictions for neutrino production. Finally, we summarise our work and present our conclusions in Sect. 5.

2. The sample

Our goal is the selection of strongly interacting SNe among all the SNe classified as IIn. An example of such objects found in the literature is SN 2010jl. These transients are the most appealing candidates for emitting HE neutrinos, although other less powerful objects cannot be ruled out. The SNe that correspond to our criteria in the literature are quite rare and often with little data, especially at very late phases. In this work, we build a sample of well-observed strongly interacting SNe to analyse their properties. The criteria we use to identify the new transients may be useful in the future to enlarge this sample.

The sample was built selecting newly discovered transients announced through services such as AstroNotes¹ and Astronomer’s Telegrams² that met our selection criteria. We evaluated both the classification spectrum and the light curve and selected SNe classified as Type IIn, with a relatively fast rising time among the luminous SNe (t_rise < 40 days), and high but not extreme luminosity (L_peak ∼ 10⁴³ erg s⁻¹). These criteria tend to exclude SLSNe, which are at least one order of magnitude more luminous (L_peak ≳ 10⁴⁴ erg s⁻¹ or <  − 21 mag) and tend to have a longer rise time to peak (t_rise ≳ 70 d) (Lunnan et al. 2015). Despite the claims in Hiramatsu et al. (2024), the debate on the mechanism powering SLSNe is not fully settled. Therefore, we decided to exclude them altogether from the sample, although this does not mean that they cannot be powered by interaction or produce HE neutrinos. Also, a preference was given to nearby events (typically within redshift z ∼ 0.1, so that a more complete follow-up could be ensured). Following our procedure, we selected 11 targets in three years, which are listed in Table 1, together with their coordinates, redshift, distance modulus, and galactic absorption. We also report some notable parameters derived from the pseudo-bolometric g, r, i light curve (for details on the computation, see Sect. 3), namely: the rise time t_rise, defined as the time between the last non-detection and the first peak; the peak luminosity L_peak; the integrated luminosity in the first ∼200 days E₂₀₀. These parameters were used as criteria to ensure that the sample is not polluted by objects that show interaction at very late times or that underwent brief bursts of interaction over short periods of time. The case of SN 2021adxl is a bit peculiar since its poorly constrained rise time is < 91 days, but given the light curve shape and luminosity, very similar to SN 2010jl (Fransson et al. 2014), it is likely that the explosion date was shortly before the peak. Therefore, we elected to keep it in the sample.

As we detail in the following, it turned out that not all SNe we selected with our criteria proved to evolve as expected and seven out of 11 candidates were excluded from more detailed analysis for not showing signs of strong interaction. Their light curve was generally dimmer than the four strongly interacting SNe that remained in the sample and faded faster. This can clearly be seen in Fig. 1, where we show their pseudo-bolometric light curve built using only optical filters (see Sect. 3 for the computation).

Fig. 1. Pseudo-bolometric g, r, i light curves of all the SNe we followed with our programme. Empty reverse triangles indicate the SNe that, after further analysis, were excluded from the final sample. All the phases are corrected for time dilation (this is true throughout the paper).

To discriminate the transients that display strong interaction, we used the measured luminosity integrated over the first 200 days and imposed that this parameter E₂₀₀ ≥ 3 × 10⁴⁹ erg. This criterion is determined a posteriori to identify the objects with prolonged, strong interaction. In fact, with longer integration times, the contamination arises from objects undergoing prolonged but mild interaction (e.g. SN 2005ip, Taddia et al. 2013, E₂₀₀ ≤ 1 × 10⁴⁹ erg). Conversely, with shorter integration times there is significant contamination from objects experiencing strong interaction for a limited period (e.g. SN 1998S). Our parameter is a compromise between the two cases and ensures that all the selected SNe are strongly interacting. The excluded SNe had a lower energy display in the first 200 days compared to the selected ones. This can be seen from Fig. 2, where we also add the values calculated for the strongly interacting SN 2010jl, the mildly interacting SN 1998S, and the normal Type IIP SN 2004et (Sahu et al. 2006) as reference for different amounts of interaction. We also added the seven SNe IIn in the sample by Taddia et al. (2013): SNe 2005ip, 2005kj, 2006aa, 2006bo, 6006jd, 2006qq, and 2008fq. These proved to be polluters since, although they show signs of interaction in the spectra, they are less luminous and their light curve has a shorter duration, thus implying the presence of less interaction. An intermediate case is that of SN 2021gci, which has significantly more energy than the others, and with a more luminous first peak, albeit still one order of magnitude below the selected transients. On the other hand, SNe 2022iaz and 2022owx have more luminous peaks but faded too fast, providing a small total amount of energy.

Fig. 2. Integrated luminosity in the first 200 days for each SN in our sample and some objects from the literature for comparison. The dashed vertical line divides the strongly interacting sub-sample from the more mildly interacting SNe.

This exercise denotes that a single spectrum and the early light curve is not enough to discern between strongly interacting and normal interacting SNe, thus causing the sample to be contaminated with objects that have to be rejected later. The very minimum information required to identify a strongly interacting SN includes at least photometric coverage of 2/3 months and a couple of good resolution spectra with enough signal-to-noise (S/N) at different epochs to check the evolution. Ideally, one should have a well-sampled light curve for the first 200 days to properly distinguish between strongly interacting SNe and SNe with low interaction based on the energy input.

The four SNe selected here that fulfill our criteria are then 2021acya, 2021adxl, 2022qml, and 2022wed. In the following, we describe each object in more detail. We note that a paper presenting detailed observations for one of the targets in our list, SN 2021adxl, was recently published (Brennan et al. 2024). We refer to their observations and results in the relevant sections.

In Fig. 3, we show the r-band absolute light curves for the four SNe in our sample along with SN 2010jl. This SN is our reference throughout the whole analysis by virtue of its spectro-photometric characteristics, which perfectly fall within our parameters. Details on the observations and data reduction are given in Appendices A and B. As the explosion date, we take the middle epoch between the last non-detection and the discovery date for all SNe but for SN 2021adxl (cf. Sect. 2.2). To compute the absolute magnitude, we corrected the apparent magnitudes for Galactic extinction and distance (see Appendix B.3 for details). The transients span a range of 2 mag in absolute magnitude at maximum light and show different light curve shapes, although all of them are in general more long-lasting than ordinary Type II SNe, which fade on faster time scales (∼2 mag in 100 days).

Fig. 3. Absolute magnitude r-band light curves for the SNe IIn in our sample as well as the well-studied Type IIn SN 2010jl (Smith et al. 2012b; Fransson et al. 2014) for comparison. Black arrows indicate upper limits. All photometry is K-corrected and the phases are corrected for time dilation.

2.1. SN 2021acya

SN 2021acya was discovered on 30 October 2021 (MJD 59518.029) in the orange band by the Asteroid Terrestrial-impact Last Alert System (ATLAS, Tonry et al. 2018) at a magnitude of 18.124, while the last non-detection was only two days prior (Tonry et al. 2021). As the explosion epoch, then, we took MJD 59517 ± 1. It was then classified as SN IIn on 25 November 2021 (Ragosta et al. 2021).

The rise to the peak is relatively slow, as it reaches the brightest absolute magnitude of −20.3 ± 0.1 in the r band on MJD 59537 ± 1, 20 days after the first detection. A slow decline follows the peak, which flattens into a plateau at ∼160 days lasting for ∼150 days. A linear luminosity decline then resumes and lasts up to 480 days after the explosion, when the SN is finally lost. The shape of the late-time r-band light curve from the starting of the plateau matches well that of SN 2010jl (Fig. 3) but it is slightly brighter at all phases, probably indicating a stronger interaction in the case of SN 2021acya.

The spectral evolution of SN 2021acya is shown in Fig. 4. At first glance, the spectra seem to show a very slow evolution, which is typical of long interacting SNe. Upon closer inspection, however, there are significant changes in the continuum shape and the width of the emission lines. The first spectrum of SN 2021acya at +24 days shows a hot (10 000 K) continuum and the only prominent features are Balmer emission lines. The line profile is almost symmetrical, with an electron scattering profile. Therefore, the width of the lines cannot be used to trace the bulk motion of the gas. The spectrum at +27 days is slightly cooler and we can detect emission from He Iλ5876. Afterwards, the continuum cools down to 8000 K and, from phase +56 days, a bump starts to emerge in the bluer part of the spectrum that becomes more and more evident. This feature is usually attributed to a plethora of Fe lines, such as Fe II, that show up when the medium is heated by a strong shock (Chevalier & Fransson 1994; Mazzali et al. 2001). At +101 days a broad emission from the Ca II NIR triplet λλλ8498, 8542, 8662 appears, and its intensity increases as the evolution proceeds. On the spectrum at +112 days we tentatively identify [Ne III] λ3869 but the low S/N makes it difficult to see the feature in other spectra. He I is visible until +259 days, while it is not detected in the spectrum at +346 days. The following spectra are all similar and dominated by Hα, Hβ, and Ca NIR.

Fig. 4. Spectral evolution of SN 2021acya as shown through the most significant spectra. Numbers close to the spectra indicate the phase from explosion corrected for time dilation. All spectra are scaled with respect to the Hα and arbitrarily shifted for better visualisation. The spectra are not corrected for extinction.

2.2. SN 2021adxl

SN 2021adxl was discovered on 3 November 2021 (MJD 59521.540) by Fremling (2021) with a magnitude of 14.41 in the r band but its last non-detection was on MJD 59431.690, almost three months before the discovery. It was then classified as SN IIn on 2 February 2022 (De et al. 2022). To aid in the comparison, we adopted as explosion date a day closer to the first observation considering the similarity of the bolometric light curve with SN 2010jl, which shows a fast (∼20 days) rise to the peak. In fact, the peak magnitude is only slightly higher while the initial decline is almost identical for the two SNe and also the colour evolution is similar (cf. Appendix B.3). For these reasons, we arbitrarily chose MJD ${59500}_{-78}^{+22}$ as the explosion epoch, to match the peak to that of SN 2010jl. Given this, the luminosity in the r band peaks at −20.4 ± 0.2 mag on MJD 59529 ± 4, 29 days after our estimated explosion date.

The spectral evolution of SN 2021adxl is shown in Fig. 5. The earliest spectrum of SN 2021adxl was taken at +89 days and the continuum is already almost flat. The main emission lines are those of the Balmer series and He Iλ5876, the latter showing a distinct P Cygni profile, as well as a small bump from Ca II NIR. The blue bump due to the Fe II forest, in this case, is less accentuated. We also clearly see the absorption component of the P Cygni profile of Fe IIλ5169. The presence of broad He and Fe P Cygni profiles when they are absent in H in a SN IIn is unusual and be further discussed in Sect. 2.6. There are also narrow emissions from [O III] λλ4959, 5007 due to the background starburst region. Notably, the Hα and Hβ lines are highly asymmetric. The profile has a flat top and electron scattering wings on the sides but with a distinct blue shoulder that was observed in SN 2013L (Taddia et al. 2020). This composite profile is attributed to the shock front being exposed, giving the boxy profile, and to electron scattering forming the wings (Taddia et al. 2020). The later spectra are almost identical until phase +197 days, when the Balmer lines have nearly symmetrical profiles with electron scattering wings. At this point, the P Cygni profiles on He I and Fe II have also disappeared. The evolution then proceeds with the shrinking of the broad emission lines until the last spectrum at +539 days, which is dominated by a faint, broad Hα emission. We do not detect high ionisation lines, probably because the spectral resolution is not very high and the intensity of these features is low. However, Brennan et al. (2024) identify [Ne V] λ3346, [O III] λ4363, [Ne V] λ3346, [Ca V] λ6086, [Fe VII] λ6087, and [Fe X] λ6365 in their spectrum at +480 days.

Fig. 5. Spectral evolution of SN 2021adxl as shown through the most significant spectra. Numbers close to the spectra indicate the phase from explosion corrected for time dilation. All spectra are scaled with respect to the Hα and arbitrarily shifted for better visualisation. The spectra are not corrected for extinction.

In Fig. 6, we show a zoom-in on the Hα region on the spectrum of SN 2021adxl at +91 days compared to a spectrum of SN 2013L (Taddia et al. 2020) at a similar phase, which showed a peculiar blue shoulder on the broad Hα emission. Taddia et al. (2020) explain the profile as the combination of a boxy profile originating in the shocked shell with the extended wings caused by electron scattering, where the red side of the boxy component would be lost due to occultation of the receding shocked shell by the inner ejecta with high optical depth. To show this, a red, dashed box is added to the plot, roughly indicating the profile of a pure shocked shell with a velocity of 3300 km s⁻¹, devoid of electron scattering. On the blue side, the box is matched to the blue shoulder, while its reflection on the red side shows the missing flux. For such strong occultation, the line-emitting region must be located just above the opaque ejecta (the similar profile of SN 2022qml are commented on in the next section).

Fig. 6. Spectral comparison between SNe 2021adxl, 2022qml, and 2013L from Taddia et al. (2020) (zoom on the Hα line). The red dashed line highlights the boxy region of the line profile due to the shock. All spectra have been redshift corrected, continuum subtracted, and rescaled for better visualisation.

2.3. SN 2022qml

The SN was discovered on 2 August 2022 (MJD 59794.030) at a magnitude 18.137 in the cyan band (Tonry et al. 2022) and classified as SN IIn on 27 August 2022 (Gutierrez et al. 2022). The peak is bright, at −19.46 ± 0.1 mag in the r band on MJD 59800 ± 1, and showing a linearly declining light curve after that. The rise was very rapid, given the last non-detection only one day before the discovery. We adopt a best estimate of the explosion epoch as MJD 59793 ± 1.

The spectral evolution of SN 2022qml is shown in Fig. 7. The spectrum at +24 days has a blue continuum, implying a high temperature, and the only notable feature is the narrow Hα on top of a very broad but shallow emission. The spectrum at +49 days starts to show a distinctive blue bump due to a forest of Fe II lines and a broad, asymmetric Hα with a boxy profile smoothed by electron scattering wings. We also identify the narrow line of [Fe X] λ6375 and, possibly, [Ne V] λ3426, [Fe XI] λ7892 and a broader emission that could be due to a blend of coronal [Fe IV] λ5303 with [Ca V] λ5309. The absorption in the P Cygni profile of the narrow Hα is also clearly visible, which allows the velocity of the progenitor wind to be measured ³ (∼100 km s⁻¹). This feature is not identifiable in the other spectra of this object due to their lower resolution. At +56 days the only prominent features are Balmer lines, [Fe X] λ6375, and HeIıλ5876. After that, the high ionisation lines disappear while the spectra remain dominated by the Hα, until phase +111 days, when the Ca II NIR λλλ8498, 8542, 8662 emerges with a broad, symmetric emission. A zoom on the Hα region of this spectrum is also plotted in Fig. 6. The shape of the blue shoulder is similar to SNe 2021adxl and 2013L but it is broader, indicating a higher shock velocity.

Fig. 7. Spectral evolution of SN 2022qml as shown through the most significant spectra. Numbers close to the spectra indicate the phase from explosion corrected for time dilation. All spectra are scaled with respect to the Hα and arbitrarily shifted for better visualisation. The spectra have not been corrected for extinction. We also plot the narrow Hα in velocity space to show the P Cygni profile due to the wind that formed the CSM.

The strong blue bump of SN 2022qml is reminiscent of a feature seen in Type Ia-CSM SNe (e.g. SN 2002ic, Kotak et al. 2004) but also in Type Ic SNe (e.g. SN 1998bw Galama et al. 1998; Kulkarni et al. 1998). As we mentioned, the presence of the blue bump is usually attributed to a forest of Fe II lines. In the case of SN 2022qml, the strong bump, which is the most extreme in our sample, suggests a high Fe abundance in the ejecta. A massive star that produces Fe must also show high abundance of O and Ca, as SN 1998bw did. However, SN 2022qml spectra do not show any O I features and only a small bump likely due to Ca NIR. Also, the shape of the light curve and its faster decline, compared to the other SNe in our sample, is suggestive of a different origin for this SN. It fact, it is possible that it was not a SN IIn but rather a Type Ia-CSM SN. This would be in line with what found by Leloudas et al. (2015), who showed that simulated thermonuclear SN spectra were consistently misclassified as SNe IIn once the underlying SN flux was a fraction ∼0.2 − 0.3 or below with respect to the continuum. If this is the case, we probably missed the emergence of Si lines because the SN ejecta were still embedded in the CSM cocoon when the SN faded.

2.4. SN 2022wed

SN 2022wed was discovered on 21 September 2022 (MJD 59843.954) (Fremling 2022) at a magnitude 20.46 in the r band, while the last non-detection was on MJD 59839.484. As explosion epoch we chose MJD 59841 ± 2. It was then classified as SN IIn on 27 February 2023 (Hiramatsu et al. 2023).

The r light curve shows a first peak at −19.20 ± 0.03 mag on MJD 59863 ± 2, 22 days after the explosion. After a short decline lasting around 40 days, the light curve shows a second, broader and brighter peak at −19.68 ± 0.01 mag +230 days after the explosion followed by a very slow decline up to +460 days.

The spectral evolution of SN 2022wed is shown in Fig. 8. The first spectrum was taken at +122 days and shows an already cool (6000 K) continuum with the distinctive blue bump of Fe II, which is however less pronounced than for SN 2022qml. The main emission lines are the Balmer series and He I at λ5876 and λ7065 Å. The Hα profile in this SN is symmetric with electron scattering wings. The spectrum at +347 days shows a small bump at the position of [Ca II] λλ7291, 7324. The Hα feature has shrunk considerably and the Lorentzian wings are almost invisible. The last spectrum shares the same features, but Hα is narrower and more symmetric, while the blue bump seems to decrease slightly. Moreover, the Ca NIR triplet λλλ8498, 8542, 8662 finally appears.

Fig. 8. Spectral evolution of SN 2022wed. Numbers close to the spectra indicate the phase from explosion corrected for time dilation. All spectra have been scaled with respect to the Hα and arbitrarily shifted for better visualisation. The spectra have not been corrected for extinction.

2.5. Host galaxies

The main characteristics of the host galaxies of the SNe in our sample are summarised in Table 2. The SN host galaxies as retrieved from the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS, PS1, Chambers et al. 2019) in the r band are shown in Fig. 9. Interestingly, all known hosts are small, probably dwarf galaxies.

Fig. 9. Pre-explosion r-band images of the host galaxies of our SNe IIn retrieved from PS1. Red circles indicate the SNe position. Top left: SN 2021acya. Top right: SN 2021adxl. Bottom left: SN 2022qml. Bottom right: SN 2022wed.

The case of SN 2022wed deserves some discussion, as the SN is close to three IR sources listed on NED (WISEA J072415.72+190456.3, WISEA J072415.05+190455.7, and WISEA J072415.95+190449.4) whose redshifts are unknown. Given the angular separation and assuming these galaxies have the same redshift of the SN, the radial distances from the centre of each galaxy would be 11, 15, and 17 kpc, respectively, which allows the SN to belong to any of them. However, in the pre-explosion image produced by stacking PS1 observations between MJD 55182 − 56638, a faint source can be seen at the location of the SN (see Fig. 9, the bottom right panel). Its position is slightly offset from the SN, 0.01 arcsec in right ascension and 0.8 arcsec in declination, which assuming it is at the same redshift of the SN, gives a radial distance of 1.7 kpc. Its full width at half maximum (FWHM) is about 1.3 arcsec and is comparable with a point-like source in the image we have. We applied a point-spread function (PSF) fit to the source and found an apparent magnitude r = 22.4 ± 0.1 mag, which if we assume the same redshift of the SN, translates to an absolute magnitude of −16 mag, similar to that of the host galaxies of the other SNe in the sample.

Motivated by the presence of narrow [O III] and by the archival image of the host of SN 2021adxl, which shows a bright spot in correspondence with the SN position (see Fig. 9, the panel on the upper right), we attempted to use the line ratios of [N II] λ6583/Hα, [O III] λ5007/Hβ, and [O II] λ3727/[O III] λ5007 to derive an estimate of the oxygen abundance (which is assumed to trace the metallicity). We performed the measurements on the spectrum at +539 days, since it is the one where the SN contribution is smaller and the narrow lines due to the H II region are more evident. Unfortunately, we are only able to measure an upper limit for [N II] λ6583 because the line cannot be resolved from the strong and broad Hα but the metallicity is consistent with a solar/subsolar composition according to the N2-versus-O3 calibrator. Based on a similar analysis but with spectra of better resolution, Brennan et al. (2024) report a subsolar (0.1 Z_⊙) metallicity. Their relative flux measurements are consistent with ours but for Hα, which is three times lower than ours. This is not surprising considering that there is probably some residual contamination from the SN. For [N II] they have an actual measure 100 times smaller than our conservative limit. The metallicity measured by Brennan et al. (2024) is lower than the ones reported for SNe 2010jl and 2013L, both around 0.3 Z_⊙ (Stoll et al. 2011; Taddia et al. 2020). This is consistent with regions of high star formation in low-mass galaxies such as the host of SN 2021adxl (Yates et al. 2012). For the other SNe, the hosts are too faint to properly distinguish the structure of the galaxy but, interestingly, the hosts of SNe 2021acya and 2022qml are also UV bright (as per the magnitudes reported in NED), which is indicative of a high specific star formation rate (SFR).

We attempted a measure of the spectral energy distributions (SEDs) of the hosts of our SNe in the template images from PS1 (see Sect. B.1). On each image, we used the code SExtractor (Bertin & Arnouts 1996) to obtain the Kron magnitude of the source. For SN 2021adxl, we repeated the measurement twice, first fitting the whole galaxy and then just the H II region. We show the fluxes as a function of wavelength in Fig. 10, the measurements retrieved from NED and, in the case of SN 2021adxl, from Brennan et al. (2024) are also reported. Also, for comparison, we added the SED of the host galaxy of SN 2015bn (Nicholl et al. 2016), a superluminous (SL) Type-I SN, which was chosen because its host has one of the most complete SEDs. From Fig. 10, it appears that all hosts are faint and that their SEDs are quite blue. In particular, there is no significant difference between the SED of the whole host of SN 2021adxl and the H II region alone aside from the total flux. The shape of the SED of our SN hosts matches well that of the host of SN 2015bn. This is interesting because the hosts of SLSNe at low redshift are often blue dwarf galaxies with high star formation (Lunnan et al. 2014; Schulze et al. 2018).

Fig. 10. Spectral energy distribution of the host galaxies. Larger filled points represent our measurements on the PS1 template images, while empty points are from NED (or Brennan et al. 2024, in the case of SN 2021adxl). The fluxes have been corrected for the luminosity distance, and they are shown at the rest-frame wavelength position.

In general, dwarf galaxies in the local Universe produce more stars per unit mass than massive galaxies and, furthermore, they seem to have a top-heavy initial mass function (IMF) that allows them to produce a higher fraction of massive stars than with standard Saltpeter IMF (Dabringhausen et al. 2012; Marks et al. 2012). Interestingly, SNe IIn in general do not tend to follow the star formation as traced by the Hα emission (Habergham et al. 2014), possibly indicating multiple progenitor scenarios (Ransome et al. 2022). However, luminous SNe IIn do occur more frequently in metal-poor environments with young stellar populations (Moriya et al. 2023). This is consistent with the notion that strongly interacting SNe are massive stars from the higher end of the IMF and appear with higher frequency in dwarf galaxies.

2.6. Analysis on the line profile and emission

During the interaction between the SN ejecta and the CSM, four main regions should be considered: from the outside, i) the unshocked CSM, ii) the CSM that was shocked by a forward shock (FS), iii) the SN ejecta shocked by the reverse shock (RS), and iv) the unshocked SN ejecta that are expanding fast (Chevalier & Fransson 1994). If the CSM is dense (≥3 × 10⁻¹⁵ g cm⁻³), a cool dense shell (CDS) forms between the FS and RS (Chevalier & Fransson 1985). An opaque CDS has the same effect of an expanding photosphere (Chugai 2001). In the case of strong interaction, the heated CSM dominates the optical emission. Since emissions from the CSM can remain bright for months or years, it is possible that the ejecta thermal energy fades before the CSM becomes transparent, and thus one may never see the characteristic P Cygni profiles indicative of an expanding photosphere (Smith et al. 2017). In general, in SNe IIn one may expect a narrow (v ∼ 100 km s⁻¹) component due to the unshocked CSM and a broader (v ∼ 5000 − 10 000 km s⁻¹) component, which depending on the optical depth of the CSM can originate either from the ionised SN ejecta or from electron-scattering in the CSM (Huang & Chevalier 2018). Furthermore, if the CSM is asymmetric (Stritzinger et al. 2012), or if dust is present (Fox et al. 2011), they can affect the shape of the line profiles, thus making the theoretical interpretation more difficult.

In Fig. 11, we plot a zoom on the Hα region of spectra taken around 100 days from the explosion for all the SNe in our sample. SNe 2021adxl and 2022qml have broad, blue, asymmetric profiles with a narrow emission centered at rest wavelength. On the other hand, SNe 2021acya and 2022wed have a more symmetric broad line. In the case of SN 2021acya, there is a small dip that may be a narrow P Cygni absorption but it is probably an instrumental artifact, since it is not visible at other epochs.

Fig. 11. Zoom on the Hα region of medium-resolution spectra at ∼100 days. All spectra have been normalised with respect to the Hα and are shown in logarithmic scale for better visualisation.

Given the diversity in shape of the Hα across our sample, it is difficult to directly compare them in terms of FWHM, line position, and total flux. Opting to apply the same method to all of them while being conscious that the procedure will not be optimal, we performed a multi-component fit for the Hα lines, including the contribution of a broad Gaussian and a narrow Lorentzian function. We chose these rather simple profiles because they can ensure a more direct comparison among all the SNe in the sample, even if they might not provide the most perfect fit, especially for complicated profiles such as those of SNe 2021adxl and 2022qml. In some spectra, adding a second, intermediate-width Gaussian component representing the post-shock gas (Smith et al. 2017) allowed for a better fit. This is a relatively faint feature that was ignored for the rest of the analysis. For each component, we fit the position, peak intensity, and width. Given the asymmetric profile of the line in some objects, the position of the peak of the broad component was allowed to vary with respect to the narrow one. All the measurements were performed on reddening-corrected spectra (tables of the Hα and Hβ fits are available online only, see the Data availability section).

For all the SNe in the sample, position, FWHM, and intensity of the narrow component are constant within the errors (calculated by summing in quadrature the root mean square and the resolution). Also, the FWHM of the narrow component is not resolved in the fit with the exception of the spectrum of SN 2022qml at +49 days, which yields a FWHM of 12.9 Å. The average position of the narrow line was taken as the rest-frame reference for each transient.

In the upper left plot of Fig. 12, the rest-frame position of the centre of the broad Hα component is shown. All our SNe show evolution in the rest-frame position. In SN 2021acya, in particular, the position of the line centre is initially at zero velocity with respect to the reference, then rapidly shifts to the blue with a maximum offset of about 2000 km s⁻¹ that later reduces to ∼1000 km s⁻¹ at 300 days. SN 2021adxl appears to show the same late-time behaviour, while for SN 2022wed the variation is less extreme. SN 2022qml is the most extreme, with a broad peak significantly shifted towards the blue, especially at early epochs. This is the effect of the blue shoulder, which has peculiar prominence in this SN. In fact, while the blue shoulder is fairly common, such extreme shifts are more rarely observed (see, e.g., SNe 1997cy, Turatto et al. 2000, and 2007rt, Trundle et al. 2009).

Fig. 12. Results of the multifit on the broad Hα and Hβ for the SNe in the sample. Top left: central wavelength of the Hα line. Top right: full width half maximum of the Hα line. Bottom left: Hα luminosity. Bottom right: ratio of $\frac{\mathrm{H}\alpha }{\mathrm{H}\beta }$ from the values measured for our spectra. The dashed line indicates the division between recombination and collision, which is considered to be 3.1, as in Osterbrock & Ferland (2006) (case B). Data from SN 2010jl taken from Fransson et al. (2014) and from SN 2013L taken from Taddia et al. (2020) are also added given their spectral similarity to the SNe in the sample.

The shift of the broad peak could have different explanations. One possibility is that the red wing of the Hα is obscured by dust, which forms after the explosion but could also be already present in the CSM (Lucy et al. 1989). However, this does not explain why the position of the peak shifts back again to the rest frame. The blueward emission could be due to a mechanism by which the photons passing through the shock front multiple times acquire energy that generates bumps in the bluer part of the emission line (Ishii et al. 2024). Our favoured interpretation, however, is that the shift is due to occultation of the receding line-forming region: if the Hα is forming in a region close to the photosphere just above a CDS, this would imply a deficit in the redward flux, which depends heavily on the density distribution (Dessart & Hillier 2005). This is depicted in model A of Dessart et al. (2016), with a massive SN ejecta ramming into a dense CSM. During the luminous phase of the light curve, the photosphere corresponds to the CDS and the emitting region is moving outward. This last scenario naturally explains the progressive shift to more symmetrical lines as the optical depth decreases and the part of the line-forming region that is receding along our line of sight is revealed. The mechanism of the line formation is hence the same for all the sample of SNe, while the optical depth of the electron scattering region varies. This same argument has been used also to explain the blueshift often observed in emission lines of SNe II (Anderson et al. 2014). However, in this case we do not observe a correlation with the light curve peak luminosity.

In the upper right panel of Fig. 12, the FWHM of the broad component is shown. All SNe in our sample start with very broad (≳80 Å) FWHMs (≳4000 km s⁻¹) that then shrink to ∼40 − 60 Å (2000 − 3000 km s⁻¹) around 300 days. However, SN 2021acya starts with a lower FWHM and then broadens significantly. Considering the gaps in the spectral follow-up, it is possible that all SNe underwent the same evolution. A viable explanation for this behaviour is that at early epochs the photosphere is not (yet) hot enough to show a significant electron scattering effect, which becomes dominant later on.

In the bottom left plot of Fig. 12, the flux of the broad Hα is shown. A progressive increase in the flux is observed in all SNe, with a peak between 200 and 300 days after the explosion. This is followed by a decrease in the case of SNe 2021adxl and 2022wed, while it appears to remain constant (albeit within very large errorbars) for SN 2021acya. A high Hα luminosity is often correlated to strong interaction (Chugai 1991). If this is true also for our SNe, the strength of interaction increases with time and peaks after the light curve peak. The epoch at which the Hα flux starts to decrease is likely related to the extent and density profile of the interacting CSM. A more extended CSM would sustain the interaction for a longer time, while a rapidly declining density would drop the strength of interaction and thus the luminosity. In this respect, the decrease for SN 2022wed happens earlier than for the others, around the maximum of the second peak. Possibly, this denotes a rapid change in the CSM density for this particular object. On the other hand, for SNe 2021acya and 2021adxl the change happens later on, perhaps indicating a less steep decline in the CSM density. Also, we note that the total Hα luminosity is higher for SN 2021acya compared to SN 2021adxl, possibly, due to a higher CSM mass for the former.

Finally, the ratio of $\frac{\mathrm{H}\alpha }{\mathrm{H}\beta }$ for both the narrow and the broad components is calculated. For the narrow line, the ratio is constant within the errors and close to the 3.1 limit value for the case B hydrogen recombination (Osterbrock & Ferland 2006). For the narrow lines, this indicates that H is excited by radiation. The broad lines ratio is shown in the bottom right panel of Fig. 12, where we consider only ratios with an error below 25%. The values are significantly higher than for the narrow lines and the ratio changes over time. This happens when electrons are pushed to higher energy levels through collisions, in regions of high density and optical depth, confirming that all SNe in our sample have a strong shock component in their light emission (Chevalier & Fransson 1994). This is also in agreement with what is found in the literature, for example, for SNe 2010jl and 2013L (Fransson et al. 2014; Taddia et al. 2020).

Brennan et al. (2024) pointed to a possible two-component absorption feature in the He I P Cygni profile of SN 2021adxl, which they explain as either a blend with the Na I D λλ5890, 5896 or the effect of an asymmetric explosion. In Fig. 13, a comparison of the spectrum of SN 2021adxl at +91 days with the one of SN 2022wed at +136 days is shown. The starting position of the blue bump excess matches the one of SN 2022wed and the presence of Fe in this SN is confirmed by Fe IIλ5169 in multiple spectra. This shows that the possible high-velocity He I component mentioned by Brennan et al. (2024) is rather better understood as the emergence of the blue bump due to Fe II emission (cf. Sect. 2.7). Instead, we agree on the interpretation of the lower velocity component as He I at ∼3000 km s⁻¹. This feature is discussed later on.

Fig. 13. Spectral comparison between SNe 2021adxl and 2022wed. All spectra have been redshift corrected, continuum subtracted, and rescaled for better visualisation. We can also clearly see the P Cygni profiles of He Iλ5876 and Fe IIλ5169 for SN 2021adxl.

When the shocked region is exposed, it is possible to measure the shock velocity v_sh from the blue shoulder of the Hα line. The CDS is confined between FS and RS and it is optically thick but radially thin. Photons emitted from such a structure will give a boxy profile to the line (Dessart et al. 2015) and the blueward limit of the emission is a direct measure of v_sh. Even in the case of a composite profile with electron scattering wings, it is still possible to recover the shock velocity as long as the blue shoulder is visible. As shock velocity, we took the Doppler shift measured on the left corner of the blue shoulder compared to the position of the narrow emission. The measurements were performed on the spectra where the boxy component is visible and are shown in Fig. 14, while the highest and lowest values are reported in Table 4. The measurements have a significant scatter due to the uncertainty of the blue shoulder position, but it is still possible to fit them to a declining power-law, which gives v_sh ∝ t^−0.61. The same was done for SN 2022qml, finding v_sh ∝ t^−0.23, more similar to the velocity evolution Taddia et al. (2020) find for SN 2013L. The v_sh measured on the spectra are reported in Table 4, along with other parameters inferred from the light curves and spectra.

Fig. 14. Photospheric expansion velocity of SN 2021adxl from the He I P Cygni absorption (green diamonds) and Fe II (blue squares) along with the shock velocity derived from the blue shoulder of the broad Hα (magenta reverse triangles). The black line represents the best exponential decline fit of the Hα shock velocity based on our measurements.

In Fig. 14, the measurements on the expansion velocity performed on the P Cygni absorption of He I and Fe II for SN 2021adxl are also added. Considering the large errorbars, they have a slightly higher velocity and a similar trend to the shock velocity measured on the Hα blue shoulder at all epochs. The P Cygni profile is indicative of a somewhat extended gas shell above a photosphere and thus of a continuum, while Hα shows no absorption components. This could imply the presence of an inner He layer, on top of which there is an optically thin H layer that allows for the formation of the boxy profile on the Hα line. Our P Cygni measurements, in fact, are comparable with the expansion velocity measured, for example, for SN 2004et (Sahu et al. 2006), which was a normal SN II. A possible scenario is that all H comes from the shocked CSM, while He and Fe are part of the actual SN ejecta. Considering the mass loss (for which we calculate the rate in Sect. 4), it is possible that SN 2021adxl is actually a SE SN (Woosley et al. 1995) and the narrow H line comes from the fraction that is unshocked.

Estimating the shock velocity is very tricky when the line profile does not clearly show a blue shoulder, as it is the case for SNe 2021acya and 2022wed. Following Fransson et al. (2014), we show in Fig. 15 a zoom of the Hα spectral region of SN 2021acya, with the blue side folded over the red one. The lines were shifted so that the broad wings could match, introducing a velocity shift of the narrow peaks. Initially, the shift is low, about 100 km s⁻¹, and grows to a maximum of 500 km s⁻¹ at +119 days. These values are similar to those determined by Fransson et al. (2014), who interpreted the shift as the sign of the bulk velocity of the scattering medium. In this context, the first measurements correspond to the wind velocity and the subsequent acceleration is possibly due to the shock catching up with the wind.

Fig. 15. Zoom on the Hα profile of selected spectra of SN 2021acya. Magenta lines represent the reflected profile. The profiles have a velocity shift of 150 (+24 days), 100 (+61 days), 500 (+112 days), and 300 km s−1 (+347 days).

Because of the higher optical depth in the scattering region, the shocked gas remains hidden in the spectra of these SNe, thus precluding a direct measurement of the shock velocity. Therefore, given the similarities between the spectra of SNe 2021acya and 2022wed with 2010jl, we adopted the same value as Fransson et al. (2014), 3000  km s⁻¹ at +320 days. To estimate its evolution, we took the trend of Eq. (9) of Taddia et al. (2020) and fit it to this point. As a second estimate, we also used the same trend of SN 2021adxl fitted to this point. Finally, we also fit the formula of Taddia et al. (2020) to SNe 2021adxl and 2022qml. This gives us an estimation of how the shock velocity could realistically vary. Both estimates were used for the calculations in Sect. 4.

When the line profile is fully attributed to electron scattering, such as for SNe 2021acya and 2022wed, it is possible to calculate the electron scattering temperature T_e from the FWHM of the broad emission using the equation from Fransson et al. (2014): $\mathrm{FWHM}\approx 647{\tau }_{\mathrm{e}}{(T_{\mathrm{e}}/{10}^4\mathrm{K})}^{\frac{1}{2}}\mathrm{km}{\mathrm{s}}^{-1}$, where τ_e is the optical depth, and inverting it to extrapolate the values of T_e for each spectrum of SNe 2021acya and 2022wed. A symmetric Lorentzian profile implies τ_e ≥ 1, and typical values found in the literature range around 4 – 5 (Chugai 2001). For our calculations, τ_e = 5 was adopted, thus giving temperatures in the range 6.0 × 10³ − 2.6 × 10⁴ K. This is in line with what is derived for SN 2010jl (Fransson et al. 2014). The ranges of T_e are reported in Table 4 and its evolution is shown in Fig. 16, where only the fits with errors < 25% were considered. The T_e evolution is better followed for SN 2021acya. It shows a rapid increase to T_e ∼ 30 000 K with the emergence of the shock followed by a similarly rapid decrease and a long tail at a roughly constant T_e ∼ 8000 K. The few measurements for SN 2022wed show a similar trend but with a higher temperature during the rapidly declining phase. This is coherent with the observed luminosity evolution that suggests a delayed phase of enhanced interaction.

Fig. 16. Electron-scattering temperature for SNe 2021acya and 2022wed as inferred from the broad Hα FWHM, assuming an optical depth of τe = 5.

2.7. Blue bump excess

As we already noticed, a peculiar feature in the spectra of many interacting SNe is a blue bump bluewards from 5500 Å. This was first observed in SN 1988Z (Turatto et al. 1993) and it is typically attributed to a forest of Fe II lines (see, e.g., Foley et al. 2007; Smith et al. 2012a). The relative strength of this feature varies depending on the objects considered and the phases. In general, the spectrum of interacting SNe is a combination of a hot continuum and emission lines, some isolated and others, such as those contributing to the blue bump, blended. To estimate the relative contribution of the two components, we fit a black body (BB) function in the red part on the redshift-corrected spectra (λ5800 Å, onward), avoiding the most prominent emission lines (in particular, Hα). However, we caution that the BB fits fails after roughly 200 days (100 days for SN 2022qml). An example of the fit region and the estimated BB is shown in Fig. 17. Here, the black line is the observed spectrum of SN 2021acya at +61 days, while the red bands show the location of the fitting region and the green line the fitted BB function. A blue line representing the BB fit to the whole spectrum is also added to show that the fit is less optimal in this case due to the blue bump. The blue bump excess is defined as the ratio between the difference in flux between the observed flux in the blue (avoiding the most prominent emission lines, in particular, Hβ) and the flux below the BB fitting on the redder part of the spectrum at the same wavelengths, all divided by the total BB flux.

Fig. 17. Spectrum of SN 2021acya at +61 days (black). The green dashed line is the estimated BB flux from the fit on the red part of the spectrum, while the blue dotted line is the estimated BB flux from the fit on the whole spectrum. We note that the blue BB does not fit the spectrum between 5000 and 6000 Å well.

Considering the gaps in the observations and the different phases of the SNe in our sample, a mean blue bump excess was computed every 100 days since the explosion to aid in the comparison. It is shown, compared to the mean luminosity at the same phases derived from the g, r, i pseudo-bolometric light curve (see Sect. 3 for details on the computation), in Fig. 18, where three reference objects are also added to represent the effect of different amounts of interaction: the strongly interacting SN 2010jl, the mildly interacting SN 1998S, and the Type IIP SN 2004et, where moderate interaction only started after ∼800 days (Kotak et al. 2009). The evolution with time is significant. The excess is small (in some cases, even negative, meaning that the fitted BB exceeds the measured flux, probably due to line blanketing) at early phases, and tends to grow after roughly 100 days, while the luminosity decreases. SNe 2010jl and 2021adxl seem to show milder blue bump excesses at early phases than SNe 2021acya and 2022wed, while SN 2022qml is the most extreme, with a stark increase in the blue bump excess after 100 days. SN 1998S exhibits a similar behaviour, albeit with lower absolute values of the blue bump excess. SN 2004et, on the other hand, shows an always negative excess. This is in line with what we would expect from a non-interacting SN⁴.

Fig. 18. Blue bump excess as a function of the g, r, i luminosity. Numbers close to the points indicate the average phase.

The Hα luminosity can be used as a tracer of the strength interaction (Chugai 1991) to verify whether it correlates with the strength of the blue bump excess; however, there appears to be none. The time at which the interaction component is dominant with respect to the BB continuum is probably due to a combination of CSM density, strength of interaction, and possibly the relative Fe abundance in the ejecta with respect to the other elements. A detailed spectral modelling may help disentangle the origin of the feature in different objects, but it appears that the presence of the blue bump excess alone is not indicative of strong interaction. A combination of parameters such as high bolometric and Hα luminosity is also required. However, there seems to be a mild correlation with the CSM optical depth that works as follows: interacting SNe with a strong blue bump excess also show a boxy Hα line profile. In turn, the boxy profile is exposed when the electron scattering optical depth of the heated CSM is smaller.

3. Bolometric light curves

The SNe in our sample have different coverage in terms of wavelength and phases. To ensure a more objective comparison among them, we computed a pseudo-bolometric light curve for all SNe integrating the flux in the g,r,i filters. The magnitudes are corrected for extinction and converted to flux densities using photometric zero points⁵ The flux is then integrated in the sampled spectral region using a trapezoidal rule and assuming zero flux below/above the limit defined by the filter equivalent width of the bluer/redder filter, respectively. The measured flux is converted into luminosity using the adopted distance modulus. The results for our four SNe is shown in Fig. 19, where the pseudo-bolometric light curves of SNe 2004et and 2010jl, calculated in the same way, are also included for comparison. A logarithmic scale is used to emphasise the change of slope between the early and late phases.

Fig. 19. Pseudo-bolometric light curves computed using only optical filters. We also include the bolometric light curves, calculated in the same way, of SNe 2010jl and 2004et for comparison.

SN 2021acya is the most energetic SN in our sample, with an extremely high E₂₀₀ (Table 1). After the very bright peak, SN 2021acya shows a slow decline followed by a plateau (better seen in Fig. 1) that lasts from +150 days until +300 days. The decline rate both before and after the plateau is similar to that of SN 2010jl. SN 2021adxl shares some similarities with SN 2010jl, in particular, the break in the light curve around +400 days. However, the decline rate is much higher and, in fact, SN 2021adxl is about 1.5 mag fainter than SN 2010jl at ∼500 days. SN 2022qml has the shortest light curve, since it is observed for less than 200 days, after which the SN is lost because of conjunction with the Sun. When it re-emerged from behind the Sun, the SN was no longer detected. The decline rate is also similar to that of SN 2010jl. Finally, SN 2022wed has the most peculiar light curve shape, with its two peaks, the second much brighter than the first. This is probably due to a CSM arranged in two shells with different mass and density. It also does not show a break at 400 days; however, the decline rate after the second peak is almost identical to that of SN 2021acya after the plateau.

All the SNe in our sample are bright and long lasting like SN 2010jl and in some cases also the decline rate is similar. The prototype of Type II plateau SN 2004et, on the other hand, is completely different, since it is 1 − 2 orders of magnitude fainter at all phases and also its plateau phase is considerably shorter. This gives a measure of the significant energy contribution from interaction for the selected SNe.

Along with the direct comparison, we tried to gauge the true bolometric luminosity taking into account UV and NIR bands. In this respect, SN 2021acya has the best photometric coverage. To estimate the contribution of UV and NIR photometry, in Fig. 20 the bolometric light curve computed with the contribution from the UV or NIR flux is compared to the one computed using only g, r, i. The contributions are calculated as the ratio between the UV or NIR flux and the total flux. While the NIR contribution at early phases is small (around 20%) and rises up to 40% at later times, the UV contribution is significant, up to 60%, but then it decreases rapidly. A similar UV excess was also measured in SN 2010jl and other interacting SNe (Fransson et al. 2014). In SN 2010jl, the early NIR light curves followed the decline of the optical ones, while after ∼200 days they flattened in J and H and even increased the flux in K. This was attributed to pre-existing dust in the CSM that was re-heated by the optical emission (Fransson et al. 2014). For SN 2010jl, the NIR bands contribute 20–50% of the bolometric flux in the phase range 200 − 600 days. This is in line with the NIR contribution that we measure 100 − 200 days after the maximum for SN 2021acya.

Fig. 20. Ultraviolet and NIR contribution for SN 2021acya. The upper panel shows the differences between the luminosities computed using UV and NIR observations with respect to the bolometric light curve computed taking into consideration only optical bands.

The same check was also performed on SNe 2021adxl (for which we have a sequence of u,z,J,H,K observations), 2022qml (for which we have a sequence of u observations), and 2022wed (u observations). The u band contribution is consistent in all SNe, adding alone between ∼5 and 10% of the total flux at early phases. The NIR observations of SN 2021adxl, on the other hand, behave like in SN 2021acya, giving a flux ∼25 − 35% higher. Also, the NIR contribution appears to increase with time.

The considerable flux difference that is found when adding the NIR contribution is in line also with Martinez et al. (2022), where they show that NIR observations are fundamental to properly reconstruct the bolometric light curve of SNe II. To compute the total luminosity, which we use in the next section, the contribution of UV and NIR bands, when available, was added to the optical ones.

In Table 3 the measured decline rates in the log-log scale for SNe 2021acya, 2021adxl, and 2022qml are reported (the measurement was not performed on SN 2022wed since its light curve shape is very different from the others). Fransson et al. (2014) interpret the break as the time when the FS emerges at the outer edge of the dense CSM, while the specific slope is a sign of the CSM density. In particular, a steeper slope is indicative of a smaller density, since the radiation diffuses earlier.

In Khatami & Kasen (2024), the authors explore the diversity of SN light curves when interaction with the CSM is involved. The key parameter in their model is the break-out parameter ξ:

$$\begin{array}{c}\hfill \xi ={\tau }_0{\beta }_0{\eta }^{-\alpha }\sim t_{\mathrm{esc}}/t_{\mathrm{sh}},\end{array}$$(1)

where t_esc and t_sh are, respectively, the time scale for photons to escape ahead of the shock and the dynamical time scale of the shock, η = M_CSM/M_ej is the ratio of CSM to ejecta mass, with M_CSM and M_ej the masses of CSM and ejecta, respectively, β₀ = v_ej/c is the ejecta velocity with respect to the speed of light, ${\tau }_0=\frac{\kappa M_{\mathrm{CSM}}}{4\pi R_{\mathrm{CSM}}^2}$ is the characteristic optical depth of the CSM, with R_CSM the radius of the CSM and κ the opacity, and α = 1/2 if η > 1, while α = 1/(n − 3) if η < 1, where n is the power-law exponent of the density profile of the ejecta. Depending on the values of ξ and η, Khatami & Kasen (2024) identify four different scenarios: i) edge-breakout with light CSM (ξ ≫ 1, η ≪ 1), ii) edge-breakout with heavy CSM (ξ ≫ 1, η ≳ 1), iii) interior-breakout with light CSM (ξ ≲ 1, η ≪ 1), and iv) interior-breakout with heavy CSM (ξ ≲ 1, η ≳ 1). Different combinations of ξ and η generate different light curves (Khatami & Kasen 2024, their Fig. 3). Judging by the light curve shape of our targets, they are well represented by case iv, with a CSM comparable to or even more massive than the ejecta mass colliding with it and a dynamical time scale comparable to or longer than the escape time (although SN 2022wed, with its long second peak, could also be part of case ii, with an escape time scale longer than the dynamical one).

4. Mass-loss rate and mass of the circumstellar medium

For the four SNe in our sample, the high luminosity, slow decline and narrow emission lines in the spectra are all indicative of CSM/ejecta interaction as source of the luminosity. In fact, their luminosity is at all phases 1 − 2 orders of magnitude brighter than typical CC SNe. In turn, if the luminosity is attributed to radioactive material, it would require a ⁵⁶Ni mass of the order of 1 M_⊙ or more.

In the contest of CSM/ejecta interaction, it is possible to have an estimate of the mass loss following the approach of Ishii et al. (2024), in particular their Eq. (15), which is derived in Kokubo et al. (2019). The formula is

$$\begin{array}{cc}\hfill \dot{M}& =0.04{{\mathrm{M}}_{\odot }\mathrm{yr}}^{-1}\left(\frac{L_{\mathrm{tot}}}{8\times {10}^{43}{\mathrm{erg}\mathrm{s}}^{-1}}\right)\left(\frac{v_{\mathrm{w}}}{40{\mathrm{km}\mathrm{s}}^{-1}}\right)\hfill \\ \hfill & ·{\left(\frac{ϵ}{0.3}\right)}^{-1}{\left(\frac{v_{\mathrm{sh}}}{4400{\mathrm{km}\mathrm{s}}^{-1}}\right)}^{-3},\hfill \end{array}$$(2)

where L_tot is the total bolometric luminosity, v_sh is the shock velocity, ϵ is the radiation conversion efficiency, which is assumed to be 0.3 as in Kokubo et al. (2019), and v_w is the wind velocity of the CSM. The wind velocity of SN 2022qml is inferred from the position of the narrow Hα P Cygni absorption in the high-resolution spectrum at + 51 days, which gives 100 km s⁻¹. For the other SNe the resolution is too poor to properly identify this feature but an upper limit is measured from the FWHM of the narrow line, which is reported in Table 4. Upper limits are consistent with the value adopted as reference by Ishii et al. (2024), that is ≤47.2 km s⁻¹. We remind that LBVs have been proposed as progenitors of strongly interacting SNe IIn. On the one hand, LBVs have faster winds, of the order of 200 km s⁻¹ (Humphreys et al. 1988) and on the other hand they are known to emit a huge amount of mass through short eruptive episodes (e.g. Eta Carinae has a calculated mass-loss rate of 0.075 M_⊙ yr⁻¹ Andriesse & Viotti 1979). For the sake of the calculation, we assumed v_w = 100 km s⁻¹ for all SNe, while the shock velocity v_sh comes from the estimates in Sect. 2.6.

From the shock velocity and the duration of the interaction phase, it is possible to calculate the radius of the CSM and, from the wind velocity, the duration of the mass-loss phase. Then, the mass of the shocked CSM M_CSM is derived by integrating the mass-loss rate and considering the duration of the interaction phase as seen from the light curve. These parameters for each SN are reported in Table 5, along with the energy calculated from the integration of the total luminosity, and the ranges are calculated using one or the other v_sh estimate described in Sect. 2.6. Several factors can affect our computation: as mentioned in Sect. 3, the bolometric luminosity is likely underestimated because of the of lack UV and NIR observations for most of the SNe in the sample. This means that the mass-loss rates and hence the CSM masses are lower limits, and the SNe could have shed even more mass. On the other hand, a steady wind velocity of 100 km s⁻¹ was assumed. This value was measured only on the spectrum of SN 2022qml, and it may be different for the other SNe. However, given the upper limit from the narrow FWHM reported Table 4, the actual value for the other SNe is likely of the same order of magnitude. The conversion efficiency is also an unknown factor, but it cannot affect the estimate of the mass loss by more than a factor of two. The parameter that most affects the calculation is the shock velocity, both because of the difficulty to estimate it and because of its weight in the equation. While we have a good estimate of the shock velocity for SNe 2021adxl and 2022qml, the evolution of v_sh for SNe 2021acya and 2022wed is a mystery due to their higher opacity. However, it cannot be much higher than our estimate, otherwise a blue shoulder would have been visible in the spectra. If, on the other hand, the values were lower than our estimate, this would imply an even larger mass-loss rate than what was derived.

Even with all the caveats due to the assumptions described above, the values of Ṁ and M_CSM are much higher than those expected for a steady wind and more in line with eruptive processes (Matsumoto & Metzger 2022), in agreement with an LBV progenitor (for reference, Fassia et al. 2001 find for SN 1998S a mass-loss rate of 2 × 10⁻⁵ M_⊙ yr⁻¹, two orders of magnitude lower than our values, which are instead more in line with the eruptive episodes in Eta Carinae). To produce a CSM mass of 5 − 20 M_⊙ requires a huge progenitor, the formation of which could be challenging for stellar evolution models. The only SN in our sample with a considerably smaller M_CSM is SN 2022qml. This, together with the different spectra and shorter light curve duration, points towards a possible thermonuclear explosion whose progenitor was embedded in a H-rich CSM.

It is also interesting to explore the mass-loss rate as a function of time. Regardless of the CSM origin (wind or eruption), we assume that it is distributed in a shell within which the velocity is constant. Assuming a steady wind for simplicity, it is possible to derive the value of Ṁ on the time scale of the wind. The results for each SN are shown in Fig. 21. Here, there is a strong dependence especially on the ratio between the shock and the wind velocity, which gives the overall shape to the trend, as well as the duration of the bolometric light curve. As shock velocity, we considered the trend derived on the Hα blue shoulder fit for SNe 2021adxl and 2022qml, while for SNe 2021acya and 2022wed we take the same trend of SN 2021adxl rescaled to the velocity derived by Fransson et al. (2014) at 320 days for SN 2010jl, as explained in Sect. 2.6. From Fig. 21, it appears that the mass-loss rate of SN 2022qml was more constant, compared to the other SNe, and overall lower, once again in agreement with a different origin for this SN. The trends of SNe 2021adxl and 2022wed are similar but shifted, with SN 2022wed starting to significantly increase later on, around 20 years after the onset, while SN 2021adxl shows a high mass-loss rate already ten years in. SN 2021acya, on the other hand, has a peculiar behaviour with possibly three distinct peaks that could indicate the onset of different mass-loss episodes. Interestingly, the timeline of the peaks seems to align with the other SNe: the second one happens around the same time as the one of SN 2021adxl and the third is almost in coincidence with the peak of SN 2022wed. With these considerations, it appears even more evident that the progenitor stars of these SNe had to be massive, peculiar objects to be able to expel such a huge amount of mass within such a limited time period, possibly, in the case of SN 2021acya, in two separate episodes. A mass-loss rate that starts high and increases for decades after the onset of the episode is also in line with what is found by Moriya et al. (2014), although it may decline closer to the explosion.

Fig. 21. Mass-loss rate as a function of the time from the onset of the mass-loss episode, assuming a steady wind and shock velocity as derived in Sect. 2.6.

In Table 4, the total radiated energy obtained by integrating the bolometric light curve is reported. It can be converted to the kinetic energy E_K produced in the explosion by assuming a conversion factor ϵ, which is not well constrained and can vary from 0.1 to 1 (Moriya & Maeda 2014). As before, we assumed ϵ = 0.3 as in Kokubo et al. (2019), while the result is reported in Table 5.

The CSM radius, its mass, and the kinetic energy are useful to compare our findings to models of neutrino production in shocked regions of interacting SNe. In particular, the semi-analytical model in Pitik et al. (2023) shows that the largest probability to produce HE neutrinos is found for SNe with E_K ≳ 10⁵¹ erg, 1 ≲ M_CSM ≲ 30 M_⊙, and R_CSM ≳ 10¹⁶ cm. These parameters are favoured in SNe that display luminous light curves (L_peak ≳ 10⁴³ − 10⁴⁴ erg s⁻¹) but average rise time (10 d ≲ t_rise ≲ 40 d). All our SNe fulfill these criteria but for SN 2022qml (see Table 4), whose total energy is smaller than the rest due to the shorter duration of the light curve, as well as its overall dimness compared to the others.

The parameters of the SNe in our sample are plotted over the contour plots of the integrated neutrino energy Pitik et al. (2023) calculated for their model in the upper panel in Fig. 22, where it is assumed that a fraction of 10% of the kinetic energy of the shock is used to accelerate protons. The plot shows the integrated neutrino energy ${\epsilon }_{\nu +\overline{\nu }}$ for energy above 1 TeV as a function of M_CSM and R_CSM at fixed kinetic energy (E_K = 10⁵¹ erg) and ejecta mass (M_ej = 10 M_⊙). To our measurements are also added those of SN 2010jl (Fransson et al. 2014), SN 2013L (Taddia et al. 2020), and a prototypical Type Ia-CSM whose parameters were inferred from an average on the sample presented in Sharma et al. (2023). While SNe 2022qml and 2013L are outside the calculated range of the model (which could be due to a different progenitor for SN 2022qml), SNe 2021acya and 2022wed have a CSM that is too massive, while the others fall within the limits for HE neutrino production. However, we caution that there is some degeneracy in the parameters. The same exercise is repeated in the bottom panel of Fig. 22, but this time plotting the kinetic energy with respect to the CSM mass at fixed ejecta mass (M_ej = 10 M_⊙) and CSM radius (R_CSM = 10¹⁶ cm). In this case, SNe 2021acya and 2022wed are within the model margins and actually SN 2021acya is the most favoured. Moreover, we remind that the integrated luminosity is a lower limit, since there is no full bolometric coverage at all epochs. In principle, this could push some SNe upwards in the plot.

Fig. 22. Total neutrino energy contour plots. Top: contour plots of the total neutrino energy ${\epsilon }_{\nu +\overline{\nu }}$ integrated for Eν ≥ 1 TeV from Pitik et al. (2023) in function of the CSM radius and mass, over which the same parameters for our SNe are plotted. The orange star indicates the benchmark of the model. Bottom: same as above but in function of kinetic energy and CSM mass.

This test shows that the SNe in our sample indeed have the characteristics to produce HE neutrinos. We should stress that the parameters have some degeneracy and, in fact, multiple SN models could lead to the same neutrino flux. For example, a higher mass of the CSM disfavours neutrino production, but it can be compensated by a higher kinetic energy (Pitik et al. 2023).

For reference, we can apply the same analysis also to SN 2020faa, a SN powered by hidden interaction, which has a derived kinetic energy and CSM mass and radius: E_K = 4.8 × 10⁵¹ erg, M_CSM = 1 M_⊙, and R_CSM = 10¹⁴ − 10¹⁵ cm (Salmaso et al. 2023). The small CSM radius disfavours SN 2020faa as a likely HE neutrino source; however, the kinetic energy and the CSM mass are in agreement with what is found for the other SNe in our sample. As mentioned, degeneracy in the parameters implies uncertainties that are difficult to disentangle but it is clear that high values of CSM mass and kinetic energy, hence requiring a massive progenitor, seem to be needed.

5. Conclusions

This paper shows that strongly interacting SNe can display a variety of light curves and spectra. These differences mainly depend on the progenitor and its pre-explosion history, since this will determine the mass, radius, density, and distribution of the CSM, apart from the final mass at the explosion. We also inferred some interesting clues from a comparison among the events and with prototypical events from the literature. Our analysis of similarities and differences shows the range of diversity in the density and profile of the CSM and in the properties of the nuclear engine. In particular, there is circumstantial evidence that SN 2022qml could hide a Type Ia SN explosion rather than a CC, which is the case for the other SNe.

The hosts of these SNe are dwarf galaxies consistent with high SFR. Top-heavy IMFs, which favour the production of massive stars, could cause these strongly interacting SNe from massive progenitors to be more frequent in smaller galaxies. This would be important for the chemical enrichment of such environments. Also, in a multimessenger context, massive SN progenitors are interesting because the compact remnant after the explosion could be a massive black hole in the range of those observed by LIGO-Virgo (Abbott et al. 2020).

Despite their differences, all the SNe analysed here are consistent with energetic explosions and strong interaction with the surrounding CSM. The CSM masses are high (extremely so in the case of SNe 2021acya and 2022wed), and the mass-loss rates are similar, around 30 years before the explosion. Instead, the trend of the mass-loss rate differs, since during this time SN 2022qml had a more constant mass loss, while SNe 2021adxl and 2022wed underwent a significant acceleration, and SN 2021acya possibly even had three separate mass-loss episodes.

The energies and CSM mass and radius derived from our analysis are also mostly consistent with HE neutrino production. Murase (2023) shows that for SNe IIn with parameters similar to those analysed here, the ideal time scale of neutrino detection is between 10⁶ and 10⁸ s, that is, between ten days and three years, which is in line with other models (Sarmah et al. 2022). Eventually, these transients appear to be located in regions of the parameter space where the production of HE neutrinos is marginal but not excluded.

Data availability

Tables with photometric measurements, spectral log of observations, and fit on the Hα and Hβ are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/695/A29. Spectra are available to download in WISeREP at the following links: https://www.wiserep.org/object/19799 (SN 2021acya), https://www.wiserep.org/object/20180 (SN 2021adxl), https://www.wiserep.org/object/21394 (SN 2022qml), https://www.wiserep.org/object/22674 (SN 2022wed).

Acknowledgments

We thank the anonymous referee for their helpful insights. The author acknowledges the help of Dr. T. Pitik in recreating Figure 22. This paper is supported by fundings from MIUR, PRIN 2017 (grant 20179ZF5KS). IS, AR, and NER are partially supported by the PRIN-INAF 2022 project “Shedding light on the nature of gap transients: from the observations to the models”. AR acknowledges the support GRAWITA Large Program Grant (PI P. D’Avanzo). YZC is supported by the National Natural Science Foundation of China (NSFC, Grant No. 12303054), the National Key Research and Development Program of China (Grant No. 2024YFA1611603), the Yunnan Fundamental Research Projects (Grant No. 202401AU070063), and the International Centre of Supernovae, Yunnan Key Laboratory (No. 202302AN360001). JPA was funded by ANID, Millennium Science Initiative, ICN12_009. TP acknowledges the support by ANID through the Beca Doctorado Nacional 202221222222. MN 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 UK Space Agency Grant No. ST/Y000692/1. LG, CPG and TEMB 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”, the European Union Next Generation EU/PRTR funds, the Horizon 2020 Research and Innovation Programme of the European Union, and by the Secretary of Universities and Research (Government of Catalonia), under the PID2020-115253GA-I00 HOSTFLOWS project, 2021-SGR-01270, the 2019 Ramón y Cajal program RYC2019-027683-I, the 2021 Juan de la Cierva program FJC2021-047124-I, the Marie Skłodowska-Curie and the Beatriu de Pinós 2021 BP 00168 programme, 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. TWC acknowledges the Yushan Young Fellow Program by the Ministry of Education, Taiwan for the financial support. MB, EC and TP acknowledge the financial support from the Slovenian Research Agency (grants I0-0033, P1-0031, J1-8136, J1-2460 and Z1-1853) and the Young Researchers program. PC acknowledges support via Research Council of Finland (grant 340613). MS is funded by the Independent Research Fund Denmark (IRFD, grant number 10.46540/2032-00022B). Based on observations collected at Copernico 1.82m telescope and Schmidt 67/92 telescope (Asiago Mount Ekar, Italy) INAF – Osservatorio Astronomico di Padova. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 101004719. This material reflects only the authors views and the Commission is not liable for any use that may be made of the information contained therein. 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. Observations from the Nordic Optical Telescope were obtained through the NUTS2 collaboration, which are supported in part by the Instrument Centre for Danish Astrophysics (IDA). The data presented here were obtained in part with ALFOSC, which is provided by the Instituto de Astrofisica de Andalucia (IAA). This work makes use of data from the Las Cumbres Observatory network. The LCO team is supported by NSF grants AST-1911225 and AST-1911151, and NASA SWIFT grant 80NSSC19K1639. The Liverpool Telescope is operated on the island of La Palma by Liverpool John Moores University in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias with financial support from the UK Science and Technology Facilities Council. This work was based in part on observations made with the Italian Telescopio Nazionale Galileo (TNG), operated on the island of La Palma by the Fundación Galileo Galilei of the INAF (Istituto Nazionale di Astrofisica) at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias. The Gran Telescopio CANARIAS (GTC) is a 10.4 m telescope with a segmented primary mirror. It is located in one of the top astronomical sites in the Northern Hemisphere: the Observatorio del Roque de los Muchachos (ORM, La Palma, Canary Islands). The GTC is a Spanish initiative led by the Instituto de Astrofísica de Canarias (IAC). The project is actively supported by the Spanish Government and the Local Government from the Canary Islands through the European Funds for Regional Development (FEDER) provided by the European Union. The project also includes the participation of Mexico (Instituto de Astronomía de la Universidad Nacional Autónoma de México (IA-UNAM) and Instituto Nacional de Astrofísica, Óptica y Electrónica (INAOE)), and the US University of Florida. Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programmes 1103.D-0328 and 105.20TF.001. We acknowledge the use of public data from the Swift data archive. Based 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 National Science Foundation under Grant No. AST-2034437 and a collaboration including Caltech, IPAC, the Weizmann Institute for Science, the Oskar Klein Center at Stockholm University, the University of Maryland, Deutsches Elektronen-Synchrotron and Humboldt University, the TANGO Consortium of Taiwan, the University of Wisconsin at Milwaukee, Trinity College Dublin, Lawrence Livermore National Laboratories, and IN2P3, France. Operations are conducted by COO, IPAC, and UW. The Pan-STARRS1 Surveys (PS1) and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation Grant No. AST-1238877, the University of Maryland, Eotvos Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation. This work presents results from the European Space Agency (ESA) space mission Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia MultiLateral Agreement (MLA). The Gaia mission website is https://www.cosmos.esa.int/gaia. The Gaia archive website is https://archives.esac.esa.int/gaia. This work has made use of data from the Asteroid Terrestrial-impact Last Alert System (ATLAS) project. The Asteroid Terrestrial-impact Last Alert System (ATLAS) project is primarily funded to search for near earth asteroids through NASA grants NN12AR55G, 80NSSC18K0284, and 80NSSC18K1575; byproducts of the NEO search include images and catalogues from the survey area. This work was partially funded by Kepler/K2 grant J1944/80NSSC19K0112 and HST GO-15889, and STFC grants ST/T000198/1 and ST/S006109/1. The ATLAS science products have been made possible through the contributions of the University of Hawaii Institute for Astronomy, the Queen’s University Belfast, the Space Telescope Science Institute, the South African Astronomical Observatory, and The Millennium Institute of Astrophysics (MAS), Chile.

References

Appendix A: Observations

Our targets were monitored using the Schmidt and Copernico telescopes of the Asiago Observatory⁶, INAF (Italy); the Nordic Optical Telescope (NOT)⁷, the Telescopio Nazionale Galileo (TNG)⁸, the Liverpool Telescope (LT)⁹, and the Gran Telescopio CANARIAS (GTC)¹⁰, all located in La Palma (Spain); the Rapid Eye Mount (REM)¹¹ of INAF at La Silla and the Very Large Telescope (VLT)¹² of ESO at Paranal (Chile). We also exploited the time allocated to other facilities to international collaborations, namely the NOT via NUTS2 (Nordic-optical-telescope Un-biased Transient Survey)¹³ and the ESO NTT¹⁴ via ePESSTO+ (advanced Public ESO Spectroscopic Survey for Transient Objects, Smartt et al. 2015). These observations were complemented with data from transient surveys available from public archives, in particular, the Zwicky Transient Facility (ZTF, Bellm et al. 2019), the Las Cumbres Observatory (LCO, Brown et al. 2013), the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS, PS1, Chambers et al. 2019), the Asteroid Terrestrial-impact Last Alert System (ATLAS, Tonry et al. 2018). Space observations obtained by the Gaia mission (Gaia Collaboration 2018), and from NASA’s Neil Gehrels Swift Observatory (SWIFT, Roming et al. 2005) were also retrieved.

Appendix B: Data reduction

B.1. Photometry

Photometric observations were reduced with standard techniques using IRAF¹⁵ recipes and a number of different python packages, in particular astropy (Astropy Collaboration 2022) and affiliated packages (astroquery, ccdproc, photutils). At first, the detector signature was removed with bias and flat-field corrections. Then, an algorithm for cosmic-ray rejection¹⁶ was applied. For astrometric and photometric calibration and for the measurement of the SN magnitudes the required recipes in the ecsnoopy package¹⁷ were implemented. In most cases (when the SN magnitude was below ∼16 − 17 mag), the SN magnitude was measured after subtracting a template image from the observed frame, to better remove the contamination from the host galaxy. To this aim, reference archival images from public surveys such as PS1, SDSS (Sloan Digital Sky Survey, Kollmeier et al. 2019), or Skymapper (Keller et al. 2007) were used. With ecsnoopy, the registration of the template image was secured to the same pixel grid of the science image and then the code hotpants (Becker 2015) was used for the convolution of the two images to the same PSF and photometric scale. Finally, the instrumental magnitudes (or upper limits) were calibrated using photometric zero points measured from local stars with photometry retrieved from the public surveys mentioned above.

For the SWIFT observations of SN 2021acya, the magnitudes were measured with aperture photometry using uvotsource from heasoft (Nasa High Energy Astrophysics Science Archive Research Center (Heasarc) 2014) adopting a circular radius of 5 arcsec and subtracting the sky background measured in an offset empty region.

Finally, for the near-infrared (NIR) observations, the same procedure as for the optical observations was used but including the preliminary subtraction of the sky background, which was obtained from the median combination of the dithered images for each filter. In this case, the nightly zero-points were computed with respect to the Two Micron All-Sky Survey (2MASS, Skrutskie et al. 2003) source catalogue photometry. The magnitudes of the targets are available online only (see the Data availability section), together with the magnitudes acquired from the public surveys mentioned above.

The available pre-supernova imaging was also checked for possible eruptive episodes in our SNe. In particular, through the available web tool,¹⁸ ATLAS forced photometry at the SN locations was obtained. In all cases, we concluded that there is no evidence of variation in the observed flux in the last decade. However, we stress that the limiting magnitude of ATLAS (20.0 mag in cyan and 19.5 mag in orange), corresponds to an absolute magnitude that ranges from −14.7 to −18.8 mag for the distance of our SNe, so the limit is not stringent and it is possible that any low-intensity precursor activity went undetected.

B.2. Spectroscopy

The spectroscopic observations were in general reduced using standard prescriptions with the package foscgui¹⁹. The spectra were corrected for bias, flat-field, and cosmics rejection, calibrated in wavelength in the 2D frame, and extracted to obtain the 1D spectrum. This was then calibrated in flux and corrected for second-order contamination (if required) and telluric absorptions, the latter with the aid of the spectrum of a hot spectrophotometric standard star. For the spectra taken using X-shooter, the publicly available data reduction pipeline EsoReflex²⁰ was used, following the same steps described above. The logs of spectroscopic observations are available online only (see the Data availability section). The spectral evolution of the targets through the most significant spectra is shown in Figs. 4, 5, 7, and 8.

B.3. Extinction and redshift corrections

Galactic extinction values were obtained from NED²¹, assuming R_V = 3.1 (Schlafly & Finkbeiner 2011). Moreover, we examined the spectra of all transients to search for evidence of Na I D λλ5890, 5896 interstellar gas in the host galaxy. This was only identified in the X-shooter spectrum of SN 2022qml at +49 days. The lines have an equivalent width (EW) of 0.1168 Å and 0.1198 Å for D1 and D2, respectively. Given the line EW vs. extinction relation in Poznanski et al. (2012), these values correspond to a total reddening E(B − V) = 0.028 ± 0.2 mag, where the error comes from the uncertainties in the relation, and give an absorption A_V = 0.09 ± 0.6 if R_V = 3.1 also in this case, which is much lower than the extinction due to the MW (see Table 1). However, Phillips et al. (2013) found that the use of this method is not reliable for dust reddening estimation in SNe Ia, while Rodríguez et al. (2023) argue that it underestimates the dust contribution in SNe II. Another estimate of intrinsic reddening can be obtained through the relation proposed in Turatto et al. (2003), which gives E(B − V) = 0.038 mag and is still small compared to the MW contribution. Therefore, we conclude that in all cases the extinction inside the host galaxy is likely negligible.

The presence of strong emission lines combined with a significant redshift for some SNe in the sample suggested checking for the effect of the K-correction on the photometry. To calculate the magnitude of the K-correction, we examined all the spectra and extrapolate the flux at g, r, i bands. The measurement was performed twice, once on the original spectra and once on the redshift-corrected spectra. The difference is the K-correction that was then linearly interpolated and applied to the light curve at all epochs. The magnitude of the K-correction is shown in Fig. B.1, where the location of each transient at different phases both with and without K-correction (empty and filled symbols, respectively) is plotted in a colour-colour diagram (r − i vs. g − r). The strongly interacting SN 2010jl (Fransson et al. 2014) is also added as reference. The shift after the K-correction appears more significant for SN 2022wed, as expected being the transient at the highest redshift, and also for SN 2021acya. In the case of SN 2022qml, the transient does not show a significant evolution and is always in agreement with an extremely hot black body (BB). This is probably due to the blue bump excess, a feature that is discussed in Sect. 2.7. On the other hand, SNe 2010jl and 2021adxl are in countertrend, showing redder colours at all phases. This is interpreted as the effect of a luminous Hα emission line along with a small blue bump excess compared to the other SNe in our sample, although a contribution from higher host extinction cannot be completely ruled out. The K-correction to all the SNe in our sample for the rest of the analysis.

Fig. B.1. Colour-colour r − i versus g − r diagram of our SNe with and without the K correction. The extinction vector E(B-V) and a BB with temperature between 25000 and 5000 K are also added.

The redshift is measured from the narrow Hα emission line. To calculate the distance modulus for the absolute magnitude computation it is assumed H₀ = 70 km s⁻¹ Mpc⁻¹, while the measured redshifts are corrected to the V3K reference frame.

All Tables

All Figures

	Fig. 1. Pseudo-bolometric g, r, i light curves of all the SNe we followed with our programme. Empty reverse triangles indicate the SNe that, after further analysis, were excluded from the final sample. All the phases are corrected for time dilation (this is true throughout the paper).
In the text

	Fig. 2. Integrated luminosity in the first 200 days for each SN in our sample and some objects from the literature for comparison. The dashed vertical line divides the strongly interacting sub-sample from the more mildly interacting SNe.
In the text

	Fig. 3. Absolute magnitude r-band light curves for the SNe IIn in our sample as well as the well-studied Type IIn SN 2010jl (Smith et al. 2012b; Fransson et al. 2014) for comparison. Black arrows indicate upper limits. All photometry is K-corrected and the phases are corrected for time dilation.
In the text

	Fig. 4. Spectral evolution of SN 2021acya as shown through the most significant spectra. Numbers close to the spectra indicate the phase from explosion corrected for time dilation. All spectra are scaled with respect to the Hα and arbitrarily shifted for better visualisation. The spectra are not corrected for extinction.
In the text

	Fig. 5. Spectral evolution of SN 2021adxl as shown through the most significant spectra. Numbers close to the spectra indicate the phase from explosion corrected for time dilation. All spectra are scaled with respect to the Hα and arbitrarily shifted for better visualisation. The spectra are not corrected for extinction.
In the text

	Fig. 6. Spectral comparison between SNe 2021adxl, 2022qml, and 2013L from Taddia et al. (2020) (zoom on the Hα line). The red dashed line highlights the boxy region of the line profile due to the shock. All spectra have been redshift corrected, continuum subtracted, and rescaled for better visualisation.
In the text

Fig. 7. Spectral evolution of SN 2022qml as shown through the most significant spectra. Numbers close to the spectra indicate the phase from explosion corrected for time dilation. All spectra are scaled with respect to the Hα and arbitrarily shifted for better visualisation. The spectra have not been corrected for extinction. We also plot the narrow Hα in velocity space to show the P Cygni profile due to the wind that formed the CSM.

In the text

	Fig. 8. Spectral evolution of SN 2022wed. Numbers close to the spectra indicate the phase from explosion corrected for time dilation. All spectra have been scaled with respect to the Hα and arbitrarily shifted for better visualisation. The spectra have not been corrected for extinction.
In the text

	Fig. 9. Pre-explosion r-band images of the host galaxies of our SNe IIn retrieved from PS1. Red circles indicate the SNe position. Top left: SN 2021acya. Top right: SN 2021adxl. Bottom left: SN 2022qml. Bottom right: SN 2022wed.
In the text

	Fig. 10. Spectral energy distribution of the host galaxies. Larger filled points represent our measurements on the PS1 template images, while empty points are from NED (or Brennan et al. 2024, in the case of SN 2021adxl). The fluxes have been corrected for the luminosity distance, and they are shown at the rest-frame wavelength position.
In the text

	Fig. 11. Zoom on the Hα region of medium-resolution spectra at ∼100 days. All spectra have been normalised with respect to the Hα and are shown in logarithmic scale for better visualisation.
In the text

Fig. 12. Results of the multifit on the broad Hα and Hβ for the SNe in the sample. Top left: central wavelength of the Hα line. Top right: full width half maximum of the Hα line. Bottom left: Hα luminosity. Bottom right: ratio of $\frac{\mathrm{H}\alpha }{\mathrm{H}\beta }$ from the values measured for our spectra. The dashed line indicates the division between recombination and collision, which is considered to be 3.1, as in Osterbrock & Ferland (2006) (case B). Data from SN 2010jl taken from Fransson et al. (2014) and from SN 2013L taken from Taddia et al. (2020) are also added given their spectral similarity to the SNe in the sample.

In the text

	Fig. 13. Spectral comparison between SNe 2021adxl and 2022wed. All spectra have been redshift corrected, continuum subtracted, and rescaled for better visualisation. We can also clearly see the P Cygni profiles of He Iλ5876 and Fe IIλ5169 for SN 2021adxl.
In the text

	Fig. 14. Photospheric expansion velocity of SN 2021adxl from the He I P Cygni absorption (green diamonds) and Fe II (blue squares) along with the shock velocity derived from the blue shoulder of the broad Hα (magenta reverse triangles). The black line represents the best exponential decline fit of the Hα shock velocity based on our measurements.
In the text

	Fig. 15. Zoom on the Hα profile of selected spectra of SN 2021acya. Magenta lines represent the reflected profile. The profiles have a velocity shift of 150 (+24 days), 100 (+61 days), 500 (+112 days), and 300 km s−1 (+347 days).
In the text

	Fig. 16. Electron-scattering temperature for SNe 2021acya and 2022wed as inferred from the broad Hα FWHM, assuming an optical depth of τe = 5.
In the text

	Fig. 17. Spectrum of SN 2021acya at +61 days (black). The green dashed line is the estimated BB flux from the fit on the red part of the spectrum, while the blue dotted line is the estimated BB flux from the fit on the whole spectrum. We note that the blue BB does not fit the spectrum between 5000 and 6000 Å well.
In the text

	Fig. 18. Blue bump excess as a function of the g, r, i luminosity. Numbers close to the points indicate the average phase.
In the text

	Fig. 19. Pseudo-bolometric light curves computed using only optical filters. We also include the bolometric light curves, calculated in the same way, of SNe 2010jl and 2004et for comparison.
In the text

	Fig. 20. Ultraviolet and NIR contribution for SN 2021acya. The upper panel shows the differences between the luminosities computed using UV and NIR observations with respect to the bolometric light curve computed taking into consideration only optical bands.
In the text

	Fig. 21. Mass-loss rate as a function of the time from the onset of the mass-loss episode, assuming a steady wind and shock velocity as derived in Sect. 2.6.
In the text

Fig. 22. Total neutrino energy contour plots. Top: contour plots of the total neutrino energy ${\epsilon }_{\nu +\overline{\nu }}$ integrated for Eν ≥ 1 TeV from Pitik et al. (2023) in function of the CSM radius and mass, over which the same parameters for our SNe are plotted. The orange star indicates the benchmark of the model. Bottom: same as above but in function of kinetic energy and CSM mass.

In the text

	Fig. B.1. Colour-colour r − i versus g − r diagram of our SNe with and without the K correction. The extinction vector E(B-V) and a BB with temperature between 25000 and 5000 K are also added.
In the text