A"blue ring nebula"surrounding a thousands of years old stellar merger

Stellar mergers are a brief-lived but common phase in the evolution of binary star systems. Among the many astrophysical implications of these events include the creation of atypical stars (e.g. magnetic stars, blue stragglers, rapid rotators), interpretation of stellar populations, and formation channels of LIGO-detected compact object mergers. Although stellar mergers are thus commonly invoked phenomena, observations of these events and details of their evolution remain elusive. While a handful of stellar mergers have been directly observed in recent years, the central remnants of these events remain shrouded by an opaque shell of dust and molecules, making it impossible to observe their final state (e.g. as a single merged star or a tighter surviving binary). Here we report observations of an unusual, ring-shaped ultraviolet nebula and the star at its center, TYC 2597-735-1. The nebula shows two opposing fronts, suggesting a bipolar outflow from TYC 2597-735-1. TYC 2597-735-1's spectrum and proximity above the Galactic plane suggest it is an old star, yet it shows abnormally low surface gravity and a detectable long-term luminosity decay, uncharacteristic for its evolutionary stage. TYC 2597-735-1 also exhibits H-alpha emission, radial velocity variations, enhanced ultraviolet radiation, and excess infrared emission, common signposts of dusty circumstellar disks, stellar activity, and accretion. The combined observations, paired with stellar evolution models, suggests TYC 2597-735-1 merged with a lower-mass companion several thousand years ago. TYC 2597-735-1 provides a look at an unobstructed stellar merger found at an evolutionary stage between its dynamic onset and the theorized final equilibrium state, directly inferring how two stars merge into a single star.

The events leading up to and following the merger of TYC 2597-735-1 and its companion shape the system we see today. As TYC 2597-735-1 and its companion approached sufficiently closely, the former overflowed its Roche lobe onto the latter, initiating the merger (Figure 3(b)).
Numerical simulations demonstrate that the earliest phase of the merger process results in the ejection of matter through the outer L 2 Lagrange point in the equatorial plane of the binary [28]. Most of this matter remains gravitationally bound around the star system, forming a circumbinary disk.
The companion, unable to accommodate the additional mass, is dragged deeper into the primary's envelope in a runaway process [11] (Figure 3(c)). This delayed dynamical phase is accompanied by the ejection of a shell of gas. A portion of the ejected material is collimated by the circumbinary disk into a bipolar outflow [29]. The balance of the mass lost during primary-companion interactions remains as a circumstellar disk, which spreads out and cools over time, eventually reaching sufficiently low temperatures to form dust.
We see evidence of this relic disk around TYC 2597-735-1 today as an infrared excess. A simple analytic model, which follows the spreading evolution of the gaseous disk due to internal viscosity over the thousands of years since the merger (see SI), is broadly consistent with both the present-day gas accretion rate (estimated from Hα; see SI) and a lower limit on the present-day disk mass obtained by fitting the infrared spectral energy distribution of TYC 2597-735-1 (M disk,dust 5 × 10 −9 M ; see SI). Accretion of disk material onto TYC 2597-735-1 could account for its observed stellar activity (e.g., Hα emission and far-ultraviolet excess). Angular momentum added to the envelope of the star by the merger, and subsequent accretion of the disk, would also increase TYC 2597-735-1's surface rotation velocity. Indeed, we find the de-projected surface rotation velocity of TYC 2597-735-1 to be ≈ 25 km/s (see SI, ED Figure 8), larger than expected for a star which has just evolved off the main sequence (v < 10 km/s) [30].
The bipolar ejecta shell expands away from the stellar merger, cooling and satisfying within weeks the conditions for molecular formation and solid condensation (Figure 3(d)), as seen directly in the observed ejecta of luminous red novae [31]. The mass ejected, inferred from merger simulations and modeling the light curves of luminous red novae, is typically ≈ 0.01-0.1 M [32], consistent with our lower mass limit of the ultraviolet nebula (0.004 M ).
As the nebula expands and sweeps up interstellar gas, a reverse shock crosses through the ejecta shell, heating electrons in its wake. These electrons excite the H 2 formed in the outflow, which fluoresces in the far-ultraviolet (Figure 3(e)). Although dust is almost always observed in the ejecta of stellar mergers [31], we find no evidence of dust in the ultraviolet nebula (e.g., ultraviolet/optical reddening). We speculate that either the dust was destroyed in the reverse shock [33] or that the bipolar ejecta shell has thinned out sufficiently such that dust is currently undetectable.
Consistent with the latter, this ultraviolet nebula marks the oldest observed stellar merger to date, being 3 − 10 times older than the previously oldest stellar merger candidate, CK Vulpeculae (1670), which is still shrouded by dust [34]. The system was caught at an opportune time -old enough to reveal the central remnant, yet young enough that the merger-generated nebula has not dissolved into the interstellar medium. 8 The discovery of an ultraviolet nebula introduces a new way of identifying otherwise hidden late-stage stellar mergers. With 1 -10 of these objects expected to be observable in the Milky Way (see SI), future far-ultraviolet telescopes may uncover more late-stage stellar mergers. TYC 2597-735-1 poses a unique opportunity to study post-merger morphology as the only known merger system not enshrouded by dust. For example, the close separation of the initial stellar binary−which shares properties broadly similar to those of protoplanetary disks−could form "second generation" planets [35]. However, given the relatively short time ( 10 8 years) until TYC 2597-735-1 reaches the end of its nuclear burning life, the potential window of habitability for such planets could be drastically reduced compared to main sequence stars.
Supplementary Information is linked to the online version of the paper at www.nature.com/nature.
Acknowledgements This research is based on observations made with the Galaxy Evolution Explorer, obtained from the MAST data archive at the Space Telescope Science Institute, which is operated by the As-

13
Use the ATLAS9 pre-set grid of synthetic stellar spectra [36] to fit the TYC 2597-735-1 spectral energy distribution to representative stellar spectra. All synthetic stellar spectra are publicly available: https://www. stsci.edu/hst/instrumentation/reference-data-for-calibration-and-tools/astronomical castelli-and-kurucz-atlas. Portions of our analysis used Python packages photutils [37] and astropy [38].  (µm)), shows enhanced near-and mid-infrared emission. Far-infrared measurements (>24 µm, IRAS; solid circles) provide upper limits only. Synthetic stellar models that best match the inferred stellar properties of TYC 2597-735-1 (gray spectra and grayed region) do not account for either the infrared excess and also suggest a far-ultraviolet excess is present. Both the far-ultraviolet and infrared are frequently observed in systems actively accreting matter from warm, gaseous and dusty disks, like T Tauri protostars [12]. Models of warm, dusty circumstellar disks reproduce the observed infrared excess (light blue dotted line: T dust ∼600 Kelvin, ∼0.2 -3 astronomical units, assuming a disk inclination angle ∼15 degrees; dark blue dashed line: T dust ∼1200 -300 Kelvin, ∼0.2 -1.5 astronomical units, assuming a disk inclination angle ∼15 degrees). b. TYC 2597-735-1 exhibits Hα emission, an unusual trait for evolved stars. The Hα line profile shows variability over short timescales. There is an enhanced blue edge to the emission, a classic signature of infalling material and traditionally interpreted as accretion flows or disk winds [14]. The Hα emission, excess infrared emission, enhanced far-ultraviolet radiation, and radial velocity variations (see ED Figure 5) all suggest TYC 2597-735-1 is actively accreting materiel from the disk creating the observed infrared excess emission. . We convolved high-resolution synthetic H 2 fluorescence spectra with the GALEX grism spectral resolution to produce the plots shown. Pumping by Lyα photons (whose source would likely come directly from TYC 2597-735-1) creates peaks in the distribution that are not seen by GALEX near 1450Å. Electron-impact fluorescence produces a spectral distribution that better matches where GALEX sees the FUV emission being produced. Extended Data Figure 3: TYC 2597-735-1 is an outlier when compared with other moderatelyevolved stars of similar mass. A large sample of moderately evolved stars is shown [39] to demonstrate that, in particular, TYC 2597-735-1's effective temperature and surface gravity are not consistent with the majority of other stars following similar evolutionary tracks. If the present-day observable properties of TYC 2597-735-1 are a consequence of a previous stellar merger, as our MESA models suggest, then we expect that TYC 2597-735-1 is currently puffed up more than usual and will continue to relax over the next thousands of year to better match the trend of evolving stars in T eff -log(g) space (see ED Figure 7).

14
Extended Data Figure  In both GALEX FUV and NUV images, starlight from TYC 2597-735-1 and field stars has been removed using a customized source removal routine in the python photutils module.
The routine masks off the point-spread function (PSF) of any star that falls above a signal-tonoise (SN) threshold (for FUV, SN>5, and for NUV, SN>2) and replaces the star with a tile of random background noise representative of the surrounding sky. The sky noise is derived from sky background maps provided with each DIS exposure and co-added to create a master sky map. This exercise was done to remove excess flux from the nebular area, to derive the flux from the diffuse nebula only.
The FUV nebula appears as a ring with a slightly-elliptical shape and a clearly-defined hole 26 at the center. The FUV emission ring of the BRN has an inner spatial extent on the sky d in = 3 .7 and an outer extent of d out = 7 .3. The geometry of the BRN is such that the FUV emission shines through where forwards and backwards outflow cones overlap (as drawn in Figure 1(f)).
The viewing angle of the nebula was derived by modeling a biconical outflow from a point source and matching the geometry of the overlapping cones to reproduce the profile of the BRN. This gives i ∼15 • (Figure 1(f)). Properties of the Ultraviolet Nebula The count rate observed by GALEX is converted to flux in standard cgs units using the GALEX conversion factors for FUV and NUV channels [42]. We Neutral Hydrogen Mass in the BRN The FUV luminosity, assuming the emission is exclusively produced by an optically-thin layer of H 2 , can be converted to a rate of H 2 fluorescence over the entire BRN. This conversion assumes that the nebular gas is very thin, such that the molecular density does not saturate the absorption cross-section for excitation.

Multi-Wavelength
To justify this assumption, we estimate the volumetric density of H 2 (n(H 2 )) in the BRN.
To do so, we estimate the rate that H 2 is collisionally excited by hot electrons in the medium, assuming that the source of H 2 excitation is associated with a reverse shock through the BRN (as is shown in the last panel of Figure 3). The nebula itself (excluding the shock filament) emits a total H 2 luminosity L H 2 = 3×10 33 erg s −1 . We convert this to a total rate of H 2 collisions required to produce this luminosity, cm 2 ), so we calculate the average expected collision rate with electrons over the area of the nebula using derived cross sections of collisions between H 2 and free electrons: is the cross-section between H 2 and warm-hot electrons (E e − > 20 eV [45]). The expected density rate of collisions between electrons and cool-warm (T ∼ 500 K) H 2 is C H 2 +e − ∼ 10 −11 cm 3 s −1 [46]. Setting the rates of collisions between warm H 2 and free electrons and the observed rate of H 2 fluorescence, we find: This density is consistent with the diffuse, electron-impact H 2 fluorescence observed in Mira's FUV tail [16].
With the assumption that the H 2 nebula is quite diffuse, we estimate how much mass of neutral hydrogen has been produced in the visible area of FUV emission exposed by this H 2 fluorescence process; each time an H 2 molecule is excited collisionally by an electron, it has a non-trivial chance (10-15%) to dissociate into neutral hydrogen (HI) [47].  Selective slit spectra along the western boundary of the nebula observe hydrogen-Balmer series emission along the shock boundary. Figure 1( shock, indicating a low Hα/Hβ ratio inconsistent with local thermodynamic equilibrium. Instead, the ratio is more consistent with those observed in solar flares (Hα/Hβ ∼ 1), which are interpreted as being non-LTE and significantly impacted by charge exchange with free electrons [49,50]. The optical and UV emission observed along the shock are therefore consistent with a non-radiative shock front with recombination of hydrogen dominating the emission from this region [51].
Constraining the Age of the Blue Ring Nebula The maximum velocity of the BRN forward shock places constraints on the age of the nebula. Assuming that, at early times after the merger, the ejecta coasts at a constant velocity matching that of the forward shock, v sh , the source age can be estimated as t = r sh /v sh . As the ejecta begins to sweep up an appreciable mass of the external interstellar medium (ISM) comparable to its own, a reverse shock will pass through the ejecta. At late times, the evolution of the forward shock will approach that of Sedov-Taylor blast wave, for which r sh ∝ t 2/5 in the case of a constant density ISM and, therefore, v sh = dr sh /dt = (2/5)r sh /t, implying a source age t = (2/5)r sh /v sh . The age of the source can therefore be constrained between 2 5 which comes out to be 2,000 years < t < 5,000 years, given the inferred size and expansion rate of the BRN. As demonstrated in our MESA models (see below), the closest-matched model to the current properties of the BRN central start, TYC 2597-735-1, is best suited to an age t ∼ 1,000 years since the stellar merger began. Therefore, we adopt an age range for the BRN between 1,000 years < t < 5,000 years throughout the course of this study.
Gaia distance The Gaia DR2 parallax of TYC 2597-735-1 is 0.518±0.029 mas, from which employing a Bayesian method with the preferred "exponentially decreasing space density" (EDSD) prior and a distance scale corresponding to a disk population of 1000 pc at a galactic latitude [52] of +39.4 • , we obtain a distance of 1935 +127 −91 pc. The fractional error on the parallax is sufficiently small that the details of the method and the adopted prior only change the distance by less than one-tenth of the quoted 1σ error bars on the distance. Excess FUV emission (both in emission lines and continuum) are well-known and reported in systems with accretion disks [57][58][59]. GALEX-FUV imaging covers a wavelength range (1350 -1750Å) where studies report excess FUV emission produced in accreting systems accounts for a small percentage (<5%) of the total measured FUV excess [58]. Using the total accretion luminosity derived from the observed Hα emission (see Stellar Hα Emission), the measured FUV excess accounts for 3% of the total accretion luminosity. This is well within reason, given that the majority of TYC 2597-735-1's FUV excess is likely not measured (i.e., shorter wavelengths than were observed with GALEX).
From 2µm out to longer wavelength, the SED exhibits excess infrared emission relative to the single temperature stellar continuum. This IR excess likely arises from a dusty, warm circumstellar disk, similar to those seen around protostars [60] and some post-AGB systems [61]. A disk-like geometry for the dust, which lies in a plane perpendicular to our line of sight (and hence approximately perpendicular to the symmetric axis of the nebula), is further supported by the lack of circumstellar reddening of TYC 2597-735-1, inferred from our Na I D analysis (see Interstellar and Circumstellar Reddening). The IR excess peaks at λ ∼8 µm, corresponding to an average dust characteristic temperature ∼600 K.
If this temperature reflects re-radiated light from thermal dust particles in orbit around TYC 2597-735-1, then the dust peak location is ∼2.2 AU (depending on the albedo of the dust, which is entirely dependent on the dust grain distribution, which ranges between 0.2 -3 AU [62][63][64]). The dust disk area is derived by modeling the dust distribution over a thin layer in the circumstellar disk to reproduce the IR SED [65]. We verify that r in ∼0.25 AU and r out ∼2.8 AU best reproduce the IR excess SED (Figure 2(a)). We empirically derive the dust dist mass from the infrared spectral energy distribution of TYC 2597-735-1, assuming all the dust is roughly the same temperature in a thin layer at the equilibrium distance and that the mean dust size is small [63,66]. We assume an average dust temperature T dust = 600 K and take the equilibrium distance D dust ∼2.2 AU. The surface density of the dust disk is taken to be small (ρ dust ≈2.5 gm/cm 3 [66]) with an average small dust particle size (a dust ∼10 −4 cm, or ∼1 µm) to ensure we are estimating a lower limit to the dust disk mass. Using the following relation [67], we derive a minimum dust mass in the present-day circumstellar disk M dust 5 × 10 −9 M . Hypothetically, if the measured Keck/HIRES RV using the iodine cell cross-correlation is linked to a close-in companion with ∼13.7 day orbital period, we argue that this companion is not sufficiently massive enough to eject the BRN we observe today. A companion with P = 13.7 days would have a semi-major axis a ∼0.1 AU. Assuming TYC 2597-735-1 has a mass 2 M and the orbital plane ranges from perpendicular to the axis of the Blue Ring Nebula (i min ∼ 15 • ) to edgeon to our vantage point (i max ∼ 90 • ), the companion mass will be between M p ≈ 3.1 − 12.1M J (ED Figure 5). The total gravitational energy liberated as the planet inspiraled from ∞ to its final semi-major axis a p is only We performed a bisector velocity span (BVS) analysis on the Keck/HIRES iodine cellcalibrated RV results, which is meant to provide a test as to whether RV signals measured from a star are tied to stellar activity, like star spots or atmospheric pulsations [13,69]. We find a statis- disk as suspected from these diagnostics, then it is reasonable to conjecture that the atmosphere is disturbed or that hot spots near the poles exist [13]. Indeed, the period of the HIRES RV is not consistent with that of the rotational velocity of TYC 2597-735-1 (6.5 km/s, corresponding to P∼20 days), so linking the activity near the stellar poles, where the period of revolution may better match the 13.8 days, seems a plausible explanation.

Keck/HIRES Observations
Moreover, we find that radial velocities from the same HIRES data measured with the telluric line method of Chubak et al. 2012 [70] do not agree with the iodine cell calibrated measurements, as shown in ED Figure 5. There are differences of up to 600 m/s, whereas the telluric method claims to be accurate to within ∼100 m/s. We do not attempt a periodic fit to the telluric derived velocities. The reason behind this significant disagreement is not known. and it raises concerns regarding the reliability of either calibration method.
To provide an ancillary interpretation to the Keck/HIRES RV results, we obtained additional RV data in the near-infrared with the Habitable-zone Planet Finder (HPF) to determine if the iodine 39 based periodicity persists at longer wavelengths.
HET/HPF Observations and RV reduction HPF is a high resolution (R ∼ 55, 000) near-infrared spectrograph [71,72] that is actively temperature stabilized to the milli-Kelvin level [73], HPF has a near-infrared (NIR) laser-frequency comb (LFC) calibrator which has been shown to enable ∼20 cm/s calibration precision in 10-minute bins, and 1.5 m/s RV precision on-sky on the bright and stable M-dwarf Barnard's Star over months [74]. Due to the faintness of the target, we elected to not use the simultaneous LFC reference calibrator for these observations to minimize any possibility of scattered LFC light in the target spectrum [75]. Instead, the drift correction was performed by extrapolating the wavelength solution from other LFC exposures on the nights of the observations [75]. This methodology has been shown to enable precise wavelength calibration and drift correction down to the ∼30 cm/s level, much smaller than the RV variability observed in this system.
To extract the RVs, we use a modified version of the SpEctrum Radial Velocity Analyzer (SER-VAL) code [78], which uses the template-matching technique to extract precise RVs, adapted for use for the HPF spectra. The stellar activity indicators used in this work -including the differential line width indicator (dLW), chromatic index (CRX), and Calcium Infrared Triplet (Ca IRT) indices -are the same as defined in the SERVAL pipeline [78].
We also see evidence of activity-induced RV variations in the HET/HPF RVs, which we interpret as features intrinsic to the star which could include stellar surface inhomogeneities (e.g., We further note that with our current sampling with HET/HPF -which covers approximately half of the phase of the 13.75-day periodicity seen in the Keck/HIRES RVs -the amplitude of the RV variations in HPF likely represents a lower limit on the NIR RV amplitude. Generally, activityinduced RV variations are expected to show lower amplitudes in the NIR than in the optical [79].
We thus speculate that given the 5 year baseline between the HIRES and HPF RVs, that the stellar surface-features could have evolved to give rise to larger activity-induced variations seen during the HPF RV observations, although continued long-term RV observations across the optical and NIR are required to further test this hypothesis.
With the combined evidence in place from two independent RV surveys, we suspect that the RV signals observed by Keck/HIRES and HET/HPF are dominated by stellar surface activity and modulations on TYC 2597-735-1. Even if a close-in binary companion does exist, we argue that it is not sufficiently massive enough to explain the mass ejection that happened at TYC 2597-735-1 >1,000 years ago that is now seen as the BRN.

Stellar Hα Emission
For what appears to be a red sub-giant star (see below), TYC 2597-735-1 displays a highly unusual characteristic: Hα line emission in its stellar spectrum. Hα emission from star systems can be produced by a variety of physical mechanisms, including accretion (e.g., T Tauri stars, Herbig Ae/Be, symbiotic stars), magnetic activity (e.g., low-mass stars), stellar rotation (e.g., active M-dwarfs), and stellar pulsations (e.g., Mira variables pulsating and heating part of their stellar atmosphere) -all systems with heightened levels of stellar surface activity [80][81][82]. In many cases, Hα emission is observed to be variable (e.g., shape, flux) over time (although the timescale of variability is not typically well characterized [82]) and may display both emission and Hα emission in young stars is widely associated with material from a circumstellar disk being funneled onto the star via magnetospheric accretion, or removed from the disk in the form of winds [14]. A similar reservoir of material could be fueling the Hα emission and variability of TYC 2597-735-1. As discussed, TYC 2597-735-1 exhibits an infrared excess that supports the presence of a circumstellar disk orbiting the star, which is an expected consequence in stellar merger simulations [29]. This disk provides the necessary material that would accrete onto TYC 2597-735-1 through magnetospheric accretion, creating signposts of stellar activity. If L(Hα) is a direct measure of the accretion luminosity, L acc [83], then L acc ∼ 10 33 ergs s −1 . Converting L acc to mass accretion rate (Ṁ ), TYC 2597-735-1 accretes 1.5×10 −7 M /yr of material. We note that converting L(Hα) to a mass accretion rate provides an upper limit for the accretion rate, as studies have shown that Hα emission is also generated by disk material itself, disk winds and other sources embedded in disks [14,84,85].
Stellar Properties of TYC 2597-735-1 We present a brief description about the determination of TYC 2597-735-1's stellar properties using Keck/HIRES spectra.
Luminosity and Radius of TYC 2597-735-1 The Johnson V-band magnitude, derived from SDSS gr photometry [86], is 11.149 mag. Using extinction laws [87] and E(B−V)=0.02 mag, the V-band magnitude corrected for extinction is V 0 = 11.087 mag. The absolute V-band magnitude of BRN, based on the Gaia parallax, is then M V = -0.34±0.12 mag. From bolometric corrections [88] for the Sun and TYC 2597-735-1, -0.193 and -0.219 respectively, we derive log(L/L ) = 2.07 for TYC 2597-735-1 [89]. Paired with the adopted T ef f = 5850 K and log(L/L ) = 2.07, we find R = 10.5 ± 0.5 R for TYC 2597-735-1. We find this de-projected surface rotation velocity to be larger than expected for a star which has just evolved off the main sequence (v < 10 km/s) [30,90,91]. However, the fact that TYC 2597-735-1 is not rotating at a much larger fraction of its break-up speed disfavors a more equalmass merger (e.g., two 1M stars).

Rotational Velocity & Macroturbulence
Effective Temperature Two spectroscopic and one photometric methods to determine the effective temperature, T ef f , of TYC 2597-735-1 were explored. The second technique employs line-by-line differential abundances [93,94]. Abundances for each Fe I line are derived relative to the same line in a standard star whose atmosphere parameters, including T ef f and [Fe/H], are accurately known. As in the absolute analysis, the model that best estimates T ef f for the star provides a zero slope solution in a differential abundance versus excitation plot, relative to the standard. Ultimately, this is relative to the Sun, whose T ef f is accu-rately known and not dependent on model atmospheres. By using this technique, the adopted gf values for each line cancel-out when the difference with the standard is computed. Therefore, any measurable clean line in both target and standard star spectra can be used, even if no gf values are known. In this differential analysis we used 42 lines from the visual Keck/HIRES CCD spectrum of TYC 2597-735-1, 16 of which are in common with the absolute analysis (above). We took line abundance differences relative to star Hip 66815 [94], giving a best fit T ef f of 5,900K, and a microturbulence velocity ξ = 1.3 km/s. We adopt a T ef f uncertainty of 51K, equal to the quadrature sum of the excitation temperature uncertainties for TYC 2597-735-1 and our standard. We speculate that the 0.2 km/s microturbulent velocity difference between absolute and differential methods is likely due to errors in the chosen standard. Combining result from all three methods, we converge on T ef f = 5,850 ± 50 K.

Stellar Properties and Abundance Analysis
The model atmosphere abundance analysis follows the methods described in various studies [93,[96][97][98]. The analysis uses the MOOG spectrum synthesis routine [99] and a grid of model atmospheres [100].

MESA Models
We use the stellar evolution code MESA [24] to explore the impact of a stellar merger on the long-term stellar properties displayed by TYC 2597-735-1 [25]. We assume a 2.17M primary star and varied the mass of the merging secondary (M c ) from those of Jupitermass planets up to low mass stars M c ≈ 0.3M . We explore the merger evolution for the upper limit primary stellar mass inferred from our stellar properties analysis (see Mass of TYC 2597-735-1), but we expect that, if we fix the stellar radius and structure of the primary and decreased the mass by a factor of 2 (to match the lower limit of the primary mass), expect to increase the companion mass by a factor of 2 (to match the energy). The lower mass primary star and higher mass companion, still safely within the low-mass star mass range, would yield roughly the same timescales of merger evolution. We expect that the luminosity of the lower-mass primary case would change in a non-trivial way, resulting in quantitative, but likely not qualitative, changes.
Future modeling efforts of the lower mass primary scenario is necessary to confirm this finding.
We consider different evolutionary stages of the primary at the time of the merger, ranging from it just leaving the main sequence, to an evolved sub-giant star on the horizontal branch.
We simulate different evolutionary stages of the primary star by changing the radius of the star, ranging from 3 -10 R . For each simulation, right after the companion star plunges into the primary, we deposit energy into the envelope following the dynamical friction-driven inspiral. We follow the star's subsequent evolution over a long timescale, ∼10,000 years. We then search for the combination of pre-merger parameters that best reproduce the present-day properties of TYC 2597-735-1, namely its luminosity, effective temperature, and surface gravity. We ED Figure 7 shows the evolution of the MESA models following the deposition of energy during the merger and assumes a 0.1 M companion. Different colored lines correspond to different sizes of the 2.17 M primary star at the time of merger, ranging from 3 − 10 R . We see that the moderately evolved sub-giant star with a radius of 5R does the best job at reproducing the present-day properties of TYC 2597-735-1 for an assumed time since merger of 1,000 years. This age is a factor of ∼ 2 smaller than that estimated for the BRN of 2,000 years.
We note that both the primary star's effective temperature and surface gravity 1,000 years after the merger energy injection adequately match the present-day derived values measured from optical spectroscopy and photometry of TYC 2597-735-1, which explains why TYC 2597-735-1's stellar properties are slightly skewed away from the bulk of moderately-evolved stars in the T ef f -logg plane [39] (ED Figure 3). As demonstrated in ED Figure 7, the primary stellar properties continue to settle back towards equilibrium after the modeled merger takes place, its surface gravity increases, shifting TYC 2597-735-1's T ef f -logg relationship in-line with other moderately-evolved stars around its effective temperature.
Our simple one-dimensional models neglect a variety of effects, which if improved upon in future work would enable a more precise comparison to data. These include, for example, the back reaction of the addition of mass (particularly unprocessed hydrogen) on the long-term stellar structure, as well as multi-D effects due to rotational mixing and the delayed accretion of mass from the remnant accretion disk. More detailed models for angular momentum transport in the remnant will also be required to make specific predictions for the present-day rotation rate of TYC  [102,103], brightening and fading of four "Hot RCB stars" [104], and the slow-then-fast fading of the 'Stingray' planetary nebula nucleus [105].
The photometric accuracy must be around 0.01 mag, for binned magnitudes, and both the old and modern magnitudes must be very carefully placed onto the identical magnitude scale. Our procedure is to perform differential photometry of TYC 2597-735-1 with respect to the average of three nearby comparison stars of closely similar color and magnitude. We chose TYC 2597-1026-1, TYC 2597-458-1, and TYC 2588-182-1, and adopted B magnitudes from the AAVSO Photometric All-Sky Survey (APASS) of 11.506, 11.513, and 11.912 mag respectively. By using the same three comparison stars with similar color and only detectors with a sensitivity close to the Johnson B system, all color terms and systematic errors will be negligibly small, despite the many plates and detectors over the last century. By averaging together many images from many nights, we can beat down measurement and systematic errors to usefully small values. images, and the formal uncertainty taken as the RMS scatter divided by the square root of the number of input images. The RMS scatter between the magnitudes for the latter four observational runs is 0.008 mag, which is substantially larger than the formal error bars. TYC 2597-735-1 is likely not varying on any fast time scales, so 0.008 mag is a measure of how accurately photometry can be compared between observers with different systems, despite having the best possible conditions and procedures. These observer-to-observer differences can be beaten down by averaging over the four observers, to get a modern 2015-2019 B magnitude of 11.767±0.004.
Our final light curve is presented in Supplemental Table 1 and ED Figure 6. We see a steady decline, from 1897.5 to 1952.5, consistent with a simple linear function with a slope of 0.127±0.016 mag/century (with a reduced-χ 2 near unity). This decline is significant at the 6-σ level, as taken from an F-Test. A linear fading of a light curve in magnitudes is the same as an exponential decline in the flux. Between 1952.5 and 2019.9, the slope has flattened out and is consistent with a zero slope. The two high points in the 1970s might represent some variability during the 1952.5-2019.9 interval, or maybe some residual errors of some sort, or perhaps a shallow decline from 1940 to 2019.9. In all cases, TYC 2597-735-1 starts with a fast decline in the first half of the 1900s, and a slow-or-zero decline from 1950 to now. A formal χ 2 fit to a bent line gives a break around 1940, with slopes of 0.131±0.014 and 0.040±0.006 mag/century before and after the break. An F-test shows the break to be significant at the 3.2-σ level. The total B-mag decay is found to be between 0.11-0.12 mag from 1895-2020, consistent with ∼0.09 -0.1 mag/century. short orbital separations 0.1 − 1 AU to interact with the star due to tidal inspiral following its post-main sequence evolution [108]. Also assume that the current state of the system is observable for a time comparable to its present age, e.g. t obs ∼ 2t age ∼ 10 4 yr; this is motivated in part by the fact that in older sources the shock-heated electrons potentially responsible for exciting the H 2 would have adiabatically cooled by further expansion of the ejecta. Given the lifetime of a 2M star of t ∼ 1 Gyr, we then expect that only a fraction

Rates of BRN Formation
of A-stars to contain a BRN. The number density of stars in the stellar neighborhood of mass 2M is n ∼ 5 × 10 −4 pc −3 [109]. Therefore, the total number of stars within a radius R = D = 1.93 kpc and scale-height H ∼ 1.5 kpc is approximately N = n πD 2 H ∼ 10 7 . Thus, the number of candidate A-stars as close as TYC 2597-735-1 that we would expect to show a BRN is i.e. of order unity.
A second independent constraint on the rate comes from the Galactic rate of stellar merger events, as inferred from observations of luminous red nova transients, which Kochanek et al.
2014 [110] estimate to be R ∼ 0.3 − 0.5 yr −1 . We would therefore expect a total numberÑ BRN ∼ Rt obs ∼ 3, 000 Galactic A-stars in a BRN-bearing present state similar to TYC 2597-735-1 out of the total numberÑ ∼ 10 9 of A-stars in the Milky Way (which we take to be a fraction ∼ 1% of the total number 4 × 10 11 of stars). If every merger produced a BRN, then the expected BRN fraction would therefore bef which is consistent (within an order of magnitude) with the equally-uncertain estimate (4). The number of candidate A-stars from the rate of Galactic stellar mergers is then Therefore, we expect from our two independent methods to find anywhere between 0.4 -10 BRNs.
Our finding of one BRN is consistent with this rate.

Analytic Processes Describing Mergers
The process of the companion falling into the stellar envelope causes the ejection of mass in a rotationally supported disk (or decretion disk) near or above the surface of the star. Mass ejected from the L 2 point during the earliest stages of a stellar merger can remain gravitationally bound [28,111], particularly for binary mass ratios q ≡ M c /M 0.06, as would be satisfied in our case for M c 0.1M . The equatorial outflow acts to shape matter which is ejected at higher velocities during the final, plunge phase of the merger, into a bipolar outflow [29,112]. Our MESA calculations demonstrate that the merger process could eject a mass M ej ∼ 0.01M from the stellar envelope.
After forming a disk of radius ∼ R ,0 and mass (of, say) M d,0 ∼ M c ∼ 0. − 1M , the disk will accrete onto the star at a characteristic rate [113] M pk ∼ M d,0 t visc ∼ 6 × 10 25 g s −1 α 0.1 ∼ days, where α is the Shakura-Sunyaev viscosity parameter. We have assumed a thick disk with a scaleheight H ∼ R d /2, which is justified because the peak accretion luminosity L acc,pk = GM Ṁ pk R ∼ 10 41 erg s −1 is several orders of magnitude larger than the Eddington luminosity of the star ∼ 10 38 erg/s.
At such highly super-Eddington accretion rates, the disk cannot cool on the accretion timescale and is therefore subject to powerful outflows [114]. These outflows may carry away an order unity fraction of the initial disk mass in the form of a wind, i.e. ∼ 0.01M . This wind would also be funneled by the disk geometry into bipolar conical geometry, perhaps contributing to the BRN ejecta on top of the material ejected dynamically during the merger itself.
What happens to the disk that remains bound to the star? For a thick disk with constant H/R d that evolves due to the viscous redistribution of angular momentum, the accretion rate onto the star evolves with time approximately as [115,116] M (t) =Ṁ pk t t visc −4/3 , t t visc (11) and therefore the disk will evolve to become sub-Eddington (L acc 2×10 38 erg/s) after a timescale t Edd ∼ 100t visc ∼ 1 year The current accretion rate would be estimated aṡ M (t age ) =Ṁ pk t age t visc −4/3 ∼ 10 19 g/s, comparable to that inferred from the present-day Hα luminosity (see Methods).
As the disk accretes, it will also viscously spread outwards in radius due to the redistribution of angular momentum. If the evolution of the disk conserves total angular momentum By the time the disk has become sub-Eddington (t ∼ t Edd ) its outer edge would reach a radius R d (t Edd ) ∼ 1 AU. The disk may spread a bit further than this after t Edd , but since the disk will become geometrically thin (H/R d 1) after becoming sub-Eddington the viscous timescale over which the spreading occurs t visc ∝ (H/R d ) −2 (eq. 9) will increase significantly. Thus, we expect R d (t age ) ∼ few 10 13 cm. The equilibrium temperature at the outer disk edge i.e. sufficiently low to allow dust formation and consistent with the inferred temperatures of the NIR excess of TYC 2597-735-1 (see Spectral Energy Distribution).
The current mass of the disk in such a scenario is approximately M d,0 (t Edd /t visc ) −1/3 ∼ 0.1M d,0 ∼ 10 −3 M , which, assuming the standard 1 : 100 value for the dust-gas mass ratio, would imply a present dust mass of ∼ 10 −5 M . This is consistent with the minimum mass of dust only 5 × 10 −9 M we derive from model-fitting the IR excess [117] (see Spectral Energy Distribution).