Dusty, Self-obscured Transients from Stellar Coalescence

We discuss the central role that dust condensation plays in shaping the observational appearance of outflows from coalescing binary systems. As binaries begin to coalesce, they shock-heat and expel material into their surroundings. Depending on the properties of the merging system, this material can expand to the point where molecules and dust form, dramatically increasing the gas opacity. We use the existing population of Luminous Red Novae (LRNe) to constrain the thermodynamics of these ejecta, then apply our findings to the progressive obscuration of merging systems in the lead in to their coalescence. Compact progenitor stars near the main sequence or in the Hertzsprung gap along with massive progenitor stars have sufficiently hot circumstellar material to remain unobscured by dust. By contrast, more extended, low-mass giants should become completely optically obscured by dust formation in the circumbinary environment. We predict that 30--50\% of stellar coalescence transients for solar-mass stars will be dusty, infrared-luminous sources. Of these, the optical transients may selectively trace complete merger outcomes while the infrared transients trace common envelope ejection outcomes.


INTRODUCTION
Common envelope phases, with their possible outcomes of stellar mergers or envelope ejection are transformational episodes in multiple star lifetimes. These phases modify the mass, luminosity, and spectral distributions of stellar populations (De Marco & Izzard 2017). They are crucial in shaping the evolutionary history of many massive stars, which are especially likely to have close companions (Sana et al. 2012;Duchêne & Kraus 2013;de Mink et al. 2014). And finally, stellar mergers and common envelope ejections are thought to be crucial in producing exotic stellar outcomes, like rapidly rotating and highly magnetized stars (e.g. Soker & Tylenda 2007;de Mink et al. 2013;Schneider et al. 2016), or compact binaries that go on to merge under the influence of gravitational radiation (e.g. Dominik et al. 2012;Ivanova et al. 2013;Postnov & Yungelson 2014;Vigna-Gómez et al. 2020;Marchant et al. 2021).
Corresponding author: Morgan MacLeod morgan.macleod@cfa.harvard.edu * NASA Einstein Fellow The growing number of transient surveys across the electromagnetic and gravitational wave spectra are offering a new channel for insight into these events. Much of traditional astronomical inference about these objects has come from comparing populations of binary and multiple stars before and after presumed interactions (e.g. Paczynski 1976;van den Heuvel 1976). With timedomain surveys, the possibility of observing these events as galactic or extragalactic transients has become feasible. Catching these events in action is poised to allow us to make direct connections between stellar binaries (and multiples), the transients they produce as they coalesce, and the outcomes of their interactions. This paper attempts a step toward those inferences on the basis of the population of Luminous Red Novae (LRNe). LRNe have been associated with stellar mergers, largely through comparison to the particularly clear case of V1309 Sco -an eclipsing binary system with a decaying orbit that went into outburst and emerged as a single, thermally expanded star (Mason et al. 2010;Tylenda et al. 2011;McCollum et al. 2014;Kamiński et al. 2015;Tylenda & Kamiński 2016;Ferreira et al. 2019). There is now a small sample of LRNe (∼ 10 systems) with known progenitor properties along with detailed observations of the outbursts themselves (for a recent tabulation, see Matsumoto & Metzger 2022, and Section 2).
As the ejecta of LRNe expands, they cool, producing progressively redder emission. As a natural consequence of this cooling, these ejecta eventually form molecules and dust (see Kamiński et al. 2018Kamiński et al. , 2021, for detailed analysis of the galactic sources V4332 Sgr, V1309 Sco, and V838 Mon). Dust formation episodes mark each of the LRNe, and these sources become highly infrared luminous in their immediate aftermath as dusty shells enshroud their ongoing evolution (e.g. Nicholls et al. 2013;McCollum et al. 2014;Chesneau et al. 2014;Tylenda & Kamiński 2016;Blagorodnova et al. 2020;Melis 2020). Further, the contribution of dust formed in common envelope ejecta is thought to be similar to that formed by single evolved stars Zhu et al. 2013;Wang et al. 2015). While dust formation appears to be ubiquitous, the shocks that heat LRNe ejecta depend on the particular system's properties. The characteristic condensation temperature (∼ 10 3 K) for dust grains largely does not. We therefore expect differences to emerge in the efficiency of dust formation and the degree of optical obscuration that results.
In this paper, we draw on scalings from the observed LRNe to constrain the uncertain thermodynamics of LRNe ejecta (Section 2). We model the phases of preoutburst mass loss into the circumbinary environment (Pejcha 2014;Pejcha et al. 2016aPejcha et al. ,b, 2017MacLeod et al. 2018a,b;MacLeod & Loeb 2020a,b), that are believed to precede outburst from stellar coalescence (Section 3). Our models predict the circumbinary mass and temperature distribution, and allow us to convert this to order-of-magnitude estimates for the degree of dust obscuration, which varies widely depending on the merging system's properties (Section 4). The qualitative picture that emerges is one of a bifurcation of transients between the optical and infrared, depending on the properties of progenitor. In Section 5, we conclude.

DUST FORMATION IN MERGER EJECTA
First, we draw on constraints from the subset of the population of LRNe with progenitor detections to understand the cooling of ejecta and dust condensation in these events.
The galactic transients OGLE-2002-BLG360 (Tylenda et al. 2013) and CK Vul (1670) (Shara et al. 1985), are both possibly associated with a stellar merger origin despite having light curves that are somewhat distinct from the other LRNe (e.g. as discussed by Tylenda et al. 2013). Unfortunately, the physical properties of the progenitors of these outbursts remain uncertain. OGLE-2002-BLG360 has a number of interesting properties in the context of dusty transients, which we discuss in Section 4.5.
It is important to highlight that each of these progenitor identifications suffers from the uncertainty of model dependence. In particular, mass transfer could potentially induce significant departures in stellar appearance relative to an isolated-star model (Blagorodnova et al. 2021). For our analysis, we adopt the tabulated properties of systems in Matsumoto & Metzger (2022).

Molecule and Dust Condensation in LRNe
The light curves of LRNe exhibit a progressive reddening after peak luminosity, which traces the cooling emission of their ejecta. This evolution was first traced in exquisite detail during in the 2002 outburst of V838 Mon (Munari et al. 2002;Kimeswenger et al. 2002;Bond et al. 2003;Tylenda 2005). The source was quickly attributed to a possible stellar coalescence origin (Soker & Tylenda 2003;Tylenda & Soker 2006;Soker & Tylenda 2006). Early spectroscopy of V838 Mon showed Hα emission that was soon accompanied by broadened P-Cygni absorption lines tracing the ejecta velocity (Munari et al. 2002). Three years after the outburst, high-resolution spectra Kamiński et al. (2009) show molecular absorption outside the remnant's photosphere, still broadened by outflow into P-Cygni profiles. Early near-infrared spectra traced the emergence of CO and other molecules (Banerjee & Ashok 2002), and an infrared excess in the object's spectral energy distribution was also noted (e.g. Tylenda 2005). Molecular CO emission in the submillimeter bands has since been a valuable tracer of the ejecta (Kamiński et al. 2007), spatially resolved approximately 5 years after the outburst (Kamiński 2008). It has now been observed with Atacama Large Millimeter Array (ALMA), where the interaction of the ejecta with a tertiary component in the system is seen in continuum dust emission and a variety of carbon and sulfur oxide lines (Kamiński et al. 2021). Further recent constraints on the long-term composition and properties of the dusty ejecta have been derived from Spitzer space telescope and SOFIA observations (Woodward et al. 2021).
Subsequent transients have also been observed to undergo epochs of significant dust formation, simultaneous with the dimming of the optical light curve. For example, Nicholls et al. (2013) and Tylenda & Kamiński (2016) showed the significance of an evolving spectral energy distribution in V1309 Sco, with a significant shift toward infrared dominance as the optical light faded (see, for example, Figure 7 of Tylenda & Kamiński 2016). Meanwhile V1309 Sco's optical and sub-millimeter spectra traced similar evolution of molecular formation and absorption seen in V838 Mon and V4332 Sgr, strengthening the connection between these transients (Kamiński et al. , 2018. V1309 Sco shows evidence for warm dust (∼ 10 3 K) even prior to the outburst, indicating that dust formation was ongoing, but at a lower level than in the aftermath of the merger (Tylenda & Kamiński 2016). As a decade-older twin of V838 Mon and V1309 Sco, V4332 Sgr's remnant has proved valuable in understanding the temporal evolution of these sources (e.g. Banerjee et al. 2003;Banerjee & Ashok 2004;Tylenda et al. 2005a;Kamiński et al. 2010;Tylenda et al. 2015;Kamiński et al. 2018). The more recentlydiscovered extragalactic LRNe all share similar cooling and molecule and dust-forming properties, which are described in detail in the references in Section 2.1 as central evidence for their association with a stellar-coalescence origin.
Theoretically, there have been some efforts to model dust formation in common envelope ejecta. Models of dust formation, usually adapted from the asymptotic giant branch (AGB) star context, have been applied to ejecta kinematics and thermodynamics thought to mimic those of merger and common envelope ejecta Zhu et al. 2013). These models estimate dust masses as a function of radius under different model parameters, and also predict the population contribution of common envelope ejecta to the total dust content of the interstellar medium Wang et al. 2015). Glanz & Perets (2018) have pointed out that the adiabatic cooling of common envelope ejecta makes dust formation inevitable, and that the increase in opacity associated with this transition may allow radiation pressure to play an important role in driving mass ejection, as in dust driven winds from pulsating AGB stars. Iaconi et al. (2019Iaconi et al. ( , 2020 perform a related, but more detailed analysis of the thermodynamic trajectories of simulated common envelope phases, and the implications for dust condensation in a dynamic model. As hinted by the V1309 Sco data (Tylenda & Kamiński 2016), they find that warm dust condenses both in early, slow outflow from the phase of Roche lobe overflow preceding the merger, and in the later common envelope ejecta, with slightly different properties of grain size distribution (Iaconi et al. 2020).
Dust formation appears to be a ubiquitous and inevitable consequence of the cooling of expanding common envelope ejecta. The composition and quantity of this dust remain areas of active observational study, and are sensitive to different parameterizations of the uncertain ejecta thermodynamics (e.g. Lü et al. 2013). In the next section we focus on the condensation of warm dust at a characteristic temperature of ∼ 10 3 K as a means of tracing the constraining the uncertain thermodynamics of LRNe ejecta.

Constraints on Ejecta Thermodynamics
The thermodynamics of cooling LRNe ejecta is uncertain. It is dictated by the equation of state and ionization state transitions of hydrogen and helium in the relevant temperature range, as well as by internal shock heating and radiative cooling. Thus, it is not clear that the ejecta temperatures should follow tracks of adiabatic cooling as they expand. By comparing the LRNe population to general power-law temperature profiles, T (r) ∝ r −β , we can constrain the normalization and slope of the temperature structure of these ejecta.
The model power-law temperature as a function of radius is where R 0 is the base of the outflow. We parameterize the temperature at the base of the outflow as where f vir = T 0 /T vir is the fraction of the Virial temperature, T vir , of the primary star at radius R 0 , where we will adopt M 0 → M * and R 0 → R * . Thus, normalizing to example parameter values. Here we note that f vir is thought to be related to the temperature to which material is shock-heated before it is expelled, but it could equivalently represent the original temperature of unshocked ejecta. This temperature profile implies that material will condense to form dust grains at a radius of the order of where T dust ∼ 10 3 K is the temperature at which significant quantities of dust form. This transition induces a major opacity increase, that changes the spectral energy distribution of the LRN from optical to infrared dominated.
During the outburst phase of a stellar coalescence transient, the luminosity and temperature of of the coalescing binary increase dramatically. The source of this increased output is the dissipation of orbital energy into the surrounding gaseous envelope. One possible effect of increasing the central object's luminosity is modifying the thermodynamic properties of the circumbinary material. In thermodynamic equilibrium, the dust formation will occur at Comparing r dust as determined by the internal shock temperature to r dust,eq determined by the equilibrium thermodynamics allows us to determine the relative importance of the increased outburst luminosity in inhibiting dust formation. We apply these considerations to the observed LRNe in Figure 1. Here we compute for each transient, where t 90 is the timescale encompassing 90% of the optical emission and v obs is the velocity inferred from doppler-broadened emission lines (L 90 is the mean luminosity over t 90 , Matsumoto & Metzger 2022). This represents the approximate radius that ejecta reach during the optical duration of the light curve. After this time, LRNe tend to form dust and become extremely infrared luminous (Matsumoto & Metzger 2022). We therefore suggest that, to order of magnitude, r 90 is the radial location where dust formation occurs, modifying the LRN spectral energy distribution. We see that for each LRN transient, r 90 r dust,eq ,  Figure 1. The characteristic radius of dust formation in LRNe as measured by r90 and as predicted by various models. The thermal equilibrium dust formation radii are much larger than r90 for all of the observed systems, which are better modeled by Virialized ejecta with fvir ∼ 10 −2 and β ∼ 1. The constant Virial fraction implied by β ∼ 1 is highly suggestive of internal shocks shaping the ejecta thermodynamics.
adopting L = L 90 . This indicates that the thermal equilibrium temperature is not established in the ejecta, and does not inhibit dust formation at radii smaller than r dust,eq . Next, we compare models of r dust based on the temperature resulting from internal shocks. We find that β ∼ 1 and f vir ∼ 10 −2 provide a good description of the data with M 0 = M * and R 0 = R * .
This indicates that ejecta are only weakly thermalized at the ejection radius, but that their temperature profile decays more gradually than the equation of state alone would indicate, most likely due to internal shock heating (e.g. Pejcha et al. 2016a;. A temperature profile of β = 1 corresponds to gas specific energy proportional to the gravitational potential, in which case the Virial fraction is constant with radius, as would be expected for an internal-shock origin for the temperature profile when the ejection velocity is similar to the local escape velocity. Compared to models of common envelope ejecta dust formation, this result is in tension with previous assumptions. Lü et al. (2013) assume shallower temperature profiles of β = 0.2-0.4, along with a shallower density profile ρ ∝ r −3/2 , noting that different assumptions about the temperature profile in that range could induce seven orders of magnitude difference in the total modeled dust mass (this is likely, in part, due to the divergent density profile that is chosen -if most of the dust is at large radii, the temperatures at those radii become very important). Iaconi et al. (2019) and Iaconi et al. (2020) adopt a very different kinematic and thermodynamic structure motivated by late-stage adiabatic expansion of unbound ejecta in their common envelope simulation. They find an expanding shell of ejecta decreases in density as ρ ∝ r −3 , and at sufficiently late times, cools adiabatically, such that T ∝ r −2 . However, their models also show evidence of the specific entropy increases at earlier times that trace shock heating. Taken together, these results could suggest that internal shocks are important over some temporal and radial regime as the ejecta expand and cool toward the dust condensation temperature but that they asymptotically approach homologous, adiabatic expansion at larger radii.

PRE-OUTBURST MASS LOSS MODEL
In this section, we describe a simplified model for the mass loss that precedes binary coalescence that is motivated by constraints from LRNe outbursts and by modeling of these phases.

Coupled Mass Loss and Orbital Decay
We define a model binary system consisting of an evolving primary star of mass M * and radius R * , which is interacting with a more compact companion of mass M a = qM * , where q is the binary mass ratio. To model the runaway orbital decay associated with dynamically unstable mass transfer, we use the RLOF python package (MacLeod 2020). RLOF adopts numerical coefficients based on hydrodynamic simulations of runaway binary coalescence (MacLeod & Loeb 2020a,b) and solves the ordinary differential equations of binary orbital evolution in the point-mass binary limit. We initialize our model system when the primary star is slightly overflowing its Roche lobe, a = 0.999a RL , where a is the orbital semi-major axis and a RL is the semi-major axis of Roche lobe overflow as estimated by the Eggleton (1983) approximation. We integrate the coupled equations of non-conservative mass loss from a binary system and orbital angular momentum conservation, which reduce tȯ where M = M * + M a , and is a dimensionless specific angular momentum of material leaving the binary. It is the ratio of the specific angular momentum, l loss , of lost material to the specific angular momentum of the binary, l bin = M * M a /M 2 √ GM a. In the expression forṀ * , α is a coefficient of order unity, the orbital period is R L is the the Roche lobe radius (Eggleton 1983), and n = 1 + 1/Γ is the polytropic index of the primary star. Equation (9) is from Paczyński & Sienkiewicz (1972), and has been numerically verified in hydrodynamic simulations by (MacLeod et al. 2018a;MacLeod & Loeb 2020a), and represents adiabatic mass loss from a polytropic donor star. We caution that in cases where the thermal adjustment of the mass-shedding layers of the donor star are important, the approximation of equation (9) will be inaccurate.
In what follows, to provide context for the qualitative trends across many binaries, we make some uniform assumptions for all of the binary systems we model. We assume all primary stars have n = 3/2 polytropic index that enters into equation (9), and we further assume that the adiabatic response of the donor star is to maintain constant radius upon mass loss. In reality, the structure of donor stars depends on their evolutionary state, with most objects ranging between effective polytropic indices of 4/3 ≤ n ≤ 5/3; the choice of n = 3/2 represents an intermediate value. The response of systems to mass loss is more complex because, in addition to the adiabatic response of a donor star (Hjellming & Webbink 1987), thermal adjust modifies the mass-radius relation dependent upon the mass transfer rate itself. More tailored parameters can be easily specified with RLOF, or similar models that account for the thermal and nuclear evolution of the primary star could be completed with the binary module of MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019, for example as recently computed in a similar context by Blagorodnova et al. (2021).

Circumbinary Material
The kinematics and thermodynamics of material in the circumbinary environment are complex and multidimensional. The continuous orbital motion launches spiral shocks through outflowing material, which provide heat and redistribute angular momentum. Early, slower mass loss trails away from the binary concentrated toward the system's orbital plane and may be either bound (Pejcha et al. 2016b;MacLeod et al. 2018b;MacLeod & Loeb 2020b), or weakly unbound (Pejcha et al. 2016aHubová & Pejcha 2019). Later, faster ejecta en-counter this torus, and are reshaped into bipolar outflows (MacLeod et al. 2018b).

Torus Model
If the pre-coalescence mass loss is bound to the binary system, as found in the hydrodynamic simulations of MacLeod et al. (2018b); MacLeod & Loeb (2020b), it forms a thick torus surrounding the binary. Within this torus, gas is shock heated by the continual spiral outflow from the binary. MacLeod & Loeb (2020b) found that these shocks also redistribute angular momentum, such that the torus has nearly constant specific angular momentum determined by the total mass loss and total angular momentum change of the decaying orbit. An analytic solution for a hydrostatic toroidal configuration with constant specific angular momentum can be derived if we assume a polytropic equation of state P = Kρ γ . Then, where R and z are cylindrical coordinates relative to the central binary mass, M . In this parameterization, l is the torus specific angular momentum and R t is the outer radius of the torus (see appendix B of MacLeod & Loeb (2020b) for the derivation of this expression).
On the basis of hydrodynamic simulation results, MacLeod & Loeb (2020b) provide fiducial torus parameters for this torus model as it is implemented in RLOF. We set l = ∆L/δm where ∆L is the cumulative binary angular momentum change, and δm is the mass change. This implies that all mass and angular momentum lost from the binary go to forming the torus. We set where q 0 is the initial binary mass ratio. This implies that the extent of the circumbinary material scales with the primary star radius. Finally, the specification of the torus polytropic constant, K, and index, γ, closes the model. These are linked to the temperature structure via an ideal gas equation of state T = (P/ρ)(µm p /k B ), where we adopt µ ≈ 1. We choose γ = 4/3, then the polytropic constant K is calculated for self-consistency between R t and the integral torus mass. The choice of γ = 4/3 is motivated by two factors. First, this is the best fit to MacLeod & Loeb (2020b)'s hydrodynamic simulation models, where internal shocks (rather than the gas' adiabatic index) are the critical factor in establishing this entropy profile. Second, this implies T (R, 0) ∝ R −1 in the bulk of the torus, just as noted in Section 2 for the Virialized ejecta from LRNe. We note that this temperature probably represents the post-shock temperature in the wake of passing internal shocks, T sh , and can be modified by radiative cooling and diffusion in the intervals between the passage of shocks (e.g. Pejcha et al. 2016a).
The Virial fraction to which circumbinary material is shock-heated is proportional to the square of the ratio of the velocity dispersion of internal shocks over to the characteristic velocity at a given radius, f vir ∼ ∆v 2 /(GM * /r) (e.g. Pejcha et al. 2016a;MacLeod & Loeb 2020b). The normalization of our models shows that f vir is a function of binary mass ratio, with f vir ∼ 0.25 for q = 1/3, and f vir ∼ 0.15 for q = 0.01 (MacLeod & Loeb 2020b). We note that this is significantly higher f vir than found for the outburst phase in Section 2; this is likely because the slow, early outflow self-intersects much more dramatically as it decelerates in the gravitation potential than does the faster, more homologous late ejecta which are being modeled in Figure 1.

Wind Model
If, rather than being bound to the binary, the circumbinary material is unbound and expands outward with an asymptotic velocity v w , we can adopt a different model to describe the circumbinary distribution. We imagine that the most likely scenario is a relatively slow, equatorially-concentrated outflow, in which internal shocks again play a central role in the thermodynamic structure (Pejcha et al. 2016a).
In this model, we specify the wind velocity to be a fraction of the stellar escape velocity, where, motivated by Pejcha et al. (2016a), we choose f v = 0.25 to represent a slow outflow solution. We use the time-evolving mass loss rateṀ * (t) from the RLOF solution, along with the substitution where ∆t = t merge −t is the time before binary merger, in our case, when a = R * . Each shell of the circumbinary distribution therefore corresponds to a certain time of emission from the binary. Then whereṀ * (r) is related toṀ * (t) by equation (15) and the sign comes from the fact thatṀ * < 0. We note that the time dependence ofṀ (t) implies a density slope steeper than r −2 surrounding the binary. Under our fiducial assumptions, we find scaling similar to ρ ∝ r −3.5 .
We use the Viral factor to set the post-shock temperature structure of the outflow, and adopt the same r −1 slope, which is motivated by internal shock heating and is equivalent to the γ = 4/3 index of the torus model. As a fiducial normalization, we adopt f vir = 0.25 for q = 1/3 (MacLeod & Loeb 2020b).

Dust Properties and Extinction
We adopt a highly simplified dust model, with the goal of estimating the order of magnitude significance of dust extinction in circumbinary material. We assume that dust forms when the post-shock temperatures are below a condensation temperature, T sh T c = 10 3 K. Below T c , we assume a constant mass fraction of material forms grains, X d = 5 × 10 −3 , typical of grain condensation in solar-composition material. We adopt a fiducial V-band opacity for these grains per gram of dust of κ d = 10 3 cm 2 g −1 . Dust opacities depend relatively sensitively on the grain size limits, distribution, and composition. Our choice is motivated by the comparatively small grains expected to initially coalesce in the ejecta (Iaconi et al. 2020).
In reality, the creation and destruction of grains takes time, and a time-dependent dust condensation, growth, and sputtering model could be used to develop much more sophisticated models of the ejecta-dust interaction (e.g. Höfner et al. 2016Höfner et al. , 2022. For now, we note that the dust-to-gas fraction and the grain size distribution will be affected by the density and temperature of the outflow (e.g. Lü et al. 2013;Glanz & Perets 2018), and that revisiting these schematic properties in the future will offer significant insight into the dust forming properties of these coalescing binary systems. Such an exercise has been recently undertaken in the limit of asymptotic expansion of simulated ejecta by Iaconi et al. (2020).
In either of our models, we integrate the total optical depth due to dust extinction along a line of sight (in the equatorial plane for the torus model), where ρ is a function of r in the wind model, and r and z in the torus model and X d is a function of the local temperature. We then specify A V = 1.086τ d .

PRE-OUTBURST OBSCURATION
In this section, we explore the dependence of dust formation and obscuration on binary system parameters. [mag] torus model wind model Figure 2. Time evolution of orbital period, circumbinary mass loss rate, and optical obscuration, in an example merging binary system. This system consists of an M * = M , R * = 30R primary star unstably transferring mass to a M /3 companion. We integrate this example until the separation equals R * and the pair merges. In both the wind and torus models of pre-outburst obscuration, the degree of obscuration increases as the mass loss intensifies and the system's merger approaches.

Time evolution
We begin, in Figure 2, by showing the temporal evolution of an example system from unstable Roche lobe overflow until the orbit shrinks to the primary star's radius, a = R * . We imagine this as a point of transition from the "lead-in" to the engulfment of the secondary object that initiates the common envelope phase itself. The system in Figure 2 consists of an M * = M , R * = 30R primary star transferring mass toward a M a = M /3 companion (q = 1/3). The system is otherwise modeled by our fiducial parameters described in Section 3.
As the binary orbital period decreases, the degree of Roche lobe overflow of the primary star increases, driv-ing the rate of mass loss from the system to increase also, as described by equation (9) and shown in the middle panel of Figure 2. This increasing mass loss rate accelerates the rate of orbital decay, as described by equation (8). The coupled mass loss and orbital decay lead to the eventual deposition of ∆m ∼ 0.25M a into the circumbinary environment (MacLeod & Loeb 2020b).
The lower panel of Figure 2 shows the time-dependent obscuration of the merging binary due to this circumbinary material. We show results from both the torus and wind models for the circumbinary distribution, which are described in Sections 3.2.1 and 3.2.2, respectively. In the torus model, we adopt a viewing angle through the torus midplane (which is equivalent to the equatorial plane of the binary system). In both models, the dust obscuration, A V , increases from very low values long before merger to higher values at later times, reflecting the increasingly mass-rich circumbinary environment of the binary. The wind model exhibits higher obscuration than the torus model, especially at early times. This difference represents the greater radial extent of the unbound wind relative to the bound torus. Within the final 10 2 d before the binary merges, the wind model levels out as newly-lost material does not have time to propagate outward to the dust condensation radius. Within the hydrostatic torus model, the circumbinary material is assumed to instantaneously assume its hydrostatic configuration. The timescale for this adjustment is not known in detail, but it is likely at least the stellar dynamical time. Therefore, the predicted obscuration in the last orbit or perhaps several orbits might not be realistically modeled in the torus case because this delay for hydrostatic adjustment is not explicitly accounted for. In Figure 2, this uncertainty would correspond to the last ∼ 30 d, for example. The immediate pre-merger behavior thus represents a factor ∼ 2 uncertainty in the torus model. Despite these differences, it is interesting to note that both models reach a similar estimate of several magnitudes of dust extinction at the time of merger. Viewed over time under this increasing obscuration, we might see an optical precursor of the outburst to fade while its infrared counterpart brightens. This evolution points to the leverage that detailed pre-outburst observations of stellar-coalescence transients can provide on these sorts of models of circumbinary mass loss Blagorodnova et al. 2017Blagorodnova et al. , 2020Blagorodnova et al. , 2021Pastorello et al. 2021a).

Dependence on binary properties
In Figure 3, we map the maximal A V measured at the time of merger (a = R * ) in our wind and torus [mag] (torus model) Figure 3. Predicted obscuration just prior to merger in the wind and torus (equatorial line of sight) models given primary star mass and radius. Each model assumes q = 1/3. The progenitor properties of known LRNe are marked with diamonds. In both models, low mass, extended radius primary stars should lead to highly dust-obscured transients, in some cases reaching AV 10. By contrast, more compact or massive primary stars should lead to low levels of preoutburst obscuration. The observed population of optical LRNe traces predicted AV 1 in either model. models. Our results show that this maximal absorption depends on both stellar mass and radius. At the simplest level this result can be traced to the characteristic temperatures of circumbinary outflows. As material is shock-heated and expelled, the post-shock temperatures depend on the compactness of the primary star. More massive, smaller radius primary stars represent deeper gravitational wells and will have hotter, faster moving circumbinary material, that is, by consequence, less dust-forming. The cool circumbinary surroundings of low mass, extended radius stars become dramatically obscuring, perhaps providing tens of magnitudes of extinction in the most extreme cases. The qualitative trends and normalization of the wind and torus model predictions are similar. However, Figure 3 does highlight that these models exhibit somewhat different parameter dependence. We empirically derive the following relations based on models in which we randomly varied model parameters. For the wind model, while for the torus model with an equatorial line of sight, (20) In this context, it is worth highlighting that our simple assumption of unchanging dust physics across the entire parameter space neglects the likely dependence of dust formation on the density and velocity of ejecta. In cases of very high predicted extinction, e.g. A V ∼ 10 2 , larger grains might instead form, impacting the opacity as a function of wavelength and the infrared appearance of obscured transients. Equation (20) only parameterizes the obscuration along the equatorial line of sight. We find that within the torus model there is viewing angle dependence: at polar angles, cos(i) 0.4, the sight line to the central binary is unobscured while at more equatorial angles, cos(i) 0.4, A V is approximately constant.
Primary star mass and radius thus clearly affect the degree of obscuration due to circumbinary mass loss. These primary star properties also correlate with different stellar evolutionary states and appearances. Figure  4 traces circumbinary obscuration in the wind and torus models in the Hertzsprung-Russell Diagram (HRD). We show evolutionary tracks of 1, 2, 4, 8, and 16M stars at solar metallicity from MIST (Choi et al. 2016). Here, we can immediately observe that merger events on the main sequence will be largely unobscured, with A V 1 mag. Similarly, mergers of Hertzsprung gap objects, with T eff 5500 K, also experience low degrees of obscuration.
This can be contrasted to the late evolutionary phases of low-mass stars. In the tracks of the 1, 2, and 4M  This figure adopts wind and torus models with q = 1/3. The primary star properties are realized from MIST stellar tracks of solar metalicity. The lower-mass extended radii primary stars that lead to the highest predicted AV occur on the giant branch of primary stars with M * 4M . These systems have low T eff 4000 K, and relatively high luminosities L 10 2 L . By contrast, the progenitors of optical LRNe have a wide range in mass, but a narrow range in effective temperature, primarily being mildly evolved yellow giants crossing the Hertzsprung gap, properties that correspond to low levels of pre-outburst obscuration.
stars, very high degrees of dust formation and obscuration A V 5 mag are observed near the tips of their respective giant branches. This location in the HRD corresponds with low effective temperatures, T eff 4000 K, and high luminosities L 10 2 L . Thus, before dust obscures these merging systems, they are preferentially old, luminous, and cool giants. Figure 5 shows these same stellar tracks in the context of a common-envelope HRD, where we show the binding energy of primary stars along with their lumi-  Figure 4. Low-binding-energy primary stars are those that are more likely to lead to successful common envelope ejection in interactions with a companion object. We find that low binding energy correlates with predictions of highly obscured mergers, while higher binding energies have low predicted AV . The population of optical LRNe traces the lowobscuration, high binding energy phase space of merger outcomes.
nosity. Common envelope ejection is believed to be related to the overall system energy balance, in which to unbind a common envelope ∆E ∼ E bind , where E bind ∼ GM 2 * /R * is the binding energy of the primary star's envelope before the coalescence (e.g. Webbink 1984). The energy source, ∆E, is largely the gravitational energy of the decaying companion star orbit, ∆E ∼ E final ∼ GM 2 * q/a final . At the crudest level, therefore, envelope ejection requires a final /R * q. In more detailed analysis, coefficients can change these scalings by an order of magnitude (e.g. Webbink 1984;Iben & Livio 1993;De Marco et al. 2011;Ivanova et al. 2013).
In Figure 5, we see that systems with lower binding energies that will tend to lead to common envelope ejection outcomes are predicted to be highly obscured. Higher binding energy systems -those closer to their main sequence radii -tend to lead to low predicted obscuration and are more likely to correspond to merger outcomes.

The progenitors of Luminous Red Novae
We compare the population of optical transient progenitors (Section 2) to our model predictions for degree of obscuration. In Figure 3, we plot the LRNe progenitors on the primary-star mass-radius plane, and in Figure 4, we plot these progenitors on the HRD. Both the mass-radius and HRD views show that LRNe progenitors lie in phase space where we expect relatively low degrees of pre-outburst dust extinction. In the torus model, systems typically exhibit A V 0.1 mag, while in the wind model A V 1 mag. In either case, these predictions are consistent with a largely unobscured transient. We note that an unobscured transient is always possible for a pole-on line of sight in the torus model, but this being a viewing-angle effect is unlikely given that V1309 Sco is known to be viewed in the equatorial plane.
In Figure 5, we see that the lower binding energy primary stars are systematically not observed in the present LRNe sample. This is highly suggestive that the unobscured, optical LRNe sample traces merger outcomes, while dusty, obscured transients will characterize common envelope ejections.
There is some evidence in precursor lightcurves of LRNe for some degree of time-variable extinction. V1309 Sco and M31 LRN 2015 each exhibit a slow rise of precursor brightening followed by a optical fade (∼ 1 mag) just before the optical outburst Blagorodnova et al. 2020). If these dips were due to time-dependent dust reddening rather than a bolometric fade that would be consistent with the preoutburst mass loss model. Indeed, the presence of infrared excess in the pre-outburst phase of V1309 Sco is suggestive that this is the case (Tylenda & Kamiński 2016). For V1309 Sco (M * ∼ 1.5M , R * ∼ 3.5R ), the predicted A V is 1 mag in the wind model and 0.1 mag in the torus model. For M31 LRN 2015 (M * ∼ 4M , R * ∼ 30R ), the model predictions are very similar, suggesting that the mild fade observed pre-outburst is entirely consistent with dust reddening.

Dusty transients in the binary population
We have argued that LRNe originate from the relatively unobscured population of binary merger events. Their more obscured counterparts will be dusty and optically faint, while perhaps being infrared luminous. We . Relative population of merger transients as a function of primary star mass. In this model, we assume all mergers have q = 1/3. We compute the predicted AV along MIST solar metalicity tracks, and show the fractions of merging systems with low, intermediate, and high AV , assuming a logarithmic merging separation distribution. Among approximately solar mass stars, on the order of half of transients will be highly dust-obscured, primarily-infrared transients. As we consider higher mass systems, a greater fraction will be primarily-optical transients.
imagine that these sources could have extremely red colors of the SPRITEs class of exclusively-infrared transients reported by Jencson et al. (2019) from the Spitzer Infrared Intensive Transients Survey (SPIRITS, Kasliwal et al. 2017).
To understand the imprint of dusty transients on the population of binary merger and common-envelope transients, we make a simple estimate of their fractional contribution to the overall merging population. To do so, we assume that binaries are distributed with a broad separation distribution with equal number of binary systems per logarithmic bin in separation (Öpik 1924). That implies that as stellar radii grow, the probability of merging with a companion per log radius is constant. In Figure  6, we trace MIST evolutionary tracks from the main sequence out the giant branch, assuming an equal number of mergers per log radius of the star. We apply equations (19) and (20) to estimate A V in these events. We assume isotropic viewing angles. In the spherical wind case of equation (19) this has no effect. However, for the torus, we assume 40% of viewing angles are sufficiently pole-on to have A V = 0, otherwise we apply equation (20). Finally, we divide these events into groups of low, intermediate, and high levels of obscuration. The most obscured systems will be those where the optical signature is nearly completely extinguished, and will manifest as infrared transients. The intermediate levels of reddening will lead to highly reddened optical and infrared transients, while the low extinction case will lead to primarily-optical transients. Figure 6 shows that at 1M , about 30-50% of mergers lead to infrared transients. This fraction decreases to nearly zero above 8M . We can understand this high fraction of infrared transients from low-mass systems in context of the observed events in two ways. First, we consider a volume-limited sample of stellar coalescence events, such as all the events occurring within the galaxy. These galactic events will be greatly dominated (by number) by low-mass events, with masses similar to 1M , because of the steeply-decreasing stellar initial mass function (IMF) dN/dM * ∝ M −2.3 * . We note that this is true despite the increasing fraction of massive stars in interacting binaries, which partially counteracts the IMF. The fraction of interacting systems is about 15% for solar-mass primary stars, but increases to nearly 100% for O-type stars, with the interacting binary fraction scaling roughly as f interact ∝ ), in the wind model, 25% of events have A V < 1, 48% have 1 < A V < 5, and 26% have A V > 5. For the same weighting in the torus model, 69% of events have A V < 1, 9% have 1 < A V < 5, and 22% have A V > 5.
Given the nature of these sources as astronomical transients, we might also consider a luminosity limited sample. Kochanek et al. (2014) has shown that em-pirically, L LRN ∝ M 3 * . This scaling more than cancels the slope of the IMF and interacting binary fraction, implying a similar number of detected LRNe across all masses, with a possible weighting toward more massive stars. While the number of sources remains small, this is at least partially born out by the broad mass distribution of the observed sample. In this luminosityweighted sample, optical events will dominate because of their higher luminosities and the greater contribution from massive stars.
Finally, it is important to consider that dustenshrouded, infrared transients are associated with different primary stars than optical counterparts. These progenitor stars tend to be cool, luminous, extended companions with lower binding energies (e.g. equations (19) and (20) and Figures 4, and 5). The correlation of these properties with possible common envelope ejection outcomes suggests that common envelope ejection transients will systematically be among the obscured, infrared population. The higher binding energy envelopes associated with merger outcomes will lead to less obscuration and optical transients.
The transient OGLE-2002-BLG360, originally thought to be a lensing event, has been interpreted as a possible LRN transient (Tylenda et al. 2013). The system gradually brightened into a nearly thousand-day outburst in 2002. The progenitor was a T eff ∼ 4300 K giant. The distance was uncertain, but adopting a distance of 8.2 kpc, Tylenda et al. (2013) derive R * ∼ 30R and L ∼ 300L . The old galactic bulge stellar population implies a low progenitor mass, perhaps M * ∼ M . Intriguingly, Tylenda et al. (2013) find that the preoutburst spectral energy distribution requires an extinction of A V ∼ 3 by warm ∼ 800 K dust, which they note is only consistent with ongoing dust formation in a continuous outflow. As the source evolved, it eventually faded in the optical as it became so dust-obscured that A V 20 is required to accommodate the constraints. Our models predict A V ∼3-10 for a 1-2M , approximately 30R primary star. The long duration of the transient is likely related to the long orbital period about the extended primary, For comparison, the expansion time of material to the dust condensation radius is, if we assume that r dust is given by equation (5) and v ej ∼ (2GM * /R * ) 1/2 . The similarity of these timescales is indicative that the ejecta expansion and cooling is ongoing even during the outburst phase of the event. Thus, the entire transient episode is similar to the wind-like pre-outburst behavior where the temporal evolution is dictated by the changing mass outflow. This is different from the other LRNe, in which the duration of the event is dictated by the expansion time of the ejecta, because t dust P orb . This is indicative that dusty transients tend to be more wind-like as opposed to impulsive, and that the duration of the events will be set by the duration of mass ejection rather than material expansion.
We conclude that OGLE-2002-BLG360 is very likely a prototype event for an intermediate level of dust obscuration, positioned in the phase space of reddened transients in Figure 6. We anticipate that more detailed modeling of the abundant data for this source will be very informative.

SUMMARY AND CONCLUSIONS
In this paper, we have analyzed the thermodynamic evolution of the ejecta of stellar coalescence transients, and of the mass lost leading up to these episodes. When dust condenses in this circumstellar material, it dramatically increases the opacity obscures optical light from the merging system while redistributing emission toward the infrared. We find that the degree of this obscuration depends on the properties of the merging system. Some key results of our study are: 1. After being shock-heated to a few percent the Virial temperature, the ejecta in LRNe appear to decline linearly with radius, proportional to the gravitational potential. This implies continuous heat injection by internal shocks are crucial in regulating the ejecta thermodynamics ( Figure 1).
2. Pre-outburst mass loss accompanies orbital decay in merging binary systems. The resulting circumstellar material can lead to time-dependent obscuration or even the optical vanishing of coalescing systems due to dust formation ( Figure 2).
3. The most obscured systems are those with lower primary-star masses and extended radii high on the giant branch, while more compact or massive systems tend to be comparatively unobscured (Figures 3 and 4).
4. Because lower envelope binding energy is correlated with higher A V , the systems most likely to lead to common envelope ejection will tend to systematically produce dusty, infrared transients ( Figure 5).
5. The population of optical LRNe are all associated with model predictions of low obscuration (Figures 3, 4, and 5). A corresponding population of dusty, optically-obscured, infrared-luminous transients arise from the ∼25% of stellar coalescence events involving lower-mass extended donor stars ( Figure 6).
6. OGLE-2002-BLG360, with properties slightly distinct from many of the optical LRNe, may be a prototype of a partially dust obscured transient intermediate between the optical and fully infrared events.
The fact that dusty, infrared transients should exist, and that their properties are correlated with the properties of the donor star in the binary system, indicates the breadth of the population of transients related to stellar coalescence. The optically-luminous LRNe appear to be just one subset of this full class of events, with properties that suggest they tend to be associated with stellar merger (rather than common envelope ejection) outcomes. More broadly exploring this parameter space, both in modeling and observationally, is a promising path toward a deeper understanding of these transformative events in binary stellar evolution. With time-domain capabilities expanding in both the optical and, more recently, infrared (e.g. De et al. 2020), we are well positioned to begin mapping the range of possible events. To this end, the anticipated near-infrared follow-up capabilities of the James Webb Space Telescope will likely be pivotal in revealing the properties of these dusty systems.

REPRODUCTION SOFTWARE AND DATA
Software and data to reproduce the results of this work are available via github at: https://github.com/ morganemacleod/dustytransients materials (MacLeod 2022).