Bridging the Gap between Luminous Red Novae and Common Envelope Evolution: The Role of Recombination Energy and Radiation Force

Luminous red novae and their connection to common envelope evolution (CEE) remain elusive in astrophysics. Here, we present a radiation hydrodynamic model capable of simulating the light curves of material ejected during a CEE. For the first time, the radiation hydrodynamic model incorporates complete recombination physics for hydrogen and helium. The radiation hydrodynamic equations are solved with Guangqi. With time-independent ejecta simulations, we show that the peaks in the light curves are attributed to radiation-dominated ejecta, while the extended plateaus are produced by matter-dominated ejecta. To showcase our model’s capability, we fit the light curve of AT 2019zhd. The central mass object of 6 M ⊙ is assumed based on observations and scaling relations. Our model demonstrates that the ejecta mass of AT 2019zhd falls within the range of 0.04–0.1 M ⊙. Additionally, we demonstrate that recombination energy and radiation force acceleration significantly impact the light curves, whereas dust formation has a limited effect during the peak and plateau phases.


INTRODUCTION
Since the discovery of V1309 Sco (Tylenda et al. 2011), a luminous red nova (LRN), also a confirmed binary merger, many more new LRNe have been discovered (Kurtenkov et al. 2015;Blagorodnova et al. 2017;Cai et al. 2019;Pastorello et al. 2019;Blagorodnova et al. 2021;Pastorello et al. 2021a,b;Cai et al. 2022a,b), to name a few.It is argued that, at least for some of the LRNe, the progenitors are binary stars that have undergone common envelope evolution (CEE) (Ivanova et al. 2013a) -one of the most mysterious events in binary evolution.Bridging the gap between the observables of LRNe and theoretical models of CEE could be crucial to getting a better understanding of binary evolution (Chen et al. 2024), including cataclysmic variables (Warner 2003), X-ray binaries (Reig 2011), gravitational wave sources (Renzo et al. 2021), and type Ia supernovae (Wang & Han 2012;Liu et al. 2023).
Corresponding author: Zhuo Chen (陈卓) chenzhuo astro@tsinghua.edu.cnObservationally, LRNe exhibit three prominent characteristic features.Firstly, they all have extended plateaus with decreasing effective temperatures in the light curves.This phenomenon suggests a substantial mass ejection from the central object, which cools as it expands outward.During this phase, hydrogen and helium may recombine and release their latent heat, contributing to the thermal energy of the plateau (Ivanova et al. 2013b(Ivanova et al. , 2015;;Matsumoto & Metzger 2022).Secondly, a notable number of LRNe present a peak preceding the plateau in their light curves.Modeling this peak is challenging due to its exponential increasing and decreasing phases.Pejcha et al. (2017) employed a hydrodynamic model with post-processing to fit the slowly increasing luminosity before the peak, attributing the peak to an eruption-like event.Although Matsumoto & Metzger (2022) successfully fitted the rapidly decreasing luminosity after the peak with their hot ejecta model in 1D, the rapidly increasing phase remains a puzzle.Thirdly, LRNe frequently exhibit strong Hα emission (Munari et al. 2002;Stritzinger et al. 2020;Blagorodnova et al. 2021;Pastorello et al. 2021a), indicating a arXiv:2402.05686v2[astro-ph.SR] 17 Feb 2024 possible collision between the ejecta and circumstellar matter (Blagorodnova et al. 2020).
A prevailing theory on the origin of LRNe implicates the progenitor binary undergoes CEE and experiences a plunge-in phase accompanied with ejection of some envelope material.Numerous 3D simulations have been conducted to study the CEE dynamics (Nandez et al. 2014(Nandez et al. , 2015;;Ivanova & Nandez 2016;Nandez & Ivanova 2016;Ohlmann et al. 2016;Chamandy et al. 2018;Iaconi et al. 2019;Prust & Chang 2019;Reichardt et al. 2019;Sand et al. 2020;González-Bolívar et al. 2023;Röpke & De Marco 2023;Chamandy et al. 2024).However, due to the complexity of physical processes and the multi-scale nature of the CEE problem, there is still no consensus on how much mass is ejected during a CEE, even if 3D simulations start with the same initial conditions.Furthermore, predicting the observational appearance of these simulations is inherently difficult (Hatfull et al. 2021).Notably, radiation hydrodynamics, a crucial factor in CEE, is often absent from these 3D models due to its difficulty and is easier to consider in 1D models (Soker et al. 2018;Bronner et al. 2023;O'Connor et al. 2023).Recent work by Matsumoto & Metzger (2022) demonstrates the effectiveness of a cooling shell model, devoid of radiation transport and radiation hydrodynamics, in producing LRNe-like light curves.
However, there is still a lack of first principle models that can relate CEE to LRNe.In CEE, the plunge-in phase marks the rapid conversion of the gravitational potential energy into the kinetic energy.Some kinetic energy would be converted to thermal energy via a shock between the plunge-in star and the envelope (MacLeod et al. 2017).We would expect the shocked gas to be close to adiabatic because of the high optical depth inside a common envelope (CE).The shocked gas may be accelerated in the radial direction due to high radiation pressure and radiation force acceleration -this has been overlooked previously -and become an eruptionlike event, eventually appearing as an LRN.CEE is a 3D problem by nature.However, solving 3D radiation hydrodynamic equations with complex physics is computationally demanding.Therefore, as the first step, and hoping to resolve the microphysics better, we start by approximating the radiation hydrodynamic problem with 1D spherical symmetry.In a follow-up work, we will revisit this problem with 2D axisymmetric models.
This Letter introduces the first radiation hydrodynamic model that comprehensively incorporates recombination energies of H and He, radiation transport, and radiation force acceleration.The model can successfully produce both the peak and plateau phases of the light curve, allowing for the estimation of ejecta mass through curve fitting.The organization of the Letter is as follows: Section 2 presents the physical model, governing equations, initial and boundary conditions, and numerical setups, including Adaptive Mesh Refinement (AMR) criteria.Simulation results are presented in Section 3. The Letter concludes in Section 4.

Physical model
We adopt the following radiation hydrodynamic equations to model the evolution of the LRNe in a 1D spherical coordinate, where ρ, v, p, and E are the density, radial velocity, pressure, and total energy of the gas.In addition, a rad (to be explained later) and g = GM ⋆ /r 2 are the radiation and gravitational force acceleration, respectively, where G is the gravitational constant and M ⋆ is the mass of the central object.The time and coordinate are denoted by t and r.The radiation-related variables E r , F r , and G are explained later in this section.We adopt a simple hydrogen and helium mixture equation of state (EoS, see Appendix A), where i is the species index, n i is the number density of species i, e g and T g are the internal energy and temperature of the gas, k b is the Boltzmann constant.Throughout this Letter, we assume the hydrogen mass ratio X = 0.74 and the helium mass ratio Y = 0.26 for simplicity, because metal does not contribute much to gas thermodynamics.The hydrodynamics, together with the complex EoS, is solved by the approximate HLLC Riemann solver in Chen et al. (2019).
To solve the radiation transport problem, we use the flux-limited diffusion (FLD) approximation, which relates the radiation flux F r to the gradient of radiation energy ∂E r /∂r with a flux limiter λ (Levermore & Pom-raning 1981), i.e., where κ R is the Rosseland mean opacity (see Appendix B).The flux limiter has the property that ∂Er ∂r optically thick, cE r optically thin. (10) The radiation and gas energy coupling is modeled by solving the following equations implicitly with subtimesteps, where G = κ P ρc(E r − a r T 4 g ) is the energy coupling strength, κ P is the Planck mean opacity, and a r is the radiation constant.Meanwhile, the radiation and gravitational acceleration are integrated explicitly through, where a rad = κ R F r /c is the radiation force acceleration and c is the speed of the light.We use Guangqi(Chen & Bai in prep) to solve Equation 1-4.Guangqi is a second-order in time and space accurate and finite volume radiation hydrodynamic code.It has HLLC Riemann solvers, realistic EoS (Chen et al. 2019), and adaptive mesh refinement (AMR).Guangqi solves the radiation transport problem with FLD approximation implicitly (similar to Kolb et al. (2013)), using iterative solvers from Petsc (Balay et al. 1997(Balay et al. , 2019)).Currently, Guangqi has spherical and Cartesian geometry in 1D and 2D.

Initial and boundary conditions
We set our initial condition to be an outflow at escape velocity with a constant mass loss rate, where ρ 0 = 10 −13 g•cm −3 , T 0 = 1000K, and r in = 10R ⊙ is the inner radius of the computational domain.The constant mass loss rate is Ṁ = 2.30 × 10 −7 M ⊙ •yr −1 , which is low compared to the mass loss rate of the ejecta of the CEE.Radiation's initial condition is assumed to be in local thermal equilibrium (LTE) with gas.The initial mass and energy in the computational domain are 8.52×10 −7 M ⊙ and 3.68×10 35 erg, respectively.They are significantly smaller than the ejecta's mass and energy, and our light curve results are insensitive to the initial condition.
The outer boundary is free, i.e., the gas and radiation can leave the computational domain freely.In our simulations, we set the out boundary at r out =4000R ⊙ .The free boundary for the radiation is, This outer boundary condition means that the radiation flux is optically thin and outward.We can confirm that the outer region of our computational domain is indeed optically thin (see Section 3.2) and always outward.We calculated the luminosity at the outer boundary by, The inner boundary is time-dependent.We specify the time-dependent density, velocity, and temperature of the ejecta at the inner boundary, i.e., [ρ(t), v ej (t), T g (t)].
where f ej (t) can be a time-dependent factor, Ṁ is the mass loss rate, and we introduce α to characterize the temperature of the ejecta.When the ejection of the CEE stops, the inner boundary of the gas is changed to free to let the fallback gas pass through the inner boundary, therefore, we do not model any fallback shocks and fallback accretion disks.Meanwhile, radiation transport is turned off at the inner boundary.Mathematically, we use the zero-gradient (∂E r /∂r = 0) radiation inner boundary condition.We assume the gas and radiation are in LTE at the inner boundary because the density and opacity are very high inside the CE, and the thermal timescale is short.

Numerical setups
We adopt a uniform base grid with spherical geometry.The computational domain is r ∈ [10, 4000]R ⊙ , and the base resolution is N = 1536.We add 5 levels of static mesh refinement where r ∈ [10, 15]R ⊙ to resolve the strong gradient of the gravitational potential.Each level of mesh refinement doubles the resolution.We also adaptively refine zones with temperature gradients up to 5 levels.The mesh refine and derefine criterion are as follows, refine where i±1 in the subscript represents the cell index (not to be confused with species index).We can capture the shock, radiative layer, and radiation-dominated zones with photon trapping with AMR.The finest cell has a length of 8.12 × 10 −2 R ⊙ .The Courant-Friedrichs-Lewy number is 0.95 in our simulations.The simulation time is 4.5 × 10 6 s=52.08 days.

RESULTS
We first show some simple simulations and get a sense of the correspondence between the properties of the ejecta and the light curves.After that, we fit the light curve of AT2019zhd with a more complex ejecta.

Simple ejecta
For simplicity, in this subsection, we set M ⋆ = 6M ⊙ , r in = 10R ⊙ .We consider seven models and list the physical properties of the seven models in Table 1.In what follows, we refer to a model as radiation-dominated if E r /e g ≫ 1, and as matter-dominated if E r /e g < 1.
In particular, we calculate the total energy E total of the ejecta by, where t ej is the duration of the ejection.
Figure 1 shows the light curves of the simple models.We can see that the light curves of m02a045v2 and m01a18 have pronounced peaks, i.e., rapid, exponential increases and decreases, while m25a18 and m25a09 have long plateaus.The light curves of other models are in between.Because the physical system is highly nonlinear, we provide a qualitative discussion here.The optical depth of an ejecta decreases as it expands, and the ejecta cools faster as the optical depth decreases.The more mass in an ejecta, the higher the optical depth; the faster the ejecta, the faster the optical depth decreases.Therefore, m25a18 and m25a09 have long plateaus because their ejecta are the most massive ones among the seven simple models; m25a18 has a longer plateau than m25a09 because m25a18 has a higher energy budget and thus a longer cooling timescale.On the other hand, m02a045v2 has a sharper peak than m01a18 because its expansion speed is faster, resulting in a more rapid decrease in the optical depth and a shorter cooling timescale.Because m02a045v2 and m01a18 have the same amount of energy budget, m02a045v2 has a higher luminosity peak.Overall, we can also relate that we observe the luminosity peak formation (exponential increase and decrease) in the case of radiation-dominating ejecta.In contrast, the plateau is the feature observed in the case of matter-dominated ejecta.

Fitting AT2019zhd
We choose AT2019zhd as a fitting example because it is a recently well-observed LRN and its light curve resembles V1309 Sco, whose progenitor is a confirmed binary (Tylenda et al. 2011;Pastorello et al. 2021a).The mass of the progenitor of AT2019zhd is unclear (Pastorello et al. 2021a).We take the following steps to estimate the total mass of the merger product.
3. The total radiation energy released during the peak and the plateau of AT2019zhd is roughly 10 times more than V1309 Sco.
If we assume that the luminosity energy mainly comes from the release of the gravitational potential energy, it may scale as, This scaling relation motivates us to set the mass of the central object to be the same as was adopted for the simple models, 6M ⊙ , and we set α = 0.018 in this subsection.
Unlike the case of simple models, the inner boundary conditions here are time-dependent, see Figure 2. As a result, the time-dependent ejecta is initially radiationdominated (denoted by the orange color in Figure 2) and then transits to matter-dominated.We anticipate that the initial high-temperature ejecta is produced by a shock between the rapidly plunge-in companion star and the envelope (MacLeod et al. 2017).The plungein speed v p is comparable to the Keplerian speed at the plunge-in radius r p , but much larger than the envelope's speed v env .On the other hand, the kinetic energy required to produce E r and e g in the ejecta, converted to speed, can be calculated by, where ρ is the ejecta's density.Due to energy conversion, The maximum v source in Figure 2 is 448km•s −1 , which means that r p < 5.69R ⊙ .More sophisticated stellar and binary evolution analysis should be considered to further narrow down r p and M ⋆ (Ge et al. 2010(Ge et al. , 2015(Ge et al. , 2020)).We carry out simulations of two sub-models: the shock model and the shock-free model.The shock model has an ejecta with a slightly increasing speed during the late stage and the shock-free model has an ejecta with a slightly decreasing speed during the late stage.Figure 2 shows the time-dependent variables of the two submodels.The functions that generate these variables can be found in Appendix C.
Figure 3 shows the light curves of the shock and shockfree models.We have shifted the time axis to let the peak be located at t = 0. We can see that the light curves of the shock and shock-free models both resemble AT2019zhd.Consequently, the time evolution of the shock-free model is similar to the shock model, therefore, we just show the time evolution of the shock model in Figure 4 to save space.In particular, we calculate the optical depth of the radiation flux by, The time axis is adjusted to be consistent with the light curves.The time evolution of H and He species can be found in Appendix D, since they are determined by ρ and T g .There are several prominent results in the time evolution plots, see Figure 5.
1.In panel 2, the temperature profile has a spike near the peak of the luminosity.This is when the radiation-dominated ejecta becomes optically thin and the radiation heats the surrounding gas.
2. The observationally important ejecta material is the one that has a temperature of about 5000K.The ejecta's matter that is colder than 5000K is closer to the observer but is mostly transparent, revealing behind it the ejecta's matter with a temperature of about 5000K.The time-evolution of the surface log 10 T g = 3.7 (see panel 2) resembles the shape of the light curve in Figure 3.
3. Panel 3 shows that the a rad exceeds the g at small radii during the first 30 days, and at large radii during the late stage.This phenomenon is consistent with the evolution of κ R .When 10 4 ≲ T g ≲ 10 5 K, most of the gas is ionized and κ R is large, radiation flux can be significant.When T g < 1400K, dust may form and provide the opacity, and the a rad becomes  26, and the temperature of the ejecta.On the right panel, from the top to the bottom, they are the mass loss rate, cumulative mass loss, and Er/eg, respectively.The shock model and shock-free model also differ in the duration of the ejecta.significant again.At this stage, the LRNe is similar to the asymptotic-giant-branch (AGB) stars (Höfner & Olofsson 2018).
4. In panel 4, the radius of τ R = 1 expands rapidly when the inner region of the ejecta cools off because dust can form at low temperatures and obscure the object.
5. The inner region (r < 2000R ⊙ ) transits from radiation-dominated to matter-dominated as the ejecta expands into the ambient and cools off.It can be seen in panel 5 of Figure 4.
Our model is comprehensive in the sense that we incorporate many physical processes.To address the impact of each specific physics on the light curve, we run some companion simulations by turning off the corresponding physics in Appendix E, F, and G.We show that our calculations are converged in Appendix H.

CONCLUSIONS
We present a 1D radiation hydrodynamic model that can obtain light curves of LRNe events.For the first time, our model incorporates recombination physics, radiation transport, and radiation force acceleration of a CEE into consideration.There are several key new physics that we determined to be intrinsic to CEE and LRNe.
First, we find that the light curve's peak and plateau are formed differently.Specifically, the peak is formed primarily by radiation-dominated (E r /e g ≫ 1), and the plateau is formed by matter-dominated ejecta (E r /e g < 1).
Secondly, radiation force plays a significant role in driving the ejecta outwards (panel 3 of Figure 4).This occurs where opacity is high, specifically when both H and He are still ionized.Radiation force may help to eject matter that initially has sub-escape velocity.At last, rapid dust formation, assuming it occurs, absorbs radiation and accelerates during the late LRNe stages, resembling environments near AGB stars.While the radiation force helps to accelerate the dusty matter, this has minimal impact on the peak and first plateau of the light curve.
There are several limitations in our model right now.For example, the temperature of the ejecta is described by Equation 21.Thus, we can not strictly distinguish the effect of changing M ⋆ /r in and α.However, one can constrain M ⋆ by a scaling relation, for example, Equation 25, and can constrain r in (which is affected by r s ) with long-term binary evolution models.Moreover, one can further constrain the thermalization efficiency α with realistic 3D radiation hydrodynamic calculations.
The effect of changing initial e g and E r in the ejecta is also not fully independent, but radiation transport can, to some extent, help to distinguish the radiationdominated cases from the matter-dominated cases.
We also find that shock and shock-free models can produce similar results.In our models, shock and shockfree models refer to whether the speed of the ejecta was increasing or decreasing with time.Indeed, a shock can generate heat and prolong the plateau phase, but adding mass and internal energy to the ejecta can also result in a prolonged plateau in the light curve.We may need other observables, such as Hα emission, to distinguish these two scenarios.We can also explore the effect of pre-CEE asymmetric matter distribution around the central binary object; this will be the subject of the follow-up 2D work.
Even with a 1D model, we can fit the AT2019zhd light curve reasonably well.Motivated by observations and scaling relations (Section 3.2), we adopted the mass of the central object to be M ⋆ = 6M ⊙ , and show that our choice of M ⋆ and r in are reasonable.With this choice, we estimate that the mass of the ejecta could be ∆M ∈ [0.04, 0.1]M ⊙ 1 .As a result, the 1D model that we proposed in this Letter can provide a crucial bridge between CEE theory and the observables of LRNe.which is calculated by, Thus, to calculate c s , we need to know (∂s/∂ρ) T and (∂s/∂ρ) ρ , which means (∂n i /∂ρ) T and (∂n i /∂ρ) ρ .They can be calculated by thermodynamic equilibrium equations, i.e., the Saha equations.To further simplify the EoS model, we divide the ρ − T diagram into pure and mixture states (Chen et al. 2019).We define pure states if the mass fraction of a single species (not including electron) exceeds 0.99999 and all other cases are mixture states.Figure 5 shows the division of H and He.The blue color represents the pure states of H and He.
In doing so, we can solve fewer Saha equations.As a first example, let us consider hydrogen's ionization and dissociation reactions (the red color zone in the Solve Equation A23 gives us n He 2+ , then n He + and n e − can be solved by Equation A20 and A22, and n He is calculated from A21.Meanwhile, we write down the differential form of Equation A19, A20, A21, and A22 and calculate (∂n i /∂ρ) T and (∂n i /∂ρ) ρ of helium species.
If only one reaction is present, the thermodynamic system is trivial to solve.disassociation+H ionization 10 8 g cm 3 10 9 g cm 3 10 10 g cm 3 10 11 g cm 3 10 12 g cm 3 10 8 g cm 3 10 9 g cm 3 10 10 g cm 3 10 11 g cm 3 10 12 g cm 3 Figure 7.The Rosseland mean opacity that includes dust opacity at different densities.Solid and dashed lines denote the gas and dust opacities, respectively.The red zone indicates the dust formation temperature, the gray zone indicates the temperature of hydrogen dissociation and H -formation, and the orange zone indicates the ionization temperature of H and He.
Finally, we can calculate all the thermodynamic variables.Figure 6 shows the specific heat capacity and adiabatic index we use in this Letter.

B. OPACITY
The opacity table combines the gas opacity table of Malygin et al. (2014) that dominates over 1500K, and the dust opacity table (1% of dust to gas ratio) of Semenov et al. (2003) that dominates below 1100-1200K, depending on the density.For temperatures in between, the maximum of the value of the two tables is taken.Figure 7 shows κ R at different densities and κ P is similar.
C. TIME-DEPENDENT BOUNDARY CONDITIONS When t < t 1 , the ejecta is radiation-dominated, where x = 20(t − t 1 )/t 1 . (C31) There are 13 parameters in this model, we list all the parameters just to show the functional form of the ejecta.
In terms of the shock and shock-free models, they have many parameters in common.The common parameters are: x 2 = 0.55.We list the different parameters in Table2.
We anticipate that there could be other fitting functions with more/less parameters.One should explore different functional forms and a thorough parameter space study is deemed appropriate in the future.

D. THE TIME EVOLUTION OF HYDROGEN AND HELIUM SPECIES
Figure 8 shows the time evolution of the mass fraction of H and He species of the shock model.For example, the mass fraction of H 2 is defined as, Although the mass fraction of each species can be derived from ρ and T g already, there are two features we would like to clarify.
1. H 2 forms in the late stage.There are two channels for H 2 formations.At high density, H can form H 2 directly (Omukai et al. 2005).When dust is present, H 2 can form on dust (Wakelam et al. 2017).In the time evolution, we find that the H 2 formation region largely coincide with the dust formation region (can be seen from the log 10 κ R panel in Figure 4) or with the high-density region (can be seen from the log 10 ρ panel in Figure 4).Therefore, Saha's EoS model is suitable for this problem.
2. There is a line showing the partial ionization of H + in the H + panel.This is due to Zeldovich's spikes in the radiation hydrodynamic calculation after a shock.The shock is caused by the fast ejecta colliding with the ambient.However, we tend not to ascribe the Hα emission to this line because the decreasing width of Hα in observations (Pastorello et al. 2021a) indicates that the shock decelerates.We believe that the Hα emission is caused by a marginally escaping ejecta colliding with the pre-existing decreation disk (Pejcha et al. 2017;Metzger & Pejcha 2017).However, the 1D model is only shell-like and cannot model disks.

E. THE IMPACT OF REALISTIC EOS
EoS can play an important role in radiation hydrodynamics.We compare our simulations with realistic EoS to the widely used in hydrodynamic codes γ-law EoS, described by where µ = 1 is the mean molecular weight.With γ-law EoS, gas energy does not contain latent heat.As a result, the cooling timescale is shorter, and all the light curves ran with γ-law EoS last for a shorter time.We can also see that the peak is less affected by the change in the EoS because the radiation-dominated ejecta governs the peak.The internal energy e g plays a minor role in this stage.

F. THE IMPACT OF RADIATION FORCE
Radiation force was largely ignored in previous hydrodynamical simulations of CEE.However, we find that the a rad is larger than g at high (> 10000K) and low temperatures (< 1400K); see panel 3 in Figure 4. To demonstrate the importance of radiation force, we compare our original simulations' light curves to those without radiation force (a rad = 0).The results are presented in Figure 10.
In the case of no radiation force in the momentum equation, the luminosity drops off faster.This faster drop is due to more material in the ejecta falling back to the central object under the influence of the gravitational force.The drop happens so early that the gas is still fully ionized and should be pushed away by the radiation force.In our model, the fallback material will pass the inner boundary and no longer contribute to the luminosity.

G. THE IMPACT OF DUST FORMATION
It is argued that molecules and dust form during the late stage of CEE (Nicholls et al. 2013;Kamiński et al. 2015;Banerjee et al. 2015;Blagorodnova et al. 2020;Iaconi et al. 2020;MacLeod et al. 2022).Numerical models also show that dust may provide an additional push away from the binary merger due to radiation acceleration (González-Bolívar et al. 2023).To test the impact of dust formation, we modify the opacity table (Figure 7) by decreasing κ R and κ P of dust by 100 times and name it the dust deficit model.We find that modification of dust opacities does now affect light curves during the peak and the plateau.This is unsurprising, as the temperature is too high during the early stage of LRNe for dust to form.The minor difference is that dust-deficit models have slightly higher luminosity during the very late stage (see Figure 11), mainly because the optical depth of the dust region is smaller.
Figure 11 also shows the time evolution of log 10 (a rad /g) and log 10 κ R of the dust-deficit model.Comparing Figure 11 to Figure 4, we notice that the difference is mainly in the right corner of the figure, caused by the lowered dust opacity.The time and radial location of the dust formation is not changed.It indicates that the dust formation may have a limited impact on the radiation hydrodynamics of LRNe events during the early stage.In the late stage, due to radiation force acceleration, dust formation may play a role in shaping the ejecta, but still not a major physical mechanism of unbinding the ejecta (Bermúdez-Bustamante et al. 2024).

H. CONVERGENCE
Convergence is very important for simulations that have iterative solvers and mesh refinement.It is also well-known that radiation hydrodynamic equations may be sensitive to choice of timesteps.We show that our simulations approach convergence by comparing the light curves of our original resolution to a lower base resolution (N = 768), which also doubles the timesteps.Figure 12 has the light curves of the high and low resolutions of the shock and shock-free models.We can see that they are similar, indicating that our results are close to convergence.

Figure 1 .
Figure 1.The light curves of simple simulations with different mass and temperature of the ejecta.The physical parameters of the ejecta are listed in Table1.

Figure 2 .
Figure2.The solid and dashed lines are the time-dependent inner boundary conditions of the shock and shock-free models, respectively.The red line indicates the escape velocity at the inner boundary.The orange region denotes the radiation-dominated ejecta.On the left panel, from top to bottom, each plot shows the velocity of the ejecta, vsource calculated from Equation26, and the temperature of the ejecta.On the right panel, from the top to the bottom, they are the mass loss rate, cumulative mass loss, and Er/eg, respectively.The shock model and shock-free model also differ in the duration of the ejecta.

Figure 3 .
Figure 3.The solid line with black dots is the observed light curve of AT 2019zhd (Pastorello et al. 2021a), and the red and blue solid lines are the light curves of the shock and the shock-free models, respectively.

Figure 5 .
Figure 5.The ρ-T diagram of H and He are divided into 6 zones to improve the computational efficiency.The pure states are colored with blue, and the rest are mixture states that require solving Saha equations.

Figure 6 .
Figure 6.The first panel shows the combined specific heat capacity of X = 0.74 and Y = 0.26, and the bottom panel shows the adiabatic index.

Figure 8 .
Figure 8. From the left to the right, the first row shows the time evolution of the mass fraction of He, He + , and He 2+ , the second row shows the time evolution of the mass fraction of H2, H, and H + .

Figure 9 .
Figure 9.The upper and lower panels show the light curves of H and He mixture EoS and perfect gas EoS with various γ of the shock (s) and shock-free (sf) models, respectively.

Figure 10 .
Figure 10.The black line with dots shows the light curve of AT2019zhd.The red and blue colors indicate the shock and shock-free (sf) models.The solid and dashed lines indicate our original models and their counterparts when we turn off the radiation force.

Figure 11 .
Figure 11.Upper panel: light curves of the last 10 days of shock (s) and shock-free (sf) models and their corresponding dust deficit models.Bottom panel: the time evolution of log 10 (κRFr/(cg)) and log 10 κR of the shock dust deficit model.

Figure 12 .
Figure 12.The black line with dots shows the light curve of AT2019zhd.The red and blue colors indicate the shock (s) and shock-free (sf) models.The solid and dashed lines indicate the original and low-resolution models, respectively.

Table 1 .
From left to right, model name, ejecta's velocity factor, ejecta duration tej, mass loss rate Ṁ , cumulative mass loss ∆M , ejecta's density, α, the temperature of the ejecta, the ratio of Er and eg, and the total energy of the ejecta of the seven simple models.

Table 2 .
) The parameters that are different in shock (s) and shock-free (sf) models.In particular, δv > 0 cases generally belong to shock models and δv < 0 cases belong to shock-free models.