We present a comprehensive analysis of the photometric and spectroscopic evolution of SN2021foa unique among the class of transitional supernovae for repeatedly changing its spectroscopic appearance from hydrogentoheliumtohydrogendominated IIntoIbntoIIn within 50 days past peak brightness. The spectra exhibit multiple narrow approx 300600kms1 absorption lines of hydrogen helium calcium and iron together with broad helium emission lines with a fullwidthathalfmaximum FWHM of sim 6000kms1. For a steady windmass loss regime light curve modeling results in an ejecta mass of sim 8 Modot and CSM mass below 1 Modot and an ejecta velocity consistent with the FWHM of the broad helium lines. We obtain a massloss rate of approx 2 Modot rm yr1. This massloss rate is three orders of magnitude larger than derived for normal Type II SNe. We estimate that the bulk of the CSM of SN2021foa must have been expelled within half a year about 15 years ago. Our analysis suggests that SN2021foa had a helium rich ejecta which swept up a dense shell of hydrogen rich CSM shortly after explosion. At about 60 days past peak brightness the photosphere recedes through the dense ejectaCSM region occulting much of the redshifted emission of the hydrogen and helium lines which results in observed blueshift sim 3000kms1. Strong mass loss activity prior to explosion such as those seen in SN2009iplike objects and SN2021foa as precursor emission are the likely origin of a complex multipleshell CSM close to the progenitor star.
['Astrophysics - High Energy Astrophysical Phenomena']
SN 2021foa The FlipFlop Type IIn  Ibn supernova
{'No Header': 'Draft version October 29, 2024  \nTypeset using L A T E X twocolumn style in AASTeX631', "SN 2021foa: The 'Flip-Flop' Type IIn / Ibn supernova": "D. Farias, 1 C. Gall, 1 G. Narayan, 2 S. Rest, 3 V. A. Villar, 4, 5 C. R. Angus, 1, 6 K. Auchettl, 7, 8 K. W. Davis, 7 R. J. Foley, 7 A. Gagliano, 4, 9, 5 J. Hjorth, 1 L. Izzo, 1, 10 C. D. Kilpatrick, 11 H .M. L. Perkins, 2 E. Ramirez-Ruiz, 7 C. L. Ransome, 4 A. Sarangi, 1 R. Yarza, 7 D. A. Coulter, 12 D. O. Jones, 13 N. Khetan, 14 A. Rest, 12, 3 M. R. Siebert, 12 J. J. Swift, 15 K. Taggart, 7 S. Tinyanont, 7, 16 P. Wrubel, 15 T. J. L. de Boer, 17 K. E. Clever, 7 A. Dhara, 7 H. Gao, 17 and C.-C. Lin 17  \n1 DARK, Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200 Copenhagen, Denmark  \n2 Department of Astronomy, University of Illinois at Urbana-Champaign, 1002 W. Green St., IL 61801, USA  \n3 Department of Physics and Astronomy, The Johns Hopkins University, Baltimore, MD 21218.  \n4 Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA  \n5 The NSF AI Institute for Artificial Intelligence and Fundamental Interactions  \n6 Astrophysics Research Centre, School of Mathematics and Physics, Queen's University Belfast, Belfast BT7 1NN, UK  \n7 Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064, USA  \n8 School of Physics, The University of Melbourne, VIC 3010, Australia  \n9 Department of Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA  \n10 INAF-Osservatorio Astronomico di Capodimonte, salita Moiariello 16, I-80121, Naples, Italy  \n11 Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA) and Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208, USA  \n12 Space Telescope Science Institute, Baltimore, MD 21218.  \n13 Institute for Astronomy, University of Hawaii, 640 N. A'ohoku Pl., Hilo, HI 96720, USA  \n14 School of Mathematics and Physics, The University of Queensland, QLD 4072, Australia  \n15 The Thacher School, 5025 Thacher Rd., Ojai, CA 93023, USA  \n16 National Astronomical Research Institute of Thailand: Chiang Mai, TH  \n17 Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822, USA", 'ABSTRACT': 'We present a comprehensive analysis of the photometric and spectroscopic evolution of SN 2021foa, unique among the class of transitional supernovae for repeatedly changing its spectroscopic appearance from hydrogen-to-helium-to-hydrogen-dominated (IIn-to-Ibn-to-IIn) within 50 days past peak brightness. The spectra exhibit multiple narrow ( ≈ 300-600 km s -1 ) absorption lines of hydrogen, helium, calcium and iron together with broad helium emission lines with a full-width-at-half-maximum (FWHM) of ∼ 6000 km s -1 . For a steady, wind-mass loss regime, light curve modeling results in an ejecta mass of ∼ 8 M ⊙ and CSM mass below 1 M ⊙ , and an ejecta velocity consistent with the FWHM of the broad helium lines. We obtain a mass-loss rate of ≈ 2 M ⊙ yr -1 . This mass-loss rate is three orders of magnitude larger than derived for normal Type II SNe. We estimate that the bulk of the CSM of SN 2021foa must have been expelled within half a year, about 12 years ago. Our analysis suggests that SN 2021foa had a helium rich ejecta which swept up a dense shell of hydrogen rich CSM shortly after explosion. At about 60 days past peak brightness, the photosphere recedes through the dense ejecta-CSM region, occulting much of the red-shifted emission of the hydrogen and helium lines, which results in observed blue-shift ( ∼ -3000 km s -1 ). Strong mass loss activity prior to explosion, such as those seen in SN 2009ip-like objects and SN 2021foa as precursor emission, are the likely origin of a complex, multiple-shell CSM close to the progenitor star.  \nKeywords: Supernovae(1668) - Stellar masss loss(1613) - Core-collapse supernovae(304)', '1. INTRODUCTION': "Massive stars ( ≳ 8 M ⊙ ) undergo different mass-loss phases, losing material from their outer layers shortly before core-collapse (CC, Smartt 2009; Smith 2017). Analysis of this expelled material, termed circumstellar material (CSM), can provide important information about the progenitor system and thus the late-stage of massive stellar evolution (see e.g., Morozova et al. 2018). The CSM surrounds the progenitor and thus, the supernova (SN) radiation and ejecta inevitably interact with the CSM. The emanating signatures arising from the interaction appear at a variety of phases during the evolution of the SN, depending primarily upon the mass distribution of the CSM (Dessart & Hillier 2022). This interaction produces SN spectra that can be dominated by narrow ( ∼ 100 -500 km s -1 ) or intermediately broad ( ∼ 1000 km s -1 ) emission lines and P-Cygni profiles (see Fraser 2020, and references therein).  \nDepending upon the progenitor system, as well as the composition and radial distribution of the CSM, different classes of core-collapse supernovae (CCSNe) have been identified.  \nCCSNe with a hydrogen (H) rich CSM and little helium (He) emission in their spectra are commonly classified as Type IIn (Schlegel 1990). Classical examples of such events are e.g., SN 1998S, SN 2005ip and SN 2010jl (Fox et al. 2009; Mauerhan & Smith 2012; Fransson et al. 2014; Gall et al. 2014). However, if the CSM is He-rich with little-to-no H emission in the SN spectra, then the CCSNe are classified as Type Ibn (Pastorello et al. 2007; Foley et al. 2007). In recent years, another class of interacting SNe has emerged, the Type Icn (Gal-Yam et al. 2022; Pellegrino et al. 2022; Davis et al. 2023). These SNe exhibit narrow oxygen (O) and carbon (C) emission lines in their spectra.  \nIn the local universe Type IIn and Ibn SNe comprise about 5% and 1% of the volumetric rate of CCSNe, respectively (Maeda & Moriya 2022; Cold & Hjorth 2023). Among the interacting CCSNe, Type Icn are the rarest, with only five members known thus far (see Davis et al. 2023). However, the classification of several interacting SNe is ambiguous, as some CCSNe appear to change their type, e.g., from Type IIn to Ibn or vice versa. Prominent examples are SN 2005la (Pastorello et al. 2008), SN 2011hw (Smith et al. 2012; Pastorello et al. 2015a), iPTF15akq (Hosseinzadeh et al. 2017) and SN 2020bqj (Kool et al. 2021). These objects form the unique group of transitional IIn/Ibn SNe.  \nDetermining the exact nature of the progenitor of interacting SNe is challenging due to the complexities of the interaction between the ejecta and the CSM. Thus, SNe IIn have diverse light curves (LCs), spanning a broad range of peak magnitudes (Nyholm et al. 2020) and LC shapes. This has led to suggest a wide range of plausible progenitors systems for Type IIn SNe, such as red supergiants (RSG) in binary systems (SN 1998S-like objects; Smith 2017), Luminous Blue Variables (LBV) for e.g., SN 2005gl (Gal-Yam et al. 2007), SN 2009ip (Foley et al. 2011; Smith et al. 2013), SN 2010jl (Smith et al. 2011) and SN 2015bh (Boian & Groh 2018), while a ∼ 20 M ⊙ star is preferred for SN 2016jbu (Kilpatrick et al. 2018; Brennan et al. 2022a). The progenitor of SN 2015bh is also proposed to be a yellow supergiant ( ∼ 50 M ⊙ ; Thone et al. 2017).  \nContrary, LCs of Type Ibn SNe show a high degree of homogeneity (Hosseinzadeh et al. 2017). Thus, the most plausible progenitor is an anevolved Wolf-Rayet star (Pastorello et al. 2008), which is consistent with the majority of Ibn SNe being found in star-forming galaxies (although see PS1-12sk; Hosseinzadeh et al. 2019). Furthermore, the mass loss rates derived from LC modeling of e.g., OGLE-2014-SN-131 and SN 2020bqj (e.g., Karamehmetoglu et al. 2017; Kool et al. 2021) favour such a progenitor. Alternatively, Type Ibn SNe may be the explosion of a low mass helium star in a binary system (Dessart et al. 2022). Unlike Type IIn SNe (see SN 2005gl, Gal-Yam et al. 2007), there are no detection of any progenitor system of Type Ibn SNe in archival data. However, late-time photometry at the location of SN 2006jc has shown a potential companion associated with the exploding star (Maund et al. 2016; Sun et al. 2020).  \nA handful of CCSNe exhibited pre-explosion activities or outbursts up to two decades prior to their terminal explosion. Precursor emission in Type IIn SNe is common (e.g., Mauerhan et al. 2013; Ofek et al. 2013b; Thone et al. 2017; Elias-Rosa et al. 2016; Pastorello et al. 2018; Hiramatsu et al. 2024), such as the case of SN 2009ip (Pastorello et al. 2013; Margutti et al. 2014). This transient was first classified as a SN impostor (e.g., Foley et al. 2011). However, after two more outbursts in 2011 and 2012, its 'final' re-brightening in 2012 reached M R ≈ -18 mag, which was proposed as the terminal explosion of a SN IIn (Mauerhan et al. 2013). On the contrary, precursor emission has only been observed for three Type Ibn supernovae such as SNe 2006jc (Fo-l  \net al. 2007), 2019uo (Strotjohann et al. 2021) and 2023fyq (Brennan et al. 2024; Dong et al. 2024).  \nInteracting SNe with observed signatures of preexplosion outbursts and a photometric and spectroscopic evolution similar to that of SN 2009ip, are termed 2009ip-like objects (Pastorello et al. 2018; Brennan et al. 2022b). Prominent examples of this class include SN 2015bh (Elias-Rosa et al. 2016; Thone et al. 2017), SN 2016jbu (Kilpatrick et al. 2018; Brennan et al. 2022b) and SN 2019zrk (Fransson et al. 2022). Typically, photometric and spectroscopic data obtained around the epochs of the outbursts suggest that these are LBV-like eruptions prior to the presumed, terminal explosion of the progenitor (Pastorello et al. 2013; Thone et al. 2017; Brennan et al. 2022a). However, whether or not the latter are indeed, stellar explosions remains unclear (Smith et al. 2022).  \nHere, we present unpublished multi-band photometry and time-series spectroscopy of the fifth transitional Type IIn/Ibn SN 2021foa. The data were collected by the Young Supernova Experiment (Jones et al. 2021; Aleo et al. 2023; Coulter et al. 2023). The SN (RA=13:17:12.29, DEC= -17:15:24.19) was discovered by ASAS-SN (Stanek & Kochanek 2021) on March 15, 2021 ( g ∼ 15 . 9) in the galaxy IC 086 (ASASSN-21dg). It was initially classified as a Type IIn due to the strong, narrow Balmer lines in the optical spectrum (Angus 2021). An analysis of its light curve and spectra until +79 days (Reguitti et al. 2022) suggests that SN 2021foa is photometrically similar to SN 2009ip-like objects such as SN 2005gl (Gal-Yam et al. 2007), SN 2009ip (Pastorello et al. 2013) and notably, SN 2016jbu (Kilpatrick et al. 2018; Brennan et al. 2022b), while spectroscopically resembles the transitional IIn/Ibn SNe. In this work, we present a comprehensive analysis of the spectroscopic and photometric evolution out to +427 days. We detail a physical picture of this unusual SN, which, amongst its peers in the 2009ip-like class, exhibits some unique characteristics.  \nWe determine the time of maximum light in r -band using a second-degree polynomial fit between MJD 59280 and 59315 to be MJD max =59302 . 35 ± 0 . 14 ( r max ). We use this as our reference time in the remainder of the paper. Throughout this work, we assume a standard ΛCDM cosmology with H 0 = 67 . 8 km s -1 Mpc -1 and Ω m = 0 . 307 (Planck Collaboration et al. 2014).  \nThe paper is structured as follows: Sect. 2 describes the data acquisition and reduction. In Sect. 3 we present our photometric and spectroscopic analysis methodology, with the results presented in Sect. 4. We discuss our results from the analysis in Sect. 5, and conclude by  \nformulating a complete physical scenario for SN 2021foa in Sect. 6.", '2. OBSERVATIONS': 'Our photometric and spectroscopic follow up observations of SN 2021foa, including archival and public data are described below.', '2.1. Photometric data': "We obtained optical photometry with the Sinistro imagers on the Las Cumbres Observatory (LCO) 1-m telescope network (Brown et al. 2013) in UBVRI,up,gp,rp,ip,z bands, starting March 15, 2021 through August 17, 2021. Initial automatic processing of the imaging data, including instrument signature removal, pixel-level corrections and astrometric calibration, was performed by the LCO BANZAI pipeline (McCully et al. 2018). Thereafter, images were processed using photpipe pipeline (Rest et al. 2005; Jones et al. 2021). We measured the flux of SN2021foa from the LCO images using an updated version of DoPhot (Schechter et al. 1993), and this photometry was calibrated using u -band SDSS (Alam et al. 2015) together with griz Pan-STARRS1 photometric standards observed in the vicinity of SN 2021foa. We calibrated the photometry of R and I bands on the Cousins photometric system.  \nObservations with the Ultraviolet and Optical Telescope (UVOT; Roming et al. 2005) on the Neil Gehrels Swift Observatory were reported in Reguitti et al. (2022). UVOT data of SN 2021foa were taken between March 16, 2021 and May 10, 2021. Following the methodolgy described in Brown et al. (2014), we use uvotsource from the HEASoft v6.26 package to perform aperture photometry within a 3 '' aperture centered on SN 2021foa. We measured the total background flux at the location of SN 2021foa from frames obtained on June 17 2022, when any residual light from the fading SN is well below the sky-background. The background emission was then subtracted from all previous observations. We detect emission of SN 2021foa at a > 3 σ level in all UVOT bands in the 2021 observations.  \nSN 2021foa was also observed by the Asteroid Terrestrial-impact Last Alert (ATLAS) System (Tonry et al. 2018) between February 22 to July 16 and March 20 to July 14 in o - and c -bands, respectively. Following the procedure described by Davis et al. (2023), we obtained the binned light curve data calculated as a 3 σ -cut weighted mean for each night. In contrast to Reguitti  \nFigure 1. Photometry of SN 2021foa. Swift UVW1, UVM2, UVW2 , Thacher griz , Nickel BVri , and LCO/Sinistro UBVRI,up,gp,rp,ip,z data (AB system) are symbol and colour coded as described in the legend. All light curves are offset by a constant, for clarity. The gray vertical lines mark our spectroscopic data at their respective epochs relative to r max at MJD max = 59302 . 35 days. Event A encompasses the time where the light-curve of SN 2021foa shows precursor emission. Uncertainties are smaller than the size of the markers. The green band delimits the phases where a short plateau is observed (from +14 to +22 days). The inset highlights that same phase region.  \n<!-- image -->  \net al. (2022, 2024), we did not find any significant detection in the c -band prior to February 22 (see Sect. B).  \nAdditional images of SN 2021foa were obtained in BVri bands with the 1-m Nickel telescope at Lick Observatory and in up, gp, rp, zp bands with the Thacher 0.7-m telescope in Ojai, CA (Swift et al. 2022). The images from Nickel telescope were calibrated using bias and sky flat-field frames following standard procedures. PSF photometry was performed, and photometry was calibrated relative to Pan-STARRS1 photometric standards (Flewelling et al. 2020). Similarly, photometry of the images from Thacher telescope was obtained using DoPhot and calibrated with the griz Pan-STARRS1 catalog.  \nFig. 1 shows our photometric data of SN 2021foa, including LCO/Sinistro ( UBVRI,up,gp,rp,ip,z ), our re-reduced and host galaxy subtracted UVOT ( UVW 1 , UV W 2 , UV M 2 , U, B and V ) data, ATLAS ( c and o ), Nickel ( BVri ) and Thacher ( griz ) data. All photometric data are summarized in Tab. A1.", '2.2. Spectroscopic data': "Fig. 2 shows our extensive optical spectroscopic follow-up observations of SN 2021foa obtained within -12 days and +427 days. Data were obtained with the Kast dual-beam spectrograph (Miller J. S. & Stone R. P. S. 1993) on the Lick Shane 3-m telescope at -7, +8, +15, +21, +35, +42 and +51 days, the Goodman spectrograph (Clemens et al. 2004) on the NOIRLab 4.1m Southern Astrophysical Research (SOAR) telescope at Cerro Pach'on at -8 day, the Alhambra Faint Object Spectrograph and Camera (ALFOSC) on the Nordic Optical Telescope (NOT) at -11, -7 and -6 days and the Wide Field Spectrograph (WiFeS) at the Australian National University (ANU) 2.3-m telescope located at Siding Spring Observatory (Dopita et al. 2007) at -11 and +51 days.  \nThe Kast observations are performed with the blue side 452/grism, 300/7500 grating, d58 dichroic and 2 '' slit. Goodman observations were carried out using the 400 lines/mm grating with the M1 wavelength setting (300 -705 nm). To reduce the Kast and Goodman spectra, we used the UCSC Spectral Pipeline 1 (Siebert et al. 2019). The ALFOSC spectra were taken with a 1.0 '' slit and grisms 4, and 8. For all reduction, extraction and calibration steps, we used standard IRAF 2 routines using PYRAF 3 .  \nThe Wide Field Spectrograph (WiFeS) is an integralfield spectrograph with a field-of-view of 38 × 25 arcsec. SN 2021foa was observed using a RT-560 beam splitter and a B3000 and R3000 diffraction gratings that covers the 3200 -5900 ˚ A and 5300 -9800 ˚ A wavelength ranges. All observations had a Y=2 binning readout mode, corresponding to a 1 × 1 sq. arcsec spaxel. Each observation was reduced using PyWiFeS (Childress et al. 2014). We extract an isolated part of the sky for background determination and subtraction.  \nFig. 3 presents our NIR spectra of SN 2021foa which were obtained with the SpeX spectrograph (Rayner et al. 2003) mounted on the 3-m NASA Infrared Telescope Facility (IRTF) at -8, +12 and +41 days. In this mode with the 0.8 '' slit, the spectral resolving power is R ≈ 1000. The SN was observed in an ABBA dithering pattern with an A0V star observed immediately before or after science observations for telluric correction. We also obtain observations of internal flat field and arcs calibration lamps at the science pointing. We reduced the data using spextool (Cushing et al. 2004), which performed flat fielding, wavelength calibration, background subtraction, and spectral extraction. We then performed telluric correction using xtellcor (Vacca et al. 2003).  \nMedium resolution spectra were obtained with the X-shooter echelle spectrograph (Vernet et al. 2011) mounted at the Very Large Telescope (VLT) at the European Southern Observatory (ESO) on Cerro Paranal, Chile. The data are presented in Figs. 2 and 3. The X-shooter instrument covers the wavelength range of 0 . 3 -2 . 5 µ m in three arms, the ultraviolet and blue (UVB), visual (VIS) and near-infrared (NIR) wavelength ranges. The slit widths (and resolving power) for UBV, VIS and NIR for these observations were 0.9 '' ( R = 5900), 1.0 '' ( R = 8900) and 1.0 '' ( R = 5600), respectively. UVB, VIS and NIR arms were reduced with esoreflex 2.11.5 4 (Freudling et al. 2013) pipeline individually. Then, using a custom Python code, UVB and VIS arms were combined in STARE mode, while the NIR arm was combined in NOD mode. Special efforts have been made to perform a detailed background subtraction around strong host galaxy emission lines such as H α (see Fig. A1).  \nTo improve the flux-calibration from the mentioned pipelines, we mangle all low resolution spectra to the interpolated photometry from EXTRABOL 5 , as described in Sect. 3.1.1. The UBV and VIS arms of the first  \nFigure 2. Low and medium (XSH;VLT/X-shooter) resolution optical spectra of SN 2021foa (Sect. 2.2, Tab. A2). Black, dashed vertical lines indicate the rest-frame wavelength of the strongest lines observed in SN 2021foa, mostly H i and He i lines. For visualization purposes, we smoothed all the spectra with a rolling Gaussian with kernel size of ∼ 15 ˚ A.  \n<!-- image -->  \nthree VLT/X-shooter spectra (+66, +95 and +129 days) were also flux-calibrated in this way. Due to the lack of NIR photometry at these epochs, we used the calibration provided by esoreflex . The same approach was applied for the last spectrum at +427 days, and as this late-time spectrum includes minimal contamina-  \ntion from the SN, we use narrow H α to determine the redshift ( z = 0 . 0086183).", '2.3. Extinction': 'The Milky Way extinction was obtained from Schlafly & Finkbeiner (2011); E ( B -V ) MW = 0 . 072 for R V = 3.1 and the host galaxy extinction, E ( B -V ) host = 0 . 129,  \nFigure 3. Low and medium (XSH;VLT/X-shooter) resolution NIR spectra of SN 2021foa (see Sect. 2.2, Tab. A2). For visualization purposes, the Xshooter spectra was trimmed from 2.2 to 2 . 5 µ m due to low signal-to-noise in that region. All spectra were smothed with a rolling Gaussian with kernel size of ∼ 15 ˚ A.  \n<!-- image -->  \nwas adopted from Reguitti et al. (2022). We confirmed this value from measurements of the equivalent width (EW) of the Na i D ( λ 5890 , 5896) absorption lines detected in our medium-resolution VLT/X-shooter spectra (see Sect. 3.2) together with the extinction-EW relation from Turatto et al. (2003). We obtain an EW of 0 . 85 ± 0 . 36 ˚ A, similar to measurements by Reguitti et al. (2022).', '3. ANALYSIS METHODOLOGY': 'The rich photometric and spectroscopic data set shows that SN 2021foa has a prolonged precursor emission similar to Type IIn SN 2009ip, but spectroscopically closer to transitional IIn/Ibn objects, with strong He i emission lines in the optical and near-infrared wavelengths. To understand the physical mechanism that drive the similarities and differences and to put together a coherent physical scenario that can explain all observables we employ a range of analysis methods. First, we model the multiband photometry to obtain information about the progenitor and the CSM associated to SN 2021foa, such as the ejecta and CSM masses, and progenitor radius. Secondly, we investigate the evolution of the most prominent lines (e.g., H i , He i , Ca ii ) of SN 2021foa by fitting their broad and narrow velocity components in emission and absorption. Such analysis provides insights  \nof the complexity of the CSM while also highlights the transitional nature of SN 2021foa among IIn and Ibn classes. Lastly, we fit two black body functions to the X-shooter spectra from optical to NIR wavelengths to estimate the temperature and amount of dust present at the respective epochs in SN 2021foa.', '3.1. Photometry': 'We begin with an analysis of the photometry of SN 2021foa, including a comparison to well-studied sources to inform the time-scales and energetics of the explosion, which in turn constrain the progenitor system. Fig. 4 shows the ATLAS o -band light curve of SN 2021foa compared to R -band-like ( o , r and R ) light curves of other interacting SNe. Our comparison sample consists of transitional Type IIn/Ibn SNe 2005la, 2011hw, iPTF15akq and 2020bqj. We also include SNe AT2016jbu, 2005gl and 2016bdu, members of the 2009ip-like class with a photometric resemblance to SN 2021foa. The prototypical Type Ibn SN 2006jc and SN 2023fyq (ATLAS photometry using ATClean , Rest et al. 2023) serve as representation of a well studied, interacting He-rich SN. Additionally, we compare SN 2021foa against the Type Ibn template from Hosseinzadeh et al. (2017), which represents the average, homogeneous photometric evolution of this SN class as well as the weighted mean ( ± σ ) of the photometry of all stripped-envelope, hydrogen poor Type Ibc SNe in Drout et al. (2011). Interestingly, there is a plateau  \nof ≈ one week duration starting from ≈ two weeks after maximum light in SN 2021foa (see inset in Fig. 1). A similar plateau lasting for about two weeks can also be observed for SN 2016jbu (Brennan et al. 2022b).  \nSN 2009ip-like SNe are characterized by outbursts of M R ∼ -11 ± 2 mag several years before the explosion, faint emission mimicking SN impostors of M R ∼ -13 ± 2 mag (event A) and a peak magnitude of M R ∼ -18 . 5 ± 0 . 5 (event B, Pastorello et al. 2018; Brennan et al. 2022a). For SN 2021foa, we set an upper limit of M o ≈ -13 . 4 mag for any outburst within ∼ 5 years before explosion (Sect. 5.2), while event A ( M r ≈ -14 mag) and event B ( M r ≈ -18 mag) agree with the values of 09ip-like class. The peak magnitude of Ibn SNe is, on average, much larger than that of SN 2021foa ( M R = -19 . 47 -0 . 32 +0 . 54 , Hosseinzadeh et al. 2017). Since rise/decay-times/slopes depend on the assumed explosion date and the interaction between the CSM and the SN ejecta, it is difficult to quantify the changing behaviour of the light curves of 09ip-like SNe. Nonetheless, the decay slopes from +0 to +30 days past r -band maximum of SN 2009ip, AT 2016jbu and SN 2021foa ( ∼ 0 . 04, ∼ 0 . 07 and ∼ 0 . 04 mag day -1 , Farias et al., in prep) are smaller than the average decay rate of SN Ibn (0 . 1 mag day -1 , Hosseinzadeh et al. 2017). We note that the light curve of SN 2021foa shows a clear re-brightening in r -bands-like after ≈ +80 days.', '3.1.1. Constructing a Bolometric light curve': 'To obtain the bolometric light curve we fit the spectral energy distribution (SED) with a blackbody (BB) model using EXTRABOL (Thornton & Villar 2022). EXTRABOL models the light curves in flux space as a function of wavelength and phase with a 2D Gaussian process (typically with zero mean and stationary kernels in wavelength and phase, though these are configurable). A model SED (a BB model in the case of SN 2021foa) is fit to match the posterior mean of the Gaussian process, and the resulting SED surface is marginalized over wavelengths to produce a quasi-bolometric light curve. We combined our optical observations in Sect. 2.1 with J , H and griz band photometry from Reguitti et al. (2022) to construct the SED from B band to near-infrared wavelengths. We computed the bolometric light curve only for epochs for which data were available in all selected passbands, as the 2D Gaussian process regresses to the mean when not constrained by observations.  \nFig. 5 shows the evolution of the photospheric BB radius ( R BB ), temperature ( T BB ), and the resulting bolometric LC (L bol ) using EXTRABOL . Additionally, we computed the bolometric light curves using SUPERBOL to show the consistency of our method, assuming a con-  \nFigure 4. R/r band light curve comparison of interacting IIn and Ibn SNe. Light curves are shown for SNe 2005la (Pastorello et al. 2008), 2011hw (Smith et al. 2012; Pastorello et al. 2015a), iPTF15akq (Hosseinzadeh et al. 2017), 2020bqj (Kool et al. 2021), 2021foa (Reguitti et al. 2022), 2009ip-like SNe 2009ip (Pastorello et al. 2013; Fraser et al. 2013), 2016jbu (Kilpatrick et al. 2018; Brennan et al. 2022b) and 2016bdu (Gal-Yam et al. 2007; Pastorello et al. 2018) including the Type Ibn SN 2006jc (Pastorello et al. 2007) and SN 2023fyq (Brennan et al. 2024; Dong et al. 2024). All data are symbol and color coded according to the legend. The gray shaded region corresponds to the Type Ibn template of (Hosseinzadeh et al. 2017, H17). The green shaded area is the weighted mean ± σ of the photometrical sample of Type Ibc SNe in Drout et al. (2011, D11). The reference epoch of SN 2016jbu is with respect to V -band maximum. Vertical lines delimit the event A of SN 2021foa in Fig. 1.  \n<!-- image -->  \nstant color extrapolation at early and late times. For comparison, we also included the bolometric light curve of SN 2009ip, using UVW2, UVM2, UVW2, U, B, R, V, I bands to construct the SED. We address that, even if our fits reproduce the overall evolution of SN 2009ip shown in Margutti et al. (2014), they fitted a twoblackbody function ( hot and cold components) to the SED up to NIR bands. As expected, the evolution of R BB , T BB and L bol of both SNe share similarities such as a monotonic increase of R BB before r -band maximum and the phase at which L bol reaches maximum ( ≈ -2 days). Notoriously, the evolution of R BB drastically changes, possibly due to inclusion of UV bands to construct the SED of SN 2009ip.  \n3.1.2. Modeling the light curve with the MOSFiT Framework  \nFigure 5. Upper , middle and bottom panels show the evolution of the radius ( R BB ), temperature ( T BB ), and bolometric luminosity ( L bol ) of SN 2021foa (magenta) and SN 2009ip (yellow) obtained through blackbody fits to the SED of SN 2021foa using EXTRABOL (squares) and SUPERBOL (circles, Nicholl 2018). The green shaded area represent the epochs of the short plateau in the light curves (see Fig. 1).  \n<!-- image -->  \nIn order to investigate the progenitor properties and its mass loss history, we adopt the radioactive decay (RD) model of 56 Ni + circumstellar interaction (CSI) from the MOSFiT (Guillochon et al. 2018) framework 6 based on Chatzopoulos et al. (2012, 2013). The RD + CSI model has been extensively used to model LCs for different SNe Ibn (Karamehmetoglu et al. 2017; Kool et al. 2021; Ben-Ami et al. 2023). MOSFiT estimates the posterior probability of a set of parameters from the RD+CSI model given the photometric data and physically-informed priors.', '3.1.3. Parameterization of the RD+CSI Model': 'The RD model accounts for the radioactive decay of 56 Ni through three parameters: the nickel fraction, f Ni , the γ -rays opacity of the ejecta, κ γ , and an optical opacity, κ . Additionally, MOSFiT includes a model for the SN ejecta-CSM interaction, with several physical quantities for the ejecta and the CSM listed as free parameters (see Villar et al. 2017, for a detailed description). This includes the total mass of the ejecta ( M ej )  \nand the inner and outer density profiles of the ejecta ( ρ in , ej ∝ r -δ , ρ out , ej ∝ r -n ). The main parameters to describe the CSM are the inner radius of the CSM ( R 0 ), the CSM density ( ρ CSM ) and a density profile of the CSM ( ρ CSM ∝ r -s ). For the latter, s = 0 corresponds to a constant density (shell-like) CSM, while s = 2 denotes a wind-like model. Additionally, MOSFiT fits for the explosion time ( t exp ) relative to the first photometric observation (Kool et al. 2021). MOSFiT also includes some nuisance parameters, such as the minimum allowed temperature of the photosphere before it recedes, T min . Finally, a white noise term σ , is included and added in quadrature to the photometric uncertainties. MOSFiT uses this white noise term to quantify the quality of the fitting. Acceptable values for this parameter are σ < 0 . 2 (Nicholl et al. 2017).', '3.1.4. RD+CSI Modeling Choices and Priors': 'While this is a high-dimensional model, the fundamental physical parameters of interest to this study are f Ni , M ej , M CSM , E kin and R 0 . We fix three parameters of the model: δ = 1, n = 12 and s . The assumed value of n = 12 is typically used to characterize RSG-like progenitors rather than for BSG/WR stars, but δ and n are not sensitive to derive the fundamental parameters (Villar et al. 2017; Kool et al. 2021). In order to investigate a shell and wind-like CSM configurations, we fixed s = 0 and s = 2, respectively. With these choices, eleven free parameters remain: f Ni , κ , M ej , M CSM , ρ CSM , R 0 , E kin , t exp (i.e. fundamental parameters) T min , κ γ and σ (the nuisance parameters). To constrain the explosion time of SN 2021foa, we only considered observations from -20 days onwards, i.e., we discard event A photometry as the RD+CSM model was not designed to model precursor emission. We used top-hat priors on the remaining parameters. We note that other physical parameters including κ γ , κ and t exp may be considered nuisance parameters as these are quantities that are not observable. These physically-informed modeling choices and priors significantly reduce the dimensionality of the problem, and help ensure that the inferred quantities are meaningful.  \nFollowing Kool et al. (2021), we elected to use the dynamic nested sampling implemented in DYNESTY (Speagle 2020). Nested samplers estimate the evidence ( Z ) i.e., the marginalized likelihood for the data given the model. Given that our prior distributions are constant for s = 0 and s = 2, we can directly compare the evidence reported by DYNESTY to check which model is preferred given our data (i.e., the model with larger evidence).  \nFigure 6. Collection of light curves of SN 2021foa in several bands reported in this work and Reguitti et al. (2022). Colour symbols are the photometric data. Light and bold colour lines are 1000 sampled MOSFiT light curves from the posterior distribution for the shell ( s = 0) and wind-like CSM configuration ( s = 2), respectively.  \n<!-- image -->  \nmax  \nIn order to reduce the computational cost without oversampling the observations per epoch, we provided MOSFiT with the following bands: UVW 2 , UV M 2 , UV W 1 , U, B, V (UVOT), u, B, g, V, r, i, z (Reguitti et al. 2022), o, c (ATLAS) and up, gp, rp, ip (LCO). Fig. 6 shows the MOSFiT model-  \nlight curves for s = 0 and s = 2 for a collection of bands considered, together with the observed data. The parameters of the model are summarized in Tab. 1 and discussed in Sect. 4.2.', '3.2. Spectroscopy': "Table 1. Fundamental parameters of the MOSFiT (CSM-Ni) models of a shell-like ( s = 0), wind-like ( s = 2).  \nNote -a Parameters of the radioative decay+CSI model are the nickel fraction ( f Ni ), optical opacity ( κ ), kinetic energy ( E kin ), ejecta mass ( M ej ), CSM inner radius ( R 0 ), CSM density ( ρ CSM ), explosion epoch with respect to r max ( t 0 ), nickel mass ( M Ni ), ejecta velocity ( v ej ) and the evidence, (log Z ), detailed in Sect. 3.1.4.  \nOur multi-epoch spectroscopy permits a detailed analysis of the evolution of the prominent emission and absorption line profiles of SN 2021foa (see Figs. 2 and 3) over a large range of epochs. We select the spectra with highest resolution, including those published by Reguitti et al. (2022) at +35, +43 days. For simplicity, we refer to any 'velocity component' of a line profile as 'component'.", '3.2.1. Spectral Line fitting': 'In order to quantify the evolution of the different components of the Balmer lines (H α , H β , H γ and H δ , panels A, B, C and D in Fig. A3), we fit the lines with multiple Gaussian and Lorentzian functions (see Fig. 7). Prior to this, we fit the continuum in the respective wavelength region of the emission lines (around ± 5000 km s -1 ) using a 1D polynomial and subtract it. Furthermore, to account for either low instrumental resolution or a low signal-to-noise ratio, we convolved the spectral region of the emission lines with a Gaussian kernel that has a standard deviation ≤ 2 × the dispersion ( ∼ resolution) of the instrument (see Tab. A2). We detail the line profile fitting method for H α as an example.  \nFor the first two epochs ( -11 and -8 days), we simultaneously fit the complex H α emission line with two  \nfunctions. Since the H α exhibits broad wings, we fit a Lorentzian profile to the broad base of the emission line and a Gaussian profile to the narrow component. The introduction of the Lorentzian profile was made to account for any electron scattering wings at early times (Fransson et al. 2014). For all epochs between -3 and +51 days we fit the H α emission line with three Gaussian functions, one each for the narrow absorption, narrow emission as well as for the broad emission profile. At all late epochs (+66 to +129 days) we add another Gaussian profile to fit the strongly blue-shifted emission of H α . For all other H i emission lines we follow a similar method with two minor differences: (i) the central wavelength used to extract the line profile (ii) the number of Gaussian functions to fit the line profiles. For example, for Pa β , a single Gaussian profile can well fit the emission line which neither exhibits narrow emission nor absorption.  \nFor the He i lines, we only quantify the evolution of the strongest emission lines in the optical (He i λ 5876 ˚ A and 7065 ˚ A). We fit these symmetric features with a single Gaussian at all epochs (see Sect. C.0.4). This is the simplest model that is consistent with the observations.  \nWe also identify narrow absorption lines of Fe ii λλ 5169 , 5276 , 5317, (likely) He i λλ 4922 , 5016, Ca ii λλ 8498 , 8542 , 8662 and O i λ 8446 (see Figs. A3 and A4). To quantify the velocity evolution of the region harbouring these elements, we fit a single Gaussian profile to the absorption lines.  \nFinally, we calculate the total line flux of all H i emission lines from the continuum subtracted line profile fits. We exclude the emission lines from the host galaxy at late epochs ( > +40 days). The total emission line flux of the He i λ 5876 ˚ A line complex is estimated from the continuum subtracted data directly. the emission line complex in the velocity range ± 6 000 km/s.', '3.3. Modeling the Late-time VLT/X-shooter Spectra': 'Asignature of either pre-existing or newly formed dust is the thermal (near) infrared excess emission of hot dust grains over the SN continuum emission. Such thermal emission has been observed in several types of core collapse SNe (see Gall et al. 2011, for a review). To test if dust is present in SN 2021foa, we fit the continuum emission from optical to the NIR wavelengths of our late-time SN 2021foa VLT/X-shooter spectra at +66, +95 and +129 days. We exclude the spectral regions containing strong emission lines such as H α and the He i λ 10830+ and Pa γ complexes prior to fitting as indicated in Fig. 8.  \nWe then simultaneously fit a BB function to the hot SN photosphere and a modified BB function (Hildebrand 1983) to the near IR excess emission to account  \nFigure 7. Modeling of the H α profiles. Zero velocity corresponds to the rest wavelength of H α line. The first two epochs were fitted using a combination of a Gaussian and Lorentzian profile (red, green), while the spectra from -3 to +51 days are fit by three Gaussians, accounting for the narrow and broader emission (green and red), and the absorption (purple). For the medium resolution spectra (+66, +95 and +129 days), an additional Gaussian component has been added to fit the profile (turquoise). Blue and red dashed vertical lines mark the velocities of the rest frame wavelength of H α and He i 6678 ˚ A, respectively.  \n<!-- image -->  \nfor any hot dust emission as:  \nF λ ( λ ) = R 2 SN D 2 L B λ ( λ, T SN ) + M d D 2 L κ abs ( λ, a ) B λ ( λ, T d ) (1)  \nwith D L the luminosity distance to the SN, R SN the radius of the photosphere and B λ ( λ, T ) the Planck function at temperatures of the photosphere, T = T SN and dust, T = T d . We assume optically thin dust, with a dust mass M d , which is distributed spherically sym-  \ntrically around the SN. For the dust mass absorption coefficient, κ abs ( λ ) we adopt the formalism of Gall et al. (2017), where κ abs ( λ ) follows a λ -β power law in the VLT/X-shooter NIR wavelength range. Thus, κ abs ( λ ) can be parameterized as A d ( λ/ 1 µ m) -β , with A d = 1 . 0 × 10 4 cm 2 g -1 and β = 1 . 5. This formalism mimics small sized ( ≲ 0 . 1 µ m) carbonaceous grains (e.g., Rouleau & Martin 1991).  \nFigure 8. SN and dust fit to VLT/X-shooter spectra at +66 days ( upper panel ), +95 days ( middle panel ) and +129 days ( lower panel ). The SN hot black-body function (green curves - dominant below ≈ 8 , 000 ˚ A) is fit simultaneously with a modified black-body function to describe the NIR dust emission (red curves) to the observed spectra (black). The blue curves show the sum of the two individual components. Cyan rectangles visualise the spectral regions excluded for the fit due to the large number of emission lines.  \n<!-- image -->', '4.1. Estimating Photosphere Expansion Velocity from the Bolometric Light Curve': 'The evolution of T BB follows the bolometric LC. The peak values of L bol and T BB , reached between -5 and 0 days, are ≈ 1 . 7 × 10 43 erg s -1 and ≈ 15000 K, respectively. The temperature evolution flattens between 14 -23 days as indicated by the green shaded region in Fig. 5. This coincides with the short LC-plateau that is observable in several optical passbands (see Fig. 1 and Sect. 3.1). Subsequently, T BB decreases to 8000 K at  \nday +35. On the other hand, R BB is largest, with about 10 15 cm at the beginning of the short LC-plateau.  \nFrom the estimated slope of the R BB evolution until this maximum at ∼ +10 days, we infer an expansion velocity of the photosphere of ≈ 2000 km s -1 . Thereafter, R BB recedes at a velocity of ∼ -3200 km s -1 to a radius of ≈ 6 × 10 14 cm at day ∼ +30, the same level as at -10 days.', '4.2. Estimating Progenitor Properties from the RD+CSI Model': 'We use the parameters of the RD+CSI model inferred with the MOSFiT framework (Sect. 3.1.2, Tab. 1) to constrain progenitor and SN properties. From the evidence values (log Z ), s = 2 model is preferred over s = 0. We obtain an ejecta mass of ≈ 2 . 0 M ⊙ for a s = 0 (CSMshell-like) model. This value is close to the maximum mass of ≈ 1 . 2 M ⊙ , as estimate for any progenitor (single or in a binary system) of Type Ibn SNe (Dessart et al. 2022). However, the ejecta mass inferred for s = 2 is ≈ 8 M ⊙ . This is lower than the average value of ≈ 16 M ⊙ , as obtained for Type Ibn SN 2019uo (Gangopadhyay et al. 2020), 2020bqj (Kool et al. 2021) and PS15dpn (Wang & Li 2020).  \nThe CSM mass is ∼ 0 . 4 and ∼ 0 . 1 M ⊙ for the s = 0 and s = 2 models, respectively. This is consistent with the CSM mass estimates of 0.016 M ⊙ and < 0 . 1 M ⊙ for SN 2016jbu (Brennan et al. 2022a, 0.016 M ⊙ ) and typical Ibn SNe (Maeda & Moriya 2022), respectively.  \nThe resulting CSM inner radius, R 0 , and the CSM density ρ CSM for the s = 0 model are ∼ 2 × 10 13 cm and ∼ 6 × 10 -13 g cm -3 . For the s = 2 model, R 0 ∼ 7 × 10 13 cm and ρ CSM ∼ 5 × 10 -11 g cm -3 .  \nWhile the MOSFiT framework is flexible, the fundamental parameters of the RD+CSI model are not always directly expressed as standard literature quantities, such as nickel mass and mass loss rate. For the s = 2 model, we can estimate the mass loss rate, ˙ M using the formalism of Ben-Ami et al. (2023) and Kool et al. (2021) as:  \nFollowing Chatzopoulos et al. (2012) approach, we can derive the value of R f , the outer radius of the CSM. For the preferred model ( s = 2), R f ≈ 1 . 5 × 10 14 cm, while for s = 0, R f ≈ 7 × 10 14 cm.  \n˙ M ( r ) = 4 πρ CSM r 2 v w , (2)  \nwhere v w is a steady CSM wind velocity. Additionally, for s = 0 and s = 2 models, we calculate the nickel mass as M Ni = M ej × f Ni and ejecta velocity v ej = √ 10 3 E kin M ej .  \nUsing Eq. 2 to estimate a mass loss rate, and assuming a wind velocity, v w ∼ 400 km s -1 as inferred from the velocity of the absorption minima of the H i lines at late times, together with values of R 0 and  \nρ CSM as resulting from the s = 2 model we obtain a mass loss rate ∼ 2 M ⊙ yr -1 . This mass loss rate is not expected for typical LBV and WR winds (Smith 2017, ∼ 10 -4 -10 -5 M ⊙ yr -1 ), but has also been estimated for the Type IIn SNe iPTF13z (Nyholm et al. 2017), 2017hcc (Smith & Andrews 2020) and the Ibn SN 2019kbj (Ben-Ami et al. 2023). From X-ray observations, mass-loss rates of SN 2009ip and 09ip-like SN 2010mc are estimated below 10 -1 M ⊙ yr -1 (Ofek et al. 2013a,b; Boian & Groh 2018), while Type Ibn SN 2022ablq has an upper limit of 0 . 5 M ⊙ yr -1 (Pellegrino et al. 2024).', '4.3. Determining CSM Structure from the Line Evolution': 'Figure 9. Evolution of the FWHM of H i lines. The broad and narrow emission and the absorption components are color-coded (red, green and magenta). The values of the FWHM are corrected for the resolution of the instrument. If the line was unresolved, an upper limit of the FWHM is estimated as the resolution of the instrument.  \n<!-- image -->  \nFig. 9 presents the evolution of the FWHM of the strongest H i lines, inferred from the spectral modeling described in Sect. 3.2.1 and summarized in Tab. A3. For the first +40 days, the FWHM of the broad and the absorption components remain constant at about 4 000 km s -1 and 800 km s -1 , respectively.  \nBetween +40 and +60 days the FWHM of all absorption components decrease by ∼ 50 %, remaining constant at this level until about 130 days. This suggests a second structural component in the CSM. The FWHM of the narrow emission component of H α and H β is about 600 km s -1 during the first 20 days. Thereafter, it rapidly declines to about 60 km s -1 at which it remains constant from +66 days onwards.  \nFig. 10 shows the evolution of the velocity of the absorption minima of He i , Ca ii , O i and Fe ii and H β , summarized in Tab. A4. It is evident that the velocities of the He i , Ca ii , O i and Fe ii absorption lines are at velocities which are ≳ 200 km s -1 lower than H β . The spectra exhibit a sudden decline of the narrow absorption velocities to about 100 km s -1 past +60 days for all transitions (see Fig. 9 and Fig. 10).  \nFigure 10. Evolution of the velocity of the absorption minima of H i (H α , H β , H γ and H δ ), He i ( λλ 4922 , 5016 , ), Ca ii ( λλ 8498 , 8542 , 8662), Fe ii ( λλ 5169 , 5276 , 5317) and O i λ 8446 lines.  \n<!-- image -->', "4.3.1. Line Fluxes and Ratios and the Curious Case of SN 2021foa's 'Flip-Flop'": "Fig. 11 shows the evolution of the calculated emission line fluxes of the strongest H i and He i lines (summarized in Tab. A5). The line fluxes of H i including Pa β and He i evolve in a similar fashion up to +50 days. They increase up to +12 days after which they stay constant until about +22 days. This agrees with the plateau phase of the r - and g -band light curves, as shown for comparison in the upper panel of Fig. 11. At later epochs ( > +60 days), all He i and H δ line fluxes decrease while most H line fluxes remain constant (or increases, as evident with H α ). Interestingly, at these same epochs, the decline in r -band halts (see Fig. 1). We find that the He i λ 5876 and H α lines reach approximately the same line flux of 10 -13 erg s -1 cm -2 around +20 days, consistent with the findings of Gangopadhyay et al. (2024). This is in contrast to Reguitti et al. (2022) who estimate that the He i λ 5876 flux is about half of that of H α at this epoch (see Fig. A6).  \nThe flux of He i λ 7065 is as strong as that of H β for most of the epochs. Fig. 12 shows the evolution of the  \nFigure 11. Upper panel: Light curves of g and r -bands. Lower panel: Evolution of the total emission flux of H i and He i lines. The shaded region encompasses the epochs of the short photometric plateau at ≈ +15 days.  \n<!-- image -->  \nH α /He i λ 5876 line ratio for SN 2021foa in comparison to other transitional IIn/Ibn SNe. These include SNe 2005la, 2011hw, 2020bqj and iPTF15akq. We also include the H α /He i ratio for the Ibn prototype SN 2006jc and the Type IIn SN 2009ip, SN 2010jl and SN 2016jbu. For consistency, we re-computed the line ratios for these objects similar to those of SN 2021foa (see Sect. 3.2.1). We find discrepancies in the line ratios of SN 2006jc and SN 2011hw a factor of two between Smith et al. (2012) and our measurements (see Sect. C.0.8). We describe our method to determine the line ratios in Sect. C.0.9, and attribute this discrepancy to the mis-estimation of the local continuum surrounding H α and He i λ 5876 lines in previous work.  \nAs shown in Fig. 12, all transitional IIn/Ibn SNe have H/He ratios larger than ≈ 1 at early epochs ( < +10 days). Thereafter, these SNe transit into a more He prominent regime (the 'flip'), with a mean H α /He i ratio of 0.93 ± 0.19 up to about 50 days past peak brightness. After day 50, the line-flux ratio of all other transitional IIn/Ibn SNe drops below one. The evolution of the H α /He i λ 5876 ratio for the transitional IIn/Ibn SNe is clearly different from the classical IIn and 09ip-like events (blue markers), as well as Ibn events such as SN 2006jc (cyan markers). Type IIn SN 2016jbu, SN 2009ip and SN 2010jl are hydrogen dominated (H α /He i > 1) at all epochs. For the Ibn SN 2006jc, the line emission of He i λ 5876 dominates at all epochs. It is the presence of a 'flip' from hydrogendominated to a line ratio of ≈ 1 that truly determines if an object is a member of the class of transitional IIn/Ibn.  \nFigure 12. Flux ratio of H α /He i λ 5876 for the transitional IIn/Ibn objects, SN Type Ibn 2006jc and SNe Type IIn SN 2016jbu, SN 2009ip and SN 2010jl. Each class is color-coded as red (IIn/Ibn), cyan (Ibn) and blue (IIn). The gray band encompasses the mean and one standard deviation of the line flux ratios of SNe 2011hw, 2020bqj and 2021foa within +10 to +55 days. Fluxes are averaged in a 3-day bin.  \n<!-- image -->  \nHowever, uniquely for SN 2021foa even amongst transitional IIn/Ibn 7 , the line-ratio flips back (the 'flop') to hydrogen-dominated after the end of the plateau of the line-flux ratio (day +50), and following the rebrightening of H line emission. This double transition 'flip-flop' behavior highlights the diversity of these transitional supernovae, which ultimately arises from the varied mass loss histories of their massive star progenitors. We discuss the implications of this observation in Sect. 5.", '4.4. Constraining the Photosphere and Dust Emission Properties from the Late-time VLT/X-shooter Spectra': 'The fits of the SN 2021foa VLT/X-shooter spectra at +66, +95 and +129 days (see Sect. 3.3), constrain the radius and temperature of the photosphere, as well as the mass and temperature of the dust. These results are presented in Tab. 2.  \nWe infer the radius of the SN photosphere at +66 and +95 days of ≈ 2 . 0 × 10 14 cm and the temperature ≈ 9000 K. We can compare this to the inferred R BB and T BB (Sect. 3.1.1) at the last epoch, ≈ +35 days, of our  \nTable 2. Fitted parameters for the BB + MBB model to the VLT/X-shooter spectra.  \nNote -Fitting uncertainties are given in 10 -1 units of each column.  \nEXTRABOL light curve modeling. The EXTRABOL results predict that R BB and T BB decline with time (see Fig. 5), and our inferred radius is within the bounds of a linear extrapolation of the EXTRABOL prediction at +35 days. We stress that we do not use any EXTRABOL extrapolation to constrain the late-time photosphere properties. As such a linear extrapolation from +35 days to +66 and +95, while the simplest possible model, is likely unphysical. Indeed, the large inferred photospheric temperature from the VLT/X-shooter spectra indicates a shallower evolution.  \nAt all epochs, the dust mass and temperature inferred for SN 2021foa are ≈ 3 × 10 -5 M ⊙ for a carbonaceous dust composition and ≈ 1 400 K, respectively. However, the inferred low dust temperature indicates that a silicate dust composition can be possible. Adopting a silicate dust composition ( A d = 0 . 2 × 10 4 cm 2 g -1 ), the inferred dust mass increases by about a factor of 5 at similar dust temperature.', '5. INTERPRETING THE OBSERVATIONS AND ANALYSIS RESULTS': "As with all members of the transitional IIn/Ibn, SN 2021foa exhibits characteristics of both Type IIn and Type Ibn supernovae, albeit with key differences to both classes. SN 2021foa exhibits a short ( ∼ 10 days) plateau in the optical light curves about two weeks past peak. However the plateau length is shorter than that of other transitional SNe such as SN 2011hw and SN 2020bqj ( ∼ 50 days, Kool et al. 2021). Similarly, while there are clear similarities between the spectra of SN 2021foa and Type Ibn supernovae after +22 days, neither SN 2021foa nor any other transitional IIn/Ibn SNe follow the Ibn template R -band light curve.  \nHowever, two aspects make SN 2021foa unique:  \n- · SN 2021foa is the first clear example of transitional IIn/Ibn that transitions back - a 'flip-flop'\n- · SN 2021foa exhibited prominent precursor emission about 50-20 days before peak brightness, as is common for SN 2009ip-like objects (Fig. 4)  \nThe photometric resemblance of SN 2021foa with 2009ip-like transients may point to a common progenitor system as already suggested by Reguitti et al. (2022). However, the mechanism of producing the precursor emission and the light curve plateau for SN 2021foa is likely different.  \nOur spectroscopic analysis of SN 2021foa shows that prominent He i λλ 5876 , 7065 and Ca ii IR emission lines have a broad ( ∼ 6000 km s -1 ) component. Contrary, the velocities of the broad components in all H i lines do not surpass 5000 km s -1 .  \nWhile such a velocity agree with the average bulk ejecta velocities of most core collapse SNe (Guti'errez et al. 2017), the velocities inferred from the Balmer lines of SN 2009ip (Pastorello et al. 2013, ∼ 10 000 km s -1 ) and SN 2016jbu (Kilpatrick et al. 2018; Brennan et al. 2022b, ∼ 7000 km s -1 ) are higher. Furthermore, all H i and the He i λ 5016 lines have narrow emission lines ( ∼ 600 -200 km s -1 ). Such narrow lines are characteristic of classical Type IIn and Ibn SNe such as SN 2010jl (Gall et al. 2014) and SN 2006jc, and originate from an extended CSM. However, most He i lines in the spectra of SN 2021foa lack a narrow emission component.  \nFigure 13. Upper panel: Evolution of the velocity of the He i λ 5876 for Type Ib (Liu et al. 2016, Liu16) and Ibn and IIn/Ibn SNe (Hosseinzadeh et al. 2017, H17). Lower panel: Same as upper panel , showing the evolution of the velocity of H α broad (bold) and narrow-absorption (light) components.  \n<!-- image -->  \nFig. 13 visualises the velocity evolution of H α and He i λ 5876 of SN 2021foa in comparison to a sample of Herich stripped-envelope Type Ib and IIn/Ibn and Ibn SNe from Liu et al. (2016) and Hosseinzadeh et al. (2017). Evidently, the velocities measured for SN 2021foa are inconsistent with those measured for these individual SN types.  \nSN 2021foa should be considered a hybrid heliumhydrogen-rich CSM interacting SN, and points to the diversity of the class of transitional IIn/Ibn supernovae. This diversity ultimately arises from the myriad of mass loss histories of massive stars. SN 2021foa's precursor emission together with our detailed photometric and spectroscopic record for SN 2021foa allows us to further constrain the properties of the progenitor environment.", '5.1. The Origin of Strongly Blue-shifted Emission Lines': 'The most intriguing spectroscopic signature of SN 2021foa is the persistent blue-shift of the peak of all H i and He i λ 5016 emission lines past +66 days (Fig. 14). There are three possible origins of these blue-shifted profiles that can either be i) an asymmetric CSM or SN ejecta, ii) an effect of dust extinction, or iii) occultation by the optically-thick photosphere of photons coming from a close line-forming region. We discuss each of these scenarios below. We note that a radiatively accelerated CSM has been proposed to explain blue-shifted asymmetries observed in SN 2010jl (Fransson et al. 2014). Nevertheless, we do not consider that scenario in this work since the blue-shifted profiles in SN 2021foa are persistent at late times, while the acceleration is stronger at peak luminosity (Smith & Andrews 2020).', '5.1.1. Asymmetric CSM': "An asymmetric geometry of either the CSM or the ejecta can lead to strongly blue/red-shifted emission line profiles, as has been suggested for e.g., Type IIn SNe 2009ip (Margutti et al. 2014), 2016jbu (Brennan et al. 2022a) and 2013L (Andrews et al. 2017), Type Ibn SNe 2006jc (Foley et al. 2007), 2019wep (Gangopadhyay et al. 2022) and 2015G (Shivvers et al. 2017). Indeed, the presence of either an asymmetric CSM or ejecta is confirmed by spectropolarimetric observations of supernovae, such as the Type IIn SN 1998S (Leonard et al. 2000), SN 2009ip (Mauerhan et al. 2014), SN 2010jl (Patat et al. 2011) and SN 2017hcc (Kumar et al. 2019), among others (see Bilinski et al. 2020, for the complete Type IIn sample). Unfortunately, we do not have spectropolarimetric observations for SN 2021foa, and cannot conclusively confirm or rule out an asymmetric CSM / ejecta. Never-  \nFigure 14. Left panel: Normalized VLT/X-shooter spectra of H α , H β and Pa β at v ≈ -1800 km s -1 (day +66). Right panel: Normalized VLT/X-shooter spectra (+95) at v ≈ -1800 km s -1 , comparing H α , H β and He i λ 5016.  \n<!-- image -->  \ntheless, a disk-like CSM configuration is possible, with a high density H-rich material moving towards the observer. The CSM from the receding side in this configuration has lower density since the red emission does not re-brights as much as blue's. In this scenario, the blue-shifted emission would be mostly dominated by the interaction of the SN ejecta with the dense CSM material. This scenario has been proposed for Type IIn PTF11iqb (Smith et al. 2015) to explain the highly asymmetric red-shifted H α profiles observed at > +500 days.  \nAlternatively, Thone et al. (2017) posit that the blueshifted profiles observed in Type IIn (SN 2009ip-like) SN 2015bh at > +126 days are explained as shocked emission from a single CSM shell expelled at -2000 km s -1 about -50 days prior to the main explosion event. Given our observations of SN 2021foa, we can compare the the scenario of Smith et al. (2015) and the single shell suggestion for SN 2015bh (Thone et al. 2017). Our observations indicate that the shell would need to be asymmetric in a specific direction towards the observer to produce the blue-shifted profile (see Fig. 14).  \nFurthermore, our observations, particularly the decrease in line-velocities at late-times, suggest that there are multiple H-rich CSM shells out to larger radii, rather than a single shell. Each of those shells can deviate from simple spherical symmetry. Thus, the composite of all these shells would likely lead to line profiles inconsistent with what is observed here for SN 2021foa.  \nWhile the disc-like scenario is more likely than the single shell scenario, a pre-existing H-rich, high-density CSM towards the observer must also show emission at early-times. However, as shown in Fig. 7, the H α line is well modeled by symmetric profiles throughout its evo-  \nution. Therefore, the high-density CSM must be placed further out of a spherical symmetric CSM.  \nThus, while we cannot unambiguously rule out an asymmetric CSM as an explanation for the blue-shifted line profiles without spectropolarimetry, the scenario is unlikely as it requires special fine-tuning of the CSM properties (e.g., density and location) to be consistent with our observations. Next, we consider newly formed dust and occultation by the photosphere as a potential origin of the blue-shifted emission.", '5.1.2. Newly formed dust': 'Newly formed dust located either close to the emission line formation region or within it causes a blueshift of the peaks of emission lines. Simultaneously, a red-blue asymmetry of the emission line profiles due to absorption of photons from the receding side of the SN is produced (Lucy et al. 1989; Gall et al. 2014; Bevan & Barlow 2016) . Thus, the blue side of the line profile remains unaffected, while the red side of the line gets extinguished (red-blue asymmetry).  \nEarly dust formation in a cool dense shell (CDS) causing a red-blue asymmetry of emission lines has been observed in, for e.g., SN 2006jc, SN 2010jl or SN 2017hcc (e.g. Smith et al. 2008; Chugai 2009; Gall et al. 2014; Chugai 2018; Smith & Andrews 2020; Bevan et al. 2020). However, the newly formed dust is composed of both large and small grains (Gall et al. 2014; Smith & Andrews 2020). For the latter case, the blueshift of the emission line profiles and red-blue asymmetry exhibits a measurable wavelength dependence, with bluer emission lines exhibiting larger blue-shifts than redder emission lines. As shown in Fig. 14, the emission line peaks of the normalised profiles of H α , H β and Pa β for SN 2021foa are nearly identical. This rules out newly formed dust in a CDS as the origin of the blueshifted H i and He i emission lines.  \nAdditionally, ejecta dust formation typically starts around one year after explosion when the temperature of the ejecta has cooled to less than about 1600 - 2000 K, which are the sublimation temperatures of silicate and carbonaceous dust, respectively (Gall et al. 2011, and references therein). For SN 2021foa we have a strong temperature constraint from our X-shooter modeling at late times. As shown in Tab. 2, the temperature of the photosphere remains at ≈ 9000 K - a factor of 4 too hot to form dust grains. Thus, ejecta dust formation can be ruled out as the origin of blue-shifted emission line profiles as well as the observed NIR excess emission in SN 2021foa.  \nConsequently, the observed thermal dust emission (see Sect. 4.4) must originate from surviving pre-existing  \ndust at large distances. The amount of dust inferred from our modified BB fits is consistent with dust masses derived at early epochs in other core collapse SNe (Gall et al. 2011; Gall & Hjorth 2018; Gan et al. 2021). Further, the surviving pre-existing dust must be at radii > 10 17 cm, while the emission line forming region is at lower radii ( < 10 15 cm). Hence, also this dust does not cause a blue-red asymmetry of the emission lines. Our observations therefore conclusively rule out newly formed dust as a source of the blue-shifted line profiles.', '5.1.3. Occultation by the photosphere': 'Occultation can be an alternative explanation for the origin of the non-wavelength dependent blue-shifts and red-blue asymmetry of the emission line profiles of SN 2021foa. In such a scenario, the emission from the line-forming region at the receding end of the CSM/SN ejecta is occulted by the optically thick continuum photosphere (Chevalier 1976; Smith et al. 2012; Dessart et al. 2015). This has been observed in non-interacting Type II SNe (Anderson et al. 2014) and suggested for some Type IIn such as SN 2010jl (Fransson et al. 2014), SN 2021adxl (Brennan et al. 2023) or SN 2013L (Taddia et al. 2020).  \nIn the case of occultation, the line-forming region producing the intrinsically symmetric emission lines needs to be very close to the photosphere, else the effect of occultation is minimal as discussed for e.g., SN 2010jl, where the wavelength dependent blue-shifts are likely due to newly formed dust in the CDS (see Sect. 5.1.2).  \nFor SN 2021foa, occultation is likely, because we neither observe a change of the red-blue asymmetry and blue-shifts with either wavelength or time ( ∼ 66 - 129 days). Furthermore, from our BB fits to the VLT/Xshooter data (Sect. 4.4) we find that the photospheric radius remains at around 2 × 10 14 cm, which is similar to the location of the outer radius of the CSM using both MOSFiT models (4.2).  \nHowever, if occultation occurs, the blue-shift and red-blue asymmetry of the emission lines should decrease over time since the photosphere continues receding, blocking less photons with time. Unfortunately, we do not have any data coverage of SN 2021foa beyond +129 days. Thus, we cannot unambiguously determine, if occultation by the photosphere is indeed, a viable explanation for the observed blue-shifts, but it is the most natural scenario that is consistent with all of our observations.', '5.2. Precursor emission of 2009ip-like objects': 'SN 2021foa has shed off most of its hydrogen envelope, as evident from the spectra (Sect. 2.2). Over its lifetime, the progenitor of SN 2021foa created the multiple  \nCSM layers with different velocities as evident from our emission and absorption line analysis (Sect. 4.3, Figs. 9 and 10).  \nThe formation of multiple discrete CSM layers requires episodes of strong mass loss as eruptions or steady winds. Indeed, precursor luminous outbursts have been observed months-to-years before terminal explosion for SN 2006jc, SN 2015bh and 2016jbu (Foley et al. 2007; Thone et al. 2017; Brennan et al. 2022b).  \nAdditionally, our high mass loss rate estimates (2 . 0 M ⊙ yr -1 ) suggest that SN 2021foa suffered from intense mass loss prior to explosion. Our analysis of the ATLAS o -band light curve data from ∼ 5 years prior to the SN 2021foa terminal explosion (Sect. B), shows that SN 2021foa had no eruption brighter than 20 mag over 5 years prior to explosion. Under the s = 2 scheme, MOSFiT gives an outer radius of the CSM of ≈ 1 . 5 × 10 14 cm. Assuming a wind velocity of 400 km s -1 , all the CSM was expelled ∼ 12 years ago, during a period of ∼ half a year. Given the lack of observations of SN 2021foa prior year 2019, we cannot confirm such an event. However, SN 2015bh showed numerous outbursts throughout 20 years before the event B in 2015 (Thone et al. 2017).  \nSeveral possibilities to explain the precursor emission of 2009ip-like objects have been suggested. For SN 2009ip, Mauerhan et al. (2013) argued that the explosion occurred at the onset of event A as a weak Type II SN, while the brighter outburst is mainly powered by SN ejecta-CSM interaction. A similar scenario was suggested for SN 2015bh (Elias-Rosa et al. 2016). In contrast, Pastorello et al. (2013) and Margutti et al. (2014) have proposed that event A of SN 2009ip is an eruption similar to those observed in the years before. Then, the event B is either due to interaction of the material expelled at event A with previous eruptions (colliding shells) or the expanding SN ejecta itself (terminal explosion). Shell-shell interaction is one of the suggested scenario to expain the precursor emission of SN 2015bh (Thone et al. 2017).  \nFor SN 2021foa, our MOSFiT calculations were only performed for event B photometry i.e., assuming that the true CC-SN occurred at the end of event A, and that interaction with one CSM (RD+CSM) is sufficient to reproduce the light curve. Based on these assumptions, the low 56 Ni mass obtained with MOSFiT , typical for SN 2009ip-like SNe, suggest that event B is most likely powered by shock breakout of the CSM + SN ejecta-CSM interaction rather than only radioactive decay.', '6. BUILDING A COMPLETE PICTURE OF SN 2021FOA': "In this section, we summarize the key features from Sect. 5 and build a complete, cohesive model for SN 2021foa.  \n- · SN 2021foa resembles 2009ip-like SNe (Reguitti et al. 2022). In particular, the precursor emission (event A) of SN 2021foa starting about ∼ -50 days prior to the peak of the light curve (event B) and the presence of a short plateau of a few days after the peak are very similar to SNe 2016jbu (Kilpatrick et al. 2018; Brennan et al. 2022b) and 2009ip (Pastorello et al. 2013; Margutti et al. 2014).\n- · For SN 2021foa, the velocity of the minimum of the narrow absorption component of H i , He i , Fe ii , Ca ii and O i lines decreases from ∼ 600 at day +15 to ≲ 300 km s -1 at day +60.\n- · Intriguing and strongly blue-shifted emission lines of H i and He i lines emerge at late times ( > +66 days) in SN 2021foa.\n- · SN 2021foa exhibit a distinctive phase where the He i λ 5876 emission line is as strong as H α . This line-ratio-plateau is observed in the transitional IIn/Ibn SNe.\n- · Unambiguously, SN 2021foa has a 'flip-flop' nature, transitioning from a IIn before peak brightness (flip) to a He-dominated (Ibn-like) SN for about +30 days, and returning to a IIn past +66 days (flop).", '6.1. The Luminosity of Event A and Explosion Date': 'As for most of the 2009ip-like events, the explosion time is uncertain and ultimately dependent on the physical mechanism employed to explain the luminosity of event A. Our spectroscopic observations strongly support a scenario where multiple CSM shells are expelled at different times prior to explosion. We sketch this scenario in Fig. 15. In this scenario, the closest CSM is created by material ejected from prior outbursts or the terminal explosion. Consequently, this material is close to the progenitor system and is rapidly overrun by the forward shock from the SN, powering the emission at event A.  \nThis scenario requires that the explosion date be at event A. For convenience, we place the assumed explosion at ≈ -25 days, consistent with our results from MOSFiT fits (Tab. 1), and the lack of prior outbursts (see Sect. 5.2). We stress that the final scenario is not  \nFigure 15. Proposed scenario for SN 2021foa, sketching the temporal evolution of the forward shock and the photospheric radius. In this scenario, SN 2021foa is embedded in a H-rich CSM composed by multiple shells. The precursor emission (event A) is originated by the interaction of the SN forward shock with an unobserved, close CSM. From day -10 to +35, a second, inner CSM with a characteristic velocity of ≈ 600 km s -1 is observed through the P-Cygni profile of the Balmer lines. The He-rich ejecta, hidden at early epochs, is visible at day +15 after the photospheric radius reaches its maximum. The strong He i lines at this epoch sets up the flip in Fig. 12. A rapid decline of the H α velocity at day +35 marks the phase where the forward shock overran most of the inner CSM. After day +50, the velocity of the H α line remains constant at ≈ 350 km s -1 . We attribute this velocity to an outer, CSM shell. The re-brightening of Balmer lines at late times (+66-+130 days) due to hydrogen recombination of a shocked SN ejecta-CSM (either the close, inner or the outer CSM), gives rise to the flop in Fig. 12. The apparent blue-shift of several lines at late times is either due to an asymmetric CSM or an effect of occultation by the optically thick photosphere close to the line forming region of these lines, located at ∼ 2 × 10 14 cm.  \n<!-- image -->  \ndependent on the precise explosion date. Even in the scenario that the outburst occurred earlier than event A, the key feature of our model for SN 2021foa is that event B itself is powered by SN ejecta-CSM interaction rather than radioactive decay.', '6.2. The Dynamics of the Forward Shock from Event B to the Plateau': 'The broad components of H i , He i and Ca ii IR emission lines suggest an ejecta velocity of ∼ 6000 km s -1 . The maximum velocity inferred from the blue wings of the broad component, indicative of the shock velocity ( v FS ), does not surpass 10 000 km s -1 at all epochs. Therefore, we assume the velocity of the forward shock as v FS = 10000 km s -1 . The first estimate of the photospheric radius is at ≈ 6 × 10 14 cm at day -10. At this epoch, the forward shock is located at ≈ 10 15 cm. Hence, the photosphere lies within a shocked CSM region. The interaction between the SN forward shock and this close CSM might be the underlying powering mechanism of event A. Before peak, we only observe emis-  \nom the H-rich CSM, thus explaining the weak He i lines in the spectra. From -10 up to +30 days, the P-Cygni profile of the Balmer lines suggest a velocity of the inner CSM layer (Fig. 15) to be 600 km s -1 . Around day +15, the photospheric radius reaches its maximum ( ≈ 10 15 cm) while the forward shock is located at ≈ 3 . 5 × 10 15 cm. Around this epoch, we observe the flip ; the He-rich ejecta is becoming visible and He i lines get stronger.', '6.3. Evolution Post-plateau to Late-time': 'After the plateau phase (+25 days), the photospheric radius starts to recede to a radius of about 6 × 10 14 cm. In the meanwhile, the forward shock continues to propagate outwards, reaching a radius of about 5 × 10 15 cm, assuming a constant forward shock velocity. At day +35, we first observe the decline of the velocity of the absorption minimum of the Balmer lines (H i in Fig.10). This decline implies that the forward shock has overrun the inner CSM entirely. Nevertheless, though weaker, the continuous interaction between the SN ejecta and  \nthe CSM still powers the emission of H i and He i lines (Fig. 11). This is evident in the plateau phase of the H α /He i λ 5876 ratio (Fig. 12).  \nPast +50 days, the line fluxes of H i lines, from H δ to Pa β , increase (or stay constant), while He i lines keeps decreasing (Fig. 11). At this epoch, we observe the flop in Fig. 12.', '6.4. Late-time Evolution': 'At +130 days, the velocities of the absorption minima of H i , He i , Ca ii , O i and Fe ii are below ∼ 350 km s -1 . This is consistent with a slow moving outer CSM layer at at a distance of > 1 . 3 × 10 16 cm and thus, has not been overrun yet by the forward shock ( v FS = 10 4 km s -1 ). The likely origin of the low velocity, outer CSM is from slow winds at early stages in the evolution of the progenitor star, while the inner CSM shells are from eruptions or faster winds closer to the explosion. In the multiple-shells scenario, the slowly receding photosphere ( ∼ 2 × 10 14 cm from +60 to +130 days) lies within the inner CSM, close to the line forming region of the intermediate/broad components ( ≈ 4000 km s -1 ). In this scenario, the blue-shifted emission lines emerging after the plateau phase are the result of occultation (Sect. 5.1.3) of the emission line region by a dense and optically thick CSM shell. Finally, at late times, the H-rich inner CSM recombines, explaining the re-brightening of the Balmer lines (Fig. 11). Alternatively, as suggested for PTFiqb (Smith et al. 2015), the interaction between the SN ejecta and a dense, outer CSM approaching to the observer at -2000 km s -1 , could also be viable option to explain the late re-brightening and asymmetry observed in the H i lines.', '6.5. Final remarks': "The early discovery and follow-up of SN 2021foa exhibits the imprints of CSM interaction on the SN's light curve and the evolution of its spectral features. From our light curve analysis (Sect. 4.2), we obtain a CSM mass ≤ 1 M ⊙ , ejecta mass ≈ 8 M ⊙ and mass loss rate of 2 M ⊙ yr -1 for a wind-like ( s = 2) scenario. This mass loss rate is higher than values typical found for Wolf-Rayet or LBV stars, favoured progenitors of Type Ibn and IIn SNe, respectively. From our detailed analysis of the evolution of the line profiles (Sect. 4.3), we conclude that SN 2021foa had a rich mass-loss history, forming multiple CSM shells before the terminal explosion. This CSM configuration, while rare, does share key similarities with other scenarios proposed for interacting SNe such as SN 2015bh (Elias-Rosa et al. 2016). We stressed out that the assumption of a one-shell CSM in MOSFiT is insufficient to explain the precursor emission  \nof SN 2021foa, and incompatible with a CSM composed by multiple-shells. Nevertheless, the overall behaviour of the LC during event B might be approximated well by this assumption. Therefore, the estimated MOSFiT physical parameters might still be valid if the interaction with the close/inner CSM in our scenario is the main contributor to the luminosity of event B.  \nSN 2021foa adds to the number of SNe with truly complex CSM structures which challenges our understanding of extreme mass-loss mechanisms in massive stars, opening up the possibilities of different progenitor scenarios for strongly interacting CC-SNe.", '7. ACKNOWLEDGMENTS': "We thank the S. Raimundo, O. Rodriguez, and the anonymous referee for constructive comments that helped to improve the manuscript.  \nThis work is supported by a VILLUM FONDEN Young Investigator Grant (project number 25501) and by research grants (VIL16599,VIL54489) from VILLUM FONDEN. C.R.A. is supported by the European 1172 Research Council (ERC) under the European Union's 1173 Horizon 2020 research and innovation programme (grant 1174 agreement No. 948381). V.A.V. acknowledges support from NSF through grant AST-2108676. This work is supported by the National Science Foundation under Cooperative Agreement PHY-2019786 (The NSF AI Institute for Artificial Intelligence and Fundamental Interactions, http://iaifi.org/). Parts of this research were supported by the Australian Research Council Discovery Early Career Researcher Award (DECRA) through project number DE230101069. The UCSC team is supported in part by NASA grant 80NSSC20K0953, NSF grant AST-1815935, the Gordon & Betty Moore Foundation, the Heising-Simons Foundation, and by a fellowship from the David and Lucile Packard Foundation to R.J.F. R.Y. received support from a Doctoral Fellowship from the University of California Institute for Mexico and the United States (UCMEXUS) and a NASA FINESST award (21-ASTRO21-0068) postdoctoral fellowship. GN gratefully acknowledges NSF support from AST-2206195 for this work. GN is also supported by NSF CAREER grant AST-2239364, supported in-part by funding from Charles Simonyi, and NSF OAC-2311355, DOE support through the Department of Physics at the University of Illinois, UrbanaChampaign (13771275), and support from the HST Guest Observer Program through HST-GO-16764. and HST-GO-17128 (PI: R. Foley). This work was performed in part at Aspen Center for Physics, which is  \nsupported by National Science Foundation grant PHY2210452.  \nThis investigation is based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme ID 107.22RH (PI C. Gall) and 109.23K3 (PI D. Farias). Research at Lick Observatory is partially supported by a generous gift from Google. Based in part on observations obtained at the Southern Astrophysical Research (SOAR) telescope, which is a joint project of the Minist'erio da Ciˆencia, Tecnologia e Inova¸c˜oes (MCTI/LNA) do Brasil, the US National Science Foundation's NOIRLab, the University of North Carolina at Chapel Hill (UNC), and Michigan State University (MSU). The data presented here were obtained in part with ALFOSC, which is provided by the Instituto de Astrofisica de Andalucia (IAA) under a joint agreement with the University of Copenhagen and NOT. We acknowledge the use of public data from the Swift data archive. This work makes use of data taken with the Las Cumbres Observatory global telescope network. The LCO group is funded by National Science Foundation (NSF) grants AST-1911151 and AST-1911225. This work has made use of data from the Asteroid Terrestrial-impact Last Alert System (ATLAS) project. ATLAS is primarily funded to search for near earth asteroids through NASA grants NN12AR55G, 80NSSC18K0284, and 80NSSC18K1575; byproducts of the NEO search include images and catalogs from the survey area. The ATLAS science products have been made possible through the contributions of the University of Hawaii Institute for Astronomy, the Queen's University Belfast, and the Space Telescope Sci-  \nInstitute. This work is based in part from data obtained at the Infrared Telescope Facility, which is operated by the University of Hawaii under contract 80HQTR24DA010 with the National Aeronautics and Space Administration.  \nYSE-PZ was developed by the UC Santa Cruz Transients Team. The UCSC team is supported in part by NASA grants NNG17PX03C, 80NSSC19K1386, and 80NSSC20K0953; NSF grants AST-1518052, AST1815935, and AST-1911206; the Gordon & Betty Moore Foundation; the Heising-Simons Foundation; a fellowship from the David and Lucile Packard Foundation to R.J. Foley; Gordon and Betty Moore Foundation postdoctoral fellowships and a NASA Einstein Fellowship, as administered through the NASA Hubble Fellowship program and grant HST-HF2-51462.001, to D.O. Jones; and an NSF Graduate Research Fellowship, administered through grant DGE-1339067, to D.A. Coulter.  \nFacilities: Swift (UVOT), LCOGT (Sinistro), Shane (Kast, Nickel), VLT (X-shooter), NOT (ALFOSC), Siding Spring (WiFeS), SOAR (Goodman), IRTF (SpeX), ATLAS, Thacher. Computational facility: High Performance Cluster (HPC) of the University of Copenhagen.  \nSoftware: astropy (Astropy Collaboration et al. 2013, 2018), esoreflex (Freudling et al. 2013), EXTRABOL (Thornton & Villar 2022), SUPERBOL (Nicholl et al. 2017), matplotlib (Hunter 2007), MOSFiT (Guillochon et al. 2018), numpy (Harris et al. 2020), specutils (Earl et al. 2023), YSE-PZ (Coulter et al. 2022).", 'A. HOST EMISSION AT H α': 'The narrow emission component of Balmer lines in the VLT/X-shooter spectra is dominated by the host galaxy emission rather than photo-ionized unshocked CSM. Fig. A1 demonstrates the challenge to perform an accurate reduction at the spectral range covering H α profile due to the contamination of the host emission. Despite special efforts to correctly address this issue, the line flux of H α (Fig. 11) at +95 days might be overestimated.', 'B. PRE-SN EMISSION': "To investigate the pre-explosion activity of SN 2021foa, we followed the analysis described in Wang et al. (2024) using ATClean (Rest et al. 2023, 2024). We obtain the ATLAS forced photometry of o - and c -bands covering ∼ 5 years up to event A at the position of SN 2021foa. Additionally, we performed forced photometry of eight control light curves within a distance of 17 '' to the SN. On average, the flux of these control light curves is expected to be zero. To emphasize the emission of a potential eruption, we defined a figure of merit (FOM) as the signal-to-noise ratio (SNR) convolved with a rolling Gaussian with a fixed kernel size determined by typical timescale of an eruption, 5 < τ G < 100 days (Ofek et al. 2013b; Strotjohann et al. 2021). The same rolling Gaussian was also applied to the control light curves to determine the FOM of the control light curves. By setting up a detection threshold FOM limit , we expect that most of the FOM of the control light curves lies below this limit. If it is not the case, then there are  \nFigure A1. VLT/X-shooter spectral region surrounding H α of SN 2021foa. Upper panel: Spectra at +95 and +427 days past peak. The black and red lines are located at the rest wavelength of H α and [N ii ] λλ 6548 , 6583, respectively. Bottom panel: the 2D spectrum at day +95 within the same spectral region.  \n<!-- image -->  \nunaccounted sources of contamination within the field of SN 2021foa. Any real detection in the pre-SN LC must have a FOM larger than FOM limit . In order to establish a magnitude limit to detect eruptions of a given peak magnitude, we added three simulated Gaussian bumps ( p 1 , p 2 and p 3 ) with increasing amplitudes to one control light curve. A non-detection of any of these peaks translates into an upper limit of a real pre-SN eruption throughout the ATLAS coverage.  \nFor the particular case of SN 2021foa, we first convolved both the SNR of the control and SN LCs with a rolling Gaussian with a kernel size of τ G = 30 days, close to the duration of the precursor emission observed in SN 2021foa. Furthermore, we added three simulated Gaussian bumps, with peak magnitudes of p 1 = 21 . 4, p 2 = 20 . 21, p 3 = 19 . 45 mag, and fixed standard deviation of 25 days to the control LC # 4. Finally, we set up the detection limit as FOM limit = 15. With this same value, Wang et al. (2024) recovered 80% of the eruptions larger than 20 mag for Type Ibn SN 2020nxt. For SN 2021foa, we recover 70% of the eruptions larger than o -band peak magnitude ≈ 20 mag.  \nThe upper panel of Fig. A2 displays the forced photometry ATLAS light curves for SN 2021foa (red) and the control LCs (cyan and blue). We find no signature of any precursor emission associated with SN 2021foa. The bottom panel shows the results of our detection analysis. On average, the control light curves are below FOM limit . For the simulated Gaussians, only p 3 was successfully detected, while p 2 lies slightly below the threshold. Therefore, we can safely conclude that, similar to Wang et al. (2024), no eruption is observed for SN 2021foa with a magnitude greater than or equal to ∼ 20 mag. This translates into a detection limit of absolute magnitude of M o ≈ -13 . 4 mag. However, we cannot discard any pre-SN activity below this magnitude limit.  \nFigure A2. Upper panel: light curve of SN 2021foa (red) and the eight control LCS (blue) before event A (red). Bottom panel: Figure of merit (FOM) of the pre-SN LC (red), eight control LCs (blue), control LC #4 (cyan) and the simulated Gaussians + control LC #4 (green) over ATLAS coverage ( ∼ 3 years) up to event A of SN 2021foa. All these FOM curves were obtained after the convolution of the SNR of each LCs with a rolling Gaussian with kernel size of τ G = 30 days. The simulated Gaussians have peak magnitudes of 21.4 mag, 20.21 mag and 19.35 respectively and a fixed standard deviation of 15 days. The detection limit for our analysis was set up as FOM limit = 15 (dashed line). Red shaded area encompasses the precursor emission of SN 2021foa.  \n<!-- image -->", 'C.0.1. H α': 'Panel D in Fig. A3 shows the evolution of H α between about one week prior r -band peak and to +427 days. The line is characterized by a narrow (FWHM ∼ 600 km s -1 ) and an intermediate component with a FWHM that increases from ∼ 1500 -4500 km s -1 . After maximum light, no major changes are observed aside from a decreasing flux at the red wing of the asymmetric profile. After +60 days, the red-blue asymmetry has flipped, i.e. the blue wing has increased in strength over the red wing. The narrow emission component has faded after about +51 days, revealing the H α host galaxy emission (FWHM ≲ 100 km s -1 ) instead. The apparent re-brightening of the blue peak of H α at +95 days is an artifact stemming from the extraction of the 2D spectrum at the location of H α (see Fig. A1).', 'C.0.2. H β': 'The evolution of H β line is shown in Panel C in Fig. A3. The H β line profile is similar to H α prior to peak magnitude. At later epochs, H β is blended with the strong He i λ 4922 emission. Both H β and He i λ 4922 exhibit narrow absorption lines with absorption velocities from -10 days onwards ( ∼ -600 and ∼ -400 km s -1 , respectively). However, the narrow absorption line of He i could be partially associated with Fe ii λ 4924 line of multiplet 42. Similar to the H α line, the line flux of the H β +He i λ 4922 complex reaches a maximum around +15 days and decreases between +28 and +66 days. After this epoch, H β shows a strong blue-shifted emission component while the velocity of the minimum of the absorption component reaches ∼ 300 km s -1 . The emission of He i λ 4922 is weak in comparison to H β , allowing us to disentangle both lines at these late epochs. Similar to H α and H β , He i λ 4922 also exhibits a blue-shifted emission component.', 'C.0.3. He i λ 5016': 'In panel G (Fig. A3), we show the evolution of He i λ 5016 line. From -10 days onwards, this line shows a P-Cygni profile, with the absorption minimum at a velocity of ∼ -400 km s -1 . Similarly to He i λ 4922, the narrow absorption component can be associated with Fe ii λ 5018 of multiplet 42. The extension of the wings of the broad emission component may indicate a maximum velocity of 4000 km s -1 (bulk velocity of ∼ 3000 km s -1 ). After maximum  \nFigure A3. Evolution of the the most prominent continuum-subtracted H i and He i lines of SN 2021foa. Dashed, vertical lines correspond to zero velocity. Solid, vertical lines in panels C, E, G and H correspond to the velocities at the wavelengths of He i λλ 4922 , 5048 , 10830 and Na i D, respectively.  \n<!-- image -->  \nlight, the shape of the profile is broad and boxy-like. Assuming one component, the FWHM of this complex is about 3600 km s -1 . This is consistent with what is seen for H α . From +95 days onward, the flux at the red portion of the He i λ 5016 emission line profile rapidly decreases. The region might be affected at early times by the emission of He i λ 5048 ˚ A.', 'C.0.4. He i λ 5876': 'Panel H in Fig. A3 displays the evolution of He i λ 5876, which is the strongest of all He i lines in the entire VLT/Xshooter spectral wavelength range. There is no indication of a narrow P-Cygni profile as in other He i λλ 4922 , 5016 lines. However, there is conspicuous, red-shifted absorption at all epochs up to +129 days, that likely is the Na i D  \ndoublet from interstellar material along the line of sight. The emission line profile appears symmetric at early epochs (two weeks past maximum) in comparison to H i lines. Thereafter, the line develops a red-shifted peak. From about +66 days onward it appears symmetric again, unlike other He i and H i lines.', 'C.0.5. He i λ 7065': 'In panel I (Fig. A3), we show the evolution of He i λ 7065. In analogy to He i λ 5876, the emission line profile of He i λ 7065 is of boxy-like shape with no evident narrow P-Cygni profile. Additionally, at all epochs the line profile exhibits a blue shoulder at around -3000 km s -1 . The origin of both emission features is unclear but may be due to another element.', 'C.0.6. NIR lines': 'The NIR spectra exhibit emission lines of He i λ 10830 + Pa γ , Pa β and He i λ 20581. H i λ 18751 (Pa α ) is detected, but it coincides with a telluric region. Panel E in Fig. A3 shows the He i λ 10830 + Pa γ line complex, which is dominated by He i λ 10830. We find that the absorption line at about -4000 km s -1 (with reference to Pa γ ) must be attributed to He i λ 10830. This, because it is unlikely that a narrow Pa γ absorption at a velocity of about -4000 km s -1 with a FWHM of only about 800 km s -1 exists. Furthermore, the velocity of the absorption minimum, if associated with He i λ 10830, remains nearly constant at about -600 km s -1 at all epochs up to +129 days.  \nPanel F in Fig. A3 displays the evolution of Pa β . At -8 days, the emission line has a FWHM of only about 2000 km s -1 and thus, is narrower than other H i lines at that epoch. On top of that is a narrow P-Cygni profile with a FWHM of about ∼ 1000 km s -1 and an absorption component with a minimum at ∼ -500 km s -1 . From +66 to +129 days, the Pa β line profile is nearly identical to the optical H i lines. However, Pa β does not show a narrow absorption component, likely because it coincides with telluric lines at that position.', 'C.0.7. Other lines profiles': 'Figure A4. Left panel : Evolution of the spectral region surrounding Ca ii H&K. The P-Cygni profile of the Balmer lines H ζ and H ϵ are prominent. Middle panel : Evolution of the spectral region surrounding the pseudo-continuum produced by the blended Fe ii lines. Right panel : Evolution of the Ca ii NIR + O i λ 8446. All the regions in each panel are continuum-subtracted.  \n<!-- image -->  \nAnother complex spectral structure is prominent in the wavelength range of Ca ii H&K ( λλ 3934 , 3968). Left panel of Fig. A4 shows the evolution of this structure at +66, +95 and +129 days. The absorption complex at the position of Ca ii H&K exhibit multiple components, some appear blue-shifted with respect to the SN redshift. Strong absorption features, potentially associated with these Ca ii lines are at about -350 km s -1 . Since there is no indication of any  \nemission or absorption lines at these wavelengths in the +427 day spectrum, it is likely that the absorption complex originates from different CSM layers around the supernova. However, the mismatch between the low velocities from the absorption trough of the Ca ii IR and the Ca ii H&K counterparts is intriguing. This apparent discrepancy could be solved considering that Ca ii H&K could be misidentified due to the strong blending with He i λλ 3936 , 3965 (vertical red lines in left panel).  \nSeveral Fe ii (multiplet 42, 48 and 49) lines are shown in middle panel in Fig. A4 at +66 to +427 days. Strong blending of Fe ii forest lines can create a pseudo-continuum as observed bluewards of 5700 ˚ A. This is not unusual for Type Ibn and other interacting SNe (Pastorello et al. 2015b). Strong P-Cygni absorption components at about < -300 km s -1 are detected for all Fe ii lines.  \nRight panel in Fig. A4 displays the evolution of the Ca ii λλ 8498 , 8542 , 8662, possibly blended with O i λ 8446. While it is difficult to disentangle these four lines we find that the bulk velocity of the Ca ii + O i complex does not surpass 6000 km s -1 . This is the maximum velocity of the red wing of Ca ii λ 8662. Furthermore, at all epochs past +15 days, narrow P-Cygni absorption is observed at a velocity of about -400 km s -1 , which continuously decreases to -200 km s -1 between +40 and +129 days.', 'C.0.8. On SN 2006jc and SN 2011hw': 'In Fig. A5 we show the decomposition of the line profiles of He i λ 5876 ( left column) and H α ( right column) of three different Type Ibn SNe: the prototype SN 2006jc, the transitional SN 2011hw and SN 2021foa two weeks after r -band maximum. It is clear that our decomposition is in very good agreement with the total line-flux of each profile. Furthermore, in contrast to Smith et al. (2008); Smith et al. (2012), we show that our decomposition correctly deblends the H α and He i λ 6678.  \nIn Fig. A6, we show the flux-calibrated H α (red) and He i λ 5876 (blue) profiles of SN 2006jc, SN 2011hw and SN 2021foa at ≈ 15 ( left panel ) and ≈ 40 ( right panel ) days r -band maximum. These two epochs encompass the flux-line ratio plateau observed in Fig. 12 for transitional IIn/Ibn SNe. For SN 2006jc, the flux-line of H α is smaller than that of He i in both epochs. This difference is not observed in either SN 2011hw and SN 2021foa, where the flux-line ratio is ≈ 1 for both epochs. Fig. A6 proves that the values obtained in Reguitti et al. (2022) of the line-flux ratio of H α /He i λ 5876, ≈ 0 . 5, are not consistent with our observations.', 'C.0.9. Line decomposition in transitional objects': 'For SN 2005la and iPTF15akq, a careful analysis was done to take into account the fact that these sources exhibit strong, broad P-Cygni absorption profiles for both H and or/He at early epochs. However, we note that this absorption was not accounted for in the line fluxes calculations.  \nFor SN 2011hw, SN 2005la, SN 2020bqj and SN 2006jc, we deblended the H α from He i λ 6678 line by fitting two different Gaussian profiles to each line. This results in line fluxes of H α and He i λ 5876 for SN 2006jc and SN 2011hw that are discrepant by a factor 2 from measurements in the literature (Smith et al. 2008; Smith et al. 2012) using different methods. In Fig. A5 we show that our decomposition of the line profiles of H α and He i λ 5876 light recovers the total line-flux at ≈ two weeks after maximum.  \nFigure A5. Decomposition of continuum-subtracted line profiles of He i λ 5876 ( left column ) and H α ( right column ) of Type Ibn SN 2006jc ( upper row ), transitional Type IIn/Ibn SN 2011hw ( middle row ) and SN 2021foa ( lower row ) at two weeks after r -band maximum. Solid and dashed lines represent the fits to the emission components of the line profiles, while dotted lines correspond to the fit to any absorption trough.  \n<!-- image -->  \nFigure A6. Comparison of continuum-subtracted line profiles of H α (red) and He i λ 5876 (blue) of Type Ibn SN 2006jc ( upper row ), transitional IIn/Ibn SN 2011hw ( middle row ) and SN 2021foa ( lower row ) at ≈ 15 ( left column ) and ≈ 40 ( right column ) days past r -band maximum.  \n<!-- image -->  \nTable A1 . Photometry of SN2021foa  \nNote -Full table available in electronic form.  \nTable A2 . Spectroscopy of SN2021foa', 'T able A3 . F ull Width at Half-maxim um of H i and He i lines of SN 2021foa.': '<!-- image -->', 'T able A4 . Fitted absorption minim um of sev eral sp ectral lines of SN 2021foa.': 'Note -V alues of the absorption minim um are rep orted in units of -100 km s -1 . Uncertain ties are estimated as c · ∆ λ D /λ 0 , with c the sp eed of ligh t, ∆ λ D the disp ersion of the instrumen t in T a b. A2 , and λ 0 the cen ter of the Gaussia n profile.  \n<!-- image -->  \nTable A5 . Line fluxes of H i and He i lines of SN 2021foa.  \nNote -Values are reported in units of 10 -17 erg s -1 cm -2 . Line fluxes of He i λ 5876 correspond to the integration of the data over ≈ ± 5000 km s -1 with respect to 5876 ˚ A.', 'REFERENCES': '- -. 2018, MNRAS, 481, 3643, doi: 10.1093/mnras/sty2386  \nVillar, V. A., Berger, E., Metzger, B. D., & Guillochon, J. 2017, ApJ, 849, 70, doi: 10.3847/1538-4357/aa8fcb Wang, Q., Goel, A., Dessart, L., et al. 2024, MNRAS, 530, 3906, doi: 10.1093/mnras/stae1038  \nWang, S.-Q., & Li, L. 2020, ApJ, 900, 83,  \ndoi: 10.3847/1538-4357/aba6e9'}