The luminous red nova AT 2018bwo in NGC 45 and its binary yellow supergiant progenitor

Luminous Red Novae (LRNe) are astrophysical transients associated with the partial ejection of a binary system's common envelope (CE) shortly before its merger. Here we present the results of our photometric and spectroscopic follow-up campaign of AT2018bwo (DLT18x), a LRN discovered in NGC45, and investigate its progenitor system using binary stellar-evolution models. The transient reached a peak magnitude of $M_r=-10.97\pm0.11$ and maintained this brightness during its optical plateau of $t_p = 41\pm5$days. During this phase, it showed a rather stable photospheric temperature of ~3300K and a luminosity of ~$10^{40}$erg/s. The photosphere of AT2018bwo at early times appeared larger and cooler than other similar LRNe, likely due to an extended mass-loss episode before the merger. Towards the end of the plateau, optical spectra showed a reddened continuum with strong molecular absorption bands. The reprocessed emission by the cooling dust was also detected in the mid-infrared bands ~1.5 years after the outburst. Archival Spitzer and Hubble Space Telescope data taken 10-14 years before the transient event suggest a progenitor star with $T_{prog}\sim 6500$K, $R_{prog}\sim 100R_{\odot}$ and $L_{prog}\sim 2\times10^4L_{\odot}$, and an upper limit for optically thin warm (1000 K) dust mass of $M_d<10^{-6}M_{\odot}$. Using stellar binary-evolution models, we determined the properties of binary systems consistent with the progenitor parameter space. For AT2018bwo, we infer a primary mass of 12-16 $M_{\odot}$, which is 9-45% larger than the ~11$M_{\odot}$ obtained using single-star evolution models. The system, consistent with a yellow-supergiant primary, was likely in a stable mass-transfer regime with -2.4


Introduction
Over half of all stars in our Universe are born as binary or multiple systems, and the fraction increases above 70% among the most massive stars (Moe & Di Stefano 2017;Sana et al. 2012). Studies suggest that 30% of them have already interacted with their companions via mass transfer and 10% have merged (de Mink et al. 2014). The merger process starts when one of the stars overfills its Roche lobe (Roche lobe overflow; RLOF), ini- Table 1 is only available in electronic form at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsweb.ustrasbg.fr/cgi-bin/qcat?J/A+A/ VENI fellow e-mail: n.blagorodnova@astro.ru.nl tiating an unstable mass transfer toward its companion. This process may culminate with both stars orbiting inside a shared non corotating gaseous layer, called the common envelope (Paczynski 1976;Ivanova et al. 2013a). The less-massive component quickly spirals into the gravitational well of the primary star, transferring the angular momentum of the binary to the envelope. At the termination of this phase, part of the envelope is ejected, leaving either a more compact binary or a fully coalesced star. This phenomenon is key to our understanding of the following: stripped-envelope or Type Ia supernova progenitors, cataclysmic variables, X-ray binaries, and gravitational-wave sources (Iben & Livio 1993;Nelemans et al. 2000;Belczynski et al. 2016;Vigna-Gómez et al. 2020). Mergers are a likely explanation for the existence of blue and red stragglers (McCrea 1964;Hills & Day 1976;Schneider et al. 2015;Britavskiy et al. 2019) and rapidly rotating spotted FK Com stars (Bopp & Stencel 1981). High-mass merger products might include magnetic stars (Ferrario et al. 2009;Schneider et al. 2019), the progenitor of supernova SN 1987A (Morris & Podsiadlowski 2007), and luminous blue variables such as η Car (Mauerhan et al. 2013).
Observationally, the merger process is thought to be accompanied by a peculiar type of outbursts called luminous red novae (LRNe). These events display peak luminosities lying between those of novae and supernovae and temperatures from 2000 − 10, 000 K during peak brightness. Unlike novae, the outburst temperature quickly cools down as the envelope expands with ∼100 km s −1 velocity, shifting into near-infrared (NIR) and mid-infrared (MIR) wavelengths and triggering dust formation.
Similar to the cold shells of Mira long-period variables, the ejected envelopes of LRNe become sites of formation of water vapor and metal oxides such as titanium oxide (TiO) and vanadium oxide (VO). For example, NIR spectroscopy of the Galactic stellar merger V838 Mon at early times showed strong absorption bands corresponding to water, AlO, and CO (see Banerjee et al. 2005). These same components were observed in emission in another merger, V4332 Sgr, ∼10 years after its outburst (Banerjee et al. 2015), revealing the slow condensation of highly under-oxidized SiO, AlO, Fe, and Mg dust grains in the outflow.
New discoveries by ongoing all-sky transient surveys and archival searches for LRN-like events have effectively doubled the number of known systems in the last couple years (Pastorello et al. 2019a,b;Cai et al. 2019;Stritzinger et al. 2020;Pastorello et al. 2021a,b). However, limited by the lack of archival high-resolution data on their host galaxies, the progenitors of LRNe remain largely unexplored. To date, reliable multiband detections of the systems in quiescence have only been reported for the following five LRNe outside of the Milky Way: M31-LRN2015 (Dong et al. 2015;Williams et al. 2015), M101 OT2015-1 (Blagorodnova et al. 2017), SNHunt 248 (Mauerhan et al. 2018), AT 2019zhd (Pastorello et al. 2021a), and AT 2020hat (Pastorello et al. 2021b).
In this work, we present results of our follow-up campaign in the optical and NIR of the LRN AT 2018bwo, and we interpret its progenitor system using binary stellar-evolution models. This study helps us better understand the evolutionary stage of the binary moments before the dynamical spiral-in phase and relate it to the outburst properties and its late-time evolution. This is the first time that a LRN progenitor has been studied in agreement with its binary nature.  (Wyatt et al. 2018a;Sand 2018). The last nondetection 6 days earlier with a limiting magnitude of Clear=19.55 shows that the source was detected soon after the outburst onset. Subsequently, the transient was also reported by two other surveys: ATLAS (ATLAS 18qgb; Tonry et al. 2018b) and Gaia Science Alerts (Gaia 18blv; Hodgkin et al. 1 https://wis-tns.weizmann.ac.il/ 2013, 2021). Although the source was classified on TNS as an Intermediate Luminosity Red Transient (ILRT) on 2018-05-23 18:00:01 by the extended Public ESO Spectroscopic Survey for Transient Objects (ePESSTO; Smartt et al. 2015), its plausible classification as a LRN was not ruled out (Clark et al. 2018).

Observations
The host galaxy NGC 45 has a spectroscopic redshift of z = 0.00156 (Springob et al. 2005). The distances to NGC 45 reported in literature span a wide range of values from 4.47 Mpc (Bottinelli et al. 1985) to 13.90 (Willick et al. 1997). Given the uncertainty, we choose the method that provides the most accurate individual distances, based on the tip of the red giant branch (TRGB). In this work, we adopt one of the most recent measurements of D L = 6.64±0.10 Mpc (Tully et al. 2013), corresponding to a distance modulus of µ = 29.11 ± 0.10.
The foreground Milky Way extinction toward NGC 45 corresponds to A V = 0.058 (Schlafly & Finkbeiner 2011), which translates into a reddening of E(B − V) = 0.0187 for a nominal extinction law. The average extinction within the host galaxy was derived by Mora et al. (2007) by modelling young stellar clusters at different metallicities. The best representation of cluster colors was reached for a host-galaxy metallicity between Z = 0.006 and 0.008 and a mean reddening of E(B−V) = 0.04. In our work we assume this value to be an upper limit and adopt an extinction value of A host V = 0.03 +0.09 −0.03 , motivated by the modelling of the progenitor SED in Sect. 4.1. Thus, the total fiducial reddening used through this work is E(B − V) = 0.029 (A V = 0.089).

Optical and IR photometry
Our follow-up campaign started soon after the transient was discovered and lasted up to 134 days in the infrared (IR). Optical photometry for AT 2018bwo was obtained using the Las Cumbres Observatory Global Telescope Network (LCOGT; Brown et al. 2013) under the program NOAO2018B-008 (PI T. Kupfer). Once the transient had faded, at 250 days after discovery, we obtained the reference images of the field under the program NOAO2019A-011 (PI N. Blagorodnova). For each science image, we ran a custom adaptation of the image-subtraction algorithm ZOGY (Zackay et al. 2016) and fit the field PSF to the residual image to obtain the photometry of the transient in the difference image. In addition to our photometry, we included measurements from the time domain-surveys ATLAS, Gaia Science Alerts, and DLT40. The measurements from the ATLAS difference-image detection pipeline (Smith et al. 2020;Tonry et al. 2018a), reported in the orange o and cyan c filters, have been added with a binning of 1 day. The Gaia data are reported in the Vega system. The photometric data are provided in Table 1 and Fig. 1 shows the light curve of all combined photometric measurements.
Follow-up NIR photometry for AT 2018bwo was obtained with the Wide-Field Infrared Camera on the Palomar 200-inch Hale telescope (WIRC; Wilson et al. 2003), and the adaptiveoptics instrument NIRC2 on Keck 2 (LGS AO; Wizinowich et al. 2006). The data were reduced using custom-developed pipelines, which first obtained individual frames using standard routines for dark subtraction and flat fielding. The individual dithered images were then combined for each filter. Due to the large field of view of the WIRC camera, the zeropoint for each observation night was estimated from field stars contained in the 2MASS (Skrutskie et al. 2006) catalog. For NIRC2, the zeropoints were estimated using two Hubble Space Telescope (HST) IR photometric standard stars, one of them observed at the beginning of the night and the other at the end.   MIR data were obtained with the Infrared Array Camera (IRAC; Fazio et al. 2004) on board the Spitzer Space Telescope (SST; Werner et al. 2004;Gehrz et al. 2007) as target of opportunity under the SPIRITS program (Kasliwal et al. 2017). We used archival Spitzer images taken on MJD 54460.228 ( 10.4 years before discovery) to run the SPIRITS difference-imaging pipeline on the photometry, taken 134 days after discovery.
To complement our MIR photometry of AT 2018bwo, we searched for archival data from the cryogenic Wide-Field Infrared Explorer mission (WISE; Wright et al. 2010) and the Reactivation survey (NEOWISE-R; Mainzer et al. 2011Mainzer et al. , 2014. Initially, the cryogenic WISE mission covered the whole sky in the 3.4, 4.6, 12, and 22 µm (W1−W4) from UTC January 2010 until UTC August 2010. Since its reactivation in December 2013, the NEOWISE-R space mission has been scanning the sky in the W1 and W2 infrared bands with an approximate cadence of half a year. The latest NEOWISE 2020 Data Release comprises observations between December 13, 2018 and December 13, 2019 UTC.
The cadence of NEOWISE provides ∼12 exposures taken within one or two consecutive days. To obtain higher signal-tonoise images, we initially coadded all the exposures obtained in the same "observing block" using the ICORE tool (Masci & Fowler 2009). A visual inspection at the location of the transient showed that no emission was present back in 2010, but a new source appeared in all the exposures obtained after the detection of the outburst in 2018 (see Fig. 2). To compute the magnitude of this source, we performed aperture photometry following the instructions detailed in the ICORE manual. Unfortunately, the underlying region has non-negligible emission from the unresolved stellar population in the host galaxy, which can impact the accuracy of our measurements. Nevertheless, the agreement with our follow-up Spitzer photometry, which was performed using a difference-imaging technique, shows that this impact is minor.
Our multiband light curve (see Fig. 1) shows that AT 2018bwo quickly faded in optical magnitudes, while its IR counterpart remained bright for longer, similar to other LRNe. In the optical, the discovery epoch corresponded to a peak magnitude of M clear = −12.9 ± 1.1 mag, as reported by the DLT40 survey (Wyatt et al. 2018b). However, there appears to be a systematic offset between DLT40 photometry and our difference imaging measurements (DLT40 team, private communication). Consequently, we decided to use our own data to derive the peak magnitudes, as we find that our photometry is consistent with the ATLAS survey and with the synthetic magnitudes obtained from the flux-calibrated ePESSTO classification spectrum taken at +1 day (see open markers in Fig. 1). We compute our peak magnitudes from the photometric measurements obtained during the initial plateau phase, lasting 20 days after discovery. We compute average absolute magnitudes of M r = −10.97 ± 0.11 mag and M i = −11.40 ± 0.04. Toward the end of the plateau, at day +34, we see a possible short secondary peak in the light curve in the i and r bands, although our cadence after the re-brightening is insufficient to characterize this feature. Our best approach for obtaining an additional photometric epoch consists of using the spectrum at +50.6 days in order to derive the synthetic g, r, iband photometry. In order to mitigate flux losses, we scale the spectrum to match the broadband Gaia G-band magnitude, obtained four days later.

Spectroscopy
Spectroscopic follow-up for AT 2018bwo was done in both optical and NIR wavelengths. Although optical spectroscopy for LRNe has been covered in detail in the literature, this is the first study to present NIR data on an extragalactic LRN. The full log of the observations included in this work is provided in Table 2, and the data are available through the online WiseRep repository (Yaron & Gal-Yam 2012).
The optical AT 2018bwo spectroscopic observations were conducted with the Low Resolution Imaging Spectrometer (LRIS; Oke et al. 1995) on the Keck I 10 m telescope and with the Double Spectrograph (DBSP; Oke & Gunn 1982) on the Palomar 200-inch Hale telescope (P200). For LRIS we used a combination of grism 400/3400 and grating 400/8500, and a long 1 -wide slit, providing a resolution of ∼7 Å. For DBSP we used a combination of blue grating 600/4000 and red grating 316/7500 with a long 1.5 -wide slit, yielding a resolution of ∼10.5 Å. Finally, the PESSTO spectrum retrieved from TNS was observed with the EFOSC2 spectrograph mounted on the New Technology Telescope (NTT) in La Silla. Grisms #11 and #16 were used, with a resolution of 15.8 and 16 Å..
The LRIS spectrum was reduced using the automated reduction pipeline in IDL lpipe 2 , developed by D. Perley. The DBSP spectrum was reduced using the custom-developed DBSP pipeline pyraf-dbsp (Bellm & Sesar 2016).
For NIR spectroscopy, we used the recently commissioned prism cross-dispersed Near-Infrared Echellette Spectrometer (NIRES) on Keck I 3 . For this object, the data correspond to a single epoch taken at 4 months after discovery. The NIRES data was reduced using NSX: NIRES Spectral eXtraction pipeline 4 . The data were corrected for telluric absorption using A0 V-type stars at similar airmass and sky position. Two different calibration stars were observed, one at the beginning and one toward the end of the exposure. The flux was calibrated using the standard flux-calibration stars.

SED modelling
The evolution of the outburst luminosity, temperature, and radius was estimated using black-body fits to the optical and infrared SEDs. The optical SEDs were obtained from LCO photometry for the epochs when all g, r, i bands were available. The infrared SEDs were obtained when all J, H, K or [3.6] and [4.5] bands were available. At +51 days the best fit was performed on the flux calibrated Keck/LRIS spectrum. The fluxes were corrected for both Galactic extinction of A V = 0.058 and extinction in the host of A V = 0.03, in agreement with our best estimate for the progenitor (see Sect. 4.1). When fitting the spectrum, we excluded the regions with stronger molecular absorptions. To estimate the fit uncertainties, the initial and final wavelength of the fitted region of the spectrum were shifted in steps of 100 Å within a 500 Å window. The results of our analysis are summarized in Table 3, and the LRN evolution is shown in Fig. 3, along with other well-studied LRNe.
The luminosity for AT 2018bwo during the plateau is ∼ 3 × 10 6 L (∼ 1.1 × 10 40 erg s −1 ), which is a factor of 100 brighter than the progenitor. The total duration of the luminosity plateau is ∼70 days, which is longer than observed in the optical wavelengths alone. Similar to other LRNe, the initial increase in luminosity during the plateau is followed by a sharp decline.
At early times, the photospheric temperature of AT 2018bwo is among the lowest in the LRN sample, at only 3300 ± 70 K. One other similar transient is M101 OT2015-1, which showed a comparable range during its second peak. However, the temperature of its first peak was likely hotter, in agreement with the values inferred for other bright LRNe (Pastorello et al. 2019a). Around ∼70 days after discovery, there is a fast increase in the temperature of the emission, along with the shrinking in the photospheric radius. The most likely explanation is that the most ex-   ternal layers become optically thin due to the expansion, and the photosphere moves inwards, where temperatures are higher.

Optical
The optical spectral evolution of AT 2018bwo is shown in Fig. 4. Once corrected for the redshift NGC 45 (z=0.00156), in all three spectra we identify a blueshift of 270 km s −1 for the major absorption lines, such as the Ca ii triplet, K i λλ 7665, 7699, Ba ii, and [Ca ii]. In addition, on top of the LRN spectrum, we also detect [O ii] λλ3726, 3729 emission at the redshift of the host galaxy. Therefore, we attribute this line to the underlying H ii region where the star is located. The sky spectrum also shows an underlying Hα emission at the rest wavelength of the host galaxy, which is offset from the peak of the Hα line for the LRN.
The presence of such background contamination may have an impact of some of our measurements (e.g., the P-Cygni profile identified at −800 km s −1 may be due to oversubtraction).  Table 3. Best fit values and their 1σ uncertainties for the SED of AT 2018bwo at different epochs post discovery. In the fits, we fixed the amount of extinction to a total value of A V =0.089.
The initial PESSTO classification spectrum, taken at +1 day after discovery, reveals a cool photosphere, resembling a ∼3900 K star. Analogous to the stellar spectra, the nova is dominated by absorption lines for low-ionization elements, such as Ba ii, [Ca ii], Ca ii IR triplet, and low-ionization states of Fe ii and Sc ii. This signature also resembles the early-time spectrum of V838 Mon, taken at +36 days after discovery (see comparison in Fig. 4). Although we cannot identify any emission coming from Balmer transitions, the absorption for Hα is not as strong as for an equivalent stellar type.
The Hα line shows signs of a possible P Cygni profile at all three epochs, as displayed in  corresponds to ∼800 km s −1 . However, we interpret this result with caution, as it could be an artifact from oversubtracting the galaxy background. The line also shows a bi-modal profile, likely associated with absorption in a colder expanding shell. At +47.1 days, we derive an instrument-corrected FWHM of 450 ± 18 km s −1 and 580 ± 25 km s −1 for the later epoch at +51.2 days. The approximate velocity of the absorption component for this last epoch is located at −150 km s −1 with a FWHM of ∼150 km s −1 . Similar profiles were observed for other Galactic LRNe, such as V4332 Sgr (Martini et al. 1999) and V838 Mon (Mason et al. 2010), and the extragalactic outburst of M101 OT2015-1 (Blagorodnova et al. 2017).
Toward the end of the plateau, the last two spectra are already marked by the early appearance of molecular TiO and VO bands. The strength of the bands increases as the temperature of the spectrum drops to ∼2800 K. As a comparison, molecular absorption was also detected in M101 OT2015-1 as early as +22 days after the secondary peak. For the brightest LRNe with longer plateaus, molecular absorption was detected at later stages, e.g., AT2018jfs (Pastorello et al. 2019a) or NGC4490-2011OT1 (Pastorello et al. 2019b).

Near-Infrared
Fig . 6 shows the NIR spectrum of AT 2018bwo at +103 days, along with a stellar-atmosphere model and the spectrum of a cool giant star. The model was obtained from the Phoenix stellar spectral atlas (Allard & Hauschildt 1995) and corresponds to a star with solar metallicity, log g = 0, and a temperature of 2000 K, which is the lowest available temperature in their grid. The observed spectrum (HV 12149) was obtained from the DR2 of the X-shooter Spectral Library (XSL; Gonneau et al. 2020) and corresponds to a very luminous asymptotic-giant-branch (AGB) star, classified as M8.5 II (González-Fernández et al. 2015). The classification of AT 2018bwo at late times agrees with a late M-type giant star. The most prominent features are strong water vapor absorption bands around 1.1, 1.4, and 1.9 µm, comparable to the ones detected in atmospheres of cool stellar types (Lancon & Rocca-Volmerange 1992). These bands appear particularly noticeable in variable stars as compared to more static giants, as pulsations contribute to extending the stellar atmosphere (Lançon & Wood 2000). Similar absorption troughs were also observed to appear at later times (∼4 months after detection) in the spectrum of V838 Mon (Banerjee & Ashok 2002), when the temperature of the emission dropped down to 2200−2400 K. Fig. 7 shows the continuum-subtracted spectrum of both the LRN and the cool AGB star. Molecular TiO and VO absorption features are detected around 1.05, 1.1, and 1.25 µm, in agreement with their identification in the optical spectrum. Due to the low gravity in the atmosphere of the LRN, the disappearance of molecular TiO and VO bands due to condensation will be pushed to extremely low temperatures, so that despite the temperature of the object being consistent with an L0-type star, we cannot classify the spectrum as such based on the strength of molecular absorption alone. However, early components of dust, such as AlO and Al 2 O 3 condense at higher temperatures than TiO. Therefore, their presence in the spectra of LRNe is not unexpected. In addition, we identified lines of CN in the J-band part of the spectrum. Similar to other cool giants, the H-and Kband spectra also show strong absorption from 12 CO and 13 CO bandheads 5 .  Fig. 7. Continuum subtracted NIR spectra for AT 2018bwo (black) and a comparison star HV12149 (blue). Common IR molecular and atomic transitions are marked. The three panels correspond to J, H, K spectra. Low signal-to-noise regions due to atmospheric extinction have been excluded. The coincidence of absorption features between the two spectra show that, albeit unidentified, most of the lines in the AT 2018bwo are real.

Progenitor SED modelling
To measure the photometry of the progenitor of AT 2018bwo, we searched the HST and SST archives. The progenitor location was imaged by the Advanced Camera for Surveys Wide Field Channel (ACS/WFC) onboard HST on 2004 June 1, nearly 14 years before the detection of the merger outburst. The field was also imaged with SST on 2007 July 06 (PI Robert Kennicutt), about three years after the HST observations. The precise location of the progenitor was identified by performing an astrometric registration of one of our follow-up NIRC2 adaptive-optics observations with the HST/ACS F814W frames. Our analysis shows an association within 1σ between the location of the transient and a bright progenitor star in the HST images, shown in Fig. 2.
The photometry for the progenitor star was obtained from the analysis of NGC 45 stars presented in Mora et al. (2007). The study performed aperture photometry with the DAOPHOT task in IRAF. Because of the greater flexibility of DAOPHOT to deal with variable background subtraction, we chose this approach over the magnitudes available in the Version 3 of the Hubble Source Catalog (HSCv3) (Whitmore et al. 2016), which are based on aperture magnitudes from SExtractor V2.4.3 (Bertin & Arnouts 1996). The DAOPHOT magnitudes (reported in Vega system) are m F435W =23.878±0.038, m F555W =23.371±0.038, and m F814W =22.678±0.041, which are consistent with the magaper2 values reported in HSCv3: m F435W =24.135, m F555W =23.608, and m F814W =23.346 (AB system), once corrected for aperture 6 . The error bars for the HSCv3 magnitudes can be assumed as the average MAD in the archive for point sources of 23−24 mag, with values between 0.04 and 0.05 mag.
We used the pysynphot STScI Development Team (2013) software to convert the HST magnitudes into fluxes (corrected for extinction) and ran a custom Markov Chain Monte Carlo (MCMC) code BBFit 7 based on emcee (Foreman-Mackey et al. 2013) to fit the black-body emission from the progenitor. Because the extinction in the host is unknown, we left it as a free parameter of the fit. In addition, we performed a second analysis using fixed values of A V for the host extinction. We applied corrective steps of 0.05 mag to the A V parameter on top of the Milky Way value to retrieve the non-reddened fluxes. The maximum A V considered is 0.20 mag, which is guided by the extinction values derived for young stellar clusters in NGC 45 (Mora et al. 2007). The best-fit results for the effective temperature, radius, and luminosity are reported in Table 4.
In addition, we fit the progenitor fluxes using Kurucz model atmospheres (Kurucz 1993). The models are available through the dusty radiative-transfer code DUSTY (Ivezic & Elitzur 1997;Ivezic et al. 1999;Elitzur & Ivezić 2001). Following ; Adams et al. (2016Adams et al. ( , 2017, we used a MCMC data was not included in the analysis due to its poor resolution and contamination from nearby sources. Following the dustmodelling approach for the progenitor of M31-LRN2015 described in Blagorodnova et al. (2020)-see their Sect. 3.3-we assumed that the dust is distributed in an optically thin shell around the system. The shell is formed by silicate or carbonaceous grains of a = 0.1 µm and radiates at a uniform blackbody temperature. Fig. 8 shows the maximum dust mass that would still be consistent with the nondetections, assuming silicate grains. Due to the wavelength range of SST observations, the best constraints are obtained for warmer dust. For example, dust at T d = 1500 K, if present, would need to have a mass of M d < 10 −8.2 M . However, for colder dust at 1000 K, 500 K, and 250 K, the upper limit increases to M d < 10 −6.0 , M d < 10 −4.6 , and M d < 10 −3.6 M , respectively. For graphite dust composition, the limits drop by almost one order of magnitude: M d < 10 −9.4 M (1500 K), M d < 10 −6.8 M (1000 K), and M d < 10 −5.2 M (500 K).
In a similar analysis, we assumed a radial distribution of optically-thin dust corresponding to a stellar wind with r −2 density profile and velocity of 50 km s −1 (Kochanek 2011). We find that the dust mass-loss rate is 10 −7 M yr −1 to satisfy the Spitzer limits. If we consider a relatively small extinction around the progenitor (as suggested by our optical SED modelling), a dusty wind along the line of sight is constrained to have dust mass-loss rate 10 −9 M yr −1 to remain optically thin. However, this constraint can be accommodated by confining the mass loss in the orbital plane of the binary and orienting the line of sight off the orbital plane.

Single stellar evolution models
The progenitor was initially studied using single-star stellar evolution tracks from the code Modules for Experiments in Stellar Astrophysics (MESA; Paxton et al. 2018), available through the MIST (Choi et al. 2016) web interface 8 . The retrieved stellar evolutionary tracks corresponded to a metallicity of Z = 0.0068 ([Fe/H] = −0.32), initial v/v crit = 0.4. Our choice of metallicity was guided by the results of Mora et al. (2007), who found that values between Z = 0.006 and 0.008 better represented the colors of young stellar clusters in NGC 45. Such a metallicity is between the values for the Small and Large Magellanic Clouds.
The location of the progenitor, depicted in Fig. 9, can be associated with massive stars with initial masses between 11.0−12.0 M . The parameter space of the progenitor is consistent with the Hertzsprung gap: a fast transitioning phase between the main sequence and the red-giant branch. As shown in the bottom half of the figure, after leaving the main sequence, the radius of these stars expands considerably with age. For a lowermass star around ∼9 M , the progenitor could also be consistent with being in a "blue-loop", that is, a star that had previously expanded to the red-giant branch and then shrunk again after the onset of core helium burning. However, provided that lowermass stars would have expanded beyond the current radius at an earlier evolutionary stage, any interaction with a close companion would have occurred by that time. Therefore, we assume that the star is most likely in its first expansion after finishing core H burning, which theoretically is the most likely scenario for a stellar merger. In this case, assuming our fiducial extinction value of A V = 0.089, the most likely primary-star progenitor has a zero-age main-sequence mass of M = 11.0 +0.6 −0.2 M . For this progenitor mass, the star would be ∼20.8 Myr old and would have developed a He core with a mass of M c 2 M , with an envelope mass of M env 8.80 M .

Binary stellar-evolution models
The use of single stellar tracks to infer physical properties of the progenitor (Sect. 4.2) is only justified if at the time of observation the mass transfer was yet to commence. This requires 8 http://waps.cfa.harvard.edu/MIST/index.html  Article number, page 9 of 21 A&A proofs: manuscript no. main that the timespan from the onset of RLOF to the onset of the CE phase is very short (<10 yr in the case of AT 2018bwo). In the case of yellow-supergiant donors such as the AT 2018bwo progenitor, such a scenario is only possible in rare cases of an extreme mass ratio (q 20, see discussion in Sect. 6.2). More likely, the mass transfer is driven by thermal-timescale expansion of the radiative envelope. As a result, it takes several hundreds of years before the mass-transfer rate rises to high values 10 −3 − 10 −2 M yr −1 (e.g., Wellstein et al. 2001) and, eventually, culminates in a runaway coalescence of the binary (MacLeod & Loeb 2020b). Therefore, the progenitor observed 10−13 yr before the LRN was most likely already transferring mass at a significant rateṀ onto its unseen companion.
Significant mass loss affects the luminosity of a donor star: instead of being radiated away from the surface, the energy from nuclear burning is captured in the expanding envelope layers and used to rearrange the internal structure (e.g., Paczyński 1967). As luminosity is the key parameter used to estimate the mass of a star based on stellar tracks, it is essential to take this effect into account. This can only be achieved with detailed stellarevolution computations.
Here, we employed the MESA stellar evolution code (Paxton et al. 2011(Paxton et al. , 2015(Paxton et al. , 2019 to constrain binary models consistent with the AT 2018bwo progenitor. The physical ingredients of our models (e.g., internal mixing, wind mass loss, chemical composition) are the same as in the reference model of Klencki et al. (2020), see their Sect. 2. As previously discussed in § 4.2, we adopt a metallicity of Z = 0.0068. Notably, we assume convective core overshooting with the overshooting length of 0.345 pressure scale heights, as guided by observations of B-type giants (∼15 M ) in the LMC (Brott et al. 2011). This is roughly two times larger than the overshooting adopted in the MIST stellar tracks (Choi et al. 2016 Fig. 9. Location of the AT 2018bwo progenitor in the HR diagram. MESA single stellar evolution tracks are shown for stars between 9 and 13M . The main sequence is marked with a thicker line. Different markers represent progenitor parameters using different models (see Table  4). The insert shows a zoomed in region around the progenitor location, along with a finer grid of models. leads to a somewhat higher luminosity-to-mass ratio of our models.
For simplicity, when computing binary models we only evolve the primary (donor) star and treat the companion as a point mass. The primary is initially nonrotating and can be spunup by tidal interactions. We follow the formalism of Kolb & Ritter (1990) to calculate the mass-transfer rate through the L1 Lagrangian point. We assume that the companion is unable to accrete any of the transferred mass. This assumption is commonly made for systems evolving through a phase of rapid case B mass transfer based on a spin-up argument (once the accretor is quickly spun up, the accretion is expected to cease, e.g., de Mink et al. 2013). We note, however, that recent indirect evidence from Be X-ray binaries in the Small Magellanic Cloud suggests a possibly higher accretion efficiency of ∼0.5 in such systems (Vinciguerra et al. 2020). We further assume that specific angular momentum of the mass ejected from the system is an average of the angular momenta of the accretor and the L2 point, guided by the gas-kinematics study of MacLeod & Loeb (2020a). This assumption is somewhat degenerate with the (unknown) mass ratio of the system: a fully isotropic mode of mass ejection (i.e., specific angular momentum of the accretor) would yield similar results to our models but for steeper mass ratios. We explore initial donor masses between 12 and 18 M , mass ratios between 3 and 10, and orbital periods between 100 and 1000 days. The MESA inlists (input files) necessary to reconstruct our work as well as all the binary models are available at https://zenodo.org/communities/mesa.
In Fig. 10, we show several tracks of donor stars in the HR diagram derived from binary models that were found consistent with the location of the AT 2018bwo progenitor. In each case, as soon as the mass transfer rate increases above 10 −3 M yr −1 , the donor's luminosity begins to decrease significantly. As a result, donors with various initial masses (ranging from 12 to 16 M in Fig. 10) can all be consistent with the locus of the AT 2018bwo progenitor, albeit at different mass-transfer rates, ranging from log(Ṁ/M yr −1 ) ≈ −2.4 for the 12 M progenitor to ≈ −1.2 for the 15 M progenitor. The ages of the progenitors range from ≈ 12 Myr for the 16 M donor to ≈ 18.5 Myr for the 12 M donor. The initial mass ratio was q = 5 for models with M don = 12, 13, 14, and 15 M , and q = 3.5 for the model with M don = 16 M . Models with M don = 17 or 18 M were found inconsistent with the progenitor. We note that models with M don = 12 − 14 M could also be consistent with the progenitor for all the other mass ratios explored here (q from 3.5 to 10). For a given donor mass, steeper mass ratios lead to higher mass transfer rates at the progenitor location ( see Sec. 6.2). Fig. 10 is a clear illustration why the usage of binary rather than single stellar models is necessary to infer the initial progenitor mass. In Sect. 6.2 we further discuss the interpretation of the AT 2018bwo progenitor in view of the above analysis.

Effect of binary mass loss on the progenitor appearance
Many LRNe show a gradual brightening, lasting hundreds of days to several years before the main outburst (Tylenda et al. 2011;Blagorodnova et al. 2017Blagorodnova et al. , 2020Pastorello et al. 2021b), which has been interpreted as an evidence of internal shocks in material leaving the binary system (Pejcha et al. 2016. As a result, it is not immediately clear whether the progenitor detection corresponds to the underlying binary or the material leaving the binary. To elucidate the situation for AT 2018bwo, we modified the smoothed-particle hydrodynamics code with radiative diffusion and cooling originally used by Pejcha et al. (2017) to model the ∼200 day long pre-peak brightening in V1309 Sco. We set up several realizations of a binary system with component masses of 13 and 2.6 M on a circular orbit with separation of 1 AU. We begin injecting mass around the outer Lagrange point with a rate of 3 × 10 −2 M yr −1 (based on the analysis in Sect. 4.3), which we increase as a power law toward a singularity at 2000 days from the start of the simulation. Within 1600-1800 days from the start of the simulation, the semi-major axis is reduced by a factor of ∼30 as the ejecta take away angular momentum of the binary, and the simulation stops due to limits on the equation of state. The luminosity of the outflow is initially around 2 × 10 3 L , which increases to about 10 4 L near the end of the simulations. At this point, the mass-loss rate is a few solar masses per year. We conclude that the outflow luminosity ∼10 years before the merger is much smaller than the luminosity of the stars, and the quantities derived in this paper represent the progenitor star.

Explosion energetics
Our photometric and spectroscopic analysis allow us to compute the main observables related to the photosphere during plateau. We estimate the ejecta velocity to be v ej = 500 ± 65 km s −1 , as a weighted average of the FWHM of the Hα line. The plateau has an average luminosity of L p = (1 ± 0.15) × 10 40 erg s −1 and luminosity duration t p = 70 ± 5 days (note that this duration is longer than in the optical bands). It is possible to invert the transient observables to estimate the physical parameters of the explosion. In this section, we take advantage of the binary-evolution models and confront the predictions of common-envelope energy formalism (Sect. 5.1) with physical parameters inferred from the transient, interpreted either as a scaled-down version of Type II-P supernovae (Sect. 5.2) or a shock-powered event (Sect. 5.3).

CE prescription
Classically, the outcome of a CE phase can be estimated based on the energy formalism, in which the binding energy of the envelope, E bind , is related to the energy input from orbital inspiral, ∆E orb , through a parameter α CE that is between 0 and 1: E bind = α CE ∆E orb (van den Heuvel 1976; Tutukov & Yungelson 1979;Webbink 1984;Ivanova et al. 2013b). In our current understanding of CE evolution, events initiated by massive radiativeenvelope donors (such as AT 2018bwo) are generally expected to lead to stellar mergers rather than successful envelope ejections (Klencki et al. 2021) 10 . For AT 2018bwo, looking at the binary progenitor models with 12−16 M donors shown in Fig. 10 at the point when they are consistent with the observations, and the CE inspiral is about to commence, we find envelope binding energies of donors in the range E bind ≈ 2.5 − 3 × 10 49 erg and energy inputs ∆E orb ≈ 0.8 − 2 × 10 49 erg, indicating a merger even for an extreme assumption of α CE = 1.0. For this estimate, E bind includes a gravitational binding term, lowered by the internal energy of gas and radiation as well as the recombination energy. Fig. 11. Estimation of the amount of mass ejected from the envelope during the AT 2018bwo LRN event based on a simple energy formalism, (see Eqn. 1−4). The envelope structure was taken from a binary MESA model with a 15 M donor shown in Fig. 10. The black line shows the energy required to unbind and bring to infinity a portion of the outer envelope with mass M ej . Coloured lines correspond to the orbital energy transferred by accretors with different masses when inspiraling deeper into the primary's envelope, up to the radius when the outer shell has mass M ej .
In the case of a merger, part of the loosely bound outer envelope becomes ejected, powering the LRN transient. Here, we estimate the mass of the ejecta, M ej , based on simple energybudget considerations. We assume that: where the following terms are: 10 Although the situation is less clear in the case of stellar accretors.
Article number, page 11 of 21 A&A proofs: manuscript no. main . The selected envelope profile corresponds to when the donor star was consistent with the location of the AT 2018bwo progenitor in the HR diagram, at which point its total mass was ∼ 13.1 M , radius ∼ 107 R , and the mass transfer rate Here E bind;ej is the binding energy of the ejected part of the envelope, M don is the donor mass, u is the internal energy per unit of mass, E ejecta is the energy of the ejecta at infinity, v ej is the velocity of the ejecta at radius R ej , and E orb;Macc is the energy input from inspiral of the accretor with mass M acc starting from the donor surface R don down to some radius R D at the bottom of the ejected part of the envelope (corresponding to M ej ). We assume v ej = 500 ± 65 km s −1 and R ej ≈ R don ≈ 110 R . The exact value of R ej does not influence the result, as the kinetic term dominates the ejecta energy.
In Fig. 11 we solve Eq. 1 for the envelope structure taken from the binary model with a 15 M donor (initial mass) shown in the HR diagram in Fig. 10. The result is plotted as a function of the mass from the surface of the envelope, i.e., the ejected amount M ej . Coloured lines assume different initial mass ratios q ranging from 3 to 20, and black dots indicate solutions to Eq. 1. For the specific case of a 15 M progenitor, we also show in Fig. 12 the envelope structure of mass, radius, and local escape velocity. We see that the mass from the surface starts to dramatically increase only when we get closer to the core and that the local escape velocity increases more gradually.
The case with q = 20, possibly related to a dynamical CE onset (see discussion in Sect. 6.2), ejects only ≈ 0.025 M of material. This amount appears to be an order of magnitude too low to explain the LRN transient energetics, as further discussed in Sects. 5.2 and 5.3.
The case with q = 3 did not lead to an unstable mass transfer and a CE evolution in our binary MESA simulations (minding the caveat that we model the accretor as a point mass) and is shown in Fig. 11 as a limiting case. Examples with q between 3.75 and 7.5 match well the binary evolution scenario for AT 2018bwo, in which a phase of thermal-timescale mass transfer is culminated in a dynamical inspiral as modeled in Sect. 4.3. In these cases, the energy-budget estimation suggests ejecta masses of M ej ≈ 0.15 − 0.5 M .
We carried out this exercise for all the donor masses that were found consistent with the progenitor in our binary model analysis (see Fig. 10 We note that between the progenitor observations (10−13 yr before the merger) and the LRN event itself, the donor is expected to still lose a substantial amount of mass of 1 M . Therefore, the envelope structure models used here are only an approximation of the envelope structure at the onset of the dynamical inspiral. As such, the inferred values of M ej should be treated as a crude estimation, mainly indicating an order of magnitude for M ej . In the following sections we constrain the ejected mass in more detail based on models of LRNe energetics.

Scaled Type II-P supernova model
Light curves of hydrogen-rich explosive transients will generally show plateaus of nearly constant luminosity, as the photosphere is located near the hydrogen recombination front, where the opacity steeply changes over a small range of temperatures. This is commonly seen in Type II-P supernovae, but a similar model was developed by Ivanova et al. (2013b) to link the observational class of LRNe with common-envelope ejection events.
Both Type II-P supernovae and LRNe have similar plateau durations of ∼100 days, but Type II-P supernovae have ejecta masses ∼ 20 times higher and expansion velocities ∼10 times higher than LRNe. As a result, LRNe ejecta are a factor of ∼50 denser than those of Type II-P supernovae. Still, we expect the equation of state of optically thick LRN ejecta to be radiationdominated with density-independent Thomson scattering as the primary opacity source. Around or below hydrogen recombination, the opacity might depend on density. For example, if the Rosseland-mean opacity of the negative hydrogen ion, which scales as ρ 0.5 , dominates in the outer LRN ejecta, we expect that effective temperatures of LRNe will be lower than those of Type II-P supernovae. However, the plateau duration or luminosity should not strongly depend on the form of low-temperature opacity. The only fundamental difference between Type II-P supernovae and LRNe could be the density profile of the ejecta, which depends on the uncertain mass-ejection mechanism in LRNe.
Given these similarities between Type II-P supernovae and LRNe, we follow Ivanova et al. (2013b) and apply the same analytic scaling relations linking plateau luminosity L p , duration t p , and expansion velocity v ej to the explosion energy E, ejecta mass M ej , and initial radius R 0 (Arnett 1980;Litvinova & Nadezhin 1985). We explore different calibrations for these relations from the literature.
Recently, Sukhbold et al. (2016) found an excellent agreement between 1D flux-limited diffusion simulations of Type II-P supernova light curves and analytic scaling relations of Popov (1993) with modified absolute coefficients. The plateau luminosity L p is where we replaced E = 1 2 M ej v 2 ej and rescaled to parameters relevant for LRNe. Ivanova et al. (2013b) and MacLeod et al. (2017) used identical equations but with the absolute term based on the original work of Popov (1993). In our units, their absolute terms are 4.2 × 10 38 erg s −1 and 133 days for their default choice of opacity of 0.32 cm 2 g −1 and recombination temperature of 4500 K. Alternatively, Kasen & Woosley (2009) performed Monte Carlo multiwavelength radiation transport in the context of Type II-P supernovae and obtained an absolute term in Eq. (5) of 6.29 × 10 38 erg s −1 and slightly different exponents in Eq. (6), which then becomes However, the calibration of Kasen & Woosley (2009) might be less relevant for LRNe, because their simulations are tuned for Type II-P supernova densities and velocities. In Eqs. (5-7), R 0 should be interpreted as a radius where any heating of the ejecta stops, and the ejecta expands homologously thereafter (Goldberg & Bildsten 2020).
To find M ej and R 0 from Eqs. (5-7), we randomly sampled L p , t p , and v ej assuming uncorrelated Gaussian distributions. For the Sukhbold et al. (2016)  It is interesting to note the substantial differences in the inferred values from the relations of Sukhbold et al. (2016) and MacLeod et al. (2017), which differ only by the absolute terms. The estimate of M ej based on Sukhbold et al. (2016) is substantially larger than predicted in Sect. 5.1 (Fig. 11), but R 0 is much smaller and located inside the progenitor. We see that M ej and R 0 are roughly compatible with the progenitor structure and that v ej corresponds to the escape velocity at R 0 within few tens of percent, as we showed in Sect. 5.1 (Fig. 12). It is natural to expect that the common-envelope ejecta velocity will be similar to the escape velocity of the binary companion at radius r, where most of the envelope ejection occurs.
The M ej based on MacLeod et al. (2017) is roughly compatible with the predictions of the CE theory shown in Sect. 5.1 (Fig. 11), however, R 0 is substantially larger than the initial progenitor radius, implying that the ejecta got reheated after ejection, probably by internal shocks.
The inferences based on Kasen & Woosley (2009) relations give very small M ej and unrealistically large R 0 , much larger than the original progenitor radius or the likely binary orbit. The likely explanation is that the Kasen & Woosley (2009) relations are too tuned for Type II-P supernovae and hence not applicable to LRNe.
To summarize our results, the inferences based on recombination models are broadly compatible with other constraints on the progenitor if the event was a stellar merger, although there are discrepancies with the standard common-envelope energy formalism applied to binary evolution models of the AT 2018bwo progenitor. Analytic scaling relations give different values of M ej and R 0 . Our preferred choice of analytic scaling relation is the one of Sukhbold et al. (2016). Analytic scaling relations are especially powerful for examining relative differences in a population of objects, where the uncertain absolute terms play less of a role.

Shock-powered model
Metzger & Pejcha (2017) proposed a model for the light curves of LRNe, in which faster spherically symmetric ejecta collide with a preexisting equatorially concentrated mass distribution. This model can naturally explain double-peaked LRN light curves, where the first peak is caused by cooling of the fast unshocked polar ejecta, while the second peak or plateau arises due to continuing shock interaction acting as an embedded power source inside the ejecta. The equatorial distribution is formed by mass loss from the binary preceding the dynamical merger, which is likely caused by nonconservative mass transfer and/or by ejection through the outer Lagrange point (Shu et al. 1979;Pejcha 2014). This long-lasting pre-dynamical mass loss can explain gradual brightenings seen in many LRNe before the primary peak, as was quantified for V1309 Sco by Pejcha et al. (2017).
A power source embedded inside hydrogen-rich ejecta will in many cases lead to a plateau-like light curve. Sukhbold & Thompson (2017) provided a particularly striking example by embedding a spinning-down magnetar inside otherwise subluminous Type II-P supernovae ejecta. This setup led to plateau luminosities higher by a factor of ∼30 and light curves photometrically indistinguishable from normal Type II-P supernovae. Another possible source of internal heating could also be jets emanating from a central engine (Shiber et al. 2019;Soker 2020;Soker & Kaplan 2021).
Within the model of spherical ejecta colliding with a preexisting equatorial mass distribution,  provided analytic relations for the characteristic luminosity of the shock-powered peak and the timescale for the shock-powered light curve to rise to its peak where we identify L p with the plateau luminosity and t p with the plateau duration. Here, M w is the mass of the preexisting wind, which occupies a fraction of the solid angle f Ω and which was formed by binary mass loss with a rate exponentially increasing on a timescale t run . The radial velocity of the preexisting wind is v w , which we assume to be 10% of the fast spherically symmetric ejecta velocity v ej with mass M ej . Following , we assume that the velocity of the shocked shell is v sh = 0.5v ej . We further assume opacity κ = 0.32 cm 2 g −1 .
Following the same procedure as in Sect. 5.2, we can use the observed properties of AT 2018bwo to infer M w and M ej from Eqs. (8-9). For t run = 1000 days, we get M w = 2.34 +0.31 −0.27 M and M ej = 0.80 +0.07 −0.06 M . To our knowledge, this is the first time that Eqs. (8-9) have been used to infer parameters of a transient, which makes it worthwhile to explore the sensitivity to changes of parameters. Increasing or decreasing t run by a factor of 10 leads to M w higher by a factor of few. Fixing v w to 100 or 10 km s −1 changes M w to 1.92 +0.21 −0.19 or 3.92 +0.70 −0.57 M , respectively. If we put a pre-factor of 2 or 0.5 in front of the brackets of Eq. (9), which represents our uncertainty in relating the observed transient duration to the theoretical timescale, M ej changes to 0.20 +0.02 −0.02 or 3.22 +0.26 −0.25 M , respectively. In all cases, the correlation coefficient between log M ej and log M w is around −0.52, which does not suggest any strong degeneracy between the two parameters.
Although the analytic model of  has not yet been calibrated with a multidimensional radiation hydrodynamic simulation, our inferred parameter values appear reasonable. We find that the amount of mass lost in the pre-dynamical phase of the merger is likely a few solar masses. Binary-evolution models predict that the mass-transfer rate steeply increases with time (Fig. 14), and we expect that most of the material leaves the binary. Our simulations of the final few thousand days of the binary evolution described in Sect. 4.4 suggest that about 2 M of material lost with specific angular momentum of the outer Lagrange point is sufficient to remove most of the binary angular momentum. Our inference of M w is thus compatible with the theoretical expectations. The value of M ej also agrees with the predictions of commonenvelope energy formalism shown in Fig. 11, especially if the plateau duration is longer than indicated by Eq. (9). There is the same discrepancy as in the recombination model that the escape velocity of the progenitor layer corresponding to M ej is by a factor of two smaller than the observed v ej .
To summarize, both recombination-only and shockinteraction models can explain the observed properties of AT 2018bwo within the limits of current theoretical uncertainties and degeneracies inherent to the inverse problem. A more involved comparison would require developing the shock-interaction model to the same level of sophistication as the recombination model achieved for Type II-P supernovae. Parameters from both models are broadly compatible with the common-envelope theory (for mass ratios q 10) and the progenitor structure of AT 2018bwo. Further study should be devoted to the discrepancy between the observed v ej , which implies mass ejection from deep inside the progenitor, and the relatively low M ej , which suggests ejection of layers close to the surface.

Discussion
This is the first observational study of an extragalactic stellar merger in which the progenitor was interpreted using binary stellar-evolution models. In this section, we first discuss AT 2018bwo in the context of other LRNe. Next, we elaborate on the interpretation, limitations, and future directions for modelling LRN progenitors using stellar binary models. Finally, we consider the binary-evolution model within the context of the star cluster where AT 2018bwo is located.

Comparison of AT 2018bwo with other LRNe
Our LRNe sample contains known transients from the literature (see Table 5), except two Galactic LRNe OGLE-2002-BLG-360 (Tylenda et al. 2013) and CK Vul (Kato 2003), and the LRN in M31: M31 RV (Rich et al. 1989;Mould et al. 1990), due to the limited available data. Our comparison, shown in Fig. 13, particularly focuses on the correlations between: a) the peak magnitude and the length of the plateau, b) the peak magnitude and the FWHM of the Hα line reported in the literature, which is a proxy for the expansion velocity of the ejecta, c) Hα FWHM and the length of the plateau, and d) the peak magnitude and the inferred mass for the primary progenitor. Here we distinguish between bona-fide LRNe that showed molecular absorption at later times and LRNe candidates, where such absorption was not identified (filled vs. empty circles). In addition, we marked those LRNe that at early times showed a cool, red continuum (red markers) from those that appeared hotter (blue markers) and showed stronger Balmer emission lines from interaction with the circumstellar medium (CSM).
One possible source of confusion for LRNe is another "gap transients" class of intermediate luminosity red/optical transients (ILOT/ILRT), such as SN 2008S (Prieto et al. 2008;Botticella et al. 2009), NGC 300 OT2008-1 (Bond et al. 2009;Humphreys et al. 2011), iPTF 10fqs (Kasliwal 2011) or M85-OT (initially proposed to be a LRN by Kulkarni et al. (2007)). The progenitors of these transients are generally massive, dust-enshrouded stars (Prieto et al. 2009), and their spectra are characterized by generally cold continuum with emission for the Balmer and Ca ii IR triplet lines, as well as the characteristic lines of [Ca ii], which are absent in LRNe spectra. The colour evolution for ILRTs is also more moderate when compared to LRNe, and their latetime spectroscopy does not show the appearance of molecular absorption bands.
Within the LRNe sample, AT 2018bwo corresponds to a transient with intermediate brightness and duration. Its peak magnitude of M r = −10.91 ± 0.03 mag is in between V838 Mon and M101 OT2015-1, similar to AT 2020hat. However, its rband plateau duration of ∼41 days is in better agreement with the less luminous M31-LRN2015. This shorter timescale may be caused by the higher ejecta velocity of ∼500 km s −1 as compared to other lower luminosity LRNe, which all had FWHM ≤300 km s −1 .
One interesting highlight of the comparison is the strong correlation between the peak magnitude of the transient and other observables, e.g., the length of the plateau, the expansion velocity, and the mass of the progenitor star. Similar correlations were explored in Kochanek et al. (2014); Mauerhan et al. (2018); Pastorello et al. (2019aPastorello et al. ( , 2021b, although with slightly different approach. For example, our analysis estimates the duration of the transient by fitting an analytical function to the optical light curve, as described in Valenti et al. (2016) for core collapse supernovae (see their eq. 1). We considered the time of the first peak as our reference epoch and adopted an uncertainty of 10%. For UGC 12307-2013OT1 the time of the first peak is unknown, so that we computed the plateau duration for the available part of the lightcurve and included the time to the last nondetection as an upper error in our measurement. For V4332 Sgr, there is no available data before its discovery, so we also cautiously adopted a larger upper value for the duration of the plateau.
Another highlight is the apparent distinction between LRNe depending on their peak magnitude. On the one hand, transients with M peak ≥ −12 mag generally show one or no distinguishable peaks, followed by a plateau in the optical bands. At early times, Open circles indicate additional LRNe candidates whose nature is still debated. Circles represent LRNe with one initial peak followed by a plateau. Squares are LRNe having a slower secondary red peak. Triangles show LRNe with 3 maxima. Red markers represents transients that already showed red continuum during the first peak. Blue markers correspond to LRNe that showed a hotter continuum and Balmer emission lines at early times. Gray markers are used for transients without early times spectroscopy. The name of some transients has been shortened for clarity. this group also mostly shows a reddened continuum and a forest of narrow absorption lines. On the other hand, brighter transients tend to show two distinct peaks: an initial fast "blue" peak, followed by a more extended "'red" peak. Spectroscopy taken during the first peak shows a rather hot continuum with strong Balmer emission lines from CSM interaction. During the second peak, the characteristic cool continuum emerges, and molecular lines form. If Thompson scattering is the main mechanism responsible for the formation of emission lines, analogous to interacting SNe, then the two groups would differ in the state of ionization of the photosphere during peak.   For higher-luminosity events, the photosphere would initially be located in a fully ionized shell of previously ejected material. This shell, preceding the shock wave generated by the mass ejection during the merger, would reprocess the shock emission until the ejecta take over, pushing the rapidly cooling photosphere outwards. This cooling would also move the peak of the continuum emission toward longer wavelengths, appearing as a secondary peak in the redder bands. While during the first peak the radius is expected to remain fairly constant, the cooling and expansion of the ejecta would increase the apparent location of photosphere while keeping the overall luminosity relatively stable.
For lower-luminosity events, the photosphere would be located in a neutral (or partially ionized) shell. Therefore, any shock emission generated by the dynamical ejection of mass would be absorbed and reprocessed by this layer, transforming the shock's kinetic energy into thermal energy that would heat and partially ionize the shell. Spectroscopically, the shell would still appear as a reddened continuum at early times, as the outer neutral layers would absorb great part of the emission at shorter wavelengths.
The ionization of the outer shell appears to depend on the peak magnitude of the outburst. Brighter events with potentially more massive progenitors are likely capable of generating more energetic outflows and supply enough radiation to fully ionize the shell before the onset of the dynamical ejection.
When comparing the length of the LRNe plateau and the expansion velocity (from Hα FWHM), there seems to be yet another dichotomy between lower-and higher-luminosity events. While fainter and redder events have generally some correlation between Hα FWHM and plateau length, the higher-luminosity group has a much larger spread. This is mainly due to transients we designated as LRNe candidates (empty markers), which seem to fall outside of the increasing trend, as they correspond to shorter events with larger expansion velocities.
The final highlight is the tight correlation between the mass of the primary progenitor and the peak magnitude of the outburst, even for LRNe with massive progenitors. This trend was first proposed by Kochanek et al. (2014) based on Galactic LRNe with masses lower than 10 M . Our sample shows that such a correlation can also be extended to extragalactic events, with generally more massive progenitors. This shows that the amount of mass ejected in stellar mergers is proportional to the initial mass of the donor, although fluctuations may be associated with the mass of the secondary companion or the geometry of the system.

Binary stellar models for LRN progenitors
In this work, we analyzed the progenitor of AT 2018bwo in the rightful context of binary evolution. We interpret the star observed 10−13 yr before the LRN event as a primary (more massive and luminous) component in a binary system in which a phase of mass transfer culminated in a common-envelope event and the LRN transient. In Sect. 4.3, with the use of a grid of binary stellar models computed with MESA, we found that the location of the AT 2018bwo progenitor in the HR diagram can be consistent with a donor star of an initial mass in the range ∼12−16 M (see Fig. 10). This mass is significantly larger than inferred from single-star models (∼11 M , see Sect. 4.2), which highlights the importance of binary models for LRN progenitors.
Ideally, one would want to narrow down the range of possible initial donor masses, M don = 12 − 16 M , obtained from binary models. This, however, is not straightforward due to challenges in the modeling of interacting binaries with high masstransfer rates (Ṁ). As shown in the upper panel of Fig. 14, donors of different initial masses are consistent with the locus of the AT 2018bwo progenitor for different values ofṀ (with log(Ṁ/M yr −1 ) ranging from −2.4 to −1.2). This is because, in general, the higher the value ofṀ, the larger the decrease in the luminosity of the mass-losing donor (see Fig. 15). The details of this relation depend also on the current structure of the envelope and hence the amount of mass that has already been lost from it, which is sensitive to the assumed mass ratio. This explains the degeneracy between M don ,Ṁ, and q among the binary models that are consistent with the progenitor, as shown in Fig. 15.
It is difficult to say which value ofṀ is consistent with the fact that only 10−13 yr later the binary coalesced and produced a LRN event. The onset of a CE phase and a runaway dynamical inspiral is thought to be preceded by a stage of extensive loss of mass and angular momentum through L2 outflows, which decrease the separation to roughly the size of the donor's envelope (e.g., Pejcha et al. 2016Pejcha et al. , 2017. This process can take a few hundred dynamical timescales or several years, during which ∼15% of the companion mass needs to be ejected through L2 to provide sufficient orbital shrinkage (MacLeod et al. 2017).
Thus, the progenitor system observed 10−13 yr ago was most likely on the verge of experiencing significant L2 mass loss. Such outflows could be triggered due to the donor star overflowing its Roche lobe up to the L2/L3 potential. The critical mass-transfer rate at which that happens is largely uncertain due to limitations of mass-transfer schemes in 1D stellar codes (Pavlovskii & Ivanova 2015). Alternatively, signif- Fig. 15. Relation between the current mass transfer rate in a binary and the decrease in the donor luminosity with respect to the moment of Roche-lobe overflow. The higher the rate with which the donor looses mass, the more energy is absorbed in an expanding envelope and, as a result, the lower is the surface luminosity. This general relation is degenerated with factors such as the mass ratio.
icant L2 outflows may be triggered because of the companion star being driven out of thermal equilibrium and expanding significantly (e.g., Benson 1970; Neo et al. 1977), which could lead to it filling its Roche-lobe, a contact-binary stage (Pols 1994;Wellstein et al. 2001), and potentially CE evolution (de Mink et al. 2007;Marchant et al. 2016). Details of this process remain largely uncertain, partly because of the complicated gas dynamics in semi-detached binaries (Lubow & Shu 1975), poorly constrained specific entropy of the accreted material (Shu & Lubow 1981), and unknown efficiency of accretion by stars rotating near their breakup limit (Popham & Narayan 1991). Because of these uncertainties, in the current study we are unable to exclude any of the binary model solutions based on the mass-transfer rate 10−13 yr prior to the LRN event. We note that once significant L2 outflows commence, the mass-transfer rate is expected to quickly increase (see also Sect. 4.4), changing the evolution ofṀ in Fig. 14 to a much steeper rise with time.
Alternatively, the mass transfer in the progenitor of AT 2018bwo was still yet to commence 10−13 yr before the LRN event. Such a scenario may be possible if the donor star is shrinking with respect to its Roche lobe on the adiabatic timescale. For this to happen, the mass ratio between the binary components has to be large enough (i.e., q > q crit;ad where q = M don /M acc ). In the case of radiative-envelope supergiant donors such as the progenitor of AT 2018bwo, the critical mass ratio q crit derived from adiabatic mass-loss sequences is q crit;ad ≈5−10 (Ge et al. 2015) but increases to possibly even q crit;ad > 20 when thermal relaxation of the outer envelope layers is taken into account (Ge et al. 2020).
The onset of the CE phase might also be very rapid when driven by Darwin (tidal) instability (Darwin 1879). The condition for tidal instability is based on the ratio between momentum of inertia of the orbit and of the donor star to be I orb /I don 3 (Eggleton & Kiseleva-Eggleton 2001). This condition can be rewritten into an expression for the critical separation below which the instability happens a crit = R don 3η don (1 + q) where the parameter η don = I don M −1 don R −2 don is related to the internal structure of the donor (MacLeod et al. 2017). Using single MESA models (with same assumptions as in the binary MESA models described in Sect. 4.3), we find that 10−12 M yellow supergiants consistent with the position of the AT 2018bwo progenitor are characterized by η don ≈ 0.013−0.014. For such values of η don , the condition of the critical separation a crit being equal to the RLOF separation requires extreme mass ratio values q crit;Darwin ≈ 45−50, which would correspond to companion stars of 0.2 M .
In summary, binary-evolution models suggest a severalhundred year long phase of essentially stable mass transfer prior to the LRN event and yield a significantly higher progenitor mass (12−16 M ) compared to single-star models (∼11 M ). Additional constraints are needed to further narrow down the progenitor mass range. Such constraints could come from clues on the stellar cluster environment of AT 2018bwo (see Sect. 6.4) or, in future work, from detailed modeling of the dust content of the progenitor system (see Sect. 6.3). Excitingly, any independent constrains on the progenitor mass or age, combined with binary-evolution models, would be a unique probe of the elusive conditions when rapid mass transfer leads to CE evolution, potentially helping to address some of the open questions outlined in this section.
An alternative scenario in which the progenitor of AT 2018bwo was still in a detached state 10 yr before the merger requires extremely steep mass ratios 20. This seems unlikely, given that such systems are relatively rare (Moe & Di Stefano 2017), and this is also disfavored based on the estimated amount of mass M ej that was ejected from the system during the dynamical inspiral (see Sect. 5.1).

On the dust content in the progenitor system
Independent constraints on the binary-evolution scenario may come from limits on the dust content of the progenitor system, obtained from the combination of archival HST photometry and SST nondetections (Sect. 4.1). Based on a simple model of an optically thin and isothermal dust configuration, we find an upper limit on the mass of the dust at log(M d /M ) = −3.0 if the dust temperature is T d = 250 K and much lower M d limits if T d is higher (e.g., log(M d /M ) = −6.0 for T d = 1000 K, see Fig. 8). This amount of dust is at least a few times lower than the mass of metals M Z (i.e., future dust) that is transferred through L1 in our binary models for the progenitor, as shown in the lower panel in Fig. 14, where M Z always reaches values 2 × 10 −3 M . It is expected that most of M Z will be ejected from the system and likely form dust (for a discussion on accretion efficiency in case B mass transfer see de Mink et al. 2013;Vinciguerra et al. 2020), seemingly in tension with the upper limits on the dust mass derived from the SST data. Similarly, the mass-transfer rate through L1 is several orders of magnitude larger than the limit on optically thin dusty wind of 10 −7 M yr −1 (Sect. 4.1).
There are two possible explanations why we did not detect dust in the progenitor, despite theoretical models suggesting copious mass loss. The first possible explanation is that the dust forms at much larger radii and lower temperatures, where it would radiate primarily in mid-IR and far-IR, where we do not have data. This could happen because of internal shocks in the outflow resetting grain growth by collisions or by UV radiation. Quantifying this possibility would require a much more sophisticated calculation beyond the scope of this work. The second possible explanation is that dusty outflow is radiatively inefficient. If we consider an equatorially concentrated outflow with velocity v, where radiation escapes perpendicular to the disk, the ratio of diffusion to expansion time is approximately where we assumed X d = 0.005 and that the main contribution to emission comes from around the dust-formation radius, which is R form ≈ 9 AU for the luminosity of AT 2018bwo progenitor (Kochanek 2011). For the theoretically expectedṀ ≈ 10 −2 M yr −1 , the grains grow to around 1 µm (Kochanek 2011), which implies that the opacity of silicates at 5 µm is κ ≈ 10 3 cm 2 g −1 . If τv/c is around unity or higher, the outflow loses most of its internal energy to adiabatic expansion before it can be radiated. Given the scalings in Eq. (10), the pre-merger binary mass loss in AT 2018bwo could plausible be radiatively inefficient. This suggests that the theoretically predictedṀ could still be compatible with the SST upper limits presented in Figure 8.

Cluster environment of AT 2018bwo
A closer analysis of the AT 2018bwo progenitor site in NGC 45 suggests that the system was a likely member of a young stellar cluster, as shown in Fig. 16. The observed group of stars has an angular radius of ∼2 , which at the distance of NGC 45 represents a radial extent of ∼64 pc, consistent with an OB association. Assuming a similar formation time, the evolved primary progenitor of AT 2018bwo would necessarily have been one of the most massive stars in the group, as other more massive components would have already ended their lives as core-collapse supernovae. However, the colour-magnitude diagram in Fig. 16 shows that potentially the cluster contains stars as massive as 16−18 M , with ages nearly half of the possible progenitors (as discussed in Sect.6.2, binary models predict ages between 12 and 18.5 Myr for 16 M and 12 M donors respectively). Simultaneously, there are also RSGs of ∼10 M , which are > 7 Myr older than the progenitor.
This apparent discrepancy can easily be explained by the limited data available on the cluster. On one hand, the lack of radial velocities or proper motions limits our ability to distinguish between real cluster members and foreground or background contamination. On the other hand, unresolved binaries and unknown extinction around each star can also play a role in the location of each source in the diagram. Finally, the cluster may also have an intrinsic spread in age, and our initial assumption of co-evolution does not hold.
Despite the mentioned caveats, there seem to be two main populations in the cluster, with ages of ∼10 and 20 Myrs. This is not unexpected, as blue stragglers from stellar mergers have been consistently observed in clusters, re-populating the main sequence above the main sequence turn-off point. In particular, in young (∼5 Myr) open clusters, their fraction near the turn-off may be as high as 30% (Schneider et al. 2015). Observations of young clusters in the S Doradus region in the LMC have also shown that the fraction for their analogous evolved counterparts, the red stragglers, can be as high as ∼55% in a coeval cluster (Britavskiy et al. 2019).
Given that the progenitor is unlikely to be a merger product itself (unless it was originally a triple system), it likely belongs to the cluster's lower-mass population. Hence, we can use the observed RSGs to narrow down its mass to be in the 10−13 M The progenitor is marked with a red star. Middle: Colour-magnitude diagram in B−I (F435W − F814W) bands for HST simulated magnitudes from MESA as compared to the progenitor magnitude. The initial mass for each model is labeled in the plot. Right: Same diagram, but showing the isochrones for different stellar ages. In both plots, the location of the AT 2018bwo progenitor is marked with a red star. The black arrow indicates the location of the progenitor when corrected for A V =0.2. The cluster components are shown in blue, and the rest of field stars in the same HST pointing are shown in faint gray. regime, which is consistent with the 11−16 M mass range derived using both single and binary stellar-evolution models.

Summary and conclusions
In this work, we have presented the results of our photometric and spectroscopic follow-up campaign of the AT 2018bwo LRN in the optical and IR wavelengths. We also modeled its progenitor system, which was observed by HST ∼10 years before the outburst and, for the first time, determined the evolutionary stage of its progenitor system using binary stellar-evolution models.
AT 2018bwo had a peak absolute magnitude of M r = −10.97 ± 0.11 mag, between those of V838 Mon and M101 OT2015-1, and comparable to AT 2020hat. The duration of its r-band plateau of 41 ± 5 days is in better agreement with lower-luminosity transients, such as M31-LRN2015 and AT 2019zhd, both discovered in M31. Similarly to AT 2020hat and other fainter LRNe, the early-time spectra of AT 2018bwo were already marked by a cool continuum of ∼3000 K and strong absorption lines of Fe and low-ionization elements. The Hα lines also show the characteristic blue-shifted emission with an average FWHM of ∼500 km s −1 and a narrow absorption component on top. Toward the end of the plateau, the appearance of strong TiO and VO molecular lines also matches previous LRNe observations. Our NIR spectra taken at 103 days after discovery show an enlarged cool star, corresponding to a M8.5 II type, analogous to AGB stars in the LMC. The progressive cooling and rapid creation of dust in the remnant is also observed in the NEOWISE MIR data, which show that, although the object is not visible at NIR wavelengths, its emission in the MIR had increased 1.5 years after the outburst.
Using MESA binary stellar-evolution models, we showed that the progenitor primary star is in the 12−16 M range, which is 9−45% more massive than determined from single stellarevolution models alone. We propose that the system was likely in a prolonged stage of semi-stable thermal-timescale mass transfer, with mass-transfer rates ofṀ ∼ 10 −2 M yr −1 , allowing the primary to lose several M before the dynamical onset of the CE. Surprisingly, this mass is not detected as an MIR dust excess in the SST data of the system in quiescence. For optically thin warm (1500 K − 250 K) dust, we place constraints of 10 −8 −10 −3 M on the maximum dust mass present around the donor 10−13 years before the outburst. This shows that the outflows from the primary were likely radiatively inefficient, or that the dust formation occurred at a larger radius, where temperatures are colder. We also suggest the onset of the dynamical instability was initiated by quick loss of angular momentum, caused by increasingly high L2/L3 mass loss, starting within the last few years before the outburst.
The analysis of the primary envelope's structure and the system's orbital energy support a partial ejection of the binary CE, likely within the range of 0.15−0.5 M . This mass is also in agreement with ejecta masses derived from modeling the energetics of the LRN outburst with both scaled supernova Type-II P (following Ivanova et al. (2013b)) and shock-powered models . Provided this mass is only a fraction of the total mass of the primary's envelope at the time of dynamical onset, we confirm that the LRN outburst AT 2018bwo is related to a stellar-merger event.
Our results show that the combination of observations of progenitor systems and binary stellar-evolution models is a powerful tool to explore the conditions that may drive binary stars to unstable mass transfer and quick coalescence, which is critical to improve our understanding of the rapid mass-transfer evolution and the CE phase in binary systems. Future surveys, such as the Rubin Observatory Legacy Survey of Space and Time (LSST; LSST Science Collaboration et al. 2009) are expected to discover 20 − 750 LRNe per year (Howitt et al. 2020). Although only the closest host galaxies are expected to have archival high resolution data -critical to identify the progenitor systems in quiescence-, the depth of the survey will be ideal to exquisitely sample the precursor emission from the binary years before the LRN event, providing unique observational constraints on the intensive mass loss that drives the system to coalescence. Increasingly large samples of LRNe with detected progenitors and precursors will provide a rich opportunity to address the masstransfer stability question from a statistical perspective, which will have an impact on binary population-synthesis models and hopefully improve our understanding on how binary stars evolve.