Light Curve Model for Luminous Red Novae and Inferences about the Ejecta of Stellar Mergers

The process of unstable mass transfer in a stellar binary can result in either a complete merger of the stars or successful removal of the donor envelope leaving a surviving more compact binary."Luminous red nova"(LRN) are the class of optical transients believed to accompany such merger/common envelope events. Past works typically model LRNe using analytic formulae for supernova light curves which make assumptions (e.g., radiation dominated ejecta, neglect of hydrogen recombination energy) not justified in stellar mergers due to the lower velocities and specific thermal energy of the ejecta. We present a one-dimensional model of LRN light curves, which accounts for these effects. Consistent with observations, we find that LRNe typically possess two light curve peaks, an early phase powered by initial thermal energy of the hot, fastest ejecta layers and a later peak powered by hydrogen recombination from the bulk of the ejecta. We apply our model to a sample of LRNe to infer their ejecta properties (mass, velocity, and launching radius) and compare them to the progenitor donor star properties from pre-transient imaging. We define a maximum luminosity achievable for a given donor star in the limit that the entire envelope is ejected, finding that several LRNe violate this limit. Shock interaction between the ejecta and pre-dynamical mass-loss, may provide an additional luminosity source to alleviate this tension. Our model can also be applied to the merger of planets with stars or stars with compact objects.


INTRODUCTION
The direct interaction between the stars in a binary system is a common occurrence (e.g., Sana et al. 2012), with implications for a number of topics in stellar evolution (e.g., Podsiadlowski et al. 2004;de Mink et al. 2014;De Marco & Izzard 2017) and energetic transients (e.g., Soker & Tylenda 2006;Chevalier 2012;Smith 2014). Such interaction often begins after the (typically more evolved) donor star overflows its Roche lobe and begins transferring mass onto the more compact accretor star. Depending on the binary mass ratio and eccentricity and the evolutionary state of the donor, this masstransfer can become unstable, resulting in a runaway growth of the mass transfer rate (e.g., Soberman et al. 1997;Pavlovskii et al. 2017;Ge et al. 2020;Klencki et al. 2021).
Following an initial period of mass-loss, potentially lasting many orbits and likely concentrated in the binary plane (e.g., from the L 2 Lagrange point; Pejcha et al. 2016aPejcha et al. ,b, 2017MacLeod & Loeb 2020) it is commonly believed that the accretor will ultimately be dragged into the envelope of the donor star (Paczynski 1976;Taam et al. 1978;Iben & Livio 1993;Ivanova et al. 2013b). The final outcome of such a "common envelope" phase can be either a complete merger of both stars (e.g., Kruckow et al. 2016), or the successful ejec-tion of the donor envelope and the survival of a tighter binary comprised of the accretor and the core of the original donor (e.g., Ivanova & Nandez 2016;Sand et al. 2020;Law-Smith et al. 2020). Stellar mergers and/or common envelope events have been invoked to give rise to number of atypical stellar and binary populations, including cataclysmic variables (Patterson 1984), magnetic A-stars (Schneider et al. 2019), magnetic white dwarfs (Tout et al. 2008;Nordhaus et al. 2011), sdB-stars (e.g., Kramer et al. 2020), carbon-deficient giants (e.g., Bond 2019), and blue stragglers (e.g., Davies et al. 2004;Ferreira et al. 2019;Wang et al. 2020), among other examples. The conditions required to remove the envelope, and the final separations of the surviving binary, are also important for predicting the population of merging binary compact objects by LIGO/Virgo (e.g., Tauris et al. 2017;Belczynski et al. 2018;Vigna-Gómez et al. 2018). However, despite increasingly sophisticated multi-dimensional simulations of the common envelope phase over the past decade (e.g., Ricker & Taam 2012;Ohlmann et al. 2016;Ivanova & Nandez 2016;MacLeod et al. 2018b;Prust & Chang 2019;Chamandy et al. 2019;Law-Smith et al. 2020), the theoretical community has yet to converge on a consistent physical picture which includes all of the relevant physical processes (Clayton et al. 2017;Soker et al. 2018;Wilson & Nordhaus 2019;Reichardt et al. arXiv:2202.10478v3 [astro-ph.SR] 28 Apr 2022 2020). Observations remain extremely important to unraveling this mystery.
The dynamical merger event may generate a shortlived transient powered by the ejection of mass, primarily from the donor's envelope (Soker & Tylenda 2006). "Luminous red novae" (LRN) are the class of optical transients with typical durations of weeks to months and peak luminosities ∼ 10 38 −10 41 erg s −1 between those of classical novae and supernovae (SNe) and characteristically red colors (e.g., Martini et al. 1999;Munari et al. 2002;Bond et al. 2003;Kulkarni et al. 2007;Tylenda et al. 2011;Kurtenkov et al. 2015;Smith et al. 2016;Blagorodnova et al. 2017;Cai et al. 2019;Stritzinger et al. 2020; see Pastorello et al. 2019a for a recent review and Fig. 1 for a sample of LRNe). The watershed event which confirmed the association between LRNe and merging binary stars was the Galactic source V1309 Sco (Mason et al. 2010), detected by OGLE as a contact binary for almost a decade prior to the outburst (Tylenda et al. 2011). The observed evolution of the binary light curve and orbital period leading up to the merger, provided evidence for high rates of mass-loss, starting hundreds of orbits prior to the final dynamical coalescence and concomitant abrupt rise to optical maximum (Pejcha 2014;Pejcha et al. 2017;MacLeod & Loeb 2020). Similar "precursor" emission, potentially powered by internal shocks in outflowing streams from the L 2 point (e.g., Pejcha et al. 2016a,b), is observed in the light curves of other LRNe, suggesting that predynamical mass loss is common if not ubiquitous in stellar mergers Blagorodnova et al. 2021).
The light curves of LRNe encode information on the merger ejecta, particularly its mass, velocity, and launching radius (e.g., Ivanova et al. 2013a;MacLeod et al. 2017;Blagorodnova et al. 2021). Accurate determinations of these properties can, in combination with information on the progenitor binary from preoutburst imaging, be used to constrain the physical processes at work leading up to and during the merger (e.g., MacLeod et al. 2017;Blagorodnova et al. 2021). Ivanova et al. (2013a) pioneered this approach by applying an analytic model for the luminosity and duration of Type II SNe (SNe II) light curves from Popov (1993) to infer the ejecta properties for a sample of LRNe. Using an updated version of the Popov (1993) formulae calibrated using SNe II radiative transfer calculations (Sukhbold et al. 2016), Blagorodnova et al. (2021) constrained the ejecta mass and launching radius for the AT2018bwo. Combined with an application of binary stellar evolution models to the observed yellow supergiant progenitor, the ejecta properties enabled them to paint a consistent picture of the event (see also MacLeod et al. 2017 for a similar analysis of M31-OT2015). Most LRN progenitor stars are moderately evolved, with radii between ∼ 1 − 10 times that of main sequence stars of their mass (Fig. 2).
Despite the advances of these past works, several physical inconsistencies can arise when applying SNe II models to LRN. SNe ejecta are dominated by radiation pressure, an assumption implicit in SN II light curve models (Popov 1993;Faran et al. 2019). By contrast, given the comparatively lower velocities, higher densities, and lower thermal energies of LRN ejecta, gas pressure can dominate over radiation pressure and sources of opacity other than electron scattering may be important. Furthermore, while hydrogen recombination could be an important source of luminosity for LRN (e.g., Ivanova et al. 2013a), recombination energy is irrelevant in SNe and is not included in models like Popov (1993). SN light curves are powered by the initial thermal energy of the recently-shocked ejecta from the explosion or the radioactive decay of 56 Ni. Although the early-time peaks seen in many LRN light curves can be powered by the initial thermal energy of hot ejecta (MacLeod et al. 2017;, the longer-lived plateau emission (sometimes manifesting as a second luminosity peak) is powered by hydrogen recombination energy or other forms of radially-distributed ejecta heating (e.g., shock interaction between the fast merger ejecta and circumbinary material from the pre-dynamical phase; Metzger & Pejcha 2017, or from an accretion-powered jet; e.g., Soker 2020; Soker & Kaplan 2021).
In this paper we present a new model for LRN light curves, which self-consistently accounts for the effects of arbitrary gas and radiation pressure, thermal energy released by recombination, and a realistic distribution of ejecta velocities. After presenting the light curve model in Section 2 we describe our results in Section 3. After walking the reader through a few example light curves (Sec. 3.1) and elucidating the physics giving rise to multiple light curve peaks, we perform a parameter study covering the full range of stellar merger properties (Sec. 3.2), apply our model to extract the ejecta properties from a sample of individual LRNe (Sec. 3.3), and compare the predictions of our light model to findings obtained using SN II models (Sec. 3.4). Section 4 discusses implications of our findings, including the need for alternative energy sources in LRN from circumbinary interaction (Sec. 4.1); the populations of binaries giving rise to LRNe (Sec. 4.2); and other applications of our models to planet-star or star-compact object mergers (Sec. 4.3). We briefly summarize our conclusions in Section 5.
2. LIGHT CURVE MODEL Our light curve model follows that developed by Metzger et al. (2021) in the different context of optical transients from the tidal disruption of a star by a white dwarf, but with added physical effects such as radiation pressure and Helium ionization (Kasen & Ramirez-Ruiz 2010). The general approach is to divide the radiallystratified ejecta into successive thin shells, calculate the emission from each such mass/velocity shell separately by means of a time-dependent one-zone model, and then obtain the total light curve by summing over the contributions from each shell. We begin by describing the one-zone model, as will ultimately be applied to each shell of velocity v and external mass M ej (> v). Immediately after the stellar merger or common envelope event, the evolution of the unbound ejecta will be complex, potentially involving self-interaction between different ejecta components (e.g., MacLeod et al. 2018b;Reichardt et al. 2019). However, we may reasonably assume that on a few dynamical times the outflow geometry will become quasi-spherical and the thermal and kinetic energies roughly equal. Hence we set the initial internal energy of the shell as E 0 = dM v 2 /2 (here dM is the shell's mass) and calculate the subsequent evolution starting from time t 0 = R 0 /v, where R 0 is the initial ("ejection") radius of the ejecta (typically equal to, or smaller than, the radius of the donor star). After the initial thermal energy has been converted into kinetic energy on a timescale ∼ t 0 , the shell will thereafter expand freely with radius R = vt.
The internal energy of the shell can be written and ε i are the mean ionization degree, abundance fraction, degree of ionization, and ionization energy of species i, respectively. We calculatex and x i by solving Saha equation taking into account singly-ionized hydrogen and helium, which we have checked capture recombination effects on the light curve to high accuracy. The other quantities are the total number of nuclei N , the Boltzmann constant k B , temperature T , radiation constant a, and shell volume V = 4πR 2 dR with the shell width dR. The internal energy evolves according to the first law of thermodynamics: where the effective adiabatic index is given by a function of density ρ = dM/V and temperature T (Krishna Swamy 1961; Kasen & Ramirez-Ruiz 2010) and is the ratio of radiation pressure P rad = aT 4 /3 to gas pressure P gas = N (1 +x)k B T /V . The radiative lumi-nosity of the shell is given by where E rad = aT 4 V and the photon diffusion time through the column associated with the external mass is and the factor t lc = R/c limits the energy escape time to the light crossing time, where c is the speed of light. It would be more accurate to express Eq. (6) as an integral over the diffusion time through the external layers > v accounting for their respective opacity κ and density distribution. However, this approach would require following energy transport between the shells via radiative diffusion; such a detailed treatment is beyond the scope of this paper. Hence we adopt this simpler single-zone expression for t d because we find the calculated light curve turn outs to be relatively insensitive to this assumption (a factor of 10 change in the pre-factor of Eq. (6) changes the light curve duration by a factor 2). We approximate the Rosseland mean opacity using the analytic formula similar to that used by   Opacity : [cm 2 /g] = 10 3 g/cm 3 10 5 10 7 10 9 10 11 10 13 Figure 3. Comparison of our assumed opacity law (dashed curve; Eq. 7) as a function of gas temperature and density with results from the OPAL opacity code (solid curves) for solar metallicity gas.
which accounts for the contribution from molecular opacity κ m 0.1 Z cm 2 g −1 , H − opacity κ H − 1.1 × 10 −25 Z 0.5 ρ 0.5 T 7.7 cm 2 g −1 , Kramers (boundbound and bound-free) opacity κ K 1.2 × 10 26 Z(1 + X)ρT −7/2 cm 2 g −1 , and electron scattering κ es 0.4xµ −1 cm 2 g −1 . Here Z is the metallicity, X is the mass fraction of hydrogen, and µ ≡ dM/m p N ( 1.25 for the solar metallicity) is the mean molecular weight neglecting the number of electrons, where m p is the proton mass. We neglect dust opacity, which will become important once the dense ejecta cools below the temperature for solid condensation; our model cannot therefore be applied to the late phases of LRN emission when dust formation is seen to occur (e.g., Wisniewski et al. 2008;Nicholls et al. 2013;Banerjee et al. 2015;Kasliwal et al. 2017). Figure 3 compares our assumed opacity law (Eq. 7) to that obtained using the OPAL opacity code for Z = Z 0.02. 1 Note that our prescription of the electron scattering opacity in Eq. (7) is different from that of . Our formula is more accurate than theirs when the recombination occurs at low densities ρ 10 −11 g cm −3 . We also adopted normalization of the Kramers opacity three times larger than that used in , to obtain better agreement with the OPAL calculations.
The total light curve is obtained by summing the luminosity contributions from each mass/velocity shell (see Metzger et al. 2021). For a given mass profile M ej (> v) of the form discussed below, we divide the ejecta into thin shells of velocity from v to v + dv, which gives the shell width dR = dvt. The mass contained in each shell 1 https://opalopacity.llnl.gov/existing.html#type1GN is given by dM = (dM ej /dv)dv and hence its density is given by ρ = dM/(4πR 2 tdv). With these quantities one can calculate the one-zone light curve for each shell dL(t; v) and the total light curve by summing over all such shells L(t) = dL(t; v).
The ejecta velocity distribution depends on the physical mechanism responsible for unbinding mass in stellar mergers and is uncertain theoretically. We consider two general functional forms (exponential and powerlaw) for the distribution of ejecta mass above a given velocity, where β ex , β pl , and v ej are parameters and the velocity range is v min < v < v max (v min = v ej for the power-law case). In what follows it will prove useful to introduce a mean ejecta velocityv E , which is defined by M ejv 2 E /2 = v 2 dM/2 and related to v ej entering Eq.
where Γ(x) is the gamma function and we have assumed v min v ej v max in the exponential case, and β pl > 2 and v max v min in the power-law case. To summarize, a given light curve model is defined by three main parameters: total ejecta mass M ej , initial (ejection) radius R 0 , and mean ejecta velocityv E (equivalently, total kinetic energy E K = M ejv 2 E /2), which is related to v ej through Eq. (9). In what follows, we will sometimes scale the ejecta mass to a fixed fraction η ≡ M ej /M of the mass of the donor star M . We shall adopt a fiducial value of η 0.25. This is somewhat larger than η 0.1 found in the common-envelope simulations by MacLeod et al. (2018b) for a fixed accretor to donor mass ratio, 1/3. The envelope ejection efficiency η may also be even smaller 0.1, in the case of main sequence stellar mergers (e.g., Lombardi et al. 2002;Glebbeek et al. 2013) or mergers of high mass-ratio binaries (e.g., Blagorodnova et al. 2021).
We likewise will frequently scale the mean ejecta velocityv E to the escape velocity of the donor at the outflow launching radius, where G is the gravitational constant. If R 0 equals the radius of the donor star R , then v esc = v esc, is the surface escape speed. However, for R 0 R then v esc v esc, . Ejecta speeds spanningv E ∼ 0.3 − 3 v esc, are supported by LRN observations (Fig. 2).
The choice of velocity distribution (Eq. 8) introduces additional parameters, {β ex (β pl ), v min , v max }. Typically, we set v min = 0.1v E and v max = 3v E but the light curve is not sensitive to these choices as long as v min v E v max . The calculation also depends on the assumed composition of the ejecta through the opacity; we assume solar metallicity (X = 0.74 and Z = 0.02) throughout this work, motivated by the moderately-evolved donor stars of observed LRNe progenitor binaries (Fig. 2).
3. RESULTS 3.1. Example Light Curves Figure 4 shows an example light curve, calculated for a total ejecta mass M ej = 1 M , initial radius R 0 = 10 R , and velocityv E = 300 km s −1 . We employ the exponential velocity distribution, M ej (> v) ∝ e −(v/vej) 3 (Eq. 9). A black line shows the total light curve, while the colored lines show the separate contributions from distinct mass shells (from v = 50 km s −1 to v = 650 km s −1 with an interval of 50 km s −1 and shell width of dv = 5 km s −1 ).
The light curve is characterized by an initial peak of high luminosity and short duration 1 week followed by longer duration plateau lasting several months. Such double-peaked light curve behavior is common in LRNe (MacLeod et al. 2017;Pastorello et al. 2021a; Fig. 1). As we now describe, these two components have distinct physical origins.

Early Cooling Peak
At early stages of evolution, the fastest, lower-mass layers of the ejecta dominate the transient luminosity. As we will confirm below, these layers are typically dominated by radiation pressure and their energy source is "cooling" emission produced by the initially hot ejecta (e.g., MacLeod et al. 2017), similar to the mechanism responsible for powering the plateau of SNe II (Arnett 1980). The duration and luminosity of this early emission phase are given by (e.g., Arnett 1980;Padmanabhan 2000) t pk κ es M 4πcv E pk L pk t pk 6 × 10 44 erg ×  Example light curve for fiducial parameters: Mej = M , R0 = 10 R ,vE = 300 km s −1 (vej 320 km s −1 ). Colored curves denote one-zone light curves for shells with v = 50 to 650 km s −1 (violet to orange) with a shell width dv = 5 km s −1 . The total light curve is characterized by two peaks, early cooling peak ( 10 d) and late recombination plateau, produced by distinct physical processes. The gray dashed curve shows an otherwise equivalent model but neglecting energy release from recombination. The approximate effective temperature of the emission is shown in the top panel (Eq. 19).
where M is the ejecta mass of the fastest layer, R pk = vt pk , and we have assumed that electron scattering of ionized hydrogen κ es 0.32 cm 2 g −1 dominates the opacity during the early peak.
The assumptions entering Eqs. (11)- (14) can now be checked. The ratio of radiation to gas pressure, on the timescale ∼ t pk , can be estimated as where (assuming P rad P gas ) we have used M v 2 /2 = (4π/3)R 3 0 aT 4 0 to derive the initial temperature and T (t pk ) = T 0 (R 0 /R pk ). In the second line of Eq. (14) we have set M = M ej and scaled the ejecta mass to donor star mass M via η = M ej /M (η 0.25 = η/0.25) and the ejecta velocity to the escape speed from the outflow launching radius (Eq. 10). Note that as long as the radiation pressure dominates initially, the pressure ratio and effective adiabatic index γ3 (right panel). Dotted and dashed contours indicate different hydrogen ionization fractionsx and ratios of gas-to-radiation pressure. Prior to hydrogen recombination, each shell decompresses and cools adiabatically with γ3 ∈ [4/3, 5/3], moving from the upper right hand corner of the plot to the lower left. However, once recombination starts at T 10 4 K, the energy released from recombination holds the temperature roughly constant at T 10 4 K (γ3 1.1), causing a convergence of the mass shell trajectories. However, as the density continues to drop, once P rad > Pgas (ρ 10 −11 g cm −3 for T ∼ Ti 5000 K), γ3 again starts to rise, causing T to decrease faster and recombination to complete. The resulting decrease in the opacity leads to the light curve peaking, thus motivating our analytic estimate of the transient duration (Eq. 16). T/10 4 K /10 11 g cm 3 x /cm 2 g 1 P rad /P gas R lum /10 15 cm We see that the radiation pressure will dominate in the high-velocity tail of the ejecta v v esc and that Eqs. (11)-(14) indeed determine the early-time peak emission from these layers (and hence the early-peak of the transient as a whole). However, for the bulk of the ejecta for which v v esc , gas pressure will dominate and the early cooling emission from these layers will be reduced compared to L pk in Eq. (12) because of the greater adiabatic losses experienced by the ejecta for γ = 5/3 (T ∝ R −2 ) than γ = 4/3. Previous models for LRNe employ analytic formula for the luminosity and duration (e.g., Popov 1993) which implicitly assume radiation pressure dominates, as is justified in SNe due to their higher ejecta velocities and specific internal energy. The ratio of Kramers to electron scattering opacity during the early peak can also be checked: where again we take κ es 0.32 cm 2 g −1 . This confirms that κ es κ K , as assumed in Eqs. (11)-(14).

Late Recombination Plateau
The longer plateau-like emission phase in Fig. 4 is not powered by the ejecta's initial thermal energy, but instead by hydrogen recombination. After forming its initial luminosity peak, the ejecta cool down and hydrogen starts to recombine at temperatures T ∼ 10 4 K. Once recombination sets in, the additional energy supply causes the adiabatic index to decrease to γ 3 −1 0.1, thus slowing the temperature evolution greatly compared to that of ordinary adiabatic expansion absent recombination (γ 3 5/3). By the time a given shell completes recombination at temperatures T i 4000−6000 K most of the radiation energy has been emitted, generating a second peak in its light curve. The sum of the second peak from all mass shells powers the ∼ 100 day plateau in Fig. 4.
The duration of the plateau is roughly given by the timescale over which the recombination completes at T T i . We find that recombination completes at roughly a constant characteristic density ρ i ∼ 10 −11 g cm −3 . As explained in Fig. 5, this density corresponds to roughly where P rad ∼ P gas for T ∼ T i . The increasing adiabatic index once P rad P gas triggers the completion of ionization of each mass shell, causing the opacity to decrease and leading to the light curve peak.
The duration of the recombination plateau is thus set by the timescale for the density to reach ρ i , where ρ i,−11 = ρ i /10 −11 g cm −3 . Since the plateau is powered by the hydrogen recombination (ε H = 13.6 eV per hydrogen atom), the radiated energy and luminosity can be written where the dimensionless parameter f ad = 0.3f ad,0.3 < 1 accounts for inefficiency in radiating the recombination energy due to adiabatic losses. Helium recombination also contributes, but since it occurs at higher temperature (at hence much higher optical depth), its effect on the light curve is small compared to that of hydrogen recombination.
In Appendix A we present a one-zone toy light curve model for a hydrogen-rich explosion, in the spirit of Popov (1993) but assuming gas-pressure instead of radiation dominates and accounting for recombination energy through an (assumed constant) effective value of the adiabatic index γ 3 ; for values of γ 3 ∼ 1 expected during recombination the analytic expressions for the plateau luminosity (Eq. A14), duration (Eq. A13), and radiated energy (Eq. A15), are quite similar to those Eqs. (16)-(18) above. This is remarkable because the conceptual motivation of the two derivations are different: the plateau duration derived in Appendix A (Eq. A13) is dictated by the photon diffusion timescale through the recombining ejecta shell, while Eq. (16) was derived assuming recombination occurs at roughly fixed density.
The top panel of Fig. 4 shows the evolution of the effective temperature, which we have defined by where σ SB is the Stefan-Boltzmann constant and R lum is the radius of the shell contributing the majority of the luminosity at any epoch. During the initial peak, T eff 10 4 K before cooling during the recombination plateau to an almost fixed valued T eff 3000 − 4000 K, consistent with the "red" colors of LRNe. This temperature stability occurs for the same reason as in SNe II, namely that the recombining shell recedes monotonically in the velocity coordinate but the expansion effect compensates and leads to the almost spatially fixed location of the shell (see Fig. 6). We emphasize however, that−unlike in SNe−the plateau is powered not by radioactivity or the thermal energy released from the initial explosion, but rather from the recombination energy itself. We confirm this by showing an otherwise identical light curve model but calculated with recombination energy artificially excluded, as a gray line in Fig. 4, demonstrating the absence in this case of an extended plateau phase. This further implies that standard relationships applied to model SN light curves (e.g., Popov 1993) should not be applied to LRNe. Figure 6 shows the time evolution of key quantities for the fiducial model, as measured at the shell dominating the luminosity at each epoch (R lum ). As expected, P rad P gas during the early peak when the highest velocity layers dominate the emission, while P rad P gas during the later recombination peak. Also note the mean ionization fractionx is low 10 −2 during the pro-longed recombination plateau, and that electron scattering opacity κ es 0.32x cm 2 g −1 becomes comparable to H − opacity. During the recombination phase, the effective adiabatic index γ 3 4/3 due to the energy released by recombination offsetting losses due to adiabatic expansion. Note that although the index γ 3 of R lum is close to 4/3, this is not related to radiation pressure. Rather, during the brief epoch when a shell dominates the total luminosity, its adiabatic index quickly increases from γ 3 1.1 to 5/3 passing through an value coincidentally close to 4/3 (as illustrated by a dash-dotted curve in Fig. 6). Figure 7 shows light curve models calculated for ejecta properties otherwise identical to the fiducial model (exponential model with β ex = 3) but varying the shape of the velocity distribution at fixedv E (Eq. 8). For the exponential velocity profile, the duration of the plateau shortens with increasing values of β ex , but the qualitative shape does not change significantly unless β ex 2. Likewise, in the case of a power-law profile, the light curve shape is not qualitatively dependent on the powerlaw index β pl . Motivated thus, we hereafter fix the velocity profile to the fiducial exponential β ex = 3. Figure 8 compares the light curve for the fiducial model (M ej = 1 M ; R 0 = 10 R ;v E = 300 km s −1 ) to otherwise similar models but for which {M ej , R 0 ,v E } are tripled from their fiducial values. As expected from Eqs. (11)-(14), larger values of R 0 andv E increase the luminosity of the early peak. As expected from Eqs. (16)-(18), the properties of the recombination plateau are roughly independent of R 0 . Likewise, the plateau gets shorter but more luminous with increasing velocity, while both the duration and luminosity increase for larger ejecta mass.

Variations About the Fiducial Model
As can be seen in Fig. 1, many LRNe light curves exhibit multi-peak light curves in which the first peak is substantially brighter than the longer plateau to follow (e.g., V838 Mon: Munari et al. 2002;M31-OT2011: Kurtenkov et al. 2015MacLeod et al. 2017;AT2019zhd: Pastorello et al. 2021aAT2020kog: Pastorello et al. 2021b) and others in which a distinct first peak is more subtle or non-existent (e.g., V1309 Sco: Mason et al. 2010;Tylenda et al. 2011;AT2020hat: Pastorello et al. 2021b. This difference may indicate the presence or absence of a high-velocity ejecta component or different progenitor's properties. To understand the latter effect, using Eqs. (12) and (18), we see that the ratio of the cooling envelope peak luminosity to that of the recombination plateau can be estimated as where we have taken ρ i,−11 = 1 and assumed that the outer layers of the ejecta responsible for powering the first peak have a velocity twice as large as that of the bulk, v 2v E . Thus, a range of values L pk /L pl ∼ 1−10 can be obtained.

Parameter Study of Stellar Mergers Across the HR Diagram
We now explore the light curve behavior across a larger parameter space of ejecta mass M ej ∼ 10 −3 − 10 M and initial radius R 0 ∼ 10 −2 − 100 R covering those expected from a wide range of stellar and even sub-stellar merger events. For each model we set the mean ejecta velocity to the escape speed from R 0 , v E = v esc (Eq. 10), where the stellar mass is related to the ejecta mass by M ej = ηM with η = 0.25. For each light curve model, we calculate the total radiated energy E rad = Ldt. We also define a characteristic transient duration, t 90 , as the interval over which 90% of the total energy is radiated, 0.9 E rad = t90 Ldt. Finally, from the radiated energy and duration we define a mean luminosity L 90 ≡ E rad /t 90 . Figure 9 shows the distribution of the characteristic luminosity, duration, and radiated energy. We also plot the ratio of the radiation energy to the total hydrogen recombination energy f ad ≡ E rad /(ε H XM ej /m p ). For comparison we show tracks of stellar radius versus terminal-age main sequence (TAMS), from the MIST database of MESA stellar evolution models (Dotter 2016;Choi et al. 2016;Paxton et al. 2011Paxton et al. , 2013Paxton et al. , 2015Paxton et al. , 2019, as well as a hypothetical state of post-main sequence evolution in which the stellar radius is 10 times larger than its TAMS value ("10TAMS"). These bracket the state of LRN progenitor stars from pre-transient imaging studies (Fig. 2). For matter ejected from the stellar surface, we expect R 0 ∼ R , while for matter ejected from deeper within the donor star R 0 < R .
Across most of the parameter space, we find that the total radiated energy is dominated by the later recombination plateau, consistent with a "by eye" integration of LRN light curves (Fig. 1). The parameter dependencies of the luminosity, duration, and radiated energy are roughly consistent with those expected from Eqs. (16)-(18). In particular, forv , and E rad ∼ E pl ∝ M , consistent with the contours in Fig. 9. These relations no longer hold when the ejecta is already partially ionized at the initial time (upper left and bottom regions in Fig. 9). Figure 10 shows the ratios of the transient luminosity and duration as predicted by our analytic formulae (Eqs. 11,12) to those obtained from the full numerical calculations (Fig. 9), demonstrating good agreement (to within a factor 3) across a band encompassing the entire TAMS sequence.
, while the right panel shows the impact of varying the power-law index β pl for Mej(> v) ∝ v −β pl . Although the light curve shape changes for different velocity distributions, the qualitative features of an early peak followed by a second peak or plateau are robust to the details.  Fig. 4) to otherwise similar models in which: the ejecta mass is increased to 3 M (blue); the initial radius is increased to R0 = 30 R (red); the mean velocity is increased tovE = 900 km s −1 (orange). The trends in the luminosity and duration of the first and second peaks follow Eqs. (11), (12) and Eqs. (16), (18), respectively.
For mergers that occur when the donor is near the end of the main sequence (R ∼ R TAMS ) or after evolving moderately up the giant branch (R ∼ 10R TAMS ) the predicted luminosities are in the range ∼ 10 39 − 10 40 erg s −1 , broadly consistent with the bulk of the LRN population. Luminosities up to ∼ 10 41 erg s −1 may be possible in extreme cases for efficient ejection of the stellar envelope η ∼ 1 from deep within the potential of the stellar core (R 0 R TAMS ). The minimum launching radius of the donor star envelope for an evolved star is determined by the size of the final He core, which we show as a white dash-dotted curve in Fig. 9 (we adopt the analytic M − M He relationship from Woosley 2019 and employ the M He − R relationship from Eq. 4 of Schaerer & Maeder 1992; note that η 0.5 − 0.9 in this case). Even under the most optimistic assumptions, the model is strained to reach luminosities L > 10 41 erg s −1 , which can therefore be taken as a theoretical upper limit (absent other energy sources not included in our fiducial model; Section 4.1).
Typical values of the radiative efficiency of the recombination energy are f ad = E rad /(ε H XM ej /m p ) 0.1 − 0.3. Across most of the parameter space, the luminosity of the initial thermal peak L pk will stick out above that of the recombination plateau L pl , even though the radiated energy is dominated by the plateau. The parameter space for which L pk > L pl is delineated by a solid black line in Fig. 9, roughly consistent with our analytic estimate in Eq. (20). Note however that having L pk > L pl does not necessarily imply the early cooling peak will be detected, because its duration is significantly shorter than the plateau, t pk /t pl 0.01 − 0.1, and could be missed depending on the discovery epoch and observational cadence of the survey (Fig. 1). . Characteristic luminosity, timescale, radiated energy, and ratio of radiated energy to recombination energy calculated for a range of stellar (donor) mass and radius for our fiducial velocity profile (Eq. 8). Ejecta mass is linked to the stellar mass via Mej = ηM (η = 0.25) and velocity is set to the escape velocity. White curves denote the mass and radius relation for TAMS, ten times larger radius of TAMS ("10TAMS"), and brown dwarf and giant planet. A white dash-dotted curve shows the radius of He core for the corresponding ZMAS mass (M ), which gives a rough lower limit of the mass-ejection radius R0 from an evolved star. A black solid line denotes the contour along which the peak luminosities of the initial cooling phase and later recombination plateau are equal (see Fig. 4), as calculated directly from the light curve models. A gray dotted curve shows the contour on which the mean ionization degree of the typical shell (v =vE) isx = 0.9 (partially ionized), at the initial time. Figure 11 shows the mean luminosity L 90 and duration t 90 for a sample of LRNe (Table 1; Fig. 1) as a function of the donor mass M inferred from pre-transient imaging observations when available. 3 Observed quan- 3 We note that the measured duration t 90 and accordingly the radiated energy E rad may be underestimated in some events because of incomplete light curve coverage (some events were discovered only at the the end of the plateau phase and others were not observed until the end of the plateau; see Fig. 1). Fig. 11 we show our light curve predictions (Fig. 9) assuming M ej = ηM andv E = v esc, = (2GM /R ) 1/2 for stars on the TAMS (red curves) and with radii 10 times their TAMS radius (10TAMS; blue curves), for a range of values η = 10 −2 (thinnest lines) to η = 1 (thickest lines). For fixed η, the theoretical relations roughly follow the trends L 90 ∝ M and t 90 ∝ M 1/2 predicted by Eqs. (16,18).  Fig. 9 but showing the ratio of the plateau luminosity and duration predicted by our analytical formulae (Eqs. 16,18; for fixed f ad,0.3 = 1 and ρi,−11 = 1) to those obtained from the full numerical model (Fig. 9). The analytic formulae accurately predict the luminosity and duration within a factor of 3 along the parameter space surrounding the TAMS sequence. Figure 11 shows that, in broad brush terms, our calculations can reproduce the observed luminosities ∼ 10 38 − 10 40 erg s −1 of most LRNe, and also the general trend of higher luminosity and longer duration for more massive stellar progenitors (e.g., Kochanek  , with luminosities that approach or exceed even the η = 1 TAMS prediction for very massive stars. We return to this tension regarding the implied ejecta masses of some LRNe below.

tities are shown by circles or dot-dashed lines (when no information on M is available). For comparison in
We can determine the ejecta properties for individual LRNe by inverting our analytic formulae for the transient luminosity and duration (Eqs. 16, 18), where we note that our calculations reveal that f ad,0.3 increases moderately ∼ 0.5 − 2 with increasing stellar mass and radius (Fig. 9). Eq. (22) provides an independent velocity estimate from that obtained from spectroscopy of the transient, v obs , the consistency between which provides a check on the model (see below).
The ejecta launching radius R 0 can likewise be inferred from the observables in two ways: (1) assumingv E = v esc = (2GM /R 0 ) 1/2 using the mea-sured M when available, (2) assuming v obs = v esc , Figure 12 compare the inferred ejecta properties {M ej ,v E , R 0 } to the stellar properties or transient properties {M , R , v obs } for events for which the latter information is available (see also Table 1). At a basic level most of our results for the baseline model (shown as circles or triangles in Fig. 12) make physical sense. Inferred ejecta masses range from 10 −2 M in V4332 Sgr to ∼ 10 M in AT2018bwo, for example. Likewise, the implied ejecta launching radii are all safely within the donor's interior, R 0 < R . We return to what inferences can be drawn about the LRN population based on these data in Section 4.2.
Upon closer inspection, several inconsistencies or tensions are also noticeable in Fig. 12. Firstly, the implied ejecta masses from the most massive binary systems (e.g., NGC4490-OT2011 and AT2017jfs) are enormous 100 M and greatly violate the theoretical limit M ej > M (η > 1) in which the entire mass of the donor is ejected, while a comparable number of events hover around η ∼ 1 (see also Soker 2020). A second point of tension is that the inferred ejecta velocities (estimated from Eq. 22) systematically exceed those observed,v E ∼ (2 − 3) v obs . Equivalently, the ejection radius R 0 one infers by equatingv E to the local escape velocity (Eq. 23) is smaller than that implied if the observed spectroscopic expansion velocity equals the escape speed at the ejection radius (squares in bottom panel of Fig. 12).
Keeping in mind that our formulae (Eqs. 21 and 22) exhibit errors of a factor of a few, the latter mentioned velocity/radius discrepancy could in principle have several origins: (1) the observed spectroscopic velocities underestimate the true energy-weighted ejecta speed, for instance due how excitation and optical-depth affect what regions of the ejecta given lines probe; (2) inaccuracies in our light curve model cause an over-estimate of the plateau duration, which results in a large velocity for a given observed duration (see Eq. 16);this could result from one or more of our simplifying assumptions, e.g. spherical symmetry; (3) our model is missing physics, such as additional sources of the ejecta heating beyond that from hydrogen recombination. The latter possibility would also have implications for the inferred ejecta mass.
To illustrate this with an example, consider a hypothetical scenario in which the ejecta is heated during the plateau phase at a specific rate ∼ 10 times larger than from hydrogen recombination (e.g., due to shock interaction with circumbinary material; Section 4.1). In such a case the true ejecta mass would be ∼ 10 times smaller than inferred from Eq. (21) from the observed transient duration and luminosity; hence, for a fixed observed plateau duration t pl (Eq. 16) the implied ejecta speed v E ∝ M 1/3 ej would be ∼ 10 1/3 ∼ 2 times smaller than previously estimated neglecting the additional heating source, bringing it into closer agreement with v obs (and likewise, the implied R 0 obtained assumingv E = v esc would be ∼ 4 times smaller).
To explore this possibility further, the top left panel of Fig. 12 includes a second estimate of the ejecta mass (shown as squares and diamonds), which makes use of our expression for the plateau duration (Eq. 16) but supplants Eq. (21) (which assumes recombination as the sole energy source) with the assumption thatv E = v obs , i.e. that the ejecta speed equals that observed from the transient spectroscopy. Comparing the circles and squares (or triangles and diamonds), we observe that the required value of M ej decreases significantly, bringing even the most luminous LRNe near or into the physically allowed region (η 1). Under the same assumption (v E = v obs ), the implied mass ejection radius increases but remains below the progenitor star surface for most events (squares in the bottom panel of Fig. 1).
To summarize, the ejecta masses one infers from the optical light curve assuming hydrogen recombination as the only energy source (as in our model) represent only upper limits. If other sources of heating not accounted in our model are present, then (1) M ej will be lower than this upper limit and; (2) the ejecta velocitiesv E implied by the observed LRN duration are brought into better agreement with those observed spectroscopically.

Comparison to SNe II Analytic Formulae
The facts that recombination energy dominates the luminosity and gas pressure can exceed radiation pressure during the plateau phase (Fig. 6), both violate implicit assumptions of analytic SNe light curve models (Popov 1993;Sukhbold et al. 2016). Because several past works use these formulae to infer the ejecta properties from observed LRNe (e.g., Ivanova et al. 2013a;MacLeod et al. 2017;Blagorodnova et al. 2021), or predict observables in a population analysis (e.g., Howitt et al. 2020), it is important to quantify the implied inaccuracies. Figure 13 shows the ratio of the transient luminosity and duration as extracted from our light curve calculation ( Fig. 9) to those predicted by Popov (1993). For the latter we employ results calibrated using SN II radiative transfer calculations (Sukhbold et al. 2016;Blagorodnova et al. 2021), Remarkably, Fig. 13 shows that the calibrated Popov formulae (Eqs. 25,26) agree with our model predictions across most of the relevant parameter space to within a factor 3. However, this agreement is largely a coincidence because the models are derived under completely different assumptions (Popov's formulae "know" nothing about hydrogen recombination energy, for instance). The Popov formulae also exhibit different parameter dependencies from our model (Eqs. 16,18), however they (again, coincidentally) become similar for v E ∼ v esc ∝ (M /R 0 ) 1/2 : L Popov ∝ M   (21) and (22)  The inferred values M ej , R 0 thus depend sensitively on the observables, particularly the ejecta velocity. This is problematic because the measured velocity (e.g., as typically inferred from spectral line widths at particular epochs) may not coincide precisely with the bulk of the ejecta. By contrast, the sensitivity of inferred ejecta quantities to LRN observables as predicted by our model (Eqs. 21,22,23), particularly if the donor star properties are independently constrained, allow for a more robust extraction of {M ej ,v E , R 0 }, compared to Eqs. (27,28). Figure 14 compares the inferences on {M ej , R 0 } made employing our model (Table 1) to those obtained by instead applying the modified Popov (1993) formulae (Eqs. 27,28). We see that the ejecta masses implied by the Popov formulae can over-or under-estimate those of our model (circles) by more than 1-2 orders of magnitude. In particular, the Popov formulae underestimate the implied M ej for the most luminous LRNe (largest M ej events), masking the M ej > M tension implied by our modeling (Sec. 3.3). In applying the Popov formulae here, we have estimated the ejecta mass and radius using the observed spectroscopic velocityv E = v obs , which typically differs from the surface escape speed (Fig. 2). The sensitive dependence of the Popov formulae velocity thus accounts for the larger deviations of the ejecta properties from our estimates in Fig. 14 than naively suggested by Fig. 13 which instead assumedv E = v esc .

DISCUSSION 4.1. Shock Interaction as a Luminosity Source
Hydrogen recombination typically dominates the radiated energy in our baseline model because the ejecta's initial thermal energy is mostly lost to adiabatic expansion before it can be radiated. However, some of this thermal energy can be regenerated if the merger ejecta were to collide with slower material at larger radii R surrounding the binary ). This pre-existing slower material could arise from an earlier episode of mass-transfer leading up to the dynamical merger event, for instance a brief thermal-timescale mass-transfer phase (e.g., Blagorodnova et al. 2021) or other physical processes generating mass-loss from the binary's L 2 point (Pejcha et al. 2016a. The presence of pre-dynamical ejecta is supported observa-  AT2018bwo AT2018hso AT2019zhd AT2020kog AT2020hat Figure 11. Luminosity L90 (left) and duration t90 (right) of stellar merger transients compared to LRNe with progenitor donor masses M inferred from pre-transient imaging (Table 1). Red and blue shaded regions show the predictions of our light curve model as a function of donor mass M for a range of η = Mej/M = 0.01 (thinnest lines) to η = 1 (thickest lines), and where we have assumedvE = vesc, = (2GM /R ) 1/2 for stars on the TAMS or 10 times their TAMS radii. In the left panel, a dashed curve denotes the stellar luminosity at TAMS (Hurley et al. 2000), while a dotted line shows the luminosity at the tip of the red giant branch (Paczyński 1970;Cummings et al. 2018).
tionally by the slow early rise in the pre-maximum light curves of LRNe such as V1309 Sco (Pejcha 2014;Pejcha et al. 2017) and the redder than expected early-time colors of V838 Mon and AT2018bwo (Blagorodnova et al. 2021). The characteristic mass M pre of the slowly-expanding (velocity v pre ) material surrounding the binary at the time of merger is uncertain, but in many cases is expected to be smaller than the ejecta mass released at later times, M ej (e.g., at most M pre ∼ 0.1M is predicted via L 2 mass-loss; MacLeod et al. 2017;Pejcha et al. 2017;MacLeod & Loeb 2020; however, see Blagorodnova et al. 2021). Assuming M pre M ej and v pre v E , and that the shock-interaction is momentum conserving, the thermal energy released by the collision is given by This exceeds the energy released by hydrogen recombination, ε H XM ej /m p , for Thus, if even a few percent of the ejecta mass is present on radial scales R pre ∼v E t pl ∼ 10 14 − 10 15 cm prior to the merger, the shock will compete with recombination in powering the light curve. A related, alternative heating mechanism for the LRN ejecta is a fast outflow or jet powered by accretion onto the accretor star or donor star core (e.g., Soker & Kashi 2016;Soker 2020;Soker & Kaplan 2021).
If the duration of the transient is still set by hydrogen recombination in this scenario, then the transient duration is still given by t 90 ∼ t pl ∝ M 1/3 ejv −1 E (Eq. 16). However, instead of Eq. (17) for the plateau luminosity, the shock scenario predicts E pl ∝ E sh ∝v 2 E and hence The value of L sh pl can now reach 10 41 erg s −1 for events like AT2017jfs which possess a relatively high velocity v obs ∼ 700 km s −1 , and which otherwise required unphysically large ejecta masses when considering recombination energy alone (Fig. 12). Depending on its opacity κ, the optical depth through the circumbinary material, Ejecta velocity : Recombination heating only vE = vobs   AT2020kog AT2020hat Figure 14. For the same sample of LRNe in Fig. 12, here we compare the ejecta masses and launching radii inferred from our light curve model, Mej/R0 on the vertical axis (using the two methods described in Fig. 12 with the same symbols notation) to those obtained using the Popov (1993) supernova light curve formulae on the horizontal axis (Eqs. 27,28).
is sufficiently large to redden the early transient emission (e.g., Blagorodnova et al. 2021). Future work is required to incorporate the time-dependent shock luminosity, given an assumed density profile of the premerger cirumbinary medium (e.g., , self-consistently into our light curve model framework (Section 2). One might also expect shock interaction to manifest through observational signatures other than just the light curve, such as optical emission lines and thermal X-rays from hot shocked gas and non-thermal radio synchrotron emission from relativistic electrons accelerated at the shock (as seen from internal shocks in classical novae, for example; Chomiuk et al. 2021). However, the equatorial geometry of the slow circumbinary medium (e.g., Pejcha et al. 2016a) implies that the polar regions of the merger ejecta and hence the optical photosphere will quickly "wrap around" the equatorial shock, blocking signatures of the latter from most viewing angles (see Andrews & Smith 2018 for an example of this in the context of shock-powered Type II supernovae). At late times after recombination completes, the opacity will drop, the photosphere recedes and residual shock signatures may become visible; however, the high opticalwavelength opacity once dust formation begins around the same time may keep these inner regions opaque and continue to obscure shock-related signatures.

Inferences about the LRN Population
A goal of measuring the ejecta properties for a sample of LRNe is to deduce information about the physics of the merger process as well as the demographics of merging binary systems.
Putting aside the ambiguities that arise in measuring individual ejecta masses in the face of additional sources of ejecta heating, the top panel of Fig. 12 reveals an interesting trend that is likely to remain robust in spite of these uncertainties: the fraction of the stellar envelope being ejected η = M ej /M systematically increases with more massive donor stars. This arises because the average luminosity of LRNe increases sharply with stellar mass, L ∝ M α with α ∼ 2 − 3 (e.g., Kochanek et al. 2014), while recombination powered models predict a sub-linear dependence of L pl ∝ M 2/3 ej (Eq. 18). What could give rise to the increasing η(M ) trend?
The fraction of the donor's star envelope unbound during the merger process is likely to increase with secondary mass and hence the mass ratio of the binary (e.g., Nordhaus et al. 2010). The fact that the average mass ratio of interacting binaries is expected to increase with stellar mass (Moe & Di Stefano 2017) could therefore contribute to η increasing with M . Another contributing factor is that more massive LRN progenitors appear to be further evolved off the main sequence (Fig. 2): more evolved stars possess a more pronounced coreenvelope structure, facilitating the ejection of a higher fraction of the envelope as the accretor releases more gravitational energy spiraling into the compact core. A smaller ejection radius also increases the ejecta velocity, which acts to increase the transient luminosity (Eq. 18). Indeed, the most luminous events tend to have larger observed velocities (see Table 1).
The current LRN sample could also be biased by a selection effect: since the intrinsic nuclear luminosity of the star at all phases of evolution increases steeply with stellar mass (dotted and dashed lines in Fig. 11), low-M ej (low-η) events could be selected against in their discovery because of the lower contrast between the transient and star luminosity. Moreover, since massive progenitors are rare, their outbursts are more typically detected as extragalactic transients; hence, it is natural to miss low-luminosity events from massive stars. These would imply the existence of a hidden population of low-M ej mergers ("minor mergers") from massive stars. This hidden population would need to be taken into account when comparing the M -dependent rate of LRN to the initial mass function of stars or to rate predictions from binary population synthesis modeling (e.g., Kochanek et al. 2014;Howitt et al. 2020).

Stellar Mergers with Planets and Compact Objects
Our model can be applied to predict the optical light curves from other types of mergers involving mass ejection from hydrogen-rich bodies.
One example is the merger of a giant planet or brown dwarf with a more massive host companion star (e.g., Metzger et al. 2012). When the average density of the planet is larger than that of the star, the energy released as the planet plunges slowly below the stellar surface can give rise to mass ejection (and a longer-term evolution of the star's luminosity due to energy released deeper in the stellar envelope; MacLeod et al. 2018a). Alternatively, if the planet is less dense than the star (e.g., "inflated" Hot Jupiters; Baraffe et al. 2010), then it can be tidally disrupted above the surface, forming an accretion disk around the star. Outflows from the resulting super-Eddington accretion flow are likely to unbind a significant fraction of its mass from the planet (e.g., Metzger et al. 2012). In both cases then, a fraction of the planet/brown dwarf mass ∼ 10 −3 − 10 −2 M is ejected from radii R 0 ∼ R at close to the surface escape speed. Figure 11 shows that for M ∼ M and η ∼ 10 −3 − 10 −2 , the predicted transient lasts days to weeks with a luminosity L ∼ 10 37 − 10 38 erg s −1 . These light curve properties are comparable to those of classical novae (e.g., Strope et al. 2010), with which planet-star mergers could easily be confused. Based on the decay rate of observed gas giant planets orbits due to tidal dissipation, Metzger et al. (2012) estimate a planet-star merger rate of ∼ 0.1 − 1 per year in the Milky Way, consistent (within an order of magnitude) of the stellar merger rate (Kochanek et al. 2014) and ∼ 100 times lower than the rate of classical novae (e.g., Kawash et al. 2021). Given their low luminosities, and low relative frequency to classical novae, it is unsurprising that no confirmed planetstar mergers have yet been discovered.
Our model can also be applied to mergers of hydrogenrich stars with compact objects, such as white dwarfs (WD; e.g., Shara & Shaviv 1977;Michaely & Shara 2021;Metzger et al. 2021), neutron stars (NS), or stellar-or intermediate-mass black holes (BH; e.g., Perets et al. 2016;Fragione et al. 2020;Kremer et al. 2021Kremer et al. , 2022, which give rise to up to several solar masses of ejecta expanding at velocities of hundreds to thousands of km s −1 . Due to the deeper potential well of the WD/NS/BH, the kinetic energies of the ejecta from these events can exceed those of ordinary star mergers; however, the small ejection radii−less than the tidal disruption radius, itself typically comparable to the stellar radius−imply that most of the initial thermal energy will be lost to adiabatic expansion, as in the stellar mergers we have studied. As a result, the total radiated optical energy from the unbound ejecta will also be limited to that from its hydrogen recombination, i.e. f ad ε H XM /m p ∼ 10 46 erg per solar mass of ejecta. The predicted luminosity from star-WD/NS/BH mergers are thus unlikely to exceed ∼ 10 39 −10 40 erg s −1 , less than estimates that do not account for adiabatic losses prior to recombination (e.g., Michaely & Shara 2021). Of course, if as may be the case also in stellar mergers (Sec. 4.1), star-WD/NS/BH mergers are accompanied by an extra central luminosity source which re-energizes the unbound ejecta (e.g., internal shocks due to ongoing outflows from the central engine; e.g., Dexter & Kasen 2013;Margalit & Metzger 2016;Gottlieb et al. 2022) then the transient luminosities could be considerably higher than from recombination power alone.

CONCLUSIONS
We have presented a new model for the light curves of LRNe from stellar mergers, which self-consistently accounts for the range of relevant opacities, an arbitrary degree of radiation and gas pressure, and the energy released from hydrogen recombination. We apply our model to a sample of 14 LRNe (Fig. 1, Table 1) to infer the ejecta properties (mass, velocity, and launching radius) and compare to the mass/radios of the progenitor donor star when available (Fig. 2). Our conclusions can be summarized as follows: • Our model naturally reproduces the observation (Fig. 1) that LRN light curves typically possess two luminosity peaks: (1) a short initial peak typically lasting days powered by the initial thermal energy ("cooling emission") from the outermost fastest layers of the ejecta; and (2) a longer peak or plateau lasting weeks or months, powered by energy released from hydrogen recombination from the bulk of the mass. While radiation pressure can dominate in the ejecta during the first peak, gas pressure typically dominates the latter plateau, at odds with the common assumption of supernova light curve models (e.g., Popov 1993).
The potential for two light curve peaks was recognized previously as a result of distinct physical stages during the merger (MacLeod et al. 2017;. Here, we emphasize (again) that such a double-peaked feature results due to the stratification of even a single spherically symmetric ejecta shell into layers with distinct velocities and initial thermal energy content. Multiple distinct mass-loss episodes (e.g., Ivanova et al. 2013a;Pastorello et al. 2019a) are not required to generate multiple light curve peaks.
• The luminosity and duration of the early and late light curve peaks follow from simple analytic estimates (Eqs. 11,12 and Eqs. 16,18, respectively). These formulae reproduce within a factor of a few the average luminosity L 90 and duration t 90 (over which 90% of the total optical energy is received) across a wide parameter space, corresponding to mergers spanning the main sequence and postmain sequence phases of stars of mass ∼ 0.1 M to ∼ 100 M (Figs. 10,11). A one-zone calculation, similar in methodology to Popov (1993) but assuming gas pressure-dominated conditions and accounting for recombination energy, reproduces similar scalings (Appendix A).
The total radiated energy E rad is dominated by the recombination-powered plateau; hence the measured value of E rad is proportional to the ejecta mass (up to an adiabatic loss factor f ad ∼ 0.1; Fig. 9).
• We determine the ejecta properties {M ej ,v E , R 0 } of individual LRN events, first under the assumption that hydrogen recombination is the sole source of ejecta heating after launch, and compare them to the masses and radii of observed stellar progenitors of LRNe (Figs. 11, 12, Table 1). The derived ejecta properties differ from those derived by applying SN II light curve models such as Popov (1993) (e.g., Figs. 13, 14). The latter also come with larger uncertainties due to the greater sensitivity of the formulae to the observables (Eqs. 26,25). Remarkably, our results roughly agree with those by Popov (1993) under the assumption that the ejecta velocity equals the stellar escape velocity (Figs. 13). However, such information is not typically available in the Popov framework, and instead using the observed spectroscopic velocity results in drastically different ejecta property estimates compared to those obtained applying our model (Fig. 14).
• In broad terms, for the range of expected ejecta masses and velocities of stars at different stages of stellar evolution, we predict merger light curve durations ∼ weeks to months and luminosities ∼ 10 38 − 10 40 erg s −1 , consistent with observed LRNe, as well as a general trend of increasing luminosity with increasing stellar mass (Fig. 11).
• However, several LRNe arising from mergers involving the most massive stars 10 M , exceed the luminosity upper limit corresponding to ejection of the entire donor star envelope (M ej /M = η > 1; circles in Figs. 11 and 12). The inferred mean ejecta speedv E of many events also exceed those inferred by transient spectroscopy by a factor of ∼ 3.
Both the luminosity and velocity tension can be alleviated by postulating the existence of an additional energy source for the ejecta exceeding that from hydrogen recombination (squares in Fig. 12). This reduces the required ejecta mass to explain the radiated energy as well as the ejecta velocity needed to explain the light curve duration t pl ∝ M 1/3 ejv −1 E . • A promising source for the required "extra" heating is shock-interaction with circumbinary material generated prior to the dynamical phase of the merger (e.g., . The strong velocity dependence of shock-generated heating ∝ v 2 and its associated plateau luminosity (Eq. 31), coupled with the observation that the ejecta velocities of LRNe increase with stellar mass (Fig. 2), would together act to steepen the L − M relationship relative to the fiducial recombination-only prediction (Fig. 11), bringing it closer into alignment with LRN observations which find L ∝ M 2−3 (Kochanek et al. 2014;Fig. 11).
LRN donor star progenitors with larger masses appear systematically more evolved off the main sequence and generate merger ejecta with velocities exceeding the surface escape speed (Fig. 2); this could arise due to the better developed coreenvelope structure of evolved stars, which enables the secondary to inspiral deeper into the gravitational potential well.
• Although the robust conclusions we are able to draw are limited by the modest statistics afforded by the current sample of stellar merger transients, this population is expected to grow considerably in the years ahead. Howitt et al. (2020) predict that the Vera Rubin Observatory will discover ∼ 500 LRNe per year in nearby galaxies, for which pre-transient imaging will enable progenitor star constraints, allowing the methodology presented here to be greatly expanded. • Our model can also be generalized to other merger events involving the ejection of hydrogen-rich material, such as planet-star mergers and starcompact object mergers. However, the relatively modest luminosities achieved by hydrogen recombination power alone, coupled with the comparatively low rates of these events, will likely render these transient signals challenging to detect.
The time evolution of temperature is determined by the first law of thermodynamics: where E, P , V , and L are the internal energy, pressure, volume, and luminosity, respectively. Since the gas-pressure dominates the radiation pressure, the internal energy and pressure are given by E ∼ 3M ej k B T /2m p and P = (γ 3 − 1)E/V with the adiabatic index γ 3 = 5/3 as long as the recombination is negligible. The radiative luminosity is given by the diffusion approximation: where ρ = M ej /V is the density. Combining Eq. (A3) with (A4), we have where A = 32π 2 acm p /(27k B κM 2 ej v). Assuming the opacity is constant, we integrate this equation to obtain where κ = 0.32 κ 0.32 cm 2 g −1 . For a typical gas-pressure-dominated ejecta, ξ 0 1 and the temperature and luminosity decrease as T ∝ t −2 and L ∝ t −4 , respectively.
When the temperature drops to a value T i ∼ 10 4 K, recombination begins in the outermost layers of the ejecta. As time proceeds, the recombination layer recedes back into the ejecta shell. We denote the location of the recombination front by a dimensionless radial coordinate x ≡ R i /R and assume that the gas temperature is fixed to T i at the layer. Considering now only the volume inside the recombination front 4π(xR) 3 /3, the diffusion luminosity (Eq. A4) becomes L 4π(xR) 2 c 3κρ d(aT 4 ) dR ∼ 4π(xR) 2 acT 4 i 3κρxR .
Then Eq. (A3) gives the time evolution of the recombination front x: where we define λ ≡ (γ 3 − 1)/γ 3 . Assuming λ and κ are constant during the recombination phase, we can integrate this equation analytically to obtain x = 1 − 2(1 − λ)ξ i 9(2λ + 5) where T i = 10 4 T i,4 K and the ejecta radius at the beginning of the recombination is given by Eq. (A6), During the recombination phase, the light curve shows a plateau similar to that predicted by Popov (1993), with its duration given by setting x = 0: t pl = R i v 9(2λ + 5) 2(1 − λ)ξ i In the last equality, we have set λ = 0 because the temperature does not evolve significantly during the plateau. A characteristic luminosity may be given by the peak of the plateau, which is obtained by setting dL/dt = 0, L pl = 16π 2 acT 4 i R 4 i 9κM ej 2λ + 5 13 1/2 9(4 − λ)(2λ + 5) 13 ( The total radiated energy is then given by For λ = 0 (γ 3 = 1) and T i,4 ∼ 1 we see that the plateau energy matches that from hydrogen recombination (Eq. 17) and the parameter dependencies for L pl ∝ M