Dust Formation in Common Envelope Binary Interactions — III. Lightcurves
June 25, 2026
Abstract
Luminous red novae are transient events thought to arise from common envelope binary interactions. In this paper, we perform post-processing light-curve calculations for two, 3D hydrodynamic simulations of common envelope events. These simulations model interactions between 1.7 \(M_\odot\) and 3.7 \(M_\odot\) asymptotic giant branch stars and a 0.6 \(M_\odot\) compact companion, including dust nucleation. In both our simulations, which are carried out for 44 years under adiabatic conditions, we observe a bright, hot peak lasting \(3-5\) years, primarily due to the expansion of the photosphere before and during inspiral. Additional peaks can be seen appearing at different times and for different viewing angles, due to the asymmetry of the interaction. Dust forms about \(1-3\) years after the beginning of the simulated interaction and shortly afterwards we witness a sharp decline in the bolometric luminosity, followed by a partial recovery and a plateau with an effective temperature of \(\sim\)400 K. The dust photosphere reaches a size of \(\sim\)250 au by the end of the simulations, but we predict that between 100 and 200 years, the dust will become optically thin at visible wavelengths, revealing an inner, warmer photosphere. The lightcurves obtained have two, well-quantified, but large uncertainties: insufficient surface resolution primarily affecting the first 1-2 years of the simulated lightcurves and the adiabatic assumption that affects primarily the later years. We finally contextualise our simulations within a group of observed luminous red nova transients, drawing particular attention to the outburst of OGLE-2002-BLG-360 and AT 2025abao, which are the closest match to our simulation.
binaries: close – stars: AGB and post-AGB – stars: winds, outflows – ISM: planetary nebulae
1 Introduction↩︎
Red novae and LRNe are stellar transients arising from low- and high-mass binary interactions, respectively [1], [2]. These events occur when two non-degenerate stars—or a growing subgiant and a compact companion—undergo a CE interaction [3], culminating in either a stellar merger or a tight binary system. While detections in our own Galaxy are rare—V 838 Mon [4], V 1309 Sco [5], V 4332 Sgr [6], OGLE-2002-BLG-360 [7]—the advent of dedicated transient surveys has led to a surge in extragalactic LRN events, resulting in a total of \(\sim\)30 studied objects [2]. With the Vera Rubin Telescope now operational, we anticipate that its population will dramatically increase during the Rubin’s 10-year Legacy Survey of Space and Time [8], tracing CE evolution for a broad range of binary configurations and remnants. Characterizing these inherently three-dimensional interactions requires 3D hydrodynamic models capable of predicting lightcurves and physical observables, ultimately enabling the “fingerprinting” of the entire interaction process.
Three-dimensional hydrodynamic simulations of the CE binary interaction have come a long way [9], but the prediction of the lightcurve has been plagued with a number of challenges [10]. Recently [11] predicted the lightcurve of the LRN V 1309 Sco, which represented the merger of a low mass subgiant 1.5M\(_{\odot}\)star with a 0.1M\(_{\odot}\)companion. They used the SPH code starsmasher [12] and post-processed the light to provide observable (photospheric) properties of the outburst, as well as multi-wavelength lightcurves in the photometric Johnson-Cousins system. Since most LRNe observations show the fast production of dust [13], [14], [11] mimicked dust production, but treated it as an instantaneous process by using opacity tables that include both gas and dust opacities tabulated as a function of quantities like density and temperature, as an indication of when dust formation would be taking place.
[15] have carried out CE simulations of 1.7 M\(_{\odot}\)and 3.7 M\(_{\odot}\)AGB stars interacting in a CE with 0.6 M\(_{\odot}\)compact companions, including the nucleation of carbon-based dust. They demonstrated that dust takes a short, but non-zero time to form in the context of a CE interaction; the formation is very efficient and all the gas that can be converted into carbon grains does so in a matter of a couple of decades. Additionally, the increasing amount of dust rapidly obscures the gas photosphere. The timing of dust formation therefore leads to completely different lightcurve properties and accurate predictions of the dust formation, properties and distribution is therefore paramount.
In this paper, we carry out a first step at bridging the gap between simulations and observations, by performing a calculation of the light emitted from the dusty CE simulations of [15]. The method is similar to what was used by [10] and [11], namely it carries out radiation transport along rays from the CE to the observer. The light contribution into the ray is then calculated with a blackbody approximation assuming LTE.
In Section 2, we explain how we calculate the lightcurve, starting with the hydrodynamic models (Section 2.1), the computation of their grey luminosity (Section 2.2), and the explanation of the opacities adopted (Section 2.3). We then present our results in Section 3 starting with the lightcurve in Section 3.1, a discussion of the properties of the photosphere in Section 3.2. In Section 4 we discuss the various sources of uncertainty. We start with the SPH simulation resolution in Section 4.1; we continue with other issues stemming from the properties of the surface layers in Section 4.2, carry out a critique of the late time lightcurve in Section 4.3 and a comparison with observations in Section 4.4. We finally present our conclusions in Section 5.
2 Method↩︎
2.1 Hydrodynamic models↩︎
We calculate the lightcurve from four 3D hydrodynamic simulations performed with the smoothed particle hydrodynamics code phantom [16], fully described in [15] and [17]. Two of the simulations are CE interactions between a 1.7 M\(_{\odot}\), 260 R\(_{\odot}\), 5180 L\(_{\odot}\), 3227 K, asymptotic giant branch star and a 0.6 M\(_{\odot}\)companion, while the other two used a 3.7 M\(_{\odot}\), 330 R\(_{\odot}\), \(1.19\times 10^4\) L\(_{\odot}\), 3317 K, AGB star with the same companion mass. All simulations are started at an orbital separation such that the giant fills its Roche lobe (550 and 637 R\(_{\odot}\), respectively). Each pair of simulations computed the evolution of the CE interaction with and without dust nucleation. We employ a tabulated EoS that includes the input of recombination energy into the gas, where it is assumed to be instantly thermalised [18]. The envelope of both stars is almost entirely unbound after a simulation time of \(\sim\)40 years. Both simulations are adiabatic with no radiation transport.
Dust forms approximately 3 years after the start of the simulations, before the in-spiral has completed at about 5 years. As illustrated in figure 11 of [15], in both simulations almost all dust that can form has formed by 30 to 40 years, limited by the C/O=2.5 and a metallicity of \(Z=0.02\). No dust destruction was implemented. The dust layer quickly becomes optically thick so that the photosphere is at the outer edge of the dusty shell. It expands over four decades, becoming larger and cooler.
2.2 Calculation of the Grey Luminosity↩︎
In this section, we apply the radiative transfer equation to determine the system’s bolometric luminosity at specific time intervals. We define our coordinate system such that the observer is positioned at \(z \rightarrow +\infty\), viewing the object (centred at \(z=0\)) along the negative \(z\)-axis. To calculate the surface brightness, we project a grid onto the object’s observable face and trace a single ray per pixel—parallel to the \(z\)-axis—to solve for the specific intensity, \(I_i\) (see Figure 13 in Appendix 7, where we also provide additional details of the computation). From these values, the isotropic-equivalent luminosity, \(L\), is derived as follows \[L \approx 4 \pi \sum_i I_i \Delta A_i, \label{eq:def:luminosity}\tag{1}\] where \(\Delta A_i\) is the area of a pixel associated with the \(i\)-th ray. As rays are evenly spaced, \(\Delta A_i = \Delta A\) are the same for all rays.
To derive the specific intensity, \(I_i\), for each ray, assuming a grey opacity \(\kappa_\lambda = \kappa\) (more on opacities in Section 2.3), we solve the radiative transfer equation, for a given ray \[\frac{dI_i}{d\tau_i} = I_i(\tau_i) - S_i(\tau_i), \label{eq:rad95trans}\tag{2}\] where \(S_{i}\) is the source function along the ray, and \(\tau_i\) is the optical depth at any given position \(z\) \[\tau_i(z) = \int_z^{+\infty} \kappa_i(z') \rho_i(z') dz', \label{eq:def:tau}\tag{3}\] where \(\kappa_{i}\) is the wavelength-integrated opacity and \(\rho_i\) is the density along the ray. Ignoring any input of light from the far side of the object and any scattering into the beam, \(I_i(\tau=\tau_{i, \mathrm{max}})\), the solution to Equation 2 is therefore \[I_i(z=+\infty) = \int_{-\infty}^{+\infty} e^{-\tau_i(z')} S_i \kappa_i \rho_i dz'. \label{eq:radiative-transfer-solution-backwards}\tag{4}\] Here we note that \(S \kappa \rho = S \frac{d\tau}{dz}\) is the source function integral per physical distance instead of per optical depth.
In the specific case of SPH simulations, the specific intensity \(I_i\) is ultimately determined by the many particles that intersect with the ray, \(i\), within the particle’s kernel cutoff radius, \(2h\). This is because all quantities on the ray are derived from the quantities from these neighbouring particles, as illustrated in Figure 14, Appendix 7. More specifically, to discretize Equation 4 in SPH, we use kernel interpolation ([19]) to interpolate the product \(S \kappa \rho\) from neighbouring particles \[I_i(z=+\infty) \approx \int_{-\infty}^{+\infty} e^{-\tau_i(z')} \sum_j \frac{m_j (S \kappa \rho)_j}{\rho_j h_j^3} w(q_{ij}(z')) dz', \label{eq:calc:kernel-interpolation:I95i}\tag{5}\] where \(m_j, S_j, \kappa_j, \rho_j, h_j\) is the mass, source function, opacity, density and smoothing radius of the \(j\)-th particle; \(q_{ij}\) is the dimensionless distance from \(j\)-th particle to \(i\)-th ray, defined as \[q_{ij} (z) \equiv \frac{\sqrt{(x_i-x_j)^2+(y_i-y_j)^2+(z-z_j)^2}}{h_j}; \label{eq:def:q95ij}\tag{6}\] and \(w\) is the dimensionless kernel, derived from the smoothing kernel, \(W\), as \[w(q_{ij}(z)) \equiv h_j^3 W(q_{ij}(z), h_j).\] The value of \(w(q_{ij}(z))\) goes to zero when \(q_{ij}(z) \geq 2\) from a central peak of \(\sim\)0.318 in 3 dimensions [16]. In other words, only neighbouring particles within \(2 h_j\) from the ray are counted when interpolating the values on the ray. When kernel-interpolating quantities deriving from the product of primary variables, we first multiply the relevant variables at each particle’s location, and then we kernel-interpolate the product to the ray. We discuss the merit of this procedure in Appendix 8.
To further stress that the SPH particles are the actual source of information that drives the calculation of \(L\), we define the (wavelength and viewing-angle-dependent) effective cross-section area \(A_{\mathrm{eff}, j}\) of each SPH particle, \(j\), to be \[A_{\mathrm{eff}, j} \equiv \frac{m_j \kappa_j}{h_j^3} \sum_i \Delta A_i \int_{-\infty}^{+\infty} e^{-\tau_i(z')} w(q_{ij} (z')) dz'. \label{eq:def:Aeffj}\tag{7}\] By combining Equation 1 , 5 , and 7 , we can now find the luminosity by \[L \approx 4 \pi \sum_j S_j A_{\mathrm{eff}, j}. \label{eq:calc:luminosity}\tag{8}\]
As we assume LTE, the source function, \(S_j\), for particle \(j\) in this case is the grey blackbody radiation \[S_j = \frac{\sigma_\mathrm{sb} T_j^4}{\pi}, \label{eq:assumption:S95j}\tag{9}\] where \(\sigma_\mathrm{sb}\) is the Stefan-Boltzmann constant2. Scattering of photons into the ray is ignored even though electron scattering opacities are included and indeed dominate the opacity source in the ionised medium.
For the optical depth in Equation 7 , we use kernel interpolation ([20]) again to derive a discretised version of \(\kappa_i(z') \rho_i(z')\) in Equation 3 , namely \[\kappa_i(z') \rho_i(z') \approx \sum_j \frac{m_j (\kappa \rho)_j }{\rho_j h^3_j} w(q_{ij}(z')).\] Equation 3 now becomes \[\tau_i(z) \approx \sum_j \frac{m_j \kappa_j}{h^3_j} \int_{z}^{+\infty} w(q_{ij}(z')) dz', \label{eq:calc:tau95z}\tag{10}\] which can be used to calculate the optical depth at any given point.
The photosphere is defined as the place where \(\tau = 2/3\) along rays that are parallel to the line of sight. From the locations of the photosphere along each ray, a mean photospheric radius can be defined as explained in Section 3.2. The value of \(\kappa, \tau\) and \(T\) as a function of position along two specific rays are shown in Figure 1 and corresponding values are reported in Appendix 9.
Figure 1: Opacity, \(\kappa\), optical depth, \(\tau\), and temperature, \(T\), along one line-of-sight as viewed from the polar \(z \rightarrow +\infty\) (blue) and the orbital \(y \rightarrow +\infty\) (orange) and \(x \rightarrow +\infty\) (green) directions, for the 1.7 M\(_{\odot}\) simulation at \(t=0\) (top left), \(t=3~\mathrm{yr}\) (top right), \(t=12~\mathrm{yr}\) (bottom left), and \(t=20~\mathrm{yr}\) (bottom right). The horizzontal dashed line in the middle panels marks \(\tau=2/3\). The photospheric location is marked with blue vertical shaded lines as seen from the \(z \rightarrow +\infty\) direction, with the line width representing the local smoothing length at the photosphere, giving a sense of uncertainty of the photospheric \(\kappa, \tau\), and \(T\). Corresponding values are listed in Table 4 in Appendix 9. Full movie version are at https://zenodo.org/records/18239846..
2.3 Gas and Dust Opacities↩︎
The dust opacity, \(\kappa_\mathrm{dust}\), for our carbon dust simulations is the same as that calculated by [15], where the Planck mean opacity of the dust-gas mixture is calculated using Mie theory (see equation 7 of that paper). Planck mean dust opacities are used instead of Rosseland mean opacities, as they are more suitable for the optically thin outer regions, near the photosphere, where dust forms and light emits.
The gas opacity, \(\kappa_\mathrm{gas}\), derives from the tabulated values adopted by the mesa code, as implemented in phantom by [18], which combine a range of different opacity sources across a wide range of temperatures from \(631~\mathrm{K}\) to \(5\times10^8~\mathrm{K}\), as shown in Figure 2. We only adopt the mesagas opacity values for \(T>1450\) K because for lower temperatures the tabulated opacities are those for oxygen-rich dust ([21], [22]), while we implement our own carbon dust opacities as explained above. For \(T<1450\) K, the gas opacity is assumed to be a low value compatible with recombined gas being transparent, or \(\kappa_\mathrm{gas}=2 \times 10^{-4} \mathrm{cm}^2 \mathrm{g}^{-1}\). The opacity discontinuity introduced by this temperature cutoff is negligible for our dusty simulations, because at the temperature of the discontinuity, typically, the opacity from carbon dust (\(\sim 10~\mathrm{cm}^2 \mathrm{g}^{-1}\)) dominates the gas opacity, as seen in Figure 3. The total opacity \(\kappa\) is the sum of dust and gas opacities: \(\kappa = \kappa_\mathrm{dust} + \kappa_\mathrm{gas}\).
[11] also assumed wavelength-independent opacities. They averaged two opacity tables that extend to different temperature regions: that of [23] in high temperature regions (\(\log T > 3.75\)), and those of [24] in low temperature regions (\(\log T < 4\)). Similarly to us they set their gas opacity to \(10^{-4} \textrm{cm}^2 \textrm{g}^{-1}\) when \(T \leq 3000\) K (instead of our lower \(1450\) K) for no-dust scenarios. For dust-forming scenarios, they did not include dust nucleation, but instead represented the overall opacity with either the Planck mean opacity or the Rosseland mean opacity version of the [24] tables, for the full low temperature region. We only consider the Planck mean opacity because we do not include radiative cooling in our simulations, and the outgoing flux is determined in the region around and above the photosphere, which is optically thin when the resolution is good, where Planck opacity applies. When the resolution is poor, the emerging luminosity is overestimated anyway and it is determined only by the first particle encountered by the ray, such that changing opacity to Rosseland will make no difference.
3 Results↩︎
3.1 The lightcurve↩︎
In Figure 4 we present a total intensity image obtained by tiling \(256 \times 256\) pixels (this resolution choice is justified in Appendix 10.1). The surface brightness at each pixel, \(i\), is equal to the specific intensity, \(I_i\). An integrated luminosity is presented in Figure 5 for both simulations with and without dust.
Figure 4: Surface brightness from each pixels for the 1.7 M\(_{\odot}\) (left two columns) and 3.7 M\(_{\odot}\) (right two columns) simulations, as viewed from the polar (\(z \rightarrow +\infty\)) direction (first and third columns) and the orbital (\(y \rightarrow +\infty\)) direction (second and fourth columns), at \(t=0, 3, 5.1, 12.1, {\rm and} 20.2~\mathrm{yr}\) for the top, second, third, fourth and bottom rows, in each panel, respectively - note the progressively increasing size of the object and variable colour scales. The origin coincides with the primary’s core. Movies corresponding to this still image can be found at https://zenodo.org/records/18239846 for both 1.7M\(_{\odot}\) (2md) and 3.7M\(_{\odot}\) (4md) simulations, viewed from \(z \rightarrow +\infty\) (xyz), \(y \rightarrow +\infty\) (xzy) and \(x \rightarrow +\infty\) (zyx)..
Figure 5: Bolometric luminosity lightcurve of the of the 1.7 M\(_{\odot}\) (left panels) and the 3.7 M\(_{\odot}\) (right panels) models for dusty (top panels) and non-dusty (bottom panels) simulations, integrated using a grid of \(256\times256\) rays. Blue, orange, and green curves are the results viewed from \(+z\), \(+y\), and \(+x\) directions, respectively. The shaded areas are the calculated uncertainties (see Section [sec:uncertainties:surface-particles]). The light grey dashed line is an illustration of the orbital separation between the donor and companion sink, superimposed here for reader’s convenience..
In Figure 6, we approximate the SED by considering a wavelength dependent blackbody emissions from each particle, while still assuming grey opacity. In other words, we assume each particle has the following source function \[S_{\lambda, j} = \frac{2 h c^2}{\lambda^5 ( e^{\frac{hc}{\lambda k_B T}} - 1 ) },\] and use this expression in Equation 8 .
In Figure 7 we derive the mean effective temperature of the model, \(T_\mathrm{eff}\), determined by finding the peak of \(f_\lambda\) in the SED (as seen in Figure 6) and applying Wien’s displacement law, \(T_\mathrm{eff} = b_\mathrm{wien} / \lambda_\mathrm{peak}\), where \(b_\mathrm{wien}\) is the Wien’s displacement constant. After dust starts to form, a second cooler component appears in the SED, resulting in the sudden drop of \(T_\mathrm{eff}\) when the cool peak becomes stronger than the hotter SED component peak from the stellar gaseous photosphere. The momentary comeback of the hot peak in \(T_\mathrm{eff}\) around \(\sim\)6 yr for the \(+y\) direction in the 1.7 M\(_{\odot}\)model is caused by the hot gas component shining through the dust and temporarily outshining the dust component again.
![Figure 6: Spectral flux density as a function of wavelength for the 1.7 M_{\odot} (left column) and the 3.7 M_{\odot} (right column) dusty simulations, integrated using a grid of 256\times256 rays, seen from 1~\mathrm{kpc}. Blue, orange, and green are the three viewports, +z, +x and +y, respectively. The highest peak of f_\lambda (which is not the same as the peak of \lambda f_\lambda) is marked with a dashed vertical blue line. The colour reddens during the simulation, with a redder peak forming between 3 yr and 5.1 yr for the 1.7 M_{\odot} and 3.7 M_{\odot} simulations, caused by dust formation happening at that time [15].](/article/2606.26495/figures/Spec_256x256.png "fig:")
Next, we divide the lightcurve in Figure 5 into three phases and describe their characteristics, guided by the \(+z\) viewport — looking down the pole; blue curves in all figures.
The rise to the peak. The exact reason for the luminosity increase that we call the “outburst" in LRNe has been a source of debate. [25] interpreted the peak and plateau emission of M31 LRN 2015 as different phases of a common envelope ejection, the first comprised by the release of a small amount of fast material, while the latter by the slower release of recombination energy from a larger amount of expanding material. [26] attributed it instead to emission from an ejected and cooling shell, which shortly thereafter impacts and shocks previous, toroidal mass loss.
In our simulations, the bolometric luminosity increases to a peak value of \(\sim4\times 10^6\) L\(_{\odot}\)(Figure 5) and reaches temperatures of \(\sim 6500\) K (Figure 7), at 2.5 and 3 years after the start of the simulations, for the 1.7 and 3.7 M\(_{\odot}\)models, respectively. A strong expansion to 15–30 au (geometric radius, solid line in Figure 8) or 7 au (observed radius; dot-dashed line in the same Figure) powers the light rise. This very early expansion is due to the sudden insertion of the companion in the computational domain as well as stellar surface instabilities [17], whereas the inspiral drives the expansion shortly thereafter. This, alongside poor surface resolution (Section 4.2) is a major source of uncertainty for early times.
The side viewports (indicated by the green and orange lines in all figures) reveal that the 1.7 M\(_{\odot}\)model peaks earlier, whereas the 3.7 M\(_{\odot}\)model exhibits a delayed peak. Furthermore, the more massive model shows a broader and notably more irregular peak preceding the decline phase discussed below, characterized by multiple sub-peaks and troughs. By 5.1 years, the surface brightness maps (Figure 4 and the accompanying movie) show that the 1.7 M\(_{\odot}\)model is already dimming due to dust obscuration (discussed in detail below). In contrast, the 3.7 M\(_{\odot}\)model still exposes its bright core, particularly from the polar perspective; this occurs because dust formation is delayed in this higher-mass, lower-\(q\) (\(q = M_2 / M_1\)) simulation.
We note that the luminosity and temperature of the models at time zero is higher than we know it to be from the time-zero mesa profile, due to the effect of low surface resolution, as will be fully explained and quantified in Section 4.2 and Appendix 11. This is why the lightcurve plot has such large error-bars (shaded areas in Figure 5).
The decline to the trough. This phase goes from the peak to the trough (again using the pole-on view as reference), which starts at \(\sim 2.5\) yrs, for both models, at which point the bolometric light declines over 2.5 – 5 years for the two models (Figure 5). While the photosphere is still expanding in both models (from \(\sim\)20 to 70 R\(_{\odot}\)for the 1.7 M\(_{\odot}\)model and from 20 to 100 R\(_{\odot}\)for the 3.7 M\(_{\odot}\)model), the dominant peak in the SED for the 1.7 M\(_{\odot}\)model goes from thousands of degrees Kelvin to hundreds, as the dust enshrouds the object and dominates its colour. The 3.7 M\(_{\odot}\)model (pole view), on the other hand, is quite different: first of all, the light decline is not as severe as it is for the lower mass model, because the SED remains dominated by a (cooling) stellar photosphere. However, the other two (side) viewports display more variability during the peak, and decline later than the pole view light, but more abruptly and to a deeper through, as the dust forms a little later, and preferentially on the equator. In all cases, the light decline is that expected for an expanding, adiabatically cooling photosphere, where the only variable is when, where and how much dust forms.
The rise to the plateau. After the trough, the light recovers to a plateau value that is approximately flat (see Figure 5 insert, where the full evolution to 44 years is seen). The 1.7 M\(_{\odot}\)model rises rapidly and then more slowly to the end of the simulation. The 3.7 M\(_{\odot}\)model rises to a new peak and then rapidly declines (the peak has a width of\(\sim2\) yrs) to a plateau, which then continues to decline gently. As is clear from the movie associated with Figure 4, the 3D nature of the structure displays an irregular surface that can reveal outer or inner layers making the light temporarily brighter or dimmer. The luminosity of the plateau, at \(1-2\times 10^5\) L\(_{\odot}\)for both models is too large, due to the lack of any radiation/cooling, as explained in Section 4.3.
In Figure 5 we also show the lightcurve obtained for the two simulations carried out with no dust nucleation (and carefully compared to the dusty ones in [15]). The peak and trough behaviour lasting a few years from the start of RLOF observed in the dusty simulations is not present in the non-dusty cases. There, the CE gradually brightens and only much later dims (after about 11 years in the case of the 1.7 M\(_{\odot}\)model, while for the 3.7 M\(_{\odot}\)model no dimming is observed within the 12.5 yrs of the simulation), as the photosphere expands and adiabatically cools, with a decrease of the effective temperature.
Figure 7: Evolution of the effective temperature, \(T_\mathrm{eff}\), for the 1.7 M\(_{\odot}\) model (left panel) and the 3.7 M\(_{\odot}\) model (right panel). The thick lines indicate the temperature of the first, hot peak in the SED for two viewing directions. The presence of a corresponding thin line at a given moment in time indicates the presence of a second, cooler peak in the SED (evident in Figure 6). The horizzontal, dashed line marks the cut-off temperature \(T_\mathrm{cutoff}\) below which dust can condensate (see Figure 3)..
3.2 Properties of the photospheric surface↩︎
Figure 8: Photosphere radius as estimated by two distinct methods: a “geometric" method, using Equations eq. ¿eq:eq:def:Aph? and eq. ¿eq:eq:def:Rph? - solid lines, and an”observer’s" method - dot-dash and dashed lines for the 1.7 M\(_{\odot}\) (left panels) and 3.7 M\(_{\odot}\) (right panels) simulations, respectively. Blue (top panels) and orange (bottom panels) lines are for measurements in the \(z \rightarrow +\infty\) (polar) and \(y\rightarrow +\infty\) (orbital) directions, respectively. The dashed and dot-dashed lines are the photospheric radii \(R_\mathrm{ph, obs} = \sqrt{L / 4 \pi \sigma T_\mathrm{eff}^4}\), where \(T_\mathrm{eff}\) is the effective temperature obtained by fitting the gas (dot-dashed) and dust (dashed) peak of the SED as in Figure 7. Faint dot-dashed lines show the size of the gas photosphere seen through an optically-thin dust photosphere..
We calculate the size of the photosphere using two methods. The first method is by using the Stefan-Boltzmann relation with the calculated bolometric luminosity and temperature from the peak(s) of the SED. Observers fit SED measurements using one or more Blackbody curves, depending on what they observe. Often, in LRNs two peaks in the SED indicate that the stellar (gas) photosphere is visible behind a shell of optically thin dust. From fitting both peaks the temperatures and sizes of both can be extracted. At late times it is often the case that the dust becomes optically thick and only one, redder SED peak remains visible from which the size of the expanding dust cocoon can be deduced. This is also the case in the simulations, where a sample of the SED at different times can be seen in Figure 6. The first method therefore uses exactly the observers’ method.
The second method finds the size of the photosphere by counting the number pixels for a given viewing angle, where \(\tau > \frac{2}{3}\) \[A_\mathrm{ph} \equiv \sum_{ i \in \{ i | \tau_i > \frac{2}{3} \} } \Delta A_i. \label{eq:def:Aph}\tag{11}\] The photospheric size is hence simply \[R_\mathrm{ph} = \sqrt{\frac{A_\mathrm{ph}}{\pi}}, \label{eq:def:Rph}\tag{12}\] assuming a circular photosphere as projected on the plane of the sky3 .
The value of \(R_\mathrm{ph}\) using the two methods is plotted in Figure 8. Using the “geometric" method to measure the photospheric radius, we see a near-constant expansion rate during the 40 years (in Figure 8 we only plot the first 20 years). Both objects’ photospheres grow to about 250 au over the simulation time of 44 years. Using a simple linear regression fit, for the 1.7 M\(_{\odot}\) simulation after 5 years, the photosphere is expanding at almost 60 km s\(^{-1}\) in all three directions. For the 3.7 M\(_{\odot}\) simulation after 5 years, the photosphere is expanding at \(\sim\)70 km s\(^{-1}\) for the three directions.
Tracking the photospheric radius using the “observers" method reveals a more complex evolution. In Figure 8, the dot-dashed and dashed lines track the radius measurements of the first (warmer) and second (hotter) SED peaks, respectively. When the cooler, dust-dominated SED peak becomes stronger than the hotter, gaseous peak, the dashed line transitions to semi-transparent, signalling that the gaseous photosphere is being viewed through an increasingly optically thick layer of dust. Because the orbiting companion ejects chaotic, expanding, and eventually dusty plumes, the side (+y) view appears significantly more complex than the polar view. By roughly 20 years, the gasseous photosphere is entirely obscured by the dusty photosphere. Ultimately, while the”geometric" photosphere consistently outlines the object’s outermost boundary on the plane of the sky, the “observer" photosphere represents the temperature components directed toward the viewer. These two definitions align early on, while the expanding cocoon is symmetric and dust free.
3.3 The dispersion of the dust shell↩︎
After the opaque dust shell has fully enshrouded the object, it experiences free expansion that will reduce its density, decreasing its optical depth, which eventually will render it transparent even to optical light. Here we linearly extrapolate our results to investigate when and how this will happen.
Using only the dust opacity, \(\kappa_\mathrm{d}\), we can obtain the optical depth of the dust shell, \(\tau_\mathrm{d} (t)\), at given time, \(t\), using Equation 3 . Using \(\rho (r, t) = dm / (4 \pi r^2(m, t) dr)\), we can express Equation 3 in mass coordinates \[\tau_\mathrm{d}(t) = \int_0^M \frac{\kappa_\mathrm{d}}{4\pi r^2} dm.\] We assume that the dust opacity scales linearly with the dust/gas equilibrium temperature, \(T\), in the later stage of the simulation as demonstrated by [15] — ignoring radiative losses as well as the fact that dust and gas temperature may decouple, giving \[\kappa_\mathrm{d}(T) = \kappa_\mathrm{d, 0} \frac{T}{T_0} = \kappa_\mathrm{d, 0} \frac{r_0^{2}}{r^{2}},\] where \(\kappa_\mathrm{d,0}\), \(T_0\) and \(r_0\) are constants, representing the quantities at a given time. Where for the second equality we have used the fact that temperature, volume and radius have the following relation for adiabatic changes \[T \propto V^{1-\gamma} \propto r^{3 - 3 \gamma},\] where \(\gamma = 5/3\) for gas, and the volume dependency on \(r^3\) as the shell is expanding homologously. Substituting we obtain \[\tau_\mathrm{d}(t) = \int_0^M \frac{\kappa_\mathrm{d, 0} r_0^{2}}{4\pi r^{4}} dm.\] Finally, as free expansion (\(r(m, t) = v_0(m) (t-t_0)\)) is assumed, we can derive \[\tau_\mathrm{d}(t) = K (t - t_0)^{-4}, \label{eq:calc:dust-dispersion:tau}\tag{13}\] where \(K \equiv \int_0^M \frac{\kappa_\mathrm{d, 0} r_0^{2}}{4\pi v^{4}}dm\).
In Figure 9 we plot the dust shell’s optical depth as a function of time to forecast when \(\tau_\mathrm{dust}(t)\) drops below \(2/3\) (\(t_{\rm gone}\)). Using least squares fitting from the scipy package, we fit \(K\) and \(t_0\) in Equation 13 and obtain \(t_{\rm gone}\) as shown in Table 1. The dust shell is dispersed \(120 - 160\) years after the end of the inspiral for the 1.7 M\(_{\odot}\) simulation, and \(200 - 230\) years for the 3.7 M\(_{\odot}\) simulation.
| Case | \(t_0\) | \(t_\mathrm{drop}\) | \(t_\mathrm{gone}\) |
| (yr) | (yr) | (yr) | |
| \(+z\) | \(8.212 \pm 0.006\) | 32 | \(128\) |
| 1.7 \(+y\) | \(8.574 \pm 0.004\) | 32 | \(158\) |
| \(+x\) | \(8.562 \pm 0.003\) | 32 | \(157\) |
| \(+z\) | \(9.194 \pm 0.023\) | 41 | \(225\) |
| 3.7 \(+y\) | \(9.252 \pm 0.006\) | 41 | \(229\) |
| \(+x\) | \(9.079 \pm 0.010\) | 41 | \(207\) |
In both models dust formation starts in an equatorial belt (figures 3, 4 and 13 in [15]) with the 1.7 M\(_{\odot}\)model starting somewhat earlier than the more massive model. In Figure 9 we see that dust formation also stops first for the 1.7 M\(_{\odot}\)model. Dust stops forming first and disperses first in the equatorial direction for the 3.7 M\(_{\odot}\)model, similarly to how it started. The opposite is the case for the 1.7 M\(_{\odot}\)model, where dust formation stops first and disperses first in the polar directions.
These behaviours are related to the lower envelope mass of the 1.7 M\(_{\odot}\)model as well as its larger mass ratio which, as we saw in [15], leads to a great deal more asymmetry in the dust shell and ejecta, such as larger radial velocities in the polar direction.
Figure 9: Grey optical depth of the dust shell, \(\tau_\mathrm{d}\), for the 1.7 M\(_{\odot}\) (left panel) and the 3.7 M\(_{\odot}\) (right panel) simulations, observed from the three viewports. \(\tau_\mathrm{d}\) starts to drop after the first few decades, as expected by the free expansion model of the dust shell. Curves in the corresponding colours (dashed lines) use Equation 13 to fit data on the right of the grey vertical line and extrapolate the time when the dust becomes transparent (see Table 1)..
4 Discussion↩︎
Here, we discuss the uncertainty deriving from finite simulation resolution, particularly at the surface, as well as what observations can tell us about our simulation efforts so far.
4.1 SPH simulation resolution test↩︎
To determine the effect of resolution (i.e., the total number of SPH particles) on the resulting lightcurve, we conducted the same analysis described in Section 3.1 but using the 3.7 M\(_{\odot}\) simulation carried out with 200 000 particles by [15]; cf. the 1.3 million SPH particles of our production simulations. The two simulations have a number of SPH particles that differ by a factor of \(\sim\) 7, just shy of the factor of 8 (\(2^3\)) that achieves a factor of 2 in linear resolution.
The results presented by [27] for exactly the same two simulations, were not significantly impacted by their different resolutions, as demonstrated by them - see, e.g., their figure A1. However, in Figure 10 we show the two lightcurves that result from the two simulations with different resolutions, showing a significant difference, with the higher resolution simulation overall dimmer by a factor of 2-3. At \(t=0\), we see that the luminosity of the high resolution simulation is lower than the low resolution (by a factor of 3.3) and closer to the input value from the 1D stellar structure. The approximate shape of the luminosity peak and trough is similar, but the specific timing and brightness of the peak, the shape of the descent to the trough, and even the presence of another secondary peak are not converged. Unfortunately, a meaningfully higher particle simulation would substantially increase the computation time, and it is out of the scope of the current work.
The lightcurve uncertainty (shaded areas; explained in Section 4.2 and Appendix 12) is smaller in the high-resolution (blue curve), as expected. However, the shaded areas do not fully account for the difference between the two curves. This suggests that the error is not only due to the low surface resolution. In Appendix 11 we demonstrate that the imperfect reproduction of the surface layer, i.e., the specific position of the relatively few SPH particles at the surface of the star do impact the light calculation at some level, while it does not impact other simulation results, such as final orbital separation, unbound mass or dust mass, as already demonstrated by [17] and [15]. This problem will be greatly ameliorated by higher surface resolution.
[28] carried out a resolution test with the number of SPH particles changed by a factor of three (99 955 vs. 299 929). They remark on some notable differences between the two lightcurves generated with the two resolutions. Ultimately a convergence test is extremely difficult to execute given prohibitive computational times. Only so can convergence be tested.
4.2 Uncertainty from low resolution at the surface↩︎
In this Section, we first find a way to formally quantify the accuracy of the lightcurve and then discuss a way to generate error bars for the lightcurve itself, which are then shown as shaded areas in Figures 5, 12, 10 and 18.
At the beginning of the simulation, the gas distribution is entirely optically thick. The photosphere is therefore located somewhere within one smoothing length (\(25-30\) R\(_{\odot}\)) of the outermost SPH particle encountered by the ray. More importantly, at this location, the temperature gradient is typically large, so that the error on the photospheric temperature is even larger, with a concomitant large error on the luminosity on a per-ray basis. As the common envelope gas expands, the outer layer’s density decreases and the gas may become optically thin (at least before dust forms). This means that the photosphere recedes into the gas distribution, the local smoothing length decreases, which increases the precision with which we can determine the photosphere’s location and temperature.
To first approximation, we can think of the precision with which we can estimate the location and properties of the photosphere by estimating the number of particles, \(N_\mathrm{res}\), used for the luminosity integration in Equation 8 . This can be determined by considering each ray. Instead of integrating across all rays to determine the contribution to \(L\) that each particle \(j\) makes (as we did in Equation 8 ), we first separate the contribution to the luminosity from each particle \(j\) and ray \(i\) (\(L_{ij} \equiv 4 \pi S_{j} A_{\mathrm{eff}, ij}\)). We then select the particle that contributes the most to the luminosity of that pixel. Finally, we calculate the ratio of that particle’s contribution to the total brightness of that pixel. In other words, we find the fraction of the highest contribution to the total effective pixel area \(A_{\mathrm{eff}, i} \equiv \sum_j{A_{\mathrm{eff}, ij}}\) for a given ray \(i\) \[f_{\mathrm{contr}, i} \equiv \frac{\max_j{A_{\mathrm{eff}, ij}}}{\sum_j{A_{\mathrm{eff}, ij}}} , \label{eq:def:F95contr}\tag{14}\] where max\(_j\) is the maximum area \(A_{{\rm eff},ij}\) contribution among all particles, \(j\), for a given ray, \(i\). This shows the extent to which one particle dominates the integration of a particular ray. The lower \(f_{\mathrm{contr}, i}\), the higher the resolution, because it denotes that information for that ray was obtained from a large number of particles.
The inverse of this fraction (\(N_{\mathrm{res}, i} \equiv f_{\mathrm{contr}, i}^{-1}\)) provides an estimate of the number of particles used to resolve the specific intensity per pixel \(i\). Averaging over pixels, using their luminosity contribution, \(L_i\), as weights, we can estimate the SPH resolution of our calculated luminosity, \(L\), in terms of the number of particles per ray \[N_\mathrm{res} \equiv \frac{ \sum_i{N_{\mathrm{res}, i} L_i} }{ \sum_i{ L_i} },\] where \(L_i \equiv 4 \pi \sum_j S_j A_{\mathrm{eff}, ij}\). The quantity \(f_{\mathrm{contr}, i}\) is plotted in Appendix 10.2, similarly to how the surface brightness was plotted in Figure 4, for easy comparison, while in Figure 11 we plot \(N_{\rm res}\) as a function of time. Note that by definition \(N_\mathrm{res}\) cannot be lower than unity. The larger the number, the better the resolution, the higher the confidence we have on the luminosity output. In Figure 11, for both dusty simulations, the resolution begins at 1 to 2, indicating an extremely low resolution where one or two particles dominate the result for each ray; however, after 2 to 4 years, \(N_\mathrm{res} \sim 5\) growing to a value of \(\sim 20\), meaning that information for the specific intensity of each ray is sourced from multiple particles and thus modelled more accurately.
Our resolution criterion, \(N_\mathrm{res}\), is similar to the SPH ray tracing relative accuracy metric, \(\eta\), defined in equation 27 of [28], but distinct from that work in that ours measures the resolution by the actual contribution to the luminosity by each particle for each ray, instead of by how much the number of neighbouring particles is lower than the average. Our method has the advantage that it summarises the result into a single number for each ray, instead of for each point on the ray; our method also has the advantage that it highlights the influence on the particles while ignoring those particles with very low optical depth and temperature (which don’t matter to the final luminosity, even if they have a very low number of neighbours). Our criterion is also similar to that of [29], where pixels with particles of optical depth \(d\tau > 1/3\) present before \(\tau\) reaches their photospheric value of 1 are flagged as bad pixels; once again, our method provides a continuous assessment of photospheric resolution instead of a boolean judgement, and it is weighted by the particle’s specific intensity, which matters directly to the luminosity error estimation.
Figure 11: An estimate of the precision of the luminosity calculation, based on the number of particles, \(N_\mathrm{res}\), used in the calculation of the specific intensity for each ray, for the 1.7 M\(_{\odot}\) simulation (the 3.7 M\(_{\odot}\) simulation is similar). Blue, orange, and green lines are averaged over a grid of \(256 \times 256\) rays observed from the three three directions..
An alternative way to quantify the surface uncertainty, which also allows us to calculate approximate error bars on the value of the calculated luminosity, uses the inherent spherical symmetry of each SPH particle’s smoothing kernel, which, when the SPH particle is located at the photospheric surface, misrepresents the steepness of the gradients in all quantities (e.g., density, pressure, temperature).
When the information on the light comes from one or a few optically thick surface particles, the kernel interpolation that assumes Gaussian profiles, overestimates the source function on the observer side of the surface layer, because the smoothing length of surface particles extends \(h\) in all directions. However, at these surface locations we have a prior, in that we know that all quantities have a sharp negative gradient.
To quantify this problem, we start by treating the current lightcurve as an upper limit of the true lightcurve, considering that the SPH kernel interpolation at the first locations on each ray where quantities are non-zero would overestimate the specific intensity contribution. This is because the surface particles are statistically distributed at the surface. However, each ray encounters the outer half of each particle first, starting the integration of the source function systematically \(2h\) towards the observer side of the SPH particle. This means that each ray picks up an extra contribution to the source function, making the luminosity too large.
We then locate the depth point below which the resolution is high enough to make the calculation of the source function and optical depth accurate, and outside of which the resolution is considered too low. This depth is defined as follows: a ray passing this point from the outside, has counted 14.5 contributing particles to the lightcurve calculations (this number is a quarter of 57.9, the average number of neighbours in our phantom simulations). A detailed explanation of how we select the boundary depth is provided in Appendix 12.
To determine a lower limit to the luminosity, we assume all the particles’ above the point do not contribute any specific intensity to the ray (i.e., we set their temperature to 0 K), while their optical depth, \(\tau\), is left unchanged. In this way, we simulate the opposite extreme with the source function at zero, while the \(\kappa \rho\) has the normal behaviour. The results are the shaded areas shown in Figure 5.
A further source of systematic uncertainty resides in the imperfect reproduction of the star in 3D. The initial surface temperature and density profile at \(t=0\) after relaxation in phantom is different from the prescription of mesa. This has little to no effect on the hydrodynamics, but plays a role in the lightcurve calculation. Together with the low surface resolution, this contributes to the initial lightcurve unreliability. In Appendix 11 we carry out several tests that quantify the error in the lightcurve, demonstrating that it reduces greatly after time zero, even if, interestingly, and only for the 1.7 M\(_{\odot}\)model, we see a change of 3-6 months in the \(\sim\)3-yr width of the luminosity peak.
4.3 The lightcurve between 20 and 45 years and the adiabatic approximation↩︎
These simulations are adiabatic. As such, the luminosity towards the end of the simulation plateaus at \(\sim\)100 000 L\(_{\odot}\), a very high value that is at odds with observations. Typical extreme carbon stars’ luminosities are in the \(\sim\)10 000 L\(_{\odot}\) [30], although some extreme cases are found with luminosities as high as \(\sim 80\,000\) L\(_{\odot}\)[31].
An object radiating 100 000 L\(_{\odot}\)over 30 years (the approximate length of the plateau) radiates \(\sim 4\times 10^{47}\) ergs. This value is an order of magnitude larger than any of the typical energies describing the star and the interaction: the envelope’s gravitational potential energy, \(\sim -6\times 10^{46}\) ergs, is overcome by the delivered orbital energy, \(GM_c M_2 / 2 a_f - GM_1 M_2 / 2 a_i \approx 10^{46}\) ergs, and some of the envelope’s internal energy (particularly the recombination energy), \(\approx 6 \times 10^{46}\) erg. We also note that the undisturbed star, with 5180 L\(_{\odot}\)would radiate \(\approx 2 \times 10^{46}\) erg in 30 years.
The point is that in Nature only \(1-2\times10^{46}\) ergs would be available to radiate over that time, lower than implied by the deduced luminosity during the late plateau. If the surface were allowed to radiate in our simulations, the likely outcome would be a smaller, cooler, less luminous surface. This discrepancy is most obvious at late times. At earlier times (within the inspiral timescales of 5-10 years) the energetics are likely dominated by the orbital energy injection and subsequent delivery of recombination energy at depth powering an adiabatic expansion.
4.4 Comparison with observations↩︎
Figure 12: Absolute AB magnitude in Johnson V (solid), Cousins R (dotted), Cousins I (dash-dotted) and J bands (dashed) for the dusty 1.7 M\(_{\odot}\) (left panels) and 3.7 M\(_{\odot}\) (right panels) simulations, as viewed from polar (\(+z\), top panels) and equatorial (\(+x\), bottom panels) directions. The shaded area are uncertainties from lower surface resolution, the same as in Figure 5—See Section 4.2 for more details..
| Observations | Observations | Our 1.7 | Our 3.7 | |
| OGLE-2002 | AT 2025abao | simulation | simulation | |
| Progenitor | ||||
| ZAMS mass () | \(<4\) | \(7\pm2\) | 2 (1.7@AGB) | 4 (3.7@AGB) |
| Radius () | 30 | 350 | 260 | 330 |
| Luminosity () | 290 | 10 000 | 5180 | 11 900 |
| \(T_{\rm eff}\) (K) | 4350 | 3100 | 3227 | 3317 |
| \(M_2\) () | Unknown | Unknown | 0.6 (\(q=0.35\)) | 0.6 (\(q=0.16\)) |
| At peak | ||||
| time (yr) | 1.6\(^a\) | 0.099\(^b\) | 1-2.5\(^c\) | 2.5-6\(^d\) |
| \(R_{\rm phot}\) (au) | 1.1\(^e\) | 4.1 | \(5-6\) | \(6-3\) |
| \(L\) () | 8 100\(^f\) | \(7.2\times10^5\) | \(2-4\times 10^6\) | \(2-4\times 10^6\) |
| \(T_{\rm eff}\) (K) | 3600 | 5740 | 6000-7000 | 6000-7000 |
| \(v_{\rm exp}\) (km s\(^{-1}\)) | – | 450 | 86 | 85 |
| At 8 years post outburst start | ||||
| \(R_{\rm phot}\) (AU) | 50 | – | 50 (dust) | 2 (stellar) 50 (dust) |
| \(L\) () | 8 200 | – | \(\sim 1\times10^5\) | \(30\,000-2\times10^5\) |
| \(T_{\rm phot}\) (K) | 550 | – | 800-1000 | 500-6000\(^g\) |
| \(v_{\rm exp}\) (km s\(^{-1}\)) | 35 | – | 60 | 68 |
| \(M_{\rm dust}\) | \(>1.4\times 10^{-5}\) | – | \(1\times10^{-4}\) | \(1\times10^{-4}\) |
| \(^a\)Time of the second peak in April 2004 (MJD 53100), after outburst start on MJD 52500. | ||||
| \(^b\)36 days between the start of the rise selected by us to be at MJD 60950 to \(g\) band maximum at | ||||
| MJD 60986.3 [32]. | ||||
| \(^c\)Peak time for both viewports. | ||||
| \(^d\)Peak time for equatorial viewport. First peak for polar viewport is \(\sim\)3yrs. | ||||
| \(^e\)Measured in October 2003 instead than April 2004 - see table 1 in [7]. | ||||
| Both models’ photospheric radii are by then dominated by dust, while in the observations the | ||||
| measured photosphere is that of the gas, where a larger, 23 au, 800 K dust shell is present but optically thin. | ||||
| \(^f\)Using the distance of 8.2 kpc of [7]. Luminosity scales with the square of the distance, | ||||
| while radii scale as the distance. Note that [33] use a distance of 4.09 kpc. | ||||
| \(^g\)Viewing angle-dependant, with polar being hotter. | ||||
A class of astrophysical transients called LRNe have been identified as mergers of non-compact stars [2]. Typical lightcurves of such transients exhibit a peak and plateau or a multi-peak behaviour with a subsequent dimming, reddening, and quick dust production [14], [33]–[36]. The transients identified in our Galaxy generally corresponded to lower mass donors, while extragalactic, brighter events likely correspond to more massive objects, as shown by the correlation between progenitor mass and peak brightness [37], [38]. Archival observations of pre-outburst progenitor systems also revealed that the majority of donor stars are subgiants, which are stars in the very early post-main-sequence, experiencing a fast, thermal-timescale, expansion phase. An often cited example is V 1309 Sco, for which observations were obtained during the entire pre-outburst and outburst phases such that this merger can reasonably be inferred to have been that of a low mass, 1.5 M\(_{\odot}\), 4 R\(_{\odot}\) star with a 0.15 M\(_{\odot}\)companion [5]. This object, along with several others such as V838 Mon (a more massive merger; [39]), took place on faster timescales compared to those with RGB/AGB progenitors, because of the shorter dynamical times of the more compact, dense envelopes.
Among the observed red nova sample, only two examples are compatible with bona fide giant star progenitors, as opposed to subgiant or main-sequence progenitors. The red nova OGLE-2002-BLG-360 [7] was characterised by an outburst that took place over years, as opposed to days and weeks. It brightened by a factor of 30 (3.7 mag) in the \(I\) band to a first peak, and then by another 2 mag by the third peak, for a total rise time of 2.2 yrs (where we take the beginning of the rise time to be on MJD 52500, corresponding to 14 August 2002; [7]). The object then steadily reddened and dimmed by \(\sim\)10 mag over the following 5 years [33].
A second object, WNTR23bzdiq/WTP19aalzlk/AT 2025abao discovered in M 31, was initially presumed to represent a full common envelope ejection from a more massive AGB star [40]. The lightcurve activity was initially interpreted as a long and shallow outburst. However, in October 2025 the object experienced another outburst, showing a quick brightening over about a month of \(\sim\)5 mag in the optical, to a peak bolometric luminosity of 870 000 L\(_{\odot}\)4, analogous to other LRNe [32], [41]. Hence, the initially observed slow brightening was, in fact, just pre-outburst activity preceding a red nova [2]. [32] analysed the lightcurve for 120 days past peak, revealing a cooling photosphere (from 5740 K at peak to 3600 K 120 days later), a declining luminosity (down to 440 L\(_{\odot}\)at 120 days past peak) and an expanding photosphere (10 au at 70 days post peak, but reducing to 8.3 by 120 days). No dust is observed around this object, likely because it is too soon.
In Table 2 we compare several parameters that could be extracted from the observations of OGLE-2002 and AT 2025abao with equivalent parameters of our models, noting that the models were not designed to fit these objects, but these objects are the closest we have to the models and can therefore be very instructive. Our two models have primary star (progenitor) masses that are in the ball park of those deduced for OGLE-2002 and AT 2025abao, but are more similar to the progenitor of AT 2025abao than the OGLE-2002’s progenitor, which is likely a smaller RGB star.
In Figure 12 we present the apparent magnitudes in the Johnson-Cousins \(V\), \(I\) [42] and 2MASS \(J\) bands [43] obtained by convolving \(L_\lambda\) with filters’ sensitivity curves . Spectra are computed using the pyphot package [44], using the filter table of ground_johnson_V for the \(V\) band, ground_cousins_I for the \(I\) band and 2mass_J for the \(J\) band. We elect not to overlay the observations on these curves because the actual similarities are misleading given the misgivings we still have with our predictions. However, we notice that, particularly for AT 2025abao, which is the more similar progenitor to our more massive simulation, the peak \(J\) band brightness of 12.6 mag, translates to an absolute \(J\) band brightness of \(M_J = -11.9\) (at a distance of 0.785 Mpc, \(A_V=0.17\) adopted by [32] and \(A_J/A_V=0.282\) of [45]), a value that is indeed very similar to our simulation over the first few years.
The time of peak in Table 2 is arbitrarily selected to correspond to the second of three peaks (at MJD 53 100) in the lightcurve of [7], and it is with reference to the start of the outburst for the OGLE-2002 observation, which is selected by us as the time when the \(I\) band brightness starts to increase (MJD 52 500 in the lightcurve of [7]). For AT 2025abao, we use as reference the base of the rise at MJD60950 and \(g\) band peak at MJD60986 - a short 36 days later. OGLE-2002 took a year and a half to rise to peak, while AT 2025abao had a rise that was quite a lot faster, only a month. For the models, time zero is defined as the time of the start of the simulations, and the timing of the luminosity peak is defined as the time elapsed form this moment until reaching maximum brightness. We note that while for the 1.7 M\(_{\odot}\)model this is quite well defined, and similar for all viewports, for the 3.7 M\(_{\odot}\)model it is more dependent on the viewport (the range is 2.5\(-\)5 yrs). This said we have little information on the actual rise as the error bar due to low surface resolution gives us little knowledge over the first year.
At peak, the stellar/gas photospheric radius is of the order of 1 au for OGLE-2002, apprximately 4 au for AT 2025abao, while for the simulations it is approximately 5 au, where we have used the "observer’s" radius for the comparison (Section 3.2), and we note that at the selected time of peak the dusty shell in the models has already formed, but it is still optically thin. The OGLE-2002 luminosity peaks at a value of 8100 L\(_{\odot}\), while for AT 2025abao it is much larger, 720 000 L\(_{\odot}\), and more similar to the simulations, with values ranging between 1 and \(2\times 10^6\) L\(_{\odot}\). Another LRN, M31-2015-LRN had a luminosity peak of 309 000 L\(_{\odot}\) (Table 3) and a progenitor mass of 3-5 M\(_{\odot}\)not dissimilar from our more massive model but likely more compact.
The next comparison is made approximately 8 years after the start of the outburst (or the start of the simulations), but only for OGLE-2002. At this time the photospheric radii of observation and 1.7 M\(_{\odot}\)simulation are dominated by dust and are similar (50 au). The 3.7 M\(_{\odot}\)model has a still dominating gas photosphere at 2 au, while a dust shell is already present with a radius of 50 au. Similarly, the photospheric, dust, temperatures are similar in the observation and 1.7 M\(_{\odot}\)simulation (500\(-\)750 K), while for the 3.7 M\(_{\odot}\)simulation the dominance of dust or gas photosphere depends on the viewport: it is similar to the other model for the \(+x\) and \(+y\) views, but for the \(+z\) view, the SED is still dominated at this time by the star (\(T_{\rm phot}\sim 6000\) K), and the dust only fully obscures it at \(\sim\)12 yrs; Figure 7). The mass of the dust measured by [7] is a lower limit (\(>1.4\times10^{-5}\) M\(_{\odot}\)) due to its optically thick nature. In the simulations, the dust mass that has formed by 8 years is still relatively low but about to increase: it will reach the level of the measurement of [33] at \(\sim\)20 yrs after the start of the interaction [15].
Finally, in Table 3 we list observationally-derived parameters for several additional red novae remnants that should be specifically pursued with simulations, such as M31-2015-LRN, which is similar to OGLE-2002, with a progenitor mass and radius of 3–5.5 M\(_{\odot}\)and 30–40 R\(_{\odot}\), respectively [25] and a typical optical lightcurve evolution over \(\lesssim 1\) yr [14], [46]. Observationally-derived quantities such as these need to be correctly interpreted. During the optical outburst, when the visible optical depth is low, the observed SED reveals the outburst’s photosphere, which then leads to a measurement of its photospheric radius, effective temperature, and bolometric luminosity. Several years after the optical peak, dust formation progressively obscures the central remnant, reprocessing the stellar luminosity into longer wavelengths. Hence, at late times, the dust radius, composition, density profile, and temperature define the photosphere, and the parameters of the merged remnant star can only be derived indirectly, via radiative transfer modelling of the dust.
| Prog. | Evol. | Time\(^a\) | \(R_\star\) | \(L_\star\) | \(T_{\rm eff}\) | \(v_{\rm exp}\) | \(R_{\rm dust}\) | \(T_{\rm dust}\) | \(M_{\rm dust}\) | \(\tau_{\rm V,dust}\) | Ref. | |
| mass | stage | |||||||||||
| () | (yr) | (au) | () | (K) | (km s\(^{-1}\)) | (au) | (K) | () | ||||
| OGLE-2002- | \(<\)4 | RGB | 1\(^b\) | 1.1 | 8100 | 3600 | – | 20 | 800\(^c\) | 5.5\(\times 10^{-7}\) | 2.5 | (1) |
| -BLG-360 | 8\(^b\) | 1.4 | 8200 | 3200 | 35 | 50 | 550 | 3.44\(\times 10^{-3}\) | \(>20\) | (1,2) | ||
| 15-21\(^b\) | 0.70 | 1646\(^d\) | 3000 | 226 | 1000\(^e\) | 50-200 | 1.2\(\times 10^{-2}\) | – | (2) | |||
| M31-2015- | 3–5 | post-MS | 0.12 | – | 3.09\(\times 10^5\) | 3390 | 100 | 23 | 1760 | 1\(\times 10^{-5}\) | 0.7 | (3) |
| -LRN | 2.5 | – | 69 000 | 9780 | 100 | 38 | 1110 | 1\(\times 10^{-2}\) | 287 | (3) | ||
| 10.0 | – | 4620 | 6000 | 100 | 180 | 385 | 7.7\(\times 10^{-4}\) | 100 | (4) | |||
| AT 2025abao | \(6-9\) | AGB | 0 | 4.1 | \(7.2\times10^5\) | 5740 | – | – | – | – | – | (5) |
| 0.19 | 10 | 1100 | 4000 | – | – | – | – | – | (5) | |||
| 0.33 | 8.3 | 440 | 3600 | – | – | – | – | – | (5) | |||
| AT 2018bwo | 12–16 | post-MS | 0 | 25.8 | 3.41\(\times 10^6\) | 3330 | 220 | – | – | \(<\)3.1\(\times 10^{-3}\) | – | (6,7) |
| 6.1 | – | 1.48\(\times 10^5\) | 2600 | 100 | 140 | 670 | 2.1\(\times 10^{-5}\) | – | (4) | |||
| AT 2021blu | 13–18 | post-MS | 0 | – | 1.698\(\times 10^7\) | 8800 | 500 | – | – | – | 22 | (6) |
| 3.1 | – | 3.30\(\times 10^5\) | 2200 | 100 | 87 | 1100 | 4.2\(\times 10^{-5}\) | 22 | (4) | |||
| AT 2021biy | 17–24 | post-MS | 0 | 7.34 | 4.180\(\times 10^7\) | 11 450 | 640 | – | – | – | – | (8) |
| 3.3 | – | 4.20\(\times 10^5\) | 2400 | 250 | 240 | 700 | 3.6\(\times 10^{-5}\) | 22 | (4) | |||
| AT 2011kp | 20–40 | post-MS | 0 | – | 2.612\(\times 10^7\) | 16 000 | 450 | – | – | – | – | (9) |
| 12.5 | 2.55 | 41 000 | 3500 | – | 291 | 450 | 1.8\(\times 10^{-3}\) | 133 | (10) | |||
| AT 1997bs | 20\(-\)60 | post-MS | 0 | – | 8.62\(\times 10^6\) | – | 585 | – | – | – | – | (10,11) |
| 7 | 3.8 | 1.52\(\times 10^5\) | 4000 | – | 301 | 600 | 1\(\times10^{-3}\) | 69 | (10) | |||
| (1) [7]; (2) [33]; (3) [14]; (4) [36]; (5) [32]; | ||||||||||||
| (6) [47]; (7) [38]; (8) [48]; (9) [49]; (10) [32]; (11) [50]. | ||||||||||||
| \(^a\)Post peak unless stated otherwise; | ||||||||||||
| \(^b\) Time post outburst’s start on MJD 52 500 or August 14, 2002. | ||||||||||||
| \(^c\)This is assumed, but the value at 3-4 yrs after outburst start was 850 K. | ||||||||||||
| \(^d\)This luminosity is for a distance of 4.09 kpc. For the distance assumed by [7], 8.2 kpc, the luminosity would be 6600 ; | ||||||||||||
| \(^e\)Outer radius of the inner shell; outer shell radius is 9500 au. | ||||||||||||
5 Conclusions↩︎
We have calculated the lightcurves over 44 years for two simulations representing 1.7 M\(_{\odot}\)and 3.7 M\(_{\odot}\)AGB stars undergoing common envelope evolution with a 0.6 M\(_{\odot}\)companion, with dust nucleation. The lighturves are obtained in post-processing, and the simulations are adiabatic. Despite a series of shortcomings, we can draw a number of (primarily qualitative) conclusions.
The lightcurve for both models (Figure 5) has the following characteristics. (a) A peak at about \(2-4\times 10^6\) L\(_{\odot}\), about \(1-5\) yrs after the start of Roche lobe overflow, depending on model and viewing angle. (b) The light rise is driven by the expansion of the photosphere. (c) The effective temperature plays a more minor role: besides a cooling tendency due to adiabatic cooling, fluctuations are due to optical depth effects and on balance the temperature remains in the range 6000-7000 K during the light rise (Figure 7). (d) The broader peak for the more massive model is likely due to its lower mass ratio (\(q=0.16\) vs. 0.35), resulting in a weaker interaction leading to a delayed dust nucleation. (e) The bolometric luminosity peak subsides as dust formation generates an optically thick, expanded, and cool photosphere. Subsequent minor variations—peaks and troughs—result from the complex, asymmetrical interplay between gas and dust layers, which dynamically exposes or hides deeper surface layers. (f) After a trough, which has a temperature of \(\sim 6000\) K, the luminosity settles on a value of \(100 \, 000 - 200\,000\) L\(_{\odot}\)and a temperature that gradually decreases to \(\sim 400\) K at 20 yrs. The luminosity is too high on account of the lack of cooling in our simulations.
The light rise to the peak is driven in our simulations by the initial phase of photospheric expansion. This expansion is caused by the early inspiral, but also by stellar surface instability as well as the abrupt insertion of the companion in the computational domain, both numerical effects. In addition, even ignoring the numerical reasons for stellar expansion, the timescale that comprise the early Roche lobe overflow, the on-set of mass transfer instability and the eventual start of the inspiral is likely faster in simulations than in Nature, as is demonstrated by the fact that the timescale of the pre-inspiral lengthens with increasing resolution and is not converged (e.g., [9], [17]).
An optically thin dust shell emerges above the warmer, gasseous photosphere between 3 and 5 years after the start of the simulations, at which point the SED shows the emergence of a second, cooler peak (Figures 6 and 7 where we see the temperature drops below 1000 K). The dust shell becomes optically thick by 20 years and it continues to cool and expand. The dust will disperse sufficiently to become optically thin between 100 and 200 years later (Figure 9), at which point another, hotter, inner photosphere would be revealed.
Equivalent models with no dust show a very different lightcurve with a gently increasing luminosity throughout the simulation—only the lower mass model shows a turn-over of the bolometric luminosity, while for the higher mass, lower \(q\) model, the calculation was carried out for too short a time to see the turn-over.
The primary source of uncertainty stems from insufficient spatial resolution at the photosphere. This issue is particularly severe at the onset of the simulation; consequently, we conclude that we lack meaningful information for the entire first year of both simulations (Figures 5 and 12). Adequately resolving the photosphere—specifically, ensuring that particles at the nominal photosphere have a resolution length, \(h\), smaller than the photospheric depth—would require approximately 340 million particles for the 1.7 M\(_{\odot}\) simulation and 3.8 billion particles for the 3.7 M\(_{\odot}\) model. Such high resolution is currently computationally prohibitive. A possible future solution is to use the novel adaptive particle refinement method of [51].
Our simulations are adiabatic and, as such, result in too bright a luminosity, particularly in the last \(\sim\)30 years of the simulation, when the bolometric luminosity is approximately constant at \(\sim\)100 000 L\(_{\odot}\). During that time the object would radiate \(\sim 4\times 10^{47}\) erg or \(\sim\)10 times the initial internal energy of the star. Such radiation, would contribute to a smaller cooler and less luminous photosphere. The solution is adopting a radiative cooling approach, e.g., [52].
We carried out a preliminary comparison with the observed lightcurve of OGLE-2002-BLG-360 and of AT 2025abao, two LRNe with progenitor masses similar to our more massive model. Both are likely derived from giant progenitors, the former from a \(\sim\)30 R\(_{\odot}\), RGB star, the latter from a \(\sim\)350 R\(_{\odot}\) AGB star quite similar to our models. The timing of the peaks measured from the start of the light rise or the start of the simulations is similar for OGLE-2002 (1.6 vs. 1-6 yrs for the simulations), but much shorter for AT 2025abao (36 days) with photospheric sizes that, for OGLE-2002 start smaller (at peak it is 1.1 R\(_{\odot}\) vs. 3-6 R\(_{\odot}\) for the simulations), but becomes quite comparable at 8 years (50 R\(_{\odot}\) for the dust photosphere), while for AT 2025abao it is within range of the simulations (4.1 au) at one year (no data after that). The peak luminosities are very discrepant for OGLE-2002 (8100 L\(_{\odot}\)vs. \(2-4\times 10^6\) L\(_{\odot}\)for the simulations), but much closer for AT 2025abao (720 000 L\(_{\odot}\)). The mass of the dust is comparable between our simulations and OGLE-2002, while AT 2025abao may not have had time to make any dust yet (at the time of writing only 8 months have passed from the start of the outburst).
Bridging the gap between simulations and observations requires an ability to model the onset of CE phases in simulations as well as resolving two primary limitations identified in this study: low photospheric resolution, which chiefly degrades the light curve during the first year, and the adiabatic assumption, which overestimates the luminosity in later phases. The first issue can be mitigated using the adaptive particle refinement method developed by [51], an approach currently under development. The second limitation can be addressed by implementing radiative cooling techniques, such as those demonstrated by [52].
Acknowledgements↩︎
This research was supported by the Commonwealth through an Australian Government Research Training Program Scholarship [DOI: 10.82133/C42F-K220]. CM acknowledge the use of the python package Sarracen [53] used to post-process phantom data dumps. OD, LS, MK and LB acknowledge support through the Australian Research Council discovery programme grant DP210101094. LS is senior research associates from F.R.S.- FNRS (Belgium). ML acknowledges funding by the European Union (ERC, ExCEED, project number 101096243). NB acknowledges being funded by the European Union (ERC, CET-3PO, 101042610). Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. The publication is part of the research project PID2024-155585NA-I00 funded by MICIU/AEI /10.13039/501100011033. NB acknowledges financial support from grant CEX2024-001451-M funded by MICIU/AEI/10.13039/501100011033.
This work was supported in part by Oracle Cloud credits and related resources provided by Oracle for Research. This work was supported by the OzSTAR national facility at Swinburne University of Technology. OzSTAR is funded by Swinburne and the Australian Government’s National Collaborative Research Infrastructure Strategy (NCRIS).
6 Data Availability↩︎
Source codes are available at https://github.com/chunliangmu/clmuphantomlib and https://github.com/chunliangmu/RTinCEE-scripts-2024. (All code and movies will be uploaded to Zenodo after submission.)
7 Light Calculation Details↩︎
In Figures 13 and 14 we give a graphical representation of the post-processing calculation.
Following the description in Section 2, the calculations of the light is parallelized across rays using the python package numba’s jit(parallel=True) and prange() functions, with the following steps:
The coordinate system was rotated so that the observer is always at \(z \rightarrow +\infty\). Currently, this is achieved by simply re-arranging the axis labels, so we only have 3 view ports \(+x, +y, +z\).
A grid of parallel rays pointing towards the \(+z\) direction were generated to sample the \(xy\)-plane. The size of the grid in the \(x\) and \(y\) directions includes all SPH particle positions \(\pm\) the kernel radius (\(2h\)) along that axis. In other words, it ranges in \(\left[ \min{\{x_j - 2 h_j | \forall j\}},\; \max{\{x_j + 2 h_j | \forall j\}} \right]\), and same for the \(y\) direction.
We sort all SPH particles by their \(z\) coordinates and go through them for each ray in a backward tracing fashion, i.e., from the ones with highest \(z\) (observer’s side) to the lowest \(z\) (the opposite side).
For each ray \(i\), and each SPH particle \(j\), we first check whether that particle lies within the range of of the ray, in other words, whether the distance from the particle \(j\) to ray \(i\) is smaller than twice its smoothing length \(2h_j\). The factor of 2 here comes from the smoothing kernel radius. We record a list of all relevant particles for each ray \(i\).
Then we compute the effective area \(A_{\mathrm{eff}, ij}\) for each particle \(j\) and ray \(i\), using Equation 7 . This is achieved by integrating along 1000 sample points (per particle per ray) across the segment of the ray that lies inside the \(2h_j\) radius of the particle. The contribution to each sample points from the neighbouring particles were also included when calculating the optical depth \(\tau_i(z)\) in Equation 10 .
If the optical depth becomes larger than the stop condition (see equation 18 below), the calculation for the \(i\)’th ray is stopped, as we know that further particles no longer make a meaningful contribution to \(I_i\).
The effective area \(A_{\mathrm{eff}, j}\) for each particle \(j\) is computed by summing up among the rays using Equation 7 .
The luminosity is then computed using Equation 8 .
Step [list:lcgen:para:s] to [list:lcgen:para:e] are parallelized across the rays. To save computational time, we implement a stop condition (expressed in Equation 18 ) once \(\tau_i(z)\) becomes sufficiently large such that further particles’ contributions to \(I_i\) are smaller than the rounding error. For simplicity, we begin by approximating the solution to the radiative transfer equations (Equation 4 ) with \[I_i = \sum_j S_j e^{-\tau_{ij}} (1 - e^{-\Delta \tau_{ij}}) , \label{eq:radiative-transfer-discretized-delta95taus}\tag{15}\] where \[\Delta \tau_{ij} \equiv \frac{m_j \kappa_j}{h_j^3} \int_{-\infty}^{+\infty} w(q_{ij}(z')) dz'. \label{eq:def:Delta95tau}\tag{16}\] is the optical depth contribution of particle \(j\) for ray \(i\), the integral in the equation above (the column kernel) is a pre-calculated, tabulated quantity, and \[\tau_{ij} \equiv \sum_{n=0}^{j-1} \Delta \tau_{in} \label{eq:def:tau95max95i}\tag{17}\] is the optical depth at the location of particle \(j\), which we are approximating by summing up the optical depth contribution from all particles in front of it. In other words, by invoking the above approximation, we are ignoring—for the purpose of deciding on the stop condition for the post-processing interpolation—how spread out the particles along the ray are; instead, we assume particles in front of \(j\) are fully in front (added to \(\tau_{ij}\)), and those behind are fully behind (ignored). In other words, here we are treating SPH particles as pre-integrated 2D circles instead of 3D spheres. This should be reasonable in the context of deciding the stop condition (which requires less precision).
To determine whether it is reasonable to stop at a given particle, \(k\), we follow the method below. Given that the source function \(S_j \geq 0\) and \(\Delta \tau_{ij} \geq 0\) remains true for all particles, we can infer that \[I_{i, k} \leq I_i \leq \sum_{j=0}^{N_i} S_j,\] where \[I_{i, k} \equiv \sum_{j=0}^k S_j e^{-\sum_{n=0}^{j-1} \Delta \tau_n} (1 - e^{-\Delta \tau_j})\] is the specific intensity for \(i\)-th ray integrated up to \(k\)-th particle. Therefore, the absolute error of using \(I_{i, k}\) to estimate \(I_i\) is \[I_i - I_{i, k} = \sum_{j=k+1}^{N_i} e^{-\tau_j} \left( 1 - e^{-\Delta \tau_j} \right) S_j\] for simulations with \(N_i\) particles relevant to ray \(i\), and that \[e^{-\tau_j} \left( 1 - e^{-\Delta \tau_j} \right) = e^{-\tau_k} e^{-\sum_{n=k}^{j-1} \Delta \tau_n} \left( 1 - e^{-\Delta \tau_j} \right) \leq e^{-\tau_k}.\] Since \(j > k\) and \(\Delta \tau_j > 0\), we can deduct that \[I_i - I_{i, k} \leq e^{-\tau_k} \sum_{j=0}^{N_i} S_j.\] Hence, the relative error of \(I_i\) is \[\frac{I_i - I_{i, k}}{I_i} \leq \frac{e^{-\tau_k} \sum_{j=0}^{N_i} S_j}{I_{i, k}}.\] When the right hand side is smaller than the machine precision \(\delta_\mathrm{tol} \approx 10^{-16}\), the relative error (represented by the left hand side of the equation) is smaller than the machine precision \(\delta_\mathrm{tol}\), therefore integrating any further will make no difference, i.e., when \[\tau_k \geq \ln{\sum_{j=0}^{N_i} S_j} - \ln{\delta_\mathrm{tol}} - \ln{I_{i, k}}, \label{eq:radiative-transfer-stop-cond}\tag{18}\] we can safely ignore all particles for ray \(i\) behind this point.
8 Integration order↩︎
We note that there are two ways of combining the quantities in Equation 5 and 10 : One first interpolates each individual component (\(S\), \(\kappa\), \(\rho\)) to the ray location and then multiplies them together; alternatively, one first multiplies the quantities for each particle (\(S \kappa \rho\)) and then interpolates the product to the ray (the method of this paper). These two approaches produce somewhat different results. For example, the lightcurve in Figure 5 produced by the multiply-then-interpolate approach is brighter by a factor of approximately 2, at time zero. In theory the two methods should approach each other when the resolution is higher, or, in other words when the SPH particles contributing to each location on the ray are closer together, with smaller smoothing lengths. It is not clear to us which approach is more correct. Our method (Section 2.2) allows us to clearly identify the times when we have little knowledge of the light and those where we have relatively good knowledge. Where the shaded areas are dominant in Figure 5, the lightcurve is virtually unknown, while where the shaded area disappears into the curve, i.e., where the upper and lower limits according to our definition converge, we have reasonable information on the light.
9 Location of the photosphere↩︎
In Table 4 we present the photospheric values that correspond to Figure 1.
| 1.7 | 3.7 | |||||||
| 0 | 3.0 | 12.1 | 44.4 | 0 | 3.0 | 12.1 | 44.4 | |
| (yr) | (yr) | (yr) | (yr) | (yr) | (yr) | (yr) | (yr) | |
| \(\tau_\mathrm{ph, +z}\) | \(0.67^{+30.22}_{-0.67}\) | \(0.67^{+1.79}_{-0.54}\) | \(0.67^{+2.96}_{-0.46}\) | \(0.67^{+0.57}_{-0.35}\) | \(0.67^{+494.09}_{-0.67}\) | \(0.67^{+21.34}_{-0.58}\) | \(0.67^{+3.25}_{-0.61}\) | \(0.67^{+0.38}_{-0.26}\) |
| \(\tau_\mathrm{ph, +y}\) | \(0.67^{+15.49}_{-0.67}\) | \(0.67^{+3.56}_{-0.54}\) | \(0.67^{+1.05}_{-0.42}\) | \(0.67^{+0.31}_{-0.20}\) | \(0.67^{+725.51}_{-0.67}\) | \(0.67^{+0.62}_{-0.41}\) | \(0.67^{+3.20}_{-0.52}\) | \(0.67^{+0.52}_{-0.33}\) |
| \(\tau_\mathrm{ph, +x}\) | \(0.67^{+5.94}_{-0.67}\) | \(0.67^{+2.55}_{-0.62}\) | \(0.67^{+1.27}_{-0.55}\) | \(0.67^{+0.32}_{-0.19}\) | \(0.67^{+106.15}_{-0.65}\) | \(0.67^{+0.92}_{-0.37}\) | \(0.67^{+2.75}_{-0.35}\) | \(0.67^{+0.61}_{-0.35}\) |
| \(\kappa_\mathrm{ph, +z}\) / \(\mathrm{cm^{2}\,g^{-1}}\) | \(0.09^{+0.12}_{-0.06}\) | \(3.4^{+2.4}_{-1.9}\) | \(5.6^{+4.0}_{-1.6}\) | \(5.3^{+0.6}_{-0.9}\) | \(2.8^{+1.4}_{-0.0}\) | \(0.02^{+0.29}_{-0.02}\) | \(6.4^{+5.0}_{-3.3}\) | \(4.4^{+0.3}_{-0.4}\) |
| \(\kappa_\mathrm{ph, +y}\) / \(\mathrm{cm^{2}\,g^{-1}}\) | \(0.12^{+0.11}_{-0.00}\) | \(7.3^{+2.5}_{-0.4}\) | \(6.2^{+1.0}_{-0.8}\) | \(4.0^{+0.6}_{-0.6}\) | \(4.4^{+0.7}_{-0.0}\) | \(4.2^{+0.0}_{-3.2}\) | \(5.4^{+4.5}_{-2.4}\) | \(4.5^{+0.6}_{-0.8}\) |
| \(\kappa_\mathrm{ph, +x}\) / \(\mathrm{cm^{2}\,g^{-1}}\) | \(0.03^{+0.05}_{-0.02}\) | \(9.9^{+0.2}_{-0.2}\) | \(5.5^{+2.0}_{-2.1}\) | \(3.9^{+0.8}_{-0.6}\) | \(2.8^{+0.0}_{+0.0}\) | \(0.01^{+0.01}_{-0.00}\) | \(4.7^{+3.8}_{-1.0}\) | \(4.9^{+0.7}_{-1.3}\) |
| \(T_\mathrm{ph, +z}\) / \(\mathrm{K}\) | \(6659^{+406}_{-437}\) | \(911^{+400}_{-136}\) | \(468^{+325}_{-124}\) | \(426^{+45}_{-67}\) | \(8675^{+169}_{-0}\) | \(5528^{+1738}_{-1166}\) | \(719^{+291}_{-272}\) | \(351^{+26}_{-27}\) |
| \(T_\mathrm{ph, +y}\) / \(\mathrm{K}\) | \(6858^{+251}_{-0}\) | \(619^{+392}_{-43}\) | \(509^{+81}_{-63}\) | \(326^{+48}_{-50}\) | \(8969^{+19}_{-2}\) | \(959^{+1229}_{-130}\) | \(517^{+288}_{-132}\) | \(366^{+44}_{-55}\) |
| \(T_\mathrm{ph, +x}\) / \(\mathrm{K}\) | \(6148^{+367}_{-345}\) | \(846^{+176}_{-17}\) | \(465^{+139}_{-183}\) | \(319^{+60}_{-47}\) | \(8660^{+5}_{-0}\) | \(5600^{+286}_{-224}\) | \(392^{+323}_{-62}\) | \(401^{+55}_{-107}\) |
| \(R_{1, \mathrm{ph, +z}}\) / \(\mathrm{AU}\) | \(2.2 \pm 0.3\) | \(28 \pm 6\) | \(104 \pm 9\) | \(525 \pm 12\) | \(2.1 \pm 0.2\) | \(5 \pm 1\) | \(136 \pm 15\) | \(530 \pm 22\) |
| \(R_{1, \mathrm{ph, +y}}\) / \(\mathrm{AU}\) | \(2.1 \pm 0.2\) | \(31 \pm 7\) | \(142 \pm 11\) | \(451 \pm 14\) | \(2.1 \pm 0.2\) | \(24 \pm 10\) | \(151 \pm 16\) | \(610 \pm 19\) |
| \(R_{1, \mathrm{ph, +x}}\) / \(\mathrm{AU}\) | \(2.2 \pm 0.3\) | \(36 \pm 5\) | \(137 \pm 10\) | \(461 \pm 15\) | \(2.3 \pm 0.2\) | \(10 \pm 1\) | \(142 \pm 16\) | \(585 \pm 18\) |
10 Convergence and resolution↩︎
10.1 Convergence with respect to ray-grid resolution↩︎
Here we give an account of the effect of changing the resolution in the grid used to construct the lightcurve. It is clear from Figure 15 that grid resolution does not greatly affect the lightcurve and that a resolution of \(256\times 256\) cells is sufficient. Also, the size of the particles’ smoothing lengths at the photosphere is two to ten times larger than the pixels’ size, denoting that the grid resolution is sufficient.
10.2 SPH resolution study↩︎
In Section 4.2 we scrutinised the resolution of the simulation and its ability to accurately reproduce the light by following a different argument. There we calculate the percentage of the light emitted by pixel (ray) \(i\) by the most contributing particle. If that fraction is close to 100% it means only one particle contributes to the light of that pixel, denoting very low resolution. In Figure 16 we show polar (\(z\)) and side (\(y\)) maps of the quantity \(f\) on the same physical scale as the brightness maps in Figure 4. This gives an immediate idea of where the resolution is poor and where it is better.
Figure 16: The fractional contribution to the specific intensity of a pixel \(i\), \(I_i\), for the \(j^{\rm th}\) SPH particle, that contributes the most (see Equation 14 ), for the 1.7 M\(_{\odot}\) (left two columns) and the 3.7 M\(_{\odot}\) (right two columns) simulations. The redder colours on the top row show that one particle dominates the calculation of the light at that location, while the dark blue colours later on show that many more particles contribute to the calculation of the light at that location, with a resulting much smaller uncertainty. Similarly to Figure 4, data is viewed from the polar (\(z \rightarrow +\infty\)) direction (left panels) and the orbital (\(y \rightarrow +\infty\)) direction (right panels), at \(t=0, 3, 12, 20\) yr for the top, second, third and bottom rows respectively. The origin is fixed at the primary’s core particle’s location. Movies of this image see https://zenodo.org/records/18239846 for both 1.7 M\(_{\odot}\) (2md) and 3.7 M\(_{\odot}\) (4md) simulations, viewed from \(z \rightarrow +\infty\) (xyz), \(y \rightarrow +\infty\) (xzy) and \(x \rightarrow +\infty\) (zyx)..
11 1D vs 3D photosphere at time zero↩︎
Another source of systematic lightcurve uncertainty comes from a mismatch of the outer layers of the initial temperature profile of the AGB star, between the 1D mesa profile, and the profile mapped and relaxed in the 3D computational domain (Figure 17). The methods used to carry out the mapping and relaxation have been thoroughly investigated [17], [54], [55] with an eye to ensuring that the result of the common envelope simulations are not dependent on the imperfect reproduction of the 1D profile in the 3D domain. The lack of sensitivity of the hydrodynamic simulations to the specifics of the stellar surface layers is mostly due to the fact that the mass that expands past the 1D photosphere is small.
In our 3D AGB stars, the few low density and relatively low temperature surface particles that protrude outside the 1D photosphere are optically thick. As we can see from Figure 17 this results in a photosphere that is both hotter and larger than the original 1D photosphere, with a resulting larger luminosity as seen in Figure 5, where the luminosities of the two models at time zero are \(5\times10^5\) and \(10^6\) L\(_{\odot}\) for the 1.7 M\(_{\odot}\) and 3.7 M\(_{\odot}\) models, respectively, compared to the value they should have as given by the mesa calculation: \(5\,180\) L\(_{\odot}\) for the 1.7 M\(_{\odot}\) simulation and \(11\,900\) L\(_{\odot}\) for the 3.7 M\(_{\odot}\) one.
Here we evaluate the effect of these particles on the lightcurve after time zero. We re-run the first 5 years of the same 1.7 M\(_{\odot}\) and 3.7 M\(_{\odot}\) simulations, but the surface particles have been forced to have the same temperature, \(T\), as prescribed by mesa. By assigning the mesa temperatures to the 3D star after relaxation (which has a different density distribution compared to the original mesa profile), leads to a different entropy distribution compared to the 1D MESA model. We do this here to achieve the correct thermal properties for this purpose.
Since phantom finds the SPH particles’ temperatures by interpolating the mesa EoS from the particles’ densities, \(\rho\), and their specific internal energy, \(u\) (see also [15]), we enforce the required \(T\) for each particle by solving for the corresponding \(u\) keeping \(\rho\) fixed in the tabulated mesa EoS via root-finding, and then setting the particles to have that value of \(u\).
Since the particles’ temperatures only deviate from the mesa prescription outside \(\sim100\) R\(_{\odot}\) (see Figure 17), only particles at radii larger than this value are treated in this way. We also removed all particles outside the mesa photosphere surface (261.8 R\(_{\odot}\) for the 1.7 M\(_{\odot}\) simulation and 330.9 R\(_{\odot}\) for the 3.7 M\(_{\odot}\) simulation). The resulting \(T\) profile at \(t=0\) is illustrated in Figure 17.
Figure 17: Stellar temperature profile of the models, before (blue symbols) and after (red symbols/curve) forcing phantom initial dump to have mesa-prescribed temperatures (see text). Left panel is for the 1.7 M\(_{\odot}\) simulation, while the right panel is for 3.7 M\(_{\odot}\)..
This new, adjusted star evolved in isolation to test its stability, shows a similar stability to the original model, as demonstrated by inspecting the “particle" and”density" radii as plotted by [17], as well as their figure A1, second row, fourth column. This reassures us that the star’s hydrodynamic behaviour would be comparable to that of the main model used here.
The resulting lightcurve, calculated for only 5 years with the new star with an adjusted surface, is compared with the original case in Figure 18. The luminosity at time zero is a factor of two lower, but still much higher than the 1D, correct luminosity, due to the surface resolution issues detailed in Section [sec:uncertainties:surface-particles].
Besides the difference in luminosity at time zero, the lightcurves (red curves) follow the same general trend as the original light curves, despite some differences in the duration of the peak (especially in the 1.7 M\(_{\odot}\) case), where it can be broader or narrower by several months. At the qualitative level, the mismatch of initial surface temperature profile is not problematic, but at the quantitative level, systematic uncertainties are very likely.
Figure 18: First 5 years of the lightcurves as viewed from three directions: +x, +y and +z, top, middle and bottom row, respectively, for the 1.7 M\(_{\odot}\) (left colummn) and the 3.7 M\(_{\odot}\) (right column) simulations. The blue/orange/green curves are for the original simulations as seen in Figure 5, as reference; red curves are for simulations where the surface particles have their temperature reset to match the mesa profile values, as described in this Appendix (see text). The shaded areas denote the error range as calculated in Section [sec:uncertainties:surface-particles] and Appendix 12..
12 Upper and lower lightcurve limits↩︎
Here we provide details to quantify the uncertainty deriving from the lack of resolution near the surface, as a direct result of too few SPH particles sampling the steep increase of density, temperature and opacity as we enter the stellar surface. Effectively, surface particles are statistically distributed around the mesa-prescribed stellar surface. However, these particles effectively protrude \(2h\) outwards and our ray integration starts as soon as the ray enters the particle volume thus defined. Simply put, the integration of density, opacity, and source functions starts too soon, likely increasing the luminosity5. This is evidenced in Figure 5, where our raw bolometric luminosity estimates are two orders of magnitude larger than the value suggested by mesa.
We next calculate what must be a lower limit to the luminosity effectively eliminating the source function coming from the outer half of surface SPH particles. In order to determine what classes as “surface" SPH particles, we divide the region travelled by each ray into two parts: the poorly resolved”surface" region, which the ray meets first as it travels from the observer into the gas distribution, and the well-resolved region, farther inside the object.
What distinguishes the well-resolved from the poorly-resolved regions (as described in Section 4.2) is that in the well resolved region the ray has already crossed approximately a quarter of the average number of neighbours per particle, \(N_\mathrm{neigh}/4\) - by "cross" we intend that the ray has passed through that number of SPH particles, each with a radius of \(2h\). We choose \(N_\mathrm{neigh}/4\), because it is expected that, in the well-resolved region, any point on the ray should have an average of \(N_\mathrm{neigh}\) particles neighbouring it. This implies that well-resolved point on ray should expect \(N_\mathrm{neigh}/2\) neighbour particles on the closer-to-observer side, while the other half is on the farther-from-observer side. At the boundary between an the well-resolved and the poorly-resolved regions, this is true as well. However, we acknowledge that one does not need the full \(N_\mathrm{neigh}\) of neighbouring particles to have a reasonable idea of the various quantities (such as \(\rho\), \(\kappa\) and \(S\)) at that point through SPH interpolation. Since an ideally less-resolved region would have zero neighbouring particles on the closer-to-observer side (instead of \(N_\mathrm{neigh}/2\) neighbouring particles), we somewhat arbitrarily pick the midway point—\(N_\mathrm{neigh}/4\) neighbouring particles on the closer-to-observer side—as our limit that separates the well-resolved region to the less-resolved region.
There is only one small addendum to this accounting. Each particle has a different impact on the values of quantities interpolated on the ray, that depends on the particle’s (dimensionless) distance \(q_{xy}\) to the ray location. It is therefore important to weigh each particle’s contribution by its distance to the ray, especially within the poorly-resolved region. As such, instead of directly counting the neighbouring particles encountered by the ray, we weigh the counting of neighbouring particles, \(j\), by their contribution to the smoothing kernel integration to the ray \(i\): \[K_i(z) \equiv \sum_j \int_{z}^{+\infty} w(q_{ij}(z')) \frac{dz'}{h_j}. \label{eq:calc:K95z}\tag{19}\] In other words, \(K_i(z)\) is a quantification of the dimensionless column kernel value - where the column kernel is itself a quantification of the particle’s influence on the ray to point \(z\). This variable is essentially the sum of the dimensionless column kernel (\(w_\mathrm{col}(q_{xy, ij}) \equiv \int_{-\infty}^{+\infty} w(q_{ij}(z')) dq_z'\)) associated with each particle \(j\) up until the point \(z\). (The column kernel is proportional to each particles’ influence on the ray: note how Equation 19 is similar to Equation 10 , just without the factor of \(m\kappa / h^2\)).
A point on ray \(i\) in an ideally resolved region should have an average of \(N_\mathrm{neigh}\) neighbouring particles. Each particle \(j\), with its distance to the ray \(i\) being \(q_{xy, ij}\), contributes to the ray a dimensionless column kernel of
\[w_\mathrm{col,ij} = \int_{-R_\mathrm{kern}}^{+R_\mathrm{kern}} w\left(\sqrt{q^2_{xy,ij} + q^2_{z,j}}\right) dq_{z,j},\] where \(R_\mathrm{kern} = 2\) is the kernel radius (Equation 6 and \(w(q)=0\) when \(q>R_\mathrm{kern}\)). Since these particles are expected to be randomly distributed, the probability of a neighbouring particle residing in the shell of \(q_{xy} \rightarrow q_{xy} + dq_{xy}\) from the ray is \[\frac{2 \pi q_{xy} dq_{xy}}{\int_0^{R_\mathrm{kern}} 2 \pi q'_{xy} dq'_{xy}} = \frac{2 \pi q_{xy} dq_{xy}}{\pi R^2_\mathrm{kern}} .\] Together, the above allows us to calculate a critical value of the dimensionless column kernel for location \(z\), \(K_{\rm crit}\) on ray \(i\), that separates the poorly resolved outer layers from the well resolved inner layers \[\begin{align} K_{\rm crit} = N_\mathrm{neigh} \int_{0}^{R_\mathrm{kern}} dq_{xy} \, \left[ \frac{2 \pi q_{xy}}{\pi R^2_\mathrm{kern}} \int_{-R_\mathrm{kern}}^{+R_\mathrm{kern}} dq_z \, w\left(\sqrt{q^2_{xy} + q^2_z}\right) \right] \\ = 4.608, \end{align}\] Following the same logic, we suggest that requiring \(K_i(z) > 4.608/4\) is the equivalent of requiring \(N_\mathrm{neigh}/4 = 14.5\) neighbouring particles before declaring entry of the well-resolved region, i.e., \[K_i(z) > 1.152\] becomes the condition for the region inside of location \(z\) on the ray, to be considered well resolved.
Comparing with the approach of neglecting the first 14.5 particles for each ray, the above method accounts instead for the lower influence of particles farther away from the ray.
References↩︎
E-mail: chunliang.mu@students.mq.edu.au↩︎
We have checked that using a full Planck spectrum and adding up intensities for each pixel at each frequency, rather than adding up pre-integrated blackbodies on a single grey ray, gives the same results. The reason is that we are using gray opacities, so all the terms in \(A_{\rm eff}\) are wavelength independent, thus, simply changing the source function to be wavelength-dependent in Equation 8 (and then integrating the result across wavelength) does not change the outcome.↩︎
Note that this method is distinct from the one adopted by [15], where the photosphere size was simply the distance between the core of the star to the photosphere integrated along one ray only in the polar or equatorial direction.↩︎
The \(g\) band magnitude maximum quoted elsewhere corresponds to 720 000 L\(_{\odot}\), but the bolometric luminosity maximum took place 67 days later.↩︎
Technically, since the integration of the opacity also starts too early, the extra source function could be self absorbed completely and even obscure the emissivity from particles behind it (closer to the centre of the star), making the luminosity actually lower.↩︎