An infrared transient from a star engulfing a planet

Planets with short orbital periods (roughly under 10 days) are common around stars like the Sun1,2. Stars expand as they evolve and thus we expect their close planetary companions to be engulfed, possibly powering luminous mass ejections from the host star3-5. However, this phase has never been directly observed. Here we report observations of ZTF SLRN-2020, a short-lived optical outburst in the Galactic disk accompanied by bright and long-lived infrared emission. The resulting light curve and spectra share striking similarities with those of red novae6,7-a class of eruptions now confirmed8 to arise from mergers of binary stars. Its exceptionally low optical luminosity (approximately 1035 erg s-1) and radiated energy (approximately 6.5 × 1041 erg) point to the engulfment of a planet of fewer than roughly ten Jupiter masses by its Sun-like host star. We estimate the Galactic rate of such subluminous red novae to be roughly between 0.1 and several per year. Future Galactic plane surveys should routinely identify these, showing the demographics of planetary engulfment and the ultimate fate of planets in the inner Solar System.

Given the crowded nature of the field as well as marginal evidence for blended sources near the transient location in the WIRC images, we obtained additional high spatial resolution Adaptive Optics (AO) assisted imaging using the Gemini South Adaptive Optics Imager (GSAOI) on the Gemini-S telescope 86,87 . The source was observed on UT 2022-04-15 as part of a Director's Discretionary Program (GS-2022A-DD-102; PI: K. De). We obtained dithered exposures of the field using Laser Guide Star (LGS) correction for a total exposure time of 300 s on source and 300 s off the source (for background subtraction) in J-band. Observations were also obtained in the H and K filters with GSAOI but were unusable due to poor atmospheric conditions. The raw images were detrended and stacked using the DRAGONS pipeline 88 Table 1. The source location was also covered in optical images acquired in the PanSTARRS1 survey 73 , but the progenitor was not detected. We obtained 5σ upper limits on the progenitor optical flux by querying the PS1 source catalog in a 2 ′ region around the position and estimating the median magnitude of sources at the 5σ detection threshold. The flux limits are reported in Extended Data Table 1.
7 Spectroscopic follow-up Following the identification of the transient, we obtained one epoch of optical spectroscopy and one epoch of near-infrared spectroscopy of the outburst using the Low Resolution Imaging Spectrometer (LRIS 95 ) on the Keck-I telescope and TripleSpec 96 on the Palomar 200-inch telescope respectively. The LRIS data were reduced using standard techniques using an automated pipeline 97 , while the TripleSpec data were reduced using the spextool pipeline 98 followed by telluric correction and flux calibration using xtellcor 99 . The spectra obtained during the outburst are summarized in Extended Data Table 2. We obtained two epochs of near-infrared spectroscopic follow-up of the infrared remnant star ≈ 2 years after the outburst using the Folded-port Infrared Echellette (FIRE 100 ) on the Magellan Baade telescope and one epoch using the Near-Infrared Echellette Spectrometer (NIRES 101 ) on the Keck-II telescope. The observing setups are summarized in Extended Data Table 2. The FIRE data were reduced using the pypeit code 102 for the echelle mode data and the firehose pipeline 103 for the prism mode data.
The NIRES data were reduced as the TripleSpec data using the spextool pipeline 98 , followed by telluric correction using the xtellcor pipeline 99 . As the late-time data were acquired close in time (Extended Data Table 2), we performed an inverse variance weighted stacking of the spectra to improve the signal-to-noise ratio. The individual spectra were flux calibrated to match the contemporaneous late-time JHK photometry obtained from WIRC, and the stacked spectrum was binned to the resolution of the Keck-NIRES data (R ≈ 2700). The final reduced and stacked spectra are shown in Figure 3. the most extended (A-array) configuration on 2022 April 14 under project 22A-464 (PI: De). The observations, of which about 53 minutes on ZTF SLRN-2020, were set up to be sensitive to any SiO (maser) emission between 42.25 and 43.5 GHz. The SiO transitions were targeted to possibly help confirm the peculiar giant merger product and to possibly measure its line-of-sight velocity to get a constraint on a (kinematic) distance to support its remarkable low luminosity.
The observations used 3C286 to calibrate the flux density scale and the instrumental bandpass response. Antenna pointing calibration in the direction of the target field was performed using J1851+0035, whereas the gain and reference calibrator J1912+0518 observed 50 cycles of 25 seconds while interspersing the 65 second observations of ZTF SLRN-2020 using a cycle time of 2 minutes (with a slew time of about 10 seconds each way). The data were inspected for bad visibilities and further processed using the general PIPEAIPS pipeline reduction procedure in AIPS. Calibration and imaging did not reveal any potential issues and the resulting synthesized beam measured about 65 by 45 mas at a position angle of −60 degrees. The calibrator J1912+0518 was measured to have a flux density of about 88 ± 2 mJy/beam.
The setup was in particular targeting the (J=0-1) v=1 and v=2 transitions using a channel width of 62.5 kHz (≈ 0.5 km/s) to an RMS level of ≈ 6 mJy/channel, although all seven possibly detectable transitions were covered using 1 MHz (≈ 7 km/s) channels with an RMS of ≈ 1.2 mJy/channel. The ∼1.2 GHz line-free continuum bandwidth at 43 GHz was imaged down to an RMS of ≈ 35µJy/beam. No significant SiO line nor continuum emission was detected in the several arcseconds surrounding the infrared position of ZTF SLRN-2020, down to the 3.0σ flux uncertainty within 0.2" of the expected position. For the line emission, we searched up to ≈ 60 MHz (±420 km/s) of the transition rest frequencies which is a conservative upper limit for objects bound to the Galaxy. The SMA data were reduced using the Millimeter Interferometer Reduction (MIR) software using the online procedure 104 . The resulting calibrated visibility data were imaged using the Common Astronomy Software Application (CASA)'s tclean task. We imaged a range of channels covering a velocity interval of ≈ ±60 km s −1 around the CO 2-1 line, producing a data cube of channel maps. The synthesized beam size in the combined observations was 1.25 × 0.85 ′′  broadly exhibits a fast rise to peak starting ≈ 10 days prior to outburst peak followed by a slow decay of ≈ 2.5 mag over ≈ 200 days. While the Galactic objects V1309 Sco and OGLE-BLG-360 exhibit a long and slow rise prior to the peak of the light curve, the photometric behavior of ZTF SLRN-2020 around peak is similar to V838 Mon. The fast rise and slow decay of ZTF SLRN-2020 is also similar to many extragalactic red novae (e.g. AT 2018hso and AT 2017jfs) but the duration is longer than the fast decay of AT 2019zhd after peak. ZTF SLRN-2020 also exhibits a rapid decay around ≈ 30 days after light curve peak, similar to AT 2019zhd but with a smaller drop in magnitude. Additionally, the pre-outburst i-band variability of ZTF SLRN-2020 is similar to that of V838 Mon as well as the nearby red novae in M31 (M31-LRN-2015 and AT 2019zhd) where the progenitor variability was detected prior to the main outburst.
In Extended Data Figure 2, we also compare the color evolution of ZTF SLRN-2020 to objects in the literature. ZTF SLRN-2020 exhibits a rapid reddening around the peak of the outburst followed by a slow reddening of the g − r color, while the r − i color remains relatively constant.
Similar reddening of the g − r color is also seen in the literature red novae, although the magnitude of the progressive reddening is larger for these objects during the ≈ 100 days after peak. The slow g − r color evolution of ZTF SLRN-2020 is most similar to AT 2018bwo, which also exhibits only a small change in the color over the first ≈ 40 days, while the slow r − i color evolution is most similar to AT 2017jfs. Accounting for the estimated foreground Galactic extinction (see Section 12), ZTF SLRN-2020 remains relatively blue compared to the literature sample. As the progressive reddening with time is attributed to the cooling of an expanding photosphere in red novae followed by dust formation 59 , the relatively slow evolution in ZTF SLRN-2020 may be related to a lower amount of photospheric expansion and dust formation relative to other objects, consistent with the modeling of the bolometric light curve and SED (see Section 12).

Spectroscopic features
The optical outburst spectrum of ZTF SLRN-2020 ( Figure 3) exhibits a relatively featureless continuum with only atomic/molecular absorption features of Na, Ba II, Hα, Mg II, TiO and VO. The TiO molecular bands are characteristic in the spectra of cold, late-type giant stars, with later type giants exhibiting deeper absorption features 106,107 . To identify the corresponding spectral type, we compared the spectrum of ZTF SLRN-2020 to a library of stellar spectra from the VLT X-shooter spectral library 108 , and found a good match of the TiO/VO line depths to the M4-III type giant HV 2255. We applied a foreground extinction of A V = 5.5 mag (consistent with the inferred total dust optical depth modeled in Supplementary Information, Section 12) to the spectrum of the comparison star using a standard Fitzpatrick extinction law 50 and show the comparison in Figure 3. The inferred spectral type corresponds to an effective photospheric temperature of ≈ 3600 K 109 .
The optical spectrum shows a strong absorption line in the Na D doublet. We measure the total equivalent width (EW ) of the line by fitting a polynomial to the absorption feature, and measure the uncertainty by creating 1000 realizations of the spectrum by adding Gaussian noise scaled to the root-mean-square noise in the adjacent part of the spectrum. We find EW = 3.65 ± 0.75Å. The strength and profile of the Na D line is known to vary with time in many types of explosive transients 110,111 , and has recently been shown to be time-varying in extragalactic red novae 112 . The line becomes stronger with time in red novae, likely due to the condensation of dust in the envelope. As a result, we are unable to use the Na D feature as an indicator of the foreground extinction as for other types of Galactic transients 60, 113 .
The spectra of ZTF SLRN-2020 show no signatures of emission lines indicative of hot gas in the eruption, as commonly seen in other types of Galactic plane transients. In Extended Data Similarly, the outbursts of young stars exhibit strong emission lines at all phases indicative of hot gas in the outflowing material. The outbursts of the FU Ori stars 115 are known to be years -decades long, while their outburst spectra show spectral types of F/G-type stars together with wind-like profiles in the H and Ca NIR triplet lines 115 , inconsistent with the properties of this source. While the pre-outburst mid-IR brightening in ZTF SLRN-2020 is similar to that seen in the young star Gaia 17bpi 62 , outbursts of FU-Ori stars are known to be consistently blue in color at all phases due to the accretion of hot gas, while ZTF SLRN-2020 shows a striking evolution from blue colors near the optical peak to a red SED dominated by infrared emission ≈ 1 year later ( Figure 2).
When compared to young stars in terms of photometric evolution, ZTF SLRN-2020 shows some similarities to that of the EXor class of young stars that exhibit few month long outbursts 117 .
However, the optical spectra of EXor outbursts are dominated by a forest of emission lines of atomic species like H and Ca II at all phases 118 , indicative of magnetospheric disk emission 119 , unlike that seen in the optical spectra of ZTF SLRN-2020 (Extended Data Figure 3). While red novae also exhibit emission lines at early phases 53, 59, 112, 120 , the lines become progressively weaker with time as the photosphere becomes dominated by molecular absorption features from the newly formed molecules and dust 55 , similar to that seen in the spectrum of ZTF SLRN-2020 around 6 months after outburst peak ( Figure 3). The position of ZTF SLRN-2020 is not coincident with any known star forming regions; a search of the WISE archival images show evidence for faint diffuse emission > 10 ′ away, arguing against an association with young stars which are exclusively found in star forming regions. Together, we find the case for ZTF SLRN-2020 to be associated with a young star outburst to be unlikely. Comparing the spectra to stars in the IRTF Spectral Library 123 , the broad H 2 O absorption bands and weak molecular absorption features are reasonably matched to a M7-III type star, and we show the comparisons in Figure 3 and Extended Data Figure 4. We note that although the NIR spectrum was obtained near the optical spectrum, the NIR spectrum suggests a later spectral type than the optical spectrum. Such differences have also been identified in the late-time spectra of previous Galactic red novae 18 and dust shell. We use the radiative dust transfer code DUSTY 125,126 to fit the multi-wavelength data. We assume a spherically symmetric distribution of the dust with a ∝ r −2 density profile around the star, which is assumed to be a point source. We assume the dust grains to be composed of warm silicates as indicated by the O-rich composition of the photospheric spectra 127 , and with a MRN grain size distribution 128 (∝ a −3.5 ) with a minimum and maximum grain size of a min = 0.005 µm and a max = 0.25 µm. We fix the thickness of the dust shell to be Y = 5× the inner radius (r in ) of the shell; the model SEDs are found to be relatively insensitive to this assumption since the 3 − 5 µm emission arises primarily from the hotter, inner part of the dust shell.
We fit the observed SED of ZTF SLRN-2020 at two epochs during the outburst that have NEOWISE mid-IR coverage (≈ 120 and ≈ 320 days after outburst peak; shown in Figure  The best-fit parameters were derived using the median of the posterior sample distribution while their confidence intervals are derived from the 16th-84th percentile (68% confidence) interval of the distributions. The derived parameters and their uncertainties are listed in Extended Data Table 3, and shown in Extended Data Figures 5 and 6. During the outburst (≈ 120 d after peak), the SED is well described by an inner photosphere with a temperature of ≈ 8900 K surrounded by a warm dust shell with T d ≈ 1020K. The optical depth of the dust (τ V ≈ 1.8) is relatively low at this epoch. The foreground extinction is A V ≈ 3.6 mag, although this parameter is degenerate with the optical depth of the dust shell. While multi-epoch mid-IR data is scarce for the literature sample of red novae, similar parameters involving a relatively hot photosphere surrounded by a low optical depth dust shell were also derived at similar phases in DUSTY modeling for the nearby red nova M31-LRN-2015 53 . Towards the end of the outburst (≈ 320 d after peak), the DUSTY fitting suggests that both the internal photosphere and the surrounding dust shell have cooled substantially, while the optical depth of the dust has increased by a factor of ≈ 6.
The estimated temperature of the inner photosphere (≈ 9000 K) at the ≈ 120 d epoch is notably higher than that suggested by the spectral type of the photosphere (≈ 3500 K) obtained ≈ 60 days later. While the unavailability of contempornaeous data makes a direct comparison difficult, similar discrepanices have been observed in previous red novae with multi-wavelength datasets. Modeling for the nearby M31-LRN-2015 suggested a similar presence of a hot (> 5000 K) photosphere > 150 d after the eruption, while the optical spectral type was later than M9 (< 3000 K 59 ). In the case of the well-studied V1309 Sco, the strengths of the molecular features were shown to be substantially enhanced (i.e. the spetroscopic temperature was lower than the corresponding continuum temperature) compared to a normal stellar photosphere 19 , similar to that reported in V4332 Sgr 130 as well as V838 Mon 131,132 . As argued for previous objects, we therefore suggest that these differences arise due to the presence of extended circumstellar gas surrounding the merger remnant, that produces enhanced absorption in the molecular features.
We estimate the total mass of the dust shell at an epoch 53 as, where M d is the dust mass and κ V ≈ 50 − 100 cm 2 g −1 is the dust opacity for a typical dust-to-gas ratio 127 . In order to derive r in , we assume a distance of 4 kpc, as indicated by the distance estimates in Supplementary Information (Section 13), and list the corresponding radius and mass values in Extended Data Table 3. The observed increase in the dust shell inner radius between the ≈ 120 and ≈ 320 days epochs suggests an expansion velocity of v ej ≈ 35 km s −1 . The low inferred velocity of the shell suggests that signatures of such an outflow would be difficult to detect in our low resolution spectra. Together with the increase in the optical depth, the increased dust mass suggests that dust formation in the ejecta shifts the SED of the transient from the optical to the IR bands. Similar abrupt dust formation has also been observed in previous red novae 8, 19, 53 . The mass estimate provides a lower limit to the total (dust + gas) ejecta mass of ≈ 10 −6 M ⊙ , noting that only a fraction of the ejecta has likely reached the dust condensation radius within a year after the outburst 53, 133 .
13 Distance constraints and bolometric light curve We use the inferred foreground dust extinction to ZTF SLRN-2020 to place constraints on the distance to the source. In Extended Data Figure 7, we show published three dimensional Galactic dust extinction maps at the location of ZTF SLRN-2020 (created using the mwdust code 134 ), the estimated A V range inferred from the SED modeling as well as the allowed distance range for each map within the estimated A V range. As shown, the different extinction maps are roughly consistent with each other within the estimated A V range.
Given the possible systematic differences between the different maps, we conservatively derive a distance range of ≈ 2−7 kpc for the transient, placing this source well within the Galactic disk. As the distance ranges suggested by the different maps overlap consistently at a range of ≈ 3 − 4 kpc, we nominally adopt 4 kpc as the best estimate for the source. We caution however, that as is the case for the foreground extinction, the estimated distance is degenerate with the optical depth of the circumstellar dust shell (Extended Data Figure 5).
We use the distance and extinction estimate to derive the bolometric luminosity light curve of the optical transient around the peak of the outburst. We use the multi-color photometry from ZTF, binned in epochs of 3-days with coverage in all three (gri) filters, and fit a blackbody function to derive the effective temperature, luminosity and radius of the star for an estimated distance of 4 kpc. The ATLAS photometric filters cover very wide bandpasses over the g, r, and i bandpasses contemporaneously with the ZTF photometry, and hence we do use the c and o-band data for constructing bolometric light curves. The fitting was performed by χ 2 minimization using the Markov Chain Monte Carlo fitting Python library emcee 90 . The results are shown in Figure 2. For comparison, we also show the evolution of the corresponding parameters for well studied red novae in the literature 7 as well as the t −4/5 luminosity decay expected for the gravitational contraction of an inflated envelope surrounding the remnant 135 . We note that a foreground extinction significantly higher than the estimated interval (A V ≳ 5 mag) would imply unphysically high blackbody temperatures (>> 10 5 K), corroborating the estimated extinction and distance range.
The bolometric luminosity of ZTF SLRN-2020 shows an initial plateau at ≈ 10 35 erg s −1 for ≈ 15 days followed by a slow decline to 2 × 10 34 erg s −1 over the next ≈ 100 days. We estimate the duration of the plateau by fitting the bolometric light curve with an analytical model developed for light curves of Type II supernovae 136 and previously used to model light curves of red novae 53 . The fitting was performed by χ 2 minimization using the emcee code, and suggests a best-fit plateau duration of t p = 25.6 +5.6 −7.1 days as measured from the assumed time of eruption. The resulting fit is shown in Figure 2 and provides a reasonable match to the data; however, we caution that the plateau characteristics are relatively poorly constrained due to the short duration and subsequent smooth decline. The inferred bolometric luminosity from the DUSTY modeling at ≈ +120 days after outburst peak is consistent with the bolometric fitting at the latest epochs, while the DUSTY modeling at ≈ +320 days shows evidence of a plateau at very late phases.
We estimate the total energy (E), luminosity (L 90 ) and duration (t 90 ) of the outburst 7 by performing trapezoidal integration on the bolometric light curve. We estimate uncertainties on these parameters by creating 1000 realizations of the bolometric luminosity light curve using the MCMC uncertainty intervals and then repeating the calculations. The corresponding estimates are E ≈ (6.5±1.7)×10 41 erg, L 90 ≈ (8.0±5.1)×10 34 erg s −1 and t 90 ≈ 103±20 days for a distance of 4 kpc and foreground extinction of A V = 3.6 mag. Although the t 90 duration is found to be similar to the plateau duration (t p ) in previous objects 7 , we find t 90 ≈ 4 × t p likely due to a relatively small amount of unbound mass that contributes to the recombination luminosity on the plateau phase compared to the late-time gravitational contraction following the plateau. Considering the 68% confidence interval for the possible distance and foreground extinction, we find the total radiated energy and luminosity to be in the range (1.0 − 12.1) × 10 41 erg s −1 and (1.1 − 13.6) × 10 34 erg s −1 .
In particular, even if the source was placed on the farthest side of the Galactic disk (d ≈ 20 kpc), the estimated extinction would suggest a luminosity of ≲ 3.4 × 10 36 erg s −1 , suggesting that the interpretation as an very low-luminosity red nova (Figure 4) is insensitive to the assumed distance, which is not accurately constrained.
14 Progenitor photometry The progenitor of ZTF SLRN-2020 is detected only in the H and K filters of archival UKIRT images. Here, we attempt to constrain the progenitor star properties using the archival NIR photometry. In Extended Data Figure 8, we show the position of the progenitor in the M K vs. H −K color magnitude diagram as a function of distance along the Galactic disk using different three dimensional dust extinction maps. To constrain the progenitor mass, we also show stellar evolutionary tracks for stars with initial masses ranging from 0.8 − 3.0 M ⊙ from the MIST database 137 . Extended Data Figure 8 shows that the progenitor colors would be too blue to be consistent with any stellar tracks for distances larger than ≈ 8 − 10 kpc as per the D03 65 and M06 66 models, consistent with our constraints on the distance based on the outburst SED modeling.
We find that the location of the progenitor within the estimated A V range (shown as shaded region) is consistent with a 0.8 − 1.5 M ⊙ star on or evolving off the main sequence. Within the 90% confidence interval for the extinction, the photometric colors intersect the 1 M ⊙ stellar track at radii of ≈ 1 − 4 R ⊙ , which we adopt as the likely initial radius of the progenitor star. However, we are unable to further constrain the progenitor properties due to the large error on the photometric color. Specifically, the allowed progenitor flux and color are also consistent with a lower mass ≈ 0.8 M ⊙ progenitor star; in such a case, the engulfment would have occurred while the star was on the main sequence (given the current age of the Galactic disk population) likely driven by tidal interactions 32, 33, 140 .
The late-time WIRC photometry together with the DUSTY modeling shows that the star has both nearly returned to its original photospheric temperature (≈ 4000 − 5000 K; as suggested by the DUSTY model) while also having faded marginally below the progenitor brightness likely due to the optically thick dust shell surrounding the remnant. Given the constraints on the distance, we note that if related, the projected physical separation to the source detected next to the progenitor in the GSAOI images would be ≈ 1000 − 2000 AU, making it an extremely wide multiple star system; however, further follow-up is necessary to confirm this scenario.
15 Constraints on pre-outburst dust ZTF SLRN-2020 exhibits a mid-IR brightening in NEOWISE data starting ≈ 7 months prior to the optical outburst. Here, we use the pre-outburst brightening to constrain the evolution of the mass loss rate before the red nova outburst. In Extended Data  Table 1 and 2 for a contracted halo model). For each transient, we assume the outburst light curve follows the shape of the normalized template scaled to a peak absolute magnitude of M r = 2.0 (as inferred for ZTF SLRN-2020) with a nominal scatter following a Gaussian distribution with standard deviation σ = 0.5 mag. Using the Galactic spatial location of each simulated transient, we estimate the apparent magnitude evolution as a function of light curve phase using three-dimensional dust distribution maps 67 .
As in previous works 60 , for each value of the Galactic event rate, we simulate a population of outbursts such that the number of events follows a Poisson distribution with a mean equal to the assumed rate. We perform 100 simulations for each value of the assumed rate, repeating the actual ZTF observation schedule to estimate the number of recovered events. Using the same selection criteria for long-lived Galactic plane outbursts as used to identify ZTF SLRN-2020, we calculate the median number of events that pass our selection criteria within this simulated sample, as well as its 68% confidence interval. In order to derive the best-fit rate and its uncertainty given the one confirmed event in our search, we construct a distribution of the fraction of simulations that produce the one observed event as a function of the global rate. We fit a skewed normal function to this distribution to estimate a Galactic rate (with 68% confidence intervals) where r 0 is the Galactic event rate. However, we caution that i) the one observed event represents a lower limit to the actual number of events since it was not possible to obtain follow-up spectroscopy for all transients, ii) rate estimates from a single observed event are subject to large uncertainties and iii) the observed rates are subject to the recovery efficiency of the ZTF subtraction pipeline for faint transients in dense Galactic plane fields which have not been quantified for this simulation. Theoretically, we can put limits on the lowest plausible mass for the companion of ZTF SLRN-2020 by requiring that the characteristic change in orbital energy is larger than the radiated energy.
We note that this is a stringent lower limit; in V1309 Sco, the radiative efficiency, estimated as ϵ rad ≈ E rad /(GM 1 M 2 /2R 1 ) ≈ 3 × 10 −3 . Estimating M 2 this way for a higher radiative efficiency, we find an object as small as M 2 ≳ 6 × 10 27 g(ϵ rad /0.1) −1 = m ⊕ (ϵ rad /0.1) −1 , satisfies M 2 ≳ E rad /(ϵ rad GM 1 /2R 1 ) for a sun-like host star. Indeed, it is worth considering the possibility of higher radiative efficiencies because one might expect less adiabatic degradation of the radiated energy in ZTF SLRN-2020 relative to V1309 Sco because the transient photosphere is closer to the ejection radius of the original star.
Comparably, we can place a constraint on the plausible minimum radius for the captured object based on its rate of orbital decay and the luminosity that it injects into the surrounding envelope 4 , where R 2 is the radius of the captured companion and η = ρ/ρ, the ratio of the local to enclosed density. This ratio is a step function of radius at the stellar surface, for context, however, in a sun-like star at 90% of the radius of the star, η ∼ 3 × 10 −3 . This formula scales with R 2 2 , because in the low-mass regime of planetary captures, the geometric (rather than gravitational) cross section of the objects determines their decay rates 4 . Thus, an object as small as Earth radius would provide energy rapidly enough to power the transient peak of ZTF SLRN-2020. is pulled from the star and expelled with specific angular momentum similar to that of the planet, l 2 ≈ √ GM 1 a, a mass similar to the planet mass needs to be expelled to tighten the planetary orbit, ∆m ≈ ∆L orb /l 2 ≈ M 2 . In the energy model, we assume that ∆E orb ≈ GM 1 M 2 /2R 1 = f mech ∆mv 2 esc /2, where f mech is an efficiency factor less than unity to account for radiative losses.
Solving for ∆m, we find ∆m ≈ f mech M 2 . Thus, the pre-outburst detections of ≈ 10 −4 M ⊙ of gas and dust point to a planet of at least similar mass, or approximately 0.1M J .
We derive more detailed models that consider the temporal evolution in the angular momentum model from the separation of stellar Roche lobe overflow to the point of contact. We use the the RLOF software 148 , which integrates the cumulative mass loss given coefficients of specific angular momentum loss motivated by hydrodynamic simulations of this process 25, 29 . These models have successfully explained the dynamics of orbital decay and the observations of increasing obscuration in V1309 Sco 22, 24, 25, 149 , but we caution that we must extrapolate them to several orders of magnitude lower companion mass than they have been systematically tested for. We must assume several parameters of the binary system. We adopt a 1M ⊙ and 2R ⊙ fiducial model, and compare to 1M ⊙ , 1R ⊙ and 1M ⊙ , 4R ⊙ examples to demonstrate the degree to which varying this within the uncertainty affects the result. We further assume that the structure of the star is approximated by a Γ s = 4/3 polytrope, that the gas adiabatic index is γ = 5/3 and that the star is not initially corotating with the companion orbit. These parameters each moderately affect the specific angular momentum loss that accompanies mass loss and orbital decay (∼ 50% differences in the predicted mass loss with varying parameters 25 ). We then convert the estimated total mass lost to a dust mass assuming that the dust-to-gas mass ratio 127 is 5 × 10 −3 .
The dust masses expected from modeled systems are shown in Figure 4. We see that the modeled dust masses and non-detections in the pre-outburst data are broadly consistent with a model that exhibits an exponentially increasing circumbinary distribution of gas and dust as it progresses toward merger. In overall normalization, our curves are most consistent with a mass ratio of q ≈ 10 −3 to 10 −1 , implying a companion mass of ∼ 10 −2 M ⊙ , or 10M J . This mass is considerably larger than the energy or angular momentum order of magnitude scalings alone indicate. This is because of the characteristic timescale of the exponential decay of the orbit in the RLOF model, which implies that most of the mass loss is concentrated in the final orbits prior to coalescence. We emphasize here that to make more robust inferences will require calibrating these energy and angular momentum based models of orbital decay to the planetary, rather than stellar, mass range in question. This mass estimate is also larger than the mass inferred from the recombination powered transient's mass and energy budget, as discussed in the section above.
However, both mass estimates are strongly suggestive of a planetary (rather than stellar or brown dwarf) mass companion object producing the ZTF SLRN-2020 transient.

Understanding a Planetary Engulfment
The planetary mass range inferred from various lines of evidence points to the fact that ZTF SLRN-2020 probes new parameter space among the red novae. Hydrodynamic simulations of stellar coalescence have suggested that as the mass ratio decreases, a smaller fraction of the shock-heated material is truly unbound from the system and ejected 28, 29 . The remaining, bound material settles into a shock-heated, rotating envelope in which spiral shocks play an important role in redistributing energy and angular momentum 29 .
An extension of this is the limit that the smallest companions should be too insignificant to eject any mass from the host star when they are engulfed. Thus, we suggest that the combination of pre-cursor and lightcurve data shown in Figure 4 are indicative of a system in which the majority of the shock-heated stellar atmosphere remains bound. Indeed, this is strongly suggested by inferred velocities in Figure 4, which are much less than the local escape velocity at 3 × 10 11 cm of ≈ 300 km s −1 .
If most of the stellar atmosphere mass remains bound following the engulfment, this would attribute the plateau luminosity to the recombination of a comparatively small unbound mass fraction (≈ 3 × 10 −5 M ⊙ ). It also informs our understanding of the brightness and colors of ZTF SLRN-2020 over the hundred-day observing window. ZTF SLRN-2020's temperature remains higher than other red novae and its photosphere radius, following its initial expansion to ≈ 3 × 10 11 cm ≈ 4.3 R ⊙ , gradually decreases rather than growing. These trends both indicate that following an initial period of mass ejection, we may be witnessing the hydrodynamic and thermal relaxation of still-bound stellar material that has been disturbed by the engulfment of the substellar companion. Indeed, the late-time light curve follows a roughly t −4/5 power-law, which is consistent with observations of other merger remnants 135 .
We note that the small ejected mass is consistent with the non-detection of molecular line emission in the SMA and VLA observations. While a detailed molecular excitation analysis is beyond the scope of the work given the non-detection, we can interpret these data with comparisons  Figure   4). Noting that ZTF SLRN-2020 was observed at a much earlier phase and therefore may have different excitation conditions, the non-detection of molecular emission is consistent with the small inferred ejecta mass in the eruption. While SiO maser emission at 43 GHz was only detected for V838 Mon 151, 152 , we use similar scaling arguments and accounting for the larger distance to V838 Mon, we estimate the maser line flux in ZTF SLRN-2020 to be ∼ 10 −4 × fainter (≲ 0.5 mJy at peak for a resolution of ≈ 0.5 km s −1 ) and hence consistent with the non-detection (≲ 6 mJy).
A particularly relevant consideration is the relative density of planets and their host stars.
Planets are tidally disrupted outside of their host stars when their density is lower than that of the star they orbit 3 . When the planetary density is higher than that of the host star, the planet can plunge into the stellar envelope, shock-heating its surroundings as it is engulfed. Given the similarity of ZTF SLRN-2020 to other, engulfment-powered red novae transients, we strongly suspect that the nears the mass range of brown dwarfs, which extends from roughly 13M J to 80M J . Could ZTF SLRN-2020's companion object be a brown dwarf rather than a planet? We suggest that this is substantially less likely because brown dwarf companions at similar separations to hot Jupiters are nearly an order of magnitude more rare 1 . The close-companion portion of this "brown dwarf desert" may be caused by the more-rapid tidal decay of these objects orbits into their host stars: even if these objects were formed at similar rates to hot Jupiters, many would coalesce with their host stars during the pre-main sequence. Since ZTF SLRN-2020 shows no evidence of being a young star there seems to be substantially more support for a planetary companion.
merger with their host stars through a combination of tidal orbital decay (due to asynchronous rotation of the host star with the companion orbit) and stellar evolution (which increases the radius of the host star while lowering its spin rate) 33 . For example, the currently observed population of hot Jupiter exoplanets all have inferred tidal decay times of less than 10 10 yr. Tidal decay has been directly detected for WASP 12b 34,35,154,155 . Even without tidal decay, each of these planets will surely merge with its host star as these stars increase in radius on the sub giant and giant branches.
Additionally, planets can be scattered or torqued into eccentric orbits that cause them to interact with their host stars at periapse [156][157][158][159] .
There are uncertainties in the propagation of this observed distribution to an engulfment event rate. As Metzger et al 3 have noted, the distribution of known planets at the closest separations implies a much higher coalescence rate than the apparent number at larger separations (see their Figure 1). This points to a need to resupply planets to these close-in orbits to maintain the currently-observed distribution. The occurrence of very close, Jupiter-like planets implies a coalescence rate of ∼ 0.1 yr −1 in the galaxy 3 , while the number of hot Jupiters at slightly larger separations (with tidal decay timescales of 10 9 to 10 10 ) appears to be fewer. Whether this relates to observational selection effects or is an intrinsic property of the population is currently under discussion 155 .
With these caveats in mind, we can produce an order of magnitude estimate of the planetary engulfment rate as follows. The total mass of stars in the Milky Way is ∼ 6.4 × 10 10 M ⊙ 160 . In a Salpeter IMF, approximately 1/3 of this mass comes from stars greater than a solar mass (which is roughly the fraction that evolve off of the main sequence in less than 10 10 yr). This implies that roughly 10 10 stars have evolved off of the main sequence in the history of the galaxy. If 1% of these stars hosted hot Jupiters 1, 161 , then the time-averaged hot Jupiter engulfment rate would be (10 8 /10 10 yr) ∼ 10 −2 yr −1 . Considering planets more broadly, the known exoplanet population is consistent with a fraction on the order of unity of all stars hosting some kind of planet. Thus, the overall planetary engulfment rate from this broader population would be ∼ 1 yr −1 .
In practice, these time-averaged estimates represents a lower limit because they do not account for the apparent resupply of planets to close in orbits or stellar impacts 156  ZTF SLRN-2020's emerging properties will be crucial to guide future search efforts.
An even larger body of effort has been devoted to understanding the long-lasting effects of substellar engulfment 5,43,44,167,[171][172][173][174][175] , with the motivation that if a large fraction of stars engulf one or more planets as they evolve, perhaps these events leave long-lasting impacts on the observable stellar characteristics. One possible property that has been invoked is secular declines in luminosity, which might trace the cooling of a star following the mechanical addition of heat from an engulfment 4, 41 .
Another is enrichment of the stellar atmosphere by lithium and other elements carried by the planet bulk composition in higher abundance than the star 44 . Lithium is thought to be a particularly useful tracer because its fragile structure is dissociated in stars, but persists in substellar objects 5,43,174,176,177 .
Observations of the remnant of ZTF SLRN-2020 at relatively high spectral resolution could search for signs of lithium enrichment and serve as a useful benchmark for these theories, as well as compare to recent studies of other red novae remnants 178 . A final possible signature that has been extensively considered is stellar rotation, with enhancements in otherwise-slow rotation rates thought to be be possible due to the deposition of angular momentum as the planet is engulfed 46, 179, 180 .
Since the total angular momentum budget of a substellar engulfment is well-known, the impact on the observable spin rate (as well as the evolution of that spin rate) becomes deeply revealing about the transport of angular momentum in the stellar interior. As future work traces each of these properties in ZTF SLRN-2020 and similar future transients, we anticipate that these events will serve as a crucial missing link in connecting the properties of observed planetary systems to the transients they produce and their effects on their host stars.