© The Authors 2025

1. Introduction

The long-standing basic picture of CE evolution consists of a dynamical-timescale, adiabatic inspiral also known as the ‘plunge-in’¹. This key assumption underlies the energy formalism (Van den Heuvel et al. 1976; Tutukov et al. 1979; Webbink 1984) that equates the orbital energy released during the inspiral with the envelope binding energy. In line with this assumption, almost all 3D hydrodynamic simulations of the CE plunge-in have neglected radiation hydrodynamics. Until recent years, these simulations have mainly targeted lower-mass donor stars on the red giant and asymptotic giant branches (e.g. Passy et al. 2012; Ricker & Taam 2012; Ohlmann et al. 2016; Chamandy et al. 2018; Reichardt et al. 2019; Prust & Chang 2019; Sand et al. 2020; González-Bolívar et al. 2022; Bermúdez-Bustamante et al. 2024). In such cases, radiative losses and thermal restructuring are assumed to be negligible because the donor envelope’s thermal timescale is many orders of magnitude longer than the inspiral duration.

However, the assumption of adiabaticity is questionable for very massive red supergiant (RSG) donors, which have short thermal timescales (Ricker et al. 2019; Vigna-Gómez et al. 2022). More broadly, the evolution leading to the plunge-in may be non-adiabatic (Pejcha et al. 2016a), while the plunge-in itself may be followed by a slower, self-regulated inspiral phase (Meyer & Meyer-Hofmeister 1979; Podsiadlowski 2001), during which the luminosity from the orbital decay matches the radiated luminosity at the photosphere. These considerations call for incorporating radiation transport into 3D simulations, especially for massive donor stars.

However, this increases the computational cost of an already challenging calculation. Simulations of CE evolution with RSG donors have only become available in recent years (Lau et al. 2022a, b; Moreno et al. 2022; Vetter et al. 2024, 2025; Landri et al. 2025). All of these simulations are adiabatic, except for the preliminary results of a radiation hydrodynamics CE simulation by Ricker et al. (2019), who modelled a 82.1 M_⊙ RSG donor and a 35 M_⊙ companion black hole. More recently, Hatfull & Ivanova (2025) presented hydrodynamic simulations with radiative emission and diffusion used to model the light curve of the V1309 Sco luminous red nova (LRN) (Tylenda et al. 2011), which is believed to be caused by the merger of a slightly evolved, ≈1.5 M_⊙ star with a much lower-mass main-sequence companion. Other calculations incorporating energy transport were performed under 1D spherical symmetry (Meyer & Meyer-Hofmeister 1979; Clayton et al. 2017; Fragos et al. 2019; Bronner et al. 2024).

As a result, our understanding of the impact of radiation transport on CE evolution, especially in CEs with massive-star donors, remains limited. Clarifying its role in the outcome of CE interactions is essential for understanding the formation of compact binary systems. In massive binary stars, this includes gravitational-wave sources, X-ray binaries, and progenitors of stripped-envelope supernovae. Simulations with radiation hydrodynamics are also critical for bridging theoretical models of CEs with observations. Energy transport is needed to accurately model photospheric evolution, which is necessary for computing light curves that can be compared with observed LRNe (Pastorello et al. 2019; Hatfull & Ivanova 2025). Radiation transport is also expected to influence the morphology and evolution of post-CE bound ejecta and bipolar lobes (e.g. Vetter et al. 2024, 2025; Gagnier & Pejcha 2023), which likely contain optically thin material. Bipolar outflows in CEs are believed to contribute to the formation of planetary nebulae (e.g. De Marco 2009; García-Segura et al. 2018; Zou et al. 2020; Jacoby et al. 2021; Ondratschek et al. 2022).

In this paper, we present 3D radiation hydrodynamics simulations of the CE phase involving a 12 M_⊙ RSG donor and a 3 M_⊙ companion. Radiation hydrodynamics is modelled using the flux-limited diffusion approximation. We compare these simulations with our previous adiabatic simulation, which used an otherwise identical setup, and find that radiation transport significantly obstructs envelope ejection, resulting in less than half the unbound mass compared to the adiabatic case. Notably, almost no material becomes unbound during the plunge-in phase itself and the sustained mass ejection previously seen during the slow spiral-in is strongly suppressed by radiative diffusion. However, we did not account for the release of recombination energy, which is known to significantly increase the amount of unbound envelope mass in adiabatic simulations.

In Sect. 2, we describe our implementation of radiative diffusion and the initial setup. We present our main results in Sect. 3, showing the effect of radiation transport on ejecta morphology (3.1), the amount of unbound ejecta (3.2), and the final orbital separation (3.3). To better understand these effects, in Sect. 3.4, we identify and analyse regions where radiation transport is dynamically important. In Sect. 3.5, we specifically examine the optical depth at the location of partially ionised hydrogen. In Sect. 4, we discuss the limitations of our method (4.1), compare our results with related studies (4.2), and discuss implications of our findings for binary stellar evolution (4.3). Finally, we summarise our findings in Sect. 5.

2. Methods

2.1. Overview

The initial binary system consists of a 12 M_⊙ RSG donor and a 3 M_⊙ companion modelled as a point mass, which could represent either a main-sequence star or a black hole. As in our past simulations (Lau et al. 2022a, b, hereinafter referred to as Papers I and II, respectively), we used the smoothed particle hydrodynamics (SPH) code PHANTOM (Price et al. 2018) with an implicit two-temperature flux-limited diffusion solver based on Whitehouse & Bate (2004) and Whitehouse et al. (2005) implemented for this work. Details of the method are explained in Sect. 2.2. The initial hydrostatic 3D stellar model is identical to that in Paper I and is resolved with 2×10⁶ SPH particles. This enables a controlled comparison to infer the effects of radiation transport. We also performed an additional simulation in which the base of the RSG envelope is heated to account for luminosity from the helium core and hydrogen-burning shell. This heating primarily offsets the slow increase in central density due to radiation transport and has an insignificant energy contribution over the course of the CE phase (see Sect. 2.6).

As in Paper I, convection is initially suppressed in the RSG envelope by using a stellar model with a flat entropy stratification. However, this setup does not prohibit convection from developing later, as shock heating and radiative diffusion modify the entropy structure during the CE plunge-in. Convective flows in the bound ejecta have been observed in Papers I and II as well as in other studies (e.g. Ohlmann et al. 2017; Vetter et al. 2024; Gagnier & Pejcha 2023).

2.2. Radiation hydrodynamics scheme

The implementation of implicit flux-limited diffusion in PHANTOM was directly ported from the SPHNG code (Benz et al. 1990; Bate et al. 1995; Bate & Keto 2015). It follows the method developed by Whitehouse & Bate (2004) and Whitehouse et al. (2005), implemented in SPHNG as described in Whitehouse & Bate (2006). The radiation hydrodynamics equations we solve are given by (Mihalas & Mihalas 1984; Turner & Stone 2001)

$$\frac{D\rho }{\mathit{Dt}}=-\rho \nabla ·\text{v},$$(1)

$$\rho \frac{D\text{v}}{\mathit{Dt}}=-\nabla p+\frac{\kappa \rho }{c}{\text{F}}_{\mathrm{rad}}+\rho {\mathrm{\Pi }}_{\mathrm{shock}}-\rho \nabla \mathrm{\Phi },$$(2)

$$\rho \frac{D\xi }{\mathit{Dt}}=-\nabla ·{\text{F}}_{\mathrm{rad}}-\nabla \text{v}:{\text{P}}_{\mathrm{rad}}+a_{\mathrm{rad}}c\kappa \rho (T_{\mathrm{gas}}^4-T_{\mathrm{rad}}^4),$$(3)

$$\rho \frac{\mathit{Du}}{\mathit{Dt}}=-p\nabla ·\text{v}+\rho {\mathrm{\Lambda }}_{\mathrm{shock}}-a_{\mathrm{rad}}c\kappa \rho (T_{\mathrm{gas}}^4-T_{\mathrm{rad}}^4),$$(4)

where D/Dt:=∂/∂t+v·∇ is the convective derivative and ρ, v, p, P_rad, F_rad, κ, Φ, ξ, u, T_gas, and T_rad represent density, velocity, scalar isotropic pressure, radiation pressure tensor, radiation flux, opacity, gravitational potential, specific radiation energy, specific gas thermal energy, gas temperature, and radiation temperature, respectively. The radiation constant is a_rad and the speed of light is c. Equations (2)–(4) have been integrated over frequency, and we do not distinguish between the Planck mean, energy mean, and flux mean opacities, all of which are denoted by κ. The symbols Π_shock and Λ_shock are shock dissipation terms (Price 2012). The latter usually includes a thermal conductivity term to correctly treat contact discontinuities (Price 2008). However, since we explicitly model physical energy conduction in the form of radiative diffusion, we disabled this term, so that Λ_shock only includes viscous shock heating. The final terms in Eqs. (3) and (4) describe energy exchange via radiative absorption and emission. These terms use the expression for specific radiation energy, $\xi =a_{\mathrm{rad}}{T}_{\mathrm{rad}}^4/\rho$. We assume local thermodynamic equilibrium (LTE), meaning the gas is in thermal equilibrium with itself even if it is decoupled from the radiation field, allowing radiative emission to be written using the Planck function (e.g. Rybicki & Lightman 1986). In the optically thick limit, there is complete thermodynamic equilibrium and T_gas=T_rad. In the optically thin limit, the gas and radiation can decouple, causing them to have different temperatures.

Equations (1)–(4) are closed by (i) an equation of state (EoS), (ii) an opacity law κ(ρ,T_gas), (iii) an expression for the radiation flux, F_rad, and (iv) an expression for the radiation pressure tensor, P_rad. We adopt the EoS for an ideal gas with a fixed adiabatic index of γ = 5/3,

$$p=(\gamma -1)\rho u.$$(5)

The specific thermal energy is

$$u=\frac{\mathcal{R}{T}_{\mathrm{gas}}}{(\gamma -1)\mu },$$(6)

where ℛ is the gas constant and μ is the mean molecular weight, for which we adopt a constant value of μ = 0.61821. We do not include the energy released by atomic and molecular recombination² and the associated changes to μ and γ.

We used opacity tables from the MESA stellar evolution code, as implemented in PHANTOM by Reichardt et al. (2020). These tables combine OPAL opacities for the high-temperature regime (Rogers et al. 1996; Rogers & Nayfonov 2002) with SCVH opacities for the low-temperature regime (Saumon et al. 1995). They capture the opacity bumps associated with partially ionised hydrogen, helium, and iron, as well as molecular transitions occurring below ≈5000 K.

The radiative flux has the standard expression,

$${\text{F}}_{\mathrm{rad}}=-D\nabla E,$$(7)

where $E=\rho \xi =a_{\mathrm{rad}}{T}_{\mathrm{rad}}^4$ is the volumetric radiation energy density. Radiation energy is transported (anti-)parallel to the temperature gradient ∇T_rad. The diffusion coefficient is given by D=cλ/(κρ), where λ is the flux limiter. We adopt the flux limiter from Levermore & Pomraning (1981), which approaches λ→1/3 in the optically thick limit. In the optically thin limit, the flux limiter ensures that |F_rad|→cE, so that radiation energy does not propagate faster than light.

The radiation pressure tensor is written as P_rad=Ef, where f is the Eddington tensor, whose components depend on the flux limiter λ. In the isotropic, optically thick limit, P_rad→E/3I. In the optically thin limit, ${\text{P}}_{\mathrm{rad}}\to E\hat{n}\hat{n}$, where $\hat{n}=\nabla E/|\nabla E|$ is a unit vector aligned with the energy density gradient.

2.3. Implicit solver for radiative diffusion

Although PHANTOM previously included an explicit solver for two-temperature radiative diffusion (Biriukov & Price 2019), an explicit scheme is only stable subject to a Courant condition on the maximum time step for flux diffusion, given by Δt_rad=Ch²ρκ/(cλ) (Whitehouse & Bate 2004), where h is the particle smoothing length and C is a Courant number of order unity. It is a known issue that Δt_rad is often many orders of magnitude shorter than the Courant time step for hydrodynamic processes in optically thin regions, where Δt_rad approaches the light-crossing time of the resolution scale. These conditions are first encountered when the companion approaches the donor star’s surface, raising tides and ejecting material in its wake. As a result, using an explicit solver for radiative diffusion in CE evolution results in prohibitively small time steps, making long-term simulations computationally infeasible.

We address this problem by implementing an implicit solver for u and ξ, based on Whitehouse et al. (2005), that we directly ported from SPHNG. This solver remains stable even with arbitrarily large time steps. It uses the backwards Euler method to evolve the gas and radiation energy equations, solving them iteratively for each pair of SPH neighbours using the Gauss-Seidel method until achieving the desired accuracy. Further details can be found in Whitehouse et al. (2005) and Bate & Keto (2015). To avert a large number of random memory accesses, the implicit method requires building large arrays that cache interaction terms between SPH neighbours. This significantly increases the memory requirement compared to the explicit scheme.

We set an accuracy tolerance of 10⁻⁶ for the solutions of u and ξ, which we find to be a good compromise between accuracy and computational speed. Optically thick regions typically converge to this accuracy within a few iterations, whereas optically thin regions typically require considerably more iterations. As a result, computations become more expensive once optically thin ejecta are produced. We imposed a maximum of 250 iterations per time step and allowed the simulation to proceed even if the u and ξ solutions had not fully converged within the tolerance. To further improve computational efficiency, we optimised each iteration by only including SPH particles that failed to converge in the previous iteration. This focuses most of the computational expense on particles in optically thin regions, which require the most iterations. We find that this optimisation does not lead to significantly different results in lower-resolution runs with 50 000 particles. In our simulations, most time steps converge within 10–20 iterations (see Sect. 3). In Appendix C, we present the results of radiation hydrodynamics test problems from Whitehouse et al. (2005).

2.4. Initial setup and stellar model

The initial conditions of our simulations are identical to those of the ‘gas + radiation’ EoS simulation in Paper I. This was ensured by starting from the same simulation file. All material is therefore initially in thermodynamic equilibrium, with T_gas=T_rad. Below, we summarise key information about our 3D stellar model and refer the reader to Paper I for further details.

The donor star is a non-rotating 12 M_⊙ RSG with a radius of 618 R_⊙. We do not drive equilibrium envelope convection in the pre-CE RSG, and so the envelope was constructed to have constant entropy for convective stability. The inner 3.84 M_⊙ of the RSG, comprising the helium core and radiative intershell, is replaced by a point mass that only interacts gravitationally with surrounding material. Its gravitational potential is softened using a cubic spline with a softening length of h₁ = 9.25 R_⊙, transitioning to the exact Newtonian potential at r_core = 2h₁. This boundary corresponds to the inner edge of the convective envelope in a MESA stellar model evolved to core helium ignition. The 3 M_⊙ companion, also modelled as a softened point mass, has a softening length of h₂ = 2.15 R_⊙. The stars are initially placed in a circular orbit with a separation of 988 R_⊙, where the RSG radius exceeds the Roche radius by 25%.

The 1D stellar profile was mapped to a 3D distribution of SPH particles using stretch mapping and an asynchronous particle shifting technique that allows different regions of the star to relax on their local dynamical timescale (Appendix A, Lau et al. 2022a). Paper I demonstrates that the 3D density structure is preserved to within ≲0.1% when evolved in isolation for 90 times the surface free-fall time, even with just 50 000 SPH particles (40 times fewer than in the production simulations).

2.5. Thermal structure of the donor

In a radiation hydrodynamics simulation, the donor star’s initial thermal structure becomes an important consideration in addition to its mechanical structure. In thermal equilibrium, the total luminosity remains constant throughout the envelope of a giant star. Convection carries the luminosity generated through core and shell burning to the photosphere, where it is radiated away with the same luminosity. This is illustrated in Fig. 1, which shows the total, radiative, and convective luminosities of the reference MESA RSG model at the onset of core helium burning.

Fig. 1. Comparison of the initial radiative luminosity profile (lrad) with that of a MESA RSG reference model. Each black dot shows one of two million SPH particles. For the MESA model, the blue line shows the total luminosity, the orange line shows the convective luminosity, and the green line shows the radiative luminosity. The dashed vertical line divides the convective envelope from the radiative interior.

Our default simulation does not explicitly model cooling from the photosphere or heating from nuclear burning (though we conduct an additional simulation that includes nuclear luminosity; see Sect. 2.6). To remain consistent with this approach, we initialised the envelope’s entropy profile to be convectively stable, as described in Sect. 2.4. At first order, the additional energy injected by the companion during the CE inspiral must be transported via increased radiative or convective luminosity. We emphasise that our CE simulations allow additional convection to develop self-consistently, as observed in Papers I and II.

With this approach, the start of the simulation only includes radiative diffusion as a means of energy transport. Figure 1 shows the radiative luminosity of the 3D stellar model, l_rad(r) = 4πr²F_rad,r, where F_rad,r is the radial component of Eq. (7). Each black dot represents a single SPH particle. The profile is qualitatively consistent with the reference MESA model (shown in green), providing a reasonable representation of an RSG’s radiative luminosity at the onset of core helium burning. The two profiles do not match exactly because the 3D model is based on the flat-entropy profile we constructed in Paper I, which only has the same total mass and radius as the MESA model. However, this structure is not in thermal equilibrium because entropy flows into or out of regions where ∂l_rad/∂r is negative or positive, respectively. Specifically, there are three dips in the radiative luminosity profile that coincide with opacity bumps from partially ionised iron, HeII, and HeI (ordered from smaller to larger radii) in the OPAL opacity tables. Also, below r_core = 18.5 R_⊙ (marked by the dashed vertical line), the softened density profile results in a much smaller radiative luminosity compared to the MESA model, which maintains a flat l_rad profile down to the hydrogen-burning shell. Consequently, the star is expected to undergo thermal adjustments at these opacity bumps and below r_core on the local radiative diffusion timescale.

Indeed, we find that when evolving the 3D RSG without a companion, entropy accumulates at the base of the convective envelope. This is associated with a gradual increase in central density. At t = 4.2 yr, corresponding to 15 times the surface free-fall time, a buoyant plume of material begins to rise, marking the onset of convection despite an initially flat entropy profile. By this point, the central density has increased by 25%. However, this rising plume is not expected to introduce major artefacts to our initial conditions, since by 4.2 yr, the companion in our CE simulation has already plunged deep into the RSG envelope (see Sect. 3). Also, the plume has a minor impact on the envelope density and pressure, as convection is highly efficient in most of the envelope.

2.6. Nuclear burning luminosity

To assess the potential impact of the slow rise in central density, we conducted an additional simulation where we heated the base of the RSG envelope to compensate for the energy transported outwards. In a real RSG, this luminosity is generated by core helium burning and hydrogen shell burning. We applied a luminosity of L_★ = 1.7×10³⁸ erg s⁻¹, matching that of the MESA model, to gas particles within a radius r_core of the point mass representing the donor’s core. This heating was applied uniformly in mass, meaning each SPH particle’s thermal energy³ increased at the rate,

$$\frac{{\mathit{du}}_i}{\mathit{dt}}=\frac{L_{★}}{m(r_{\mathrm{core}})},$$(8)

where m(r_core) is the total gas mass interior to r_core. We find that the central density of the heated RSG only increases by 6% over the same 4.2 yr duration. The total energy injected by this luminosity over the course of the CE inspiral is only a few percent of the envelope’s initial gravitational binding energy, E_bind,

$$\frac{L_{★}\mathrm{\Delta }{t}_{\mathrm{inspiral}}}{E_{\mathrm{bind}}}=0.022(\frac{L_{★}}{1.7\times {10}^{38}\mathrm{erg}{\mathrm{s}}^{-1}})(\frac{\mathrm{\Delta }{t}_{\mathrm{inspiral}}}{5\mathrm{yr}})$$(9)

$$(\frac{E_{\mathrm{bind}}}{1.2\times {10}^{48}\mathrm{erg}}{)}^{-1}.$$

Here, we assumed a duration of Δt_inspiral = 5 yr, roughly the time at which the rapid plunge-in terminates (see Sect. 3). As such, we do not expect significant qualitative differences between the outcomes of the heated and non-heated simulations. This is consistent with Fragos et al. (2019), who performed a 1D CE simulation with a 12 M_⊙ RSG and found that the cumulative nuclear energy contributes only a few percent of the change in gravitational potential energy. The ‘heated’ simulation primarily serves as a control, where the central density and entropy of the envelope are better maintained. We do not consider the heated simulation to be our default calculation, as our setup also suppresses the equilibrium convection that would normally transport the nuclear luminosity, resulting in an inconsistency. For future CE simulations, it would be valuable to explore the effects of a more realistic setup where initial envelope convection is driven via central heating and cooling from a resolved photosphere.

3. Results

We report the results of three different CE simulations: (i) radiative (radiative diffusion without central heating), (ii) radiative + heating (radiative diffusion with central heating), and (iii) adiabatic (the ‘gas + radiation’ EoS simulation in Paper I). The first two cases are also collectively referred to as the radiative simulations. Below, we provide an overview of the simulations that is common to all three cases.

Figures 2 and 3 display density slices from the three simulations at different stages of the CE inspiral. Row 1 shows the initial setup, where the orbital separation is 988 R_⊙ and the RSG overfills its Roche lobe by 25%. Initially, the implicit radiation solver converges to the specified tolerance within a few to ten iterations. After one orbit (row 2), the companion tidally deforms the RSG and ejects material in its wake. As optically thin ejecta are produced, the solver requires around 30 iterations to converge. The companion enters the RSG envelope at around 3 yr (row 3) and spirals inwards rapidly due to gravitational drag. Row 4 shows the spiral shocks generated by the motion of the two stellar cores in the CE, causing some of the envelope material to expand and eventually be ejected. The orbital separation stabilises as the surrounding envelope material drops in density and reaches a high degree of co-rotation with the stellar cores. The inspiral timescale may be defined as

Fig. 2. Density in the initial orbital plane (z = 0) in 4320-by-3000 R⊙ panels. The crosses mark the locations of the donor’s core and the companion. Each column corresponds to a different simulation while each row corresponds to a different stage in the evolution. Row 5 shows the slow inspiral taking place several years after the plunge-in. The inspiral timescale, $t_{\text{inspiral}}=a/(-\dot{a}$), is several hundred to a thousand times longer than the orbital period. This figure was made with the SARRACEN package (Harris & Tricco 2023). Online movies of our simulations are listed in Table D.1.

Fig. 3. Same as Fig. 2, except each 4320-by-3000 R⊙ panel shows a slice viewed edge-on that contains the stellar cores. Online movies of our simulations are listed in Table D.1.

$$t_{\mathrm{inspiral}}=\frac{a}{-\dot{a}},$$(10)

where a is the separation between the donor’s core and the companion. This timescale becomes several orders of magnitude longer than the orbital period (P_orb/t_inspiral≲10⁻³−10⁻²). Most of the envelope remains bound at this stage, forming an extended, diffuse structure (row 5). During this phase, some time steps required an even larger number of iterations to reach convergence, and ∼0.1% of time steps did not converge after 250 iterations, the maximum number allowed. This represents a very small fraction of time steps and likely only reflects non-convergence in distant, optically thin particles, where the diffusion approximation breaks down in any case. During the plunge-in, a small number of distant, unbound SPH particles (at least 30 000 R_⊙ away from the centre of mass) were removed because these low-density (ρ≲10⁻¹⁵ g cm⁻³) particles had smoothing kernels that overlapped with a significant fraction of all particles in the simulation, causing an unbalanced load in the density calculation. This removal resulted in negligible mass loss, less than 4×10⁻³% of the total.

Energy non-conservation⁴ exceeded 10% near t = 8 yr for the radiative simulations (see Appendix A), and the simulations were terminated shortly after. In contrast, the adiabatic simulation conserved total energy to within 0.05% up to t = 11.6 yr. The much larger energy accumulation in the radiative cases occurs because the implicit method does not conserve energy exactly.

3.1. Ejecta morphology

The addition of radiation transport leads to the formation of more condensed density structures during and after the plunge-in. For instance, rows 4 and 5 of Figs. 2 and 3 show that the outer regions of the CE form dense outer shells in the radiative simulations. In the adiabatic simulation, thermal energy deposited by the bow shock is fully utilised to expand the downstream envelope gas. However, in the non-adiabatic cases, radiative diffusion can remove thermal pressure, reducing the work done by the gas and even causing compression.

The last rows of Figs. 2 and 3 show the ejecta morphology after the plunge-in. In all three simulations, the equatorial region is primarily occupied by bound turbulent gas, while the regions above and below the orbital plane are more evacuated. In the adiabatic case, the bound equatorial material exhibits finer and more strongly stratified flows. This is consistent with the expectation that any energy conduction or diffusion (in this case, radiative diffusion) smooths out density and temperature contrasts. Only larger-scale turbulent flows, which are less strongly stratified, are preserved in the radiative simulations⁵.

In fact, in the radiative simulations, the dense outer shells of the inflated CE are susceptible to the Rayleigh-Taylor instability. We observe the formation of inward-extending fingers that later evolve into the characteristic mushroom-shaped downflows near the end of the plunge-in, continuing to the end of the simulation. However, these flow patterns are disrupted by outwardly moving spiral shocks and become mixed with the bound turbulent gas. This is illustrated in the last rows of Figs. 2 and 3⁶.

To show the full extent of the CE ejecta near the end of the simulations, we present Fig. 4, showing face-on and edge-on density slices on a scale that is more than four times larger and with a colour map that displays lower densities. In the adiabatic case, low-density bipolar outflows are launched at t≈6 yr and expand homologously at speeds of ≈200 km s⁻¹. These outflows are hydrodynamically collimated by the equatorially concentrated material that is still bound, inflating bipolar lobes (a phenomenon also described in Paper II). They lead to a persistent increase in the amount of unbound envelope mass following the plunge-in. Bipolar outflows in CE evolution have also been recently discussed by Vetter et al. (2024, 2025), who performed magnetohydrodynamic CE simulations resulting in faster and more strongly collimated jets. In their case, the jetted expulsion only accounts for ∼10% of the total ejected envelope mass. They similarly report bound material that is equatorially concentrated in a centrifugally supported structure.

Fig. 4. Same as Figs. 2 and 3, but showing the late-time ejecta (t = 7.71 yr) with a larger scale (18 200-by-14 000 R⊙) and a lower-density colour map. Top row: Orbital-plane slices. Bottom row: Edge-on slices.

In the ‘radiative’ and ‘radiative + heating’ simulations, poloidal outflows are also observed a few years after the plunge-in. However, these outflows are weaker and less distinctly separated from the equatorial ejecta. Figure 5 shows the angular distribution of newly unbound ejecta in all three simulations as a function of time. In the ‘radiative’ simulation, two downward plumes of energetically unbound material are ejected from the central binary starting at t≈6.4 yr. The single-sided and episodic nature of this ejection reflects obstruction of the outflow by surrounding envelope gas that has fallen back. Similar outflows were seen in the simulations of the post-dynamical inspiral evolution by Gagnier & Pejcha (2025). As we show in Sect. 3.2, the radiative simulations contain more fallback material because less of the CE is ejected. Furthermore, radiation transport is able to concentrate this material by removing pressure support. In the ‘radiative + heating’ simulation, an upward plume and a weaker downward plume emerge instead at a slightly later time of t≈7.4 yr. These differences are not expected to be systematic effects of central heating, as the exact launch time and direction of a weak outflow are sensitive to the surrounding turbulent ejecta, particularly the equatorial material that primarily contributes to collimation.

Fig. 5. Angular distribution of particles that become unbound between consecutive simulation snapshots, meaning the sums of their kinetic and potential energies become positive. The angle θ is the polar angle measured from $\hat{\text{z}}$, which is the direction of the initial orbital angular momentum vector. Material is initially ejected near the orbital plane, then closer to the poles in single or double-sided outflows after the plunge-in.

3.2. Unbound mass

Figure 6 shows the evolution in the fraction of unbound envelope mass. Material is considered unbound if the sum of its kinetic and potential energy is positive. For comparison, we also plot results assuming an alternative criterion that includes gas and radiation energy (dashed lines).

Fig. 6. Fraction of unbound envelope mass as a function of time for the adiabatic simulation (red), the ‘radiative’ simulation (blue), and the ‘radiative + heating’ simulation (green). Solid lines assume material is unbound when the total mechanical energy is positive (ek+ep>0), while dashed lines are less restrictive, considering material to be unbound when the sum of the mechanical and internal energy is positive (ek+ep+egas+erad>0).

As with past 3D simulations of CE evolution that do not include recombination energy, only a small fraction of the total envelope is unbound within a few years to a decade after the plunge-in, which is often the end of such simulations. At t = 9 yr, the adiabatic simulation records only 20.8% (1.7 M_⊙) of the RSG envelope as unbound, increasing at a rate of 0.20 M_⊙ yr⁻¹ (see Paper I). Including radiative diffusion reduces the amount of unbound ejecta further. By a similar time, only 11.0% (0.90 M_⊙) of the envelope is ejected in the ‘radiative’ simulation, while 8.9% (0.72 M_⊙) of the envelope is ejected in the ‘radiative + heating’ simulation. In both radiative cases, the amount of unbound mass is less than half of that in the adiabatic case. However, note that the radiative simulations display ongoing mass ejection at late times that is linked to the intermittent plumes discussed in Sect. 3.1.

In the adiabatic simulation, we identify three phases of mass ejection in Fig. 6. The first phase occurs prior to the plunge-in (t≲3 yr, during which material is ejected in the companion’s wake. In this phase, the evolution in the fraction of unbound envelope mass is nearly identical across all three simulations. In each case, 4–5% of the envelope mass (0.32–0.41 M_⊙) becomes unbound, with the ejection rate rising as the companion approaches the RSG surface. This growth is abruptly halted when the companion enters the distorted RSG envelope, cutting off the outflow of high-velocity ejecta through the L₂ point.

A second ejection event occurs at t≈4 yr in the adiabatic simulation, associated with energy injected by spiral shocks during the CE plunge-in. As shown in Fig. 5a, the ejecta span a much larger angular range due to the CE homogenising energy deposited deep inside. However, the total amount of unbound ejecta remains a small fraction (<10%) of the total envelope mass, consistent with previous adiabatic CE simulations⁷. This finding has long contrasted with the widely adopted picture of dynamical CE ejection, and the radiative simulations exacerbate this disagreement because this second ejection event is almost entirely absent. Instead, the ejection rate remains small (a few times 10⁻² M_⊙ yr⁻¹) for several years after the plunge-in until intermittent outflows set in. Shortly after the plunge-in begins, the total unbound envelope mass in the radiative simulation is ≈1% higher in absolute fraction than in the centrally heated simulation. As shown in Eq. (9), this variation is unlikely to be a direct effect of energy injection. Instead, it likely results from the difference in the RSG’s ability to maintain its central density depending on whether or not luminosity is injected in the centre, as discussed in Sects. 2.5 and 2.6.

In the adiabatic simulation, a third phase of mass ejection proceeds from plunge-in, beginning around t≈6 yr. Material is ejected through the bipolar outflows described in Sect. 3.1. A similar process occurs in the simulations with radiation transport, but a clear bipolar morphology is absent. Instead, ejection happens through single-sided and intermittent plumes. The onset of this ejection is delayed by 2–3 years compared to the adiabatic case, resulting in a much lower total fraction of unbound envelope mass at the same time. Due to energy non-conservation at late times, we have not followed the simulations long enough to determine whether this form of envelope ejection can be sustained. The episodic nature of these outflows is reflected in the uneven rise in the amount of unbound mass, particularly in the ‘radiative’ case (blue curve). In the radiative simulations, the amount of unbound ejecta computed including thermal and radiation energies (dashed lines) varies irregularly with time. Dips appear at certain points, which indicate instances where ejected material has lost energy upon becoming optically thin and is subsequently considered energetically bound again.

3.3. Final separation

Figure 7 shows the evolution of the separation between the RSG core and the companion. Despite differences in the amount of ejected mass, radiative diffusion has little impact on the orbital separation up to the end of the plunge-in phase. The orbital separations remain within ≈10% of each other until t = 4 yr, even though noticeable differences in the amount of unbound mass emerge around t≈3.8 yr. By t = 8 yr, near the end of the radiative simulations, the orbital semi-major axes in the adiabatic, ‘radiative’, and ‘radiative + heating’ simulations are 40.9, 39.4, and 42.1 R_⊙, respectively. These differences are minor and are comparable to the orbital eccentricities of ≈0.05 recorded at the same time. In contrast, our previous works found that variations in envelope thermal energy (Paper I) and different sources of recombination energy (Paper II) resulted in more discernable changes to the semi-major axis near the end of the plunge-in.

Fig. 7. Separation between the RSG core and the companion as a function of time for the adiabatic simulation (red), the ‘radiative’ simulation (blue), and the ‘radiative + heating’ simulation (green).

At the end of the simulations, the semi-major axis continues to decrease at a rate of a few times 0.01–0.1% per orbit. In Papers I & II, we recorded the post-plunge-in orbital separation when this rate dropped to 0.05% per orbit (P_orb/t_inspiral<5×10⁻⁴). However, in the radiative simulations, the rate of spiral-in does not slow steadily, making it unclear whether the orbital semi-major axis is approaching a stable value. In particular, the ‘radiative’ simulation shows a pronounced decrease in orbital separation around t≈8.5 yr, coinciding with a relatively steep increase in unbound mass. The cause of this decrease is unclear and is plagued by numerical difficulties associated with simulating the late stages of the slow spiral-in, including a lack of spatial resolution, the softening length of the point masses becoming a large fraction of their Roche radii, and a potentially problematic buildup of (non-conserved) total energy, reaching ≈10% at this point. Given these issues, it would be more instructive to study this behaviour in simulations targeting the post-plunge-in evolution (e.g. Gagnier & Pejcha 2025), if it indeed occurs.

Despite these uncertainties, it is still possible to conclude that the orbital separation during and shortly after the plunge-in is not significantly affected by the inclusion of radiative diffusion. This is likely because the final separation and orbital evolution inside the CE are primarily influenced by deeper, optically thick layers of the envelope, where radiative diffusion is not dynamically important. This aligns with analogous behaviour studied in Paper II, where we found that energy released at large radii from hydrogen recombination did not affect the post-plunge-in orbital separation despite contributing to the ejection of substantially more mass, whereas helium recombination energy, released much deeper, increased the final separation by ≈16%.

3.4. Dynamical effect of radiation transport

To better understand the differences between the adiabatic and radiative simulations, we identify regions within the CE and ejecta where radiation transport is efficient. Figure 8 shows various quantities at three representative snapshots in the ‘radiative’ simulation, including density (ρ), a local diffusion timescale (t_diff), specific energy (e_kp, the sum of kinetic and potential energy), opacity (κ), and hydrogen ionisation fraction ([HII]). The e_kp plots depict bound material in blue, unbound material in red, and the transition region in white. Similarly, the [HII] plots show ionised hydrogen as green, neutral hydrogen as pink, and the narrow partial ionisation zone as white. The local diffusion timescale is defined by

Fig. 8. Top to bottom: Various quantities in orbital-plane slices from the ‘radiative’ simulation: density (ρ), local diffusion time (tdiff, Eq. (11)), kinetic plus potential energy per unit mass (ekp), opacity (κ), and the density-weighted hydrogen ionisation fraction ([HII]). Left to right: Snapshots representative of the pre-plunge-in phase, the end of the rapid plunge-in, and the late-time slow spiral-in. In rows (d) and (e), we overlay contours where the optical depth equals 1, 3, and 10, from largest to smallest size. The contours in row (d) include the contribution from all material, whereas those in row (e) include only the contribution from bound material (ekp>0). The fact that these contours enclose the hydrogen partial ionisation zone (marked by the sharp transition between green and pink) suggests that a significant fraction of hydrogen recombination energy could be absorbed by the bound material.

$$t_{\mathrm{diff}}:=\frac{\kappa \rho }{c\lambda }(\frac{E}{|\nabla E|}{)}^2.$$(11)

In optically thick regions, this is the timescale for energy diffusion across the scale height of radiation energy density. In the optically thin limit, the flux limiter, λ, varies such that t_diff→E/(c|∇E|), becoming an advection timescale. If t_diff is shorter than or comparable to the dynamical timescale of the CE ejecta, radiation transport may alter the local pressure gradient and produce a dynamical effect⁸.

In the very beginning, only the outermost ≈10% (3%) of the RSG envelope in mass has t_diff less than ten years (one year). The first column in Fig. 8 shows that, prior to the plunge-in, low-opacity and unbound material exits through the companion’s wake with t_diff≲0.01 yr. This material is directly thrown out by the companion, acting as a gravity slingshot (MacLeod et al. 2018). Even though radiation energy diffusion is efficient in this material, it does not directly influence the energy injection process. This explains why the amount of unbound ejecta closely matches between the adiabatic and radiative simulations until the companion enters the RSG envelope (Sect. 3.2). In Fig. 8b(i), although t_diff increases again farther from the binary, this actually indicates efficient radiation transport resulting in nearly uniform ejecta temperature. The large temperature scale length leads to a long t_diff.

The middle column shows a snapshot taken near the end of the plunge-in. The transition between energetically bound and unbound ejecta, shown as the white region in Fig. 8c(ii), occurs within optically thin material. Accordingly, the evolution of the amount of unbound mass begins to diverge from that of the adiabatic case, with little to no material ejected in the radiative simulations, as noted in Sect. 3.2. The outermost tidal arms that are launched near the beginning of the dynamical plunge-in exhibit t_diff values less than 10⁻³ yr despite still being energetically bound. This is consistent with our observation in Sect. 3.1 that these tidal arms form dense shells that are not seen in the adiabatic simulation. In contrast, most of the envelope has t_diff∼10 yr, which is long relative to the inspiral timescale. The spiral shock pattern is also visible in panel b(i), due to the steep change in E across a shock front.

The third column of Fig. 8 corresponds to approximately three years after the end of the plunge-in. By this time, nearly all material in the simulation has a t_diff that is comparable to or shorter than the dynamical timescale, indicating a major departure from an adiabatic calculation. This is consistent with the description of the unbound plumes in Sect. 3.1, which are smaller in extent and have denser outer shells compared to the bipolar lobes observed in the adiabatic simulation (Fig. 4). Radiation transport is therefore essential for accurately capturing the large-scale morphology of CE ejecta, which could provide circumstellar medium in interacting supernovae. If such behaviour extends to lower-mass systems, incorporating radiation transport may also be necessary to correctly model the morphology of planetary nebulae.

3.5. Location of partially ionised hydrogen

The last row of Fig. 8 shows the location of partially ionised hydrogen. Comparing rows 2, 4, and 5 reveals that the opacity drops sharply by a factor of ∼10³ across the hydrogen recombination front, resulting in a steep boundary beyond which t_diff decreases by a similar factor. The location of the hydrogen recombination front is significant because 3D adiabatic simulations of CE evolution show that recombination energy, especially hydrogen recombination energy (Lau et al. 2022b), drastically increases the fraction of unbound material and even leads to complete envelope ejection, assuming this energy does not escape (e.g. Nandez et al. 2015; Sand et al. 2020; Reichardt et al. 2020; Lau et al. 2022a, b; Vetter et al. 2024). While our simulations do not include recombination energy release, they lend support to this key assumption since we find that the partial ionisation zone is contained within the photosphere, meaning recombination energy could be thermalised in the surrounding neutral CE ejecta. To be converted into mechanical work, it further needs to drive expansion on a timescale that is shorter than the local diffusion timescale.

We show surfaces of constant optical depth, τ, in Fig. 8d. These contours correspond to τ = 1,3,10 in order of largest to smallest size. The optical depth is integrated along radial rays from infinity to the origin using the same opacity calculated for radiative diffusion. In Fig. 8e, the contours show optical depths that only include the contribution from bound material, τ_bound, to assess the ability for recombination radiation to assist with ejection. The τ_bound contours are larger and have different shapes compared to the τ contours in Fig. 8d, implying that the unbound ejecta at larger radii have non-negligible optical thicknesses and exhibit asymmetry⁹. In Fig. 8e (ii) & (iii), the ionisation zone is entirely contained within the τ_bound = 10 contour. The result that hydrogen recombines inside the photosphere agrees with Paper I and has also been reported in other studies (e.g. Sand et al. 2020; Bronner et al. 2024). It is important to verify whether the photosphere still contains the hydrogen partial ionisation zone when incorporating both radiation transport and recombination energy release. Recombination energy absorbed in optically thick regions drives additional expansion, which could cause the material to become optically thin again. Our previous adiabatic simulations that include recombination energy show that the hydrogen partial ionisation zone is much wider, of the order of thousands of solar radii, after the plunge-in (see Fig. 13 of Paper I). Simulations incorporating radiation transport and recombination energy release will also reveal the efficiency with which thermalised energy can be converted to mechanical work.

4. Discussion

4.1. Limitations

There are important limitations to our implementation of energy transport, particularly the neglect of pre-existing envelope convection, the lack of explicit cooling at the photosphere, and shortcomings of the flux-limited diffusion approximation. Pre-existing envelope convection was deliberately suppressed through adopting a flat entropy profile in the RSG envelope (Sect. 2.5). This choice reflects the assumption that such convection transports the RSG’s nuclear luminosity, which is not included in the default ‘radiative’ simulation. At first order, the energy injected during the spiral-in must be transported by additional convective flows. Pre-existing convection may be important if it interferes non-linearly with these flows, and if they are not disrupted by spiral shocks that develop during the plunge-in. Several studies have argued that additional transonic or supersonic convection could transport away orbital energy and hinder envelope ejection (e.g. Sabach et al. 2017; Wilson & Nordhaus 2020, 2022; Noughani et al. 2024). However, many of these works rely on simplifying assumptions such as mixing-length theory, spherical symmetry, and static stellar models. Their validity should be verified with self-consistent 3D simulations.

Additionally, the diffusion approximation could lead to unphysical effects when applied to an optically thin medium. The flux limiter ensures that the radiative flux approaches the correct asymptotic limit, |F_rad|→cE, but radiation energy can nonetheless only be transported anti-parallel to the temperature gradient, prohibiting the formation of shadows behind opaque obstacles. We also neglected photon scattering, assumed LTE, and adopted gray opacities. More sophisticated methods such as Monte Carlo radiative transfer (Whitney 2011), variable Eddington tensor schemes (Stone et al. 1992), and discrete ordinates methods (e.g. Jiang 2021; Ma et al. 2025) mitigate some of these problems, but are significantly more computationally expensive and challenging to implement. It would be valuable to verify our preliminary results with these approaches in the future.

Another challenge is accurately capturing cooling near the start of the simulation, which has two main difficulties. Firstly, resolving the steep temperature stratification at the RSG photosphere requires prohibitively high resolution. This long-standing issue is known to severely overestimate the effective temperature and luminosity in post-processed light curves from 3D CE simulations (Galaviz et al. 2017, but, see Hatfull et al. 2021 and Hatfull & Ivanova 2025). Secondly, the surface layer of SPH particles cannot radiate into the background ‘vacuum’ because the radiative diffusion method only allows energy exchange between particles. Modelling photospheric cooling is necessary for driving pre-existing convection in the RSG envelope. Potential solutions include implementing a separate treatment for surface cooling or introducing boundary particles to absorb energy radiated by the surface layer. However, these issues mainly affect the pre-plunge-in evolution in a CE simulation, as CE ejecta eventually provides material into which the photosphere can cool.

A final caveat is that recombination energy has not been modelled together with radiation transport. As discussed in Sect. 3.5, adiabatic simulations show that recombination energy enhances envelope ejection and increases the post-plunge-in orbital separation. In Paper II, we demonstrated that helium and hydrogen recombination energy increase the ejected mass by 48% and the final separation by 16%. Hydrogen recombination energy also leads to forming more extended and spherically symmetric ejecta. It is therefore important to check whether radiation transport interferes with recombination-assisted envelope ejection.

4.2. Comparison with related works

At the time of writing, few studies have incorporated radiative diffusion into detailed CE simulations and systematically analysed its effects. In a conference proceedings contribution, Ricker et al. (2019) presented preliminary results from a CE simulation involving an 82.1 M_⊙ RSG donor and a 35 M_⊙ black hole companion. This simulation, conducted using the 3D adaptive mesh refinement code FLASH with flux-limited diffusion, showed that vigorous envelope convection observed in an adiabatic setup is suppressed when radiative diffusion is included. They attribute this suppression to radiation carrying part of the energy flux that would otherwise be transported entirely by convection. In our setup, envelope convection does not develop before the plunge-in, and so we are unable to make the same comparison. However, as discussed in Sect. 3.1, we observe fewer small-scale flows in the turbulent ejecta after the plunge-in phase in the radiative simulation compared to the adiabatic case.

Recently, Hatfull & Ivanova (2025) implemented a ‘flux-limited emission-diffusion’ technique into the SPH code StarSmasher. This was applied to simulating the merger of V1309 Sco in 3D and computing its light curve. Unlike standard flux-limited diffusion, their method does not include radiation momentum in the momentum equation (Eq. (2)). However, this term is not dynamically relevant for the typical timescales and luminosities associated with stellar mergers and dynamical CE evolution (Glanz & Perets 2018; Bermúdez-Bustamante et al. 2024). Their approach to radiative diffusion also differs from the standard SPH discretisation used by Whitehouse & Bate (2004) and Whitehouse et al. (2005), instead relying on energy exchange between neighbouring particles that overlap with discrete rays traced from each particle. This allows cooling into empty regions devoid of particle neighbours. Although this is not possible with our method, as noted in Sect. 4.1, the photosphere in our simulations can still cool to surrounding optically thin material once the CE expands significantly. Before such expansion, the steep temperature gradient at the RSG photosphere is unresolved in any case, preventing calculation of the correct cooling rate.

Pejcha et al. (2016a, b) conducted 3D SPH simulations of L₂ mass outflows that incorporated radiative diffusion and a prescription for radiative cooling. They identified cases where bound gas falls back towards the binary and undergoes shock heating. In the limit of efficient cooling, this gas forms a circumbinary disk and a fraction of the heated gas escapes along the unobstructed poles. Despite our dissimilar setups, these findings are qualitatively similar to the bound CE ejecta and poloidal outflows discussed in Sect. 3.1.

A number of 1D CE simulations have included energy transport (Meyer & Meyer-Hofmeister 1979; Fragos et al. 2019; Bronner et al. 2024), though they often do not isolate the effects of radiative diffusion by directly comparing with an adiabatic simulation. Fragos et al. (2019) performed 1D CE hydrodynamic simulations involving a 12 M_⊙ RSG donor and a 1.4 M_⊙ neutron star companion. They incorporated radiative diffusion and convective energy transport using mixing length theory. Their model also included various energy sources and sinks, most notably, radiative losses from the photosphere, nuclear energy, and recombination energy. They found that the system transitioned to self-regulated evolution several years after the end of the plunge-in.

Both Fragos et al. (2019) and Bronner et al. (2024) showed that the majority of the envelope mass may be ejected, which starkly contrasts with the present work. However, the outcomes of 1D simulations are sensitive to assumptions regarding the size of the energy deposition region and the manner in which energy is deposited into the 1D stellar model. For instance, one extreme case in Bronner et al. (2024) required depositing energy in a shell four times thicker than the companion’s Bondi-Hoyle-Lyttleton sphere of influence to reproduce the orbital evolution seen in a 3D simulation. The inherent instantaneous energy homogenisation over spherical shells in 1D models may also overestimate mass ejection (though, see a recent approach to address this by Everson et al. 2025). Taken together, these uncertainties could plausibly explain the differences between 1D and 3D simulation results.

4.3. Implications for binary stellar evolution

Today, CE evolution remains a major bottleneck in understanding the formation of compact binaries almost half a century since its initial proposal (Paczyński 1976). Particularly, it is questionable whether 3D simulations have self-consistently demonstrated CE ejection. Simulations that do not model recombination energy release typically find that most of the envelope mass (70–90%) remains bound upon termination of the dynamical plunge-in (e.g. Lau et al. 2022a). In contrast, simulations that include recombination energy show that most if not all of the envelope could be fully ejected over many years after the plunge-in (e.g. Nandez et al. 2015; Sand et al. 2020; Reichardt et al. 2020; Lau et al. 2022a, b; Vetter et al. 2024). However, these simulations are all adiabatic, and whether complete ejection could be demonstrated with radiation transport remains an open question. One key debate centres on whether hydrogen recombination energy contributes to envelope ejection or is mostly radiated away (e.g. Sabach et al. 2017; Grichener et al. 2018; Ivanova 2018; Soker et al. 2018). Our results indicate that non-adiabatic effects influence both envelope ejection and ejecta morphology shortly after the plunge-in. The still-bound portion of the CE may be transitioning towards self-regulation (Meyer & Meyer-Hofmeister 1979; Podsiadlowski 2001) and it is unknown whether this still allows recombination-driven ejection, such as the radial disk winds characterised by Vetter et al. (2024).

On the other hand, 3D adiabatic simulations generally demonstrate that an immediate merger of the stellar cores may be averted even if most of the envelope remains bound. The bound material instead forms an extended, centrifugally supported structure around the central binary (Kuruwita et al. 2016; Gagnier & Pejcha 2023; Tuna & Metzger 2023; Vetter et al. 2024). Embedded in low-density, co-rotating gas, the central binary experiences slow inspiral, with an orbital decay timescale that is a thousand times longer than the orbital period (Lau et al. 2022a). Our radiative simulations show even more bound material persisting after the spiral-in compared to the adiabatic case, suggesting an even more massive disk. These results further challenge the prevalent picture of CE evolution involving a dynamical-timescale phase that rapidly shrinks the orbital separation while ejecting the entire envelope. Instead, most of the envelope is likely removed through a long-timescale, non-adiabatic process that may further influence the central binary. Observations of post-CE systems with distant tertiary companions support mass-loss timescales of 10²−10⁴ yr (Michaely & Perets 2019; Igoshev et al. 2020). This also suggests that the widely used ‘energy formalism’ or ‘α-formalism’ should not be applied to the binding energy of the entire hydrogen-rich envelope.

Alternative models for envelope ejection have been proposed. Recently, Hirai & Mandel (2022) suggested dividing the envelope ejection process into an adiabatic and non-adiabatic phase, with the convective-radiative boundary in the envelope serving as a possible division point. Their formalism translates into a CE efficiency parameter that is strongly dependent on system properties, ranging from α∼0.1 to >10⁴. Another line of research explores how the post-CE binary orbit evolves under the influence of a leftover circumbinary disk (e.g. Kashi & Soker 2011; Tuna & Metzger 2023; Izzard & Jermyn 2023; Wei et al. 2024; Valli et al. 2024; Unger et al. 2024). If a substantial fraction of the donor star’s envelope (≳10% of the post-CE binary mass) forms a circumbinary disk, its angular momentum could significantly alter the post-CE orbital separation and induce mild orbital eccentricities (e.g. D’Orazio & Duffell 2021; Siwek et al. 2023b, a). This has far-reaching implications for binary evolution. For example, in gravitational wave sources like binary black holes and double neutron stars, circumbinary disk torques may be necessary in certain cases to sufficiently shrink the post-CE orbit, allowing a merger to take place within a Hubble time (Moreno et al. 2022; Wei et al. 2024; Vetter et al. 2024; Landri et al. 2025). Beyond the CE phase, circumbinary disks could also drive a subsequent case BB mass transfer episode to become unstable and provide circumstellar medium in interacting supernovae (Tuna & Metzger 2023; Ercolino et al. 2024; Wei et al. 2024). For example, SN 2014c is consistent with interaction with toroidal circumstellar material that may have been ejected through CE evolution (Brethauer et al. 2022; Orlando et al. 2024).

5. Summary and conclusions

We performed 3D radiation hydrodynamics simulations of the CE phase involving a 12 M_⊙ RSG and a 3 M_⊙ companion, using an implicit two-temperature flux-limited diffusion solver implemented in the SPH code PHANTOM. By comparing these simulations to our previous adiabatic simulations, we assessed the impact of radiation transport on envelope ejection, ejecta morphology, and the final orbital separation. Our key conclusions are listed below:

• (i)

  Radiative diffusion significantly obstructs envelope ejection, reducing the unbound mass to less than half of that in the adiabatic simulation.

• (ii)

  The immediate post-plunge-in orbital separation is largely unaffected by radiative diffusion.

• (iii)

  The envelope mass ejected in the companion’s wake before the plunge-in is not affected by radiative diffusion. However, radiative diffusion causes little to no mass to be ejected during the plunge-in itself. Late-time ejection is also suppressed due to ejecta falling back towards the central binary, obstructing further mass outflow. During these phases, the outermost bound ejecta have short radiation transport timescales due to the opacity drop associated with hydrogen recombination.

• (iv)

  While poloidal outflows inflate clear bipolar lobes in the adiabatic simulation, the radiative case instead produces intermittent and sometimes single-sided plumes. As a result, the bipolar morphology is either absent or significantly less pronounced.

• (v)

  Although our simulations do not include recombination energy release. The hydrogen (and helium) partial ionisation zones have optical depths exceeding ten. This supports the notion that recombination radiation could thermalise in bound ejecta. Whether it can contribute to the expansion of the CE remains an open question that should be addressed in future simulations incorporating both recombination energy release and radiation transport.

This study represents a preliminary step towards understanding the role of radiation transport in CE evolution. As discussed in Sect. 4.1, the main limitations of our simulations include the breakdown of the flux-limited diffusion approximation in optically thin regions, difficulty capturing photospheric cooling near the start of the simulation, and the omission of recombination energy release. Addressing these limitations would be a crucial next step. It would also be valuable to validate our results using different hydrodynamics codes and radiation transport methods. Since our findings are based on a specific binary system, it would be instructive to explore the impact of radiative diffusion with different mass ratios and evolutionary stages of the donor star. Going forward towards the paradigm of partial dynamical CE ejection, the orbital separation and bound envelope mass recorded at the end of 3D simulations can serve as initial conditions for dedicated simulations of the subsequent non-adiabatic evolution.

Data availability

Movies (see Appendix D) – in ZIP file are available at https://www.aanda.org

Acknowledgments

We thank Orsola De Marco, Miguel González-Bolívar, Luis Bermúdez-Bustamante, and members of the PSO and SET groups at the Heidelberg Institute for Theoretical Studies for useful discussions. M. Y. M. L. is supported by a Croucher Fellowship. Parts of this research were supported by the Australian Research Council Centre of Excellence for Gravitational Wave Discovery (OzGrav), through project number CE230100016. The simulations presented in this work were performed on the OzSTAR national facility at the Swinburne University of Technology. The OzSTAR program receives funding in part from the Astronomy National Collaborative Research Infrastructure Strategy (NCRIS) allocation provided by the Australian Government, and from the Victorian Higher Education State Investment Fund (VHESIF) provided by the Victorian Government.

References

Appendix A: Energy and angular momentum conservation

Total energy is conserved within 10.7% for the ‘radiative’ simulation and 6.9% for the ‘radiative + heating’ simulation up to t = 8 yr. In the latter case, this does not include energy explicitly added to model the RSG’s nuclear luminosity, which totals 4.2×10⁴⁶ erg over 8 yr. The evolution of different energy components is shown in Fig. A.1, with most of the total energy drift occurring after the plunge-in. Energy conservation is significantly worse than in the adiabatic simulation from Paper I, which maintained total energy to within 0.05%. This difference arises from the implicit radiation solver, which ensures numerical stability with arbitrarily large time steps but does not conserve energy exactly. The relatively strong drift in total energy after the plunge-in could be associated with radiative diffusion in optically thin ejecta. Lowering the accuracy tolerance of the implicit solver could improve energy conservation. In contrast, angular momentum is conserved within 0.020% for the ‘radiative’ simulation and 0.022% for the ‘radiative + heating’ simulation, comparable to the adiabatic case.

Fig. A.1. Time evolution of kinetic energy (red), gravitational potential energy (blue), gas thermal energy (green), radiation energy (purple), and total energy (orange). Top panel: ‘Radiative’ simulation. Bottom panel: ‘Radiative + heating’ simulation.

Appendix B: Resolution test

To assess the sensitivity of our results to resolution, we compare the evolution of orbital separation (Fig. B.1) and unbound mass (Fig. B.2) using a lower-resolution run with ten times fewer SPH particles (N = 2×10⁵). For the purely adiabatic case, we display results from Lau et al. (2022b), corresponding to the simulation labelled ‘None, fixed μ’.

Fig. B.1. Separation of the stellar cores as a function of time, similar to Fig. 7, but including results obtained with ten times fewer SPH particles (dash-dotted lines).

Fig. B.2. Amount of unbound mass as a function of time, similar to Fig. 6, but including results obtained with ten times fewer SPH particles (dashed lines). Material is considered to be unbound if the sum of its kinetic and potential energy is positive (ek+ep>0).

In Fig. B.1, all simulations show similar orbital separation evolution until the pronounced drop in the ‘radiative’ case near 8.5 yr, which was discussed in Sect. 3.3. This drop in the orbital separation is absent in the lower-resolution run, suggesting it may be a stochastic effect caused by turbulent ejecta in the slow spiral-in or a numerical artefact. Despite this, there is little resolution dependence in the orbital separation evolution for at least ≈1−2 years after the plunge-in, preserving the key result that radiative diffusion has minimal impact on the post-plunge-in separation. At lower resolution, the onset of the plunge-in is delayed by ≈0.1 yr in all cases. This may be due to a small amount of unresolved surface mass on the RSG, which delays the onset of mass transfer that initially contracts the orbit.

Figure B.2 compares the evolution of unbound mass across the two resolutions. In the ‘radiative’ case, the results agree within 7% in the first six years. Larger deviations arise during intermittent plume ejections (Sect. 3.1), when spatial resolution in the binary’s vicinity decreases due to lower density. The adiabatic case also shows good agreement until ≈4 yr, around the transition towards a slower inspiral. At both resolutions, the adiabatic simulation exhibits sustained mass ejection after the plunge-in that is much weaker with radiative diffusion, suggesting that this qualitative difference is robust.

Appendix C: Test problems for radiation hydrodynamics

We present results of selected radiation hydrodynamics test problems following Whitehouse & Bate (2004) and Whitehouse et al. (2005). An important distinction is that our tests are performed in 3D with SPH particles arranged in a hexagonal close-packed lattice. For tests involving a discontinuity (Appendices C.3 and C.4), the setup and boundary conditions follow the description in Price et al. (2018) for the shock tube test. In these tests, the discontinuity is defined on a constant-x plane, with periodic boundaries imposed in the y- and z-directions. The first and last layers of SPH particles along the x-direction are tagged as boundary particles so that their properties are fixed. The results presented here can be reproduced using the PHANTOM test suite and are subject to regular unit testing. In all tests, we assume an ideal gas with an adiabatic index of γ = 5/3 and a mean molecular weight of μ = 2. We set an accuracy tolerance of 10⁻⁶ for the implicit solver.

C.1. Gas-radiation energy exchange

We test the gas and radiation energy exchange terms in Eqs. (3) and (4). In the first test, we initialise a static, uniform-density, and optically thick gas with T_gas≫T_rad, allowing radiative emission from the gas to heat the radiation until thermodynamic equilibrium is reached (T_gas=T_rad). The gas is assumed to have a fixed opacity of κ = 0.4 cm² g⁻¹. The initial conditions are u = 10¹⁰ erg g⁻¹, ρ = 10⁻⁷ g cm⁻³, and ξ = 10¹² erg g⁻¹. Unlike other tests, SPH particles are arranged in a cubic lattice within in a cube of side length 1 AU with periodic boundaries.

A second test is set up identically, but with T_gas≪T_rad, so that the gas is instead heated towards the radiation temperature. In this case, the initial gas thermal energy is set to be u = 10² erg g⁻¹. Figure C.1 shows the evolution of thermal energy for both tests. The results obtained with the implicit solver using 40 time steps (plotted as markers) closely match those obtained with an explicit solver (solid lines).

Fig. C.1. Evolution of thermal energy in the tests for the gas-radiation exchange terms. Markers show results obtained with the implicit radiation solver with forty logarithmically spaced time steps, which are in close agreement with results obtained from an explicit solver (solid lines). The case of gas cooling (Tgas≫Trad) is plotted in red while the case of gas heating (Tgas≪Trad) is plotted in blue.

C.2. Diffusion of a sinusoid

The radiative diffusion term in Eq. (3) is tested by simulating diffusion in a periodic particle lattice with a sinusoidally varying ξ. This test was implemented in PHANTOM by Biriukov (2020). We disable the gas-radiation exchange terms tested previously and set the flux limiter to the optically thick value, λ = 1/3. A fixed time step of Δt = 3.37×10⁻¹⁷ s is used. Particles are placed in a uniform lattice of size x∈[−0.5,0.5] AU, y,z∈[−0.1,0.1] AU, and mass density of ρ = 1.49×10⁻³⁰ g cm⁻³. The opacity is fixed at the value κ = 1 cm² g⁻¹. The initial specific radiation energy varies sinusoidally along the x-direction as ξ(t = 0) = ξ₀[1+0.1sin(2πx)], where ${\xi }_0=a_{\mathrm{rad}}{T}_{\mathrm{rad},0}^4/\rho$ with T_rad,0 = 100 K. As there is no fluid motion, Eq. (3) reduces to a heat equation, for which the analytical solution follows an exponential decay towards ξ₀, ξ(t) = ξ₀[1+0.1sin(2πx)exp(−4π²Dt)], where D=c/(3κρ) is the diffusion coefficient. Figure C.2 shows close agreement between the analytical solution and the results obtained from the implicit radiation solver, with a maximum relative deviation of 0.2%.

Fig. C.2. Diffusion of sinusoidally varying specific radiation energy. The markers represent results obtained from the implicit radiation solver, while solid lines show the analytical solution to the heat equation. The solutions at seven equally spaced times are plotted, with darker lines indicating earlier times and lighter lines representing later times.

C.3. Radiation pulse in an optically thin medium

We simulate a radiation pulse propagating in optically thin material to verify that the radiation flux is correctly limited by the speed of light. The setup consists of gas and radiation initially in thermal equilibrium in a box with dimensions x∈[−0.1,0.9] cm and y,z∈[−0.1,0.1] cm. The gas density is 0.025 g cm⁻³ and the opacity is fixed at 0.4 cm² g⁻¹. The radiation energy is initially ξ = 0.4 erg g⁻¹. At t = 0, radiation energy in the region x<0 is instantaneously increased by a factor of ten to ξ = 4 erg g⁻¹. Layers of boundary particles are then used to fix ξ at the x-boundaries. Figure C.3 shows the ξ-profile along the x-direction at t = 10⁻¹¹ s, along with the expected pulse position of ct = 0.3 cm (dotted line). Results are displayed for implicit time steps Δt equal to 1, 10, 100, and 1000 times the explicit time step, Δt_explicit. The ξ-profiles closely agree, with maximum errors relative to the case of Δt/Δt_explicit = 1 of 1.3%, 1.6%, and 4.3% for Δt/Δt_explicit = 10, 100, and 1000, respectively.

Fig. C.3. Propagation of a radiation pulse in an optically thin medium computed with different time steps ranging from the explicit time step to one thousand times the explicit time step. The discontinuity in ξ propagates from x = 0 at t = 0 to the expected position of ct = 0.3 cm (dotted line) at t = 10−11 s.

C.4. Adiabatic and isothermal shocks

We set up a shock tube with varying opacities to explore different shock regimes, ranging from adiabatic to nearly isothermal. Similar to Appendix C.3, we place particles in an elongated 3D box with a length of x∈[−10¹⁵,10¹⁵] cm. The gas has a density of ρ₀ = 10⁻¹⁰ g cm⁻³ and is initially in thermal equilibrium with a radiation field at the temperature T₀ = 1500 K. Material in the region x<0 is given an initial velocity of v₀ = 3.2×10⁵ cm s⁻¹, while material in x>0 is given a velocity of −v₀. The flows collide at x = 0 and form a shock. We simulate the shock tube with varying gas opacities, κ = 40,0.4,4×10⁻³,4×10⁻⁵ cm² g⁻¹, progressively transitioning from an adiabatic to a nearly isothermal shock. Figure C.4 shows the density and temperature distribution at t = 10⁹ s for each opacity. The analytical solution for the adiabatic (dotted line) and isothermal (dashed line) shocks are also shown for comparison. The shock speed may be obtained from the Rankine-Hugoniot jump conditions as

Fig. C.4. Solutions for the radiative shock at t = 109 s for different opacities. The top row corresponds to the optically thick (adiabatic) case, gradually transitioning to the optically thin (isothermal) case in the bottom row. The dotted line shows the analytical solution for the adiabatic case while the dashed line corresponds to the isothermal solution. The left column shows density profiles, and the right column shows gas temperature profiles.

$$D=\frac{1}{4}[({\gamma }_{\mathrm{eff}}-3)+\sqrt{\frac{16{\gamma }_{\mathrm{eff}}{p}_0}{{\rho }_0{v}_0^2}+({\gamma }_{\mathrm{eff}}+1{)}^2}]v_0,$$(C.1)

where p₀ is the pre-shock pressure and γ_eff is the effective adiabatic index, assuming the value of 5/3 in the adiabatic limit and 1 in the isothermal limit. The post-shock density, ρ₁, is given by

$$\frac{{\rho }_1}{{\rho }_0}=1-\frac{v_0}{D},$$(C.2)

while the post-shock temperature in the optically thick case is given by

$$\frac{T_1}{T_0}=\frac{{\rho }_0}{{\rho }_1}-\frac{{\rho }_0{v}_{0}D}{p_0}$$(C.3)

in the limit where gas pressure dominates. Figure C.4 closely reproduces the results of Fig. 3 from Whitehouse et al. (2005) and demonstrates agreement with both the adiabatic and isothermal solutions. In the top panels, the temperature bump and associated density dip near x = 0 are due to wall heating (Noh 1987; Rider 2000) and are also seen in Whitehouse et al. (2005). They are smoothed out in the other panels due to more efficient radiative diffusion when approaching the isothermal limit.

Appendix D: Movies of simulations

All Tables

All Figures

Fig. 1. Comparison of the initial radiative luminosity profile (lrad) with that of a MESA RSG reference model. Each black dot shows one of two million SPH particles. For the MESA model, the blue line shows the total luminosity, the orange line shows the convective luminosity, and the green line shows the radiative luminosity. The dashed vertical line divides the convective envelope from the radiative interior.

In the text

Fig. 2. Density in the initial orbital plane (z = 0) in 4320-by-3000 R⊙ panels. The crosses mark the locations of the donor’s core and the companion. Each column corresponds to a different simulation while each row corresponds to a different stage in the evolution. Row 5 shows the slow inspiral taking place several years after the plunge-in. The inspiral timescale, $t_{\text{inspiral}}=a/(-\dot{a}$), is several hundred to a thousand times longer than the orbital period. This figure was made with the SARRACEN package (Harris & Tricco 2023). Online movies of our simulations are listed in Table D.1.

In the text

	Fig. 3. Same as Fig. 2, except each 4320-by-3000 R⊙ panel shows a slice viewed edge-on that contains the stellar cores. Online movies of our simulations are listed in Table D.1.
In the text

	Fig. 4. Same as Figs. 2 and 3, but showing the late-time ejecta (t = 7.71 yr) with a larger scale (18 200-by-14 000 R⊙) and a lower-density colour map. Top row: Orbital-plane slices. Bottom row: Edge-on slices.
In the text

Fig. 5. Angular distribution of particles that become unbound between consecutive simulation snapshots, meaning the sums of their kinetic and potential energies become positive. The angle θ is the polar angle measured from $\hat{\text{z}}$, which is the direction of the initial orbital angular momentum vector. Material is initially ejected near the orbital plane, then closer to the poles in single or double-sided outflows after the plunge-in.

In the text

Fig. 6. Fraction of unbound envelope mass as a function of time for the adiabatic simulation (red), the ‘radiative’ simulation (blue), and the ‘radiative + heating’ simulation (green). Solid lines assume material is unbound when the total mechanical energy is positive (ek+ep>0), while dashed lines are less restrictive, considering material to be unbound when the sum of the mechanical and internal energy is positive (ek+ep+egas+erad>0).

In the text

	Fig. 7. Separation between the RSG core and the companion as a function of time for the adiabatic simulation (red), the ‘radiative’ simulation (blue), and the ‘radiative + heating’ simulation (green).
In the text

Fig. 8. Top to bottom: Various quantities in orbital-plane slices from the ‘radiative’ simulation: density (ρ), local diffusion time (tdiff, Eq. (11)), kinetic plus potential energy per unit mass (ekp), opacity (κ), and the density-weighted hydrogen ionisation fraction ([HII]). Left to right: Snapshots representative of the pre-plunge-in phase, the end of the rapid plunge-in, and the late-time slow spiral-in. In rows (d) and (e), we overlay contours where the optical depth equals 1, 3, and 10, from largest to smallest size. The contours in row (d) include the contribution from all material, whereas those in row (e) include only the contribution from bound material (ekp>0). The fact that these contours enclose the hydrogen partial ionisation zone (marked by the sharp transition between green and pink) suggests that a significant fraction of hydrogen recombination energy could be absorbed by the bound material.

In the text

	Fig. A.1. Time evolution of kinetic energy (red), gravitational potential energy (blue), gas thermal energy (green), radiation energy (purple), and total energy (orange). Top panel: ‘Radiative’ simulation. Bottom panel: ‘Radiative + heating’ simulation.
In the text

	Fig. B.1. Separation of the stellar cores as a function of time, similar to Fig. 7, but including results obtained with ten times fewer SPH particles (dash-dotted lines).
In the text

	Fig. B.2. Amount of unbound mass as a function of time, similar to Fig. 6, but including results obtained with ten times fewer SPH particles (dashed lines). Material is considered to be unbound if the sum of its kinetic and potential energy is positive (ek+ep>0).
In the text

Fig. C.1. Evolution of thermal energy in the tests for the gas-radiation exchange terms. Markers show results obtained with the implicit radiation solver with forty logarithmically spaced time steps, which are in close agreement with results obtained from an explicit solver (solid lines). The case of gas cooling (Tgas≫Trad) is plotted in red while the case of gas heating (Tgas≪Trad) is plotted in blue.

In the text

	Fig. C.2. Diffusion of sinusoidally varying specific radiation energy. The markers represent results obtained from the implicit radiation solver, while solid lines show the analytical solution to the heat equation. The solutions at seven equally spaced times are plotted, with darker lines indicating earlier times and lighter lines representing later times.
In the text

	Fig. C.3. Propagation of a radiation pulse in an optically thin medium computed with different time steps ranging from the explicit time step to one thousand times the explicit time step. The discontinuity in ξ propagates from x = 0 at t = 0 to the expected position of ct = 0.3 cm (dotted line) at t = 10−11 s.
In the text

Fig. C.4. Solutions for the radiative shock at t = 109 s for different opacities. The top row corresponds to the optically thick (adiabatic) case, gradually transitioning to the optically thin (isothermal) case in the bottom row. The dotted line shows the analytical solution for the adiabatic case while the dashed line corresponds to the isothermal solution. The left column shows density profiles, and the right column shows gas temperature profiles.

In the text