Gravitational wave emission from dynamical stellar interactions

We are witnessing the dawn of gravitational wave (GW) astronomy. With currently available detectors, observations are restricted to GW frequencies in the range between ${\sim} 10\,\mathrm{Hz}$ and $10\,\mathrm{kHz}$, which covers the signals from mergers of compact objects. The launch of the space observatory LISA will open up a new frequency band for the detection of stellar interactions at lower frequencies. In this work, we predict the shape and strength of the GW signals associated with common-envelope interaction and merger events in binary stars, and we discuss their detectability. Previous studies estimated these characteristics based on semi-analytical models. In contrast, we used detailed three-dimensional magnetohydrodynamic simulations to compute the GW signals. We show that for the studied models, the dynamical phase of common-envelope events and mergers between main-sequence stars lies outside of the detectability band of the LISA mission. We find, however, that the final stages of common-envelope interactions leading to mergers of the stellar cores fall into the frequency band in which the sensitivity of LISA peaks, making them promising candidates for detection. These detections can constrain the enigmatic common-envelope dynamics. Furthermore, future decihertz observatories such as DECIGO or BBO would also be able to observe this final stage and the post-merger signal, through which we might be able to detect the formation of Thorne-\.Zytkow objects.


Introduction
Mass transfer between stars in binary systems can lead to many different outcomes.Some of these outcomes can dramatically change the systems (Langer 2012;De Marco & Izzard 2017).Whenever the stars that form a binary get close enough to each other so that one of them fills its Roche lobe (because of either a change in the evolution stage of the stars or in the orbital parameters), mass transfer will start.Unstable mass transfer can become dynamical and lead to a common-envelope (CE) event or even a full stellar merger.In a CE event, the accretor is engulfed in the envelope of the giant donor (Paczynski 1976) and the stellar cores transfer orbital energy to the envelope.This shrinks the orbital separation.The energy injected in the envelope can lead to its ejection (successful CE event).As a result, the orbit will cease shrinking, and the outcome of the interaction is a binary system in a tighter orbit (Ivanova et al. 2013).In contrast, if the companion fails to eject the envelope of the giant star, it will continue to inspiral and ultimately result in a merger of the core of the giant star and the companion (CE merger).
With the upcoming launch of the Laser Interferometer Space Antenna (LISA; Baker et al. 2019), we may have a new tool for the study of these dynamical phases of CE events and stel-lar mergers, as they are promising sources for gravitational wave (GW) signals (Ginat et al. 2020;Renzo et al. 2021).Compared to the detectors of the LIGO1 -Virgo-KAGRA2 network (Aasi et al. 2015;Acernese et al. 2014;Akutsu et al. 2021), which operates in a frequency band centered on hundreds of Hz, LISA will have a peak sensitivity centered at ∼4 mHz.The gap in between those frequency bands is proposed to be covered in the future with decihertz observatories such as the DECi-hertz Interferometer Gravitational wave Observatory (DECIGO; Seto et al. 2001) or the Big Bang Observer (BBO; Phinney 2003).The frequency range of LISA will allow it to detect signals from wider binary systems (the frequency of the GWs released by a binary system is twice the orbital frequency: f GW = 2 × f orb ).Similarly to the case of compact object mergers, the GWs produced by a stellar binary system shift to higher frequencies as the orbit shrinks.In the case of stellar mergers and CE events, the rate at which the orbital separation decreases is generally not governed by the emission of GW radiation, but is mainly driven by mass transfer and drag forces.Although not dynamically

Model
M 1 (total) M 1 (core) MS-merger 9 -8 6.2 5.3 × 10 −5 5.0 7.3 × 10 −5 CE-merger 0.63 -0.42 0.047 2.0 × 10 −2 0.039 2.5 × 10 −2 Successful CE 2 0.37 0.99 49 6.8 × 10 −7 3.9 3.0 × 10 −5 Notes. a 0 is the initial orbital separation, a final is defined as the separation between the stars right before disruption in the case of the MS star merger and the CE merger.In the case of the successful CE, as the orbit is eccentric, a final is defined as the semi-major axis of the orbital separation between the core and the companion at the end of the simulation.f GW0 and f GWfinal are computed as twice the orbital frequency assuming a Keplerian orbit and separation equal to a 0 and a final , respectively.In the case of the successful CE ejection, only the mass of the core and the companion are considered for the calculation of the Keplerian frequency.
dominant, GWs provide an interesting opportunity for observing such short-lived and elusive events.Three-dimensional hydrodynamic (3D HD) or magnetohydrodynamic (3D MHD) simulations are one of the most reliable (albeit expensive) ways of modeling CE events (see Ricker & Taam 2012;Ohlmann et al. 2016;Staff et al. 2016;Prust & Chang 2019;Reichardt et al. 2019;Sand et al. 2020;Chamandy et al. 2020;Reichardt et al. 2020;Moreno et al. 2022;Glanz & Perets 2021;Lau et al. 2022a,b;Zou et al. 2022;Ondratschek et al. 2022, for some recent examples) and mergers between different types of stars, such as main-sequence (MS) stars (Schneider et al. 2019) and white dwarfs (Munson et al. 2021;Dan et al. 2012;Zhu et al. 2015).In setups in which general relativity (GR) is dynamically not relevant, these methods provide the most consistent predictions of the merger process, and GWs can be derived from these simulations, but this has rarely been done.In the case of white dwarf mergers, predictions of the GWs that are released were obtained using smoothed particle hydrodynamics (Rasio & Shapiro 1994;Lorén-Aguilar et al. 2005), as they are a known target for LISA (Baker et al. 2019).Most of the research done has focused on the GWs that are emitted during the inspiral phase rather than during the merger, however.There are predictions of the GWs emitted during CE events (Ginat et al. 2020;Renzo et al. 2021), but because it is difficult to fully model them, all of the research so far has been based on analytic or semi-analytical models.The reliability of these predictions is limited by the fact that the interaction between the inspiraling companion and the envelope is not modeled in a self-consistent 3D MHD approach.Consequently, they ignore the contribution of the envelope itself to the production of GWs and employ an approximate drag force that drives the inspiral.
For neutron-star mergers, however, 3D HD simulations accounting for GR effects are common (Baiotti & Rezzolla 2017;Bauswein & Stergioulas 2019;Dietrich et al. 2021;Shibata 2015).This has broadly been applied for the study of GW signals that are produced by neutron-star mergers as they can be detected by the LIGO-Virgo-KAGRA network.
In this work, we fill the gap and derive GW signals from up to date 3D MHD simulations of binary interactions for which GR is not dynamically relevant.We revisit an already published simulation of a main-sequence star merger (Schneider et al. 2019) and run two new simulations of a successful CE ejection and a CE merger.
In Sect. 2 we briefly describe the methods we used for the different simulations and how we computed the GW signals.In Sect. 3 we explain the outcome of the simulations and discuss the GW signals.We discuss the detectability of these events in Sect.4, and we summarize our main findings in Sect. 5.

Methods
All simulations were carried out with the code arepo (Springel 2010;Pakmor et al. 2011), which is a 3D MHD code that uses an unstructured moving Voronoi mesh with a second-order finitevolume approach.Gravity is Newtonian and was computed using a tree-based algorithm.The energy loss by GWs is negligible and was therefore not taken into account.For the successful CE ejection simulation, we included the impact of recombination energy by using the OPAL (Rogers et al. 1996(Rogers et al. , 2002) ) equation of state (EoS).We also used the OPAL EoS for the main-sequence star merger, whereas for the CE merger, in which we consider only the stellar cores, for which ionization effects are irrelevant, we used the Helmholtz EoS (Timmes et al. 2000).The initial models of the main-sequence stars as well as the giant star were created using the one-dimensional stellar evolution code MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015)).The setups of the different simulations are summarised in Table 1.

Main-sequence star merger
For the merger between main-sequence stars, we used the simulation described in Schneider et al. (2019), involving two early main-sequence stars with masses of 8 M and 9 M with radii of 1.7 R and 1.9 R , respectively.It is computationally not feasible to simulate the merger from the point of Roche-lobe overflow until the merger, so the merging process was sped up artificially for 1.5 orbits, and only then was the evolution of the merger followed by the simulation.

Successful CE ejection
In giant stars, the dynamical timescale of the outer parts of the envelope is orders of magnitude longer than in the core.Accurately resolving the core requires that the simulation uses time steps on the order of milliseconds, whereas dynamical processes involving the envelope will occur on a timescale of days.For this reason, the computational cost of simulating the common-envelope phase with a fully resolved core and companion becomes too expensive, and therefore, approximations are necessary.
For the successful CE ejection, we used one of the models described in Ohlmann (2016) as primary star.We took a model with a zero-age main-sequence (ZAMS) mass of 2 M , evolved with MESA until the red giant (RG) stage, at which point it had a radius of 49 R .This model was mapped onto an arepo grid.For this simulation, we used a version of arepo with the A9, page 2 of 11 same methods and the OPAL EoS as in Sand et al. (2020).We replaced the core of the giant star and the companion with a gravitation-only particle to limit the range in timescales, resolving their softening lengths with a minimum of 20 cells following Ohlmann (2016).The softening length of both particles was initially set to 2.47 R , and this value was reduced during the simulation to ensure that it remained below 40% of the separation between the core and the companion.We did not include magnetic fields in this simulation.The masses of the core and companion were 0.37 M and 0.99 M , respectively.This stellar model was then relaxed to re-establish hydrostatic equilibrium (Ohlmann et al. 2017).For the relaxation process, we placed the star in a box with a size of 7×10 13 cm and a uniform background density of ρ = 10 −13 g cm −3 and evolved it for ten dynamical timescales (8.2 × 10 6 s).During the first half of the relaxation, we applied a damping force to remove spurious velocities; during the second half, the star was let free to evolve on its own, and we verified that the model remained hydrostatically stable.We then simulated the common-envelope phase by placing the star and its companion in a box with a side length of 1.6 × 10 15 cm (2.35 × 10 4 R ) at an initial separation of 49 R , which corresponds to the surface of the RG, with a corotation fraction of 95%.

CE merger
As we cannot simulate the entire common-envelope phase with a fully resolved core and companion, we set up an entirely new simulation to study the GW signals released by a CE merger in which we considered only the part that generates the loudest highest-frequency GWs during the event: the merger of the red giant core and its companion.To ensure that the companion was not disrupted too far from the RG core, we chose a compact object; in particular, we took a 0.6 M white dwarf (WD) as a companion for which 50% of the mass was carbon and 50% oxygen, and whose initial radius was 0.012 R .As first-order approximation, we modeled the RG core as a 0.4 M WD made entirely of helium with a radius of 0.016 R .The methods and setup of the white dwarfs in this simulation are similar to those of Pakmor et al. (2021).The stellar models were constructed as simple hydrostatic equilibrium configurations assuming a uniform temperature of 2 × 10 7 K for the RG core and 10 6 K for the companion.We obtained their density profiles following the Helmholtz equation of state, and then we gave the arepo cells the appropriate masses.Both stars were given a seed magnetic field set up in a dipole configuration with a polar surface field strength of 10 G.The stars were relaxed separately to ensure hydrostatic equilibrium because mapping the initial 1D models to 3D can lead to discretization errors in the hydrostatic equilibrium.The relaxation process was the same for both stars: They were placed in a box of 10 10 cm with a static, uniform background grid density of ρ = 10 −5 g cm −3 .During the first half of the relaxation, we applied a damping force that reduced potential spurious velocities from discretization errors.During the second half of the relaxation, the stars were left to evolve on their own (without the damping force), and we verified that the density profiles were stable.The duration of the relaxation was chosen to be ten times the dynamical timescale of each star, that is, 20 s for the 0.6 M WD and 63 s for the 0.4 M RG core.After the relaxation, we placed the two relaxed stars into the same box with a box size of 7.3 × 10 11 cm with a static, uniform background grid with ρ = 10 −5 g cm −3 .This background density was chosen to resemble the denser regions of a stellar envelope, but it has very little effect on the dynamics of the merger because at this stage, it will be driven by the mass transfer and tidal interaction between the core and the companion.The mass resolution of the simulation was m cell ≈ 5 × 10 −7 M .The initial distance between the WDs was chosen to be three times the separation at the point of Roche-lobe overflow, and then we artificially shrank the orbit until the onset of a steady mass transfer in a similar fashion as for the main-sequence star merger.This occurred at an orbital separation of 3.25 × 10 9 cm.At this point, we started the unmodified simulation of the merger.

Computing the gravitational waves
To compute the gravitational wave signal, we calculated the approximate quadrupole radiation from Newtonian gravity.We followed the same approach as in Seitenzahl et al. (2015) and used the method by Blanchet et al. (1990) and Nakamura & Oohara (1989) to compute the second derivative of the quadrupole moment.In this way, the gravitational quadrupole radiation field in the transverse-traceless gauge h TT takes the form Here, P is the transverse-traceless projection operator, with the normalized position vector n = x/R and the distance to the source R = |x|.As usual, G is the gravitational constant, c is the speed of light, ρ is the density, u is the velocity, ∂ i represents a partial derivative with respect to the spatial coordinate i, and Φ is the Newtonian gravitational potential.
The amplitude of the GW radiation field can be written in terms of the two unit linear polarization tensors e + and e × , and when a line of sight is chosen, the amplitude of the GWs takes the form where A + and A × are the amplitudes of the two polarizations as a function of time.We worked in Cartesian coordinates and considered as line of sight the z-direction, corresponding to the initial axis of rotation of our binary systems.In this line of sight, the GW emission is strongest and therefore provides the most optimistic scenario in terms of detectability of the signal.Along this axis, the polarization amplitudes take the form with A i j defined as We implemented the computation of the terms A xx , A xy , A xz , A yy , A yz , and A zz in arepo following Eq.( 5) in order to output the emitted GW signal on the fly during the simulations.
The characteristic strain of a GW signal is defined as h c ( f ) = 2 f h( f ), where h is the Fourier transform of the time-domain amplitude of the GW, with where A(t) is the amplitude of each polarization over time.We calculated these values using a discrete fast Fourier transform and a Tukey window (Abbott et al. 2016).To capture the main features of the signal at its start and end, we used an asymmetric Tukey window in the case of the successful CE ejection.The signal-to-noise ratio (SNR) of a signal then takes the form (e.g., Moore et al. 2015) where S n is the noise power spectrum density of the detector.We took the approximate values for LISA from Robson et al. (2019).

Main-sequence star merger
A full description of the hydrodynamic outcome of the merger between the main-sequence stars can be found in Schneider et al.  2019) is shown by the dashed orange line.The region shaded in gray represents the frequency range spanning between twice the orbital frequency at the beginning of the simulation and twice the orbital frequency at the moment the star is disrupted (assuming the orbits to be perfectly circular and Keplerian).
(2019).The initially more massive (primary) star transfers mass to the initially less massive (secondary) star for a few days, draining angular momentum and shrinking the orbit.By the time the orbital separation reaches the tidal disruption radius of the primary star, commonly approximated as where M 1 and M 2 are the masses of the primary and the secondary stars and R 1 is the radius of the disrupted star (primary), it is disrupted, and its central material merges with that of the secondary star.The resulting remnant consists of a ∼14 M central core with a 3 M torus surrounding it.The gravitational waves determined from the 3D MHD simulation of the merger event are shown in Fig. 1.The amplitude of the strain increases as the orbital separation shrinks, and it peaks at the moment of the disruption of the primary star.Immediately after the merger takes place, the emission ceases.This decline in signal strength is expected as the remnant rapidly evolves to a mostly axisymmetric configuration.The characteristic strain of the GWs is shown in Fig. 2 assuming a distance of 1 kpc between source and detector.For this distance, we obtain an SNR of 0.062 with LISA.The peak of the characteristic strain is found at a frequency of about 5.4 × 10 −5 Hz, which corresponds to an orbital frequency of 2.7 × 10 −5 Hz.
The Keplerian orbital separation for a system orbiting at this frequency is 4.3 × 10 11 cm.This value is very similar to the separation between the stars at the beginning of our computations.The characteristic strain integrates the signal over the observation time, so that the longer the system spends at a certain orbital frequency, the higher the associated characteristic strain for that frequency.In our simulation, the binary system only evolves through a small range of orbital separations before the disruption between 4.3 × 10 11 cm and 3.5 × 10 11 cm, which corresponds to a narrow GW frequency range from 5.4 × 10 −5 Hz to 7.3 × 10 −5 Hz (shaded in gray in Fig. 2).Theoretically, we A9, page 4 of 11 expect most of the GWs produced in our simulation to be within this frequency range.We lack simulation data at larger separations, which explains the cutoff at frequencies below this range.We could also have lower frequencies from ejected material after the disruption, so the characteristic strain is not zero.The signal at higher frequencies (smaller separations) originates from the phase after disruption of the primary star, and hence, we see a drop-off in the characteristic strain at frequencies higher than the frequenciy corresponding to the separation equal to the tidal disruption radius of the star.The characteristic strain at these higher frequencies ( f GW > 10 −4 Hz) comes from the weak post-merger signal.
It is difficult to compute the tidal disruption radii in noncompact interacting systems.The usual approximation for the tidal disruption radius, shown in Eq. ( 9), is obtained by equating the self-gravity force of the disrupted star with the gravitational pull that it experiences from the companion object and assuming that the radius of the star is much smaller than the separation between the stars.Therefore, this approximation does not hold true when the tidal disruption radius is on the order of magnitude of the radius of the star.
To obtain a better approximation, we calculated the tidal force across the disrupted star.First, we computed the gravitational acceleration at the points of the star that are closest and farthest from the companion: where d is the distance to the companion.To obtain the tidal acceleration, we subtracted the contribution from the acceleration on the near point from the contribution from the point far out, The acceleration due to self-gravity is At the point of tidal disruption, d = r t and a tidal = a self .Therefore, equating Eqs. ( 11) and ( 12), we obtain Solving Eq. ( 13) numerically for the tidal disruption radius, we obtain a better approximation of its value.In the following, we define the tidal disruption radius as the solution of Eq. ( 13) unless stated otherwise.
Another source of uncertainty is the fact that the stellar radius is difficult to define during ongoing mass transfer.For example, using the original radius of the star from the MESA models, we obtain a tidal disruption radius of 5.5 × 10 11 cm, which is larger than the initial separation of our system.When we consider as the stellar radius the distance from the center in which 99% of the stellar mass is contained, we obtain a more sensible value for the tidal disruption radius of 3.5 × 10 11 cm. Figure 1 shows that the latter value matches the actual point of disruption well.With this, we aim to show that the tidal disruption radius of a star in a binary system can give us a rough estimate of the highest frequency of the GWs that the merger would release during the inspiral phase.

Successful common-envelope ejection
The initial separation between the core of a giant star and its companion is generally too large for the detection of the system with LISA.The frequencies of the associated GWs are too low because the radii of giant stars are on the order of 10 to 100 R .In order for these events to be detectable, the orbit must shrink to a separation at which the orbital frequency is in the appropriate frequency band for LISA.The evolution of our simulation of a successful CE ejection is very similar to that of Sand et al. (2020), but with an RG instead of an asymptotic giant branch (AGB) star, and the whole interaction occurs on a shorter timescale.
As soon as the simulation starts, the companion plunges into the envelope and its orbit shrinks for a few hundred days, until it ejects most of the envelope mass.From this point on, there is little orbital evolution on the simulated timescale, and the binary remains on an eccentric orbit with a semi-major axis of 3.6 R after 200 days.The evolution of the orbital separation is shown in Fig. 3.At this final separation, the companion orbits the core of the RG with a frequency of 0.017 mHz.In contrast to the other two simulations, the GW signal is not shut down.Instead, the main component of the signal evolves from low to higher frequencies.The characteristic strain of this simulation is shown in Fig. 4. It peaks at a frequency of 3.4 × 10 −5 Hz, which corresponds to exactly twice the orbital frequency after the envelope has been ejected.The final orbit is not perfectly circular; it has an eccentricity of 0.1.Keplerian elliptic orbits emit GWs at multiple harmonics (frequencies that are integer multiples) of the orbital frequency (Maggiore 2008).We observe this in the characteristic strain of the signal as further peaks located at frequencies equal to the harmonics of the final orbital frequency.The frequency of maximum strain is also still too far from the peak sensitivity of LISA to be detectable; we obtain an SNR of 4.28 × 10 −4 .
A common approximation used for calculating the GWs generated by CEs is to only consider the core of the RG and the companion as point-mass particles and neglect the effect of the envelope on the GW signal (Ginat et al. 2020;Renzo et al. 2021).With our simulation, we can test how accurate this treatment is. Figure 5 shows the difference in the characteristic strain A9, page 5 of 11 of the GWs released by the entire system and the strain resulting from the cores alone.The characteristic strain differs slightly at the lowest frequencies, but is essentially the same at frequencies higher than 3 × 10 −6 Hz.As most of the envelope mass is ejected early on in the simulation when the orbital frequency of the system increases rapidly, the envelope is not expected to contribute significantly at GW frequencies much higher than twice the initial orbital frequency.As the system moves to higher orbital frequencies by ejecting the envelope, the contribution to the GW signal from the envelope decreases at higher frequencies.The difference between the two characteristic strains should therefore mostly be present at GW frequencies close to the one corresponding to the initial orbital frequency.In any case, the magnitude of this discrepancy is very small, therefore we can safely state that the approximation holds true for these systems in regard to their detectability with LISA.It is not clear whether this also applies to a CE evolution that leads to a merger.If the envelope is not completely ejected, it may contribute to the GW signal at higher frequencies.Furthermore, while there is no significant direct contribution of the envelope to the formation of the GW signal, an accurate representation of the core-envelope interaction in the simulations is still necessary to reliably determine the rate of inspiral and the final orbital separation between the cores that ultimately shape the GW signal.

Common-envelope merger
The hydrodynamical evolution of the CE merger is summarized in Fig. 6.From the beginning of the simulation, the less massive, less dense core of the RG primary star transfers mass onto the more massive and more compact companion.The orbit shrinks for a few hundred seconds due to this mass transfer until it reaches the point of tidal disruption and the RG core is disrupted around the inspiraling WD companion.The companion accretes a fraction of the mass of the former RG core, approxi- mately 0.11 M .The remaining mass forms a disk-like structure around the remnant.This structure is not axisymmetric immediately after disruption, but it becomes so after a few hundred seconds.
The strain of the GWs (Fig. 6, bottom) is quite similar to the emission of the merger between MS stars, with a few differences.For the WDs, the merger occurs on a shorter timescale (hundreds of seconds instead of days).The amplitude of the GWs does not noticeably increase immediately before the merger, but its maximum value rather remains mostly constant until the disruption.Moreover, the post-merger signal is initially more pronounced than in the case of the MS star merger.The explanation for this stronger post-merger signal is the initial lack of axial symmetry in the remnant.As the disk homogenizes, the signal declines.
There is again a narrow range of orbital separations before the binary is disrupted.Our initial separation is 3.25 × 10 9 cm, which corresponds to a GW frequency of 2 × 10 −2 Hz.No significant contribution to the signal is expected at frequencies lower than this value in our simulation.The computation of the tidal disruption radius for this simulation is even more challenging than in the case of the MS star merger because the disrupted star is a WD.This means that its radius will increase as it loses mass.Eq. ( 13) shows that the tidal disruption radius increases with the stellar radius.Therefore, the tidal disruption radius is underestimated when the initial value of the radius of the RG core is used.When the distance from the center in which 99% of the stellar mass is contained is used as the radius, we find r t = 2.4 × 10 9 cm following Eq.( 13), which corresponds to a GW frequency of 3.17 × 10 −2 Hz, assuming the orbit to be perfectly circular and Keplerian.The top panel of Fig. 6 shows that this value matches the point of disruption well.
Computing the characteristic strain (Fig. 7), we find that it peaks at a frequency of 2 × 10 −2 Hz, which matches twice the orbital frequency at the beginning of the simulation.This is expected because this is the largest orbital separation that we A9, page 6 of 11 simulated, therefore we do not predict any significant GW signal at lower frequency.In this frequency range, LISA is far more sensitive than in the case of the MS star merger.For this reason, the characteristic strain of the CE merger ejection peaks above the sensitivity curve of LISA for a signal that lasts only a few minutes, generating an SNR of 6.6 at a distance of 1 kpc.The post-merger signal also contributes some high-frequency components to the characteristic strain.There is some relatively high strain up to frequencies close to 0.1 Hz.These higher frequencies are in the correct range for future decihertz GW observatories such as DECIGO or BBO.We computed the SNR for the CE merger simulation with DECIGO and BBO using the sensitivity curves described in Yagi & Seto (2011), and we obtained values of 896 and 1865, respectively, at a distance of 1 kpc.If designed according to current plans, these instruments would be able to detect the final disruption of CE-mergers with a much higher SNR due to the increase in sensitivity and the higher frequency band.
The computed value of the tidal disruption radius has been underestimated and thus overestimates the highest frequencies released by the binary before the merger (end of the region shaded in gray in Fig. 7).It is still within a factor of two from the peak frequency, however, and can be used as an upper bound.

Toy model for CE inspiral
Based on our successful CE ejection simulation, we modeled the plunge-in of the companion inside the red giant envelope with simple analytical formulae in order to determine how the characteristic strain obtained from them compares with the strain of our simulation.We used three functional forms to model the orbital decay: an exponential, a Gaussian, and a linear function, defined as Here a 0 and a final have the values specified in Table 1, and τ in each formula was chosen to be the timescale of the orbital decay such that it fits the rate observed in our successful CE simulation.We find τ = 1.47 × 10 6 s for the exponential and Gaussian models, and τ = 2.94 × 10 6 s for the linear model.We followed these models for a time equal to the longest time in our simulation.Figure 3 shows that the Gaussian model captures the early part of the plunge-in best, whereas the exponential model matches it better at a later stage.In these simple models, we assumed no eccentricity, such that there is a smooth decrease of the orbital separation.We can then compute the GWs that would be released by two point-masses following the orbital evolution as described in Creighton & Anderson (2011), A where µ is the reduced mass, ν(t) is defined as GM total a(t) and φ(t) is the phase of the GWs, which we obtained by numerically solving φ = 2π f GW (t) = 4π f orb (t).For these toy models, we assumed that the point-masses are always on a Keplerian circular orbit in order to obtain the instantaneous value of the orbital frequency.We calculated the characteristic strain of these toy models and compare it with the data from our simulations in Fig. 5.As we explicitly specified a final orbit, we have a sharp peak at twice the final orbital frequency of the system, and no features at its harmonics as there is no eccentricity.In general, the shape of the strain of the toy model GWs reasonably matches the strain of the arepo simulation.For the lowest frequencies produced by our simulation, the three toy models give the same value for the GW strain as in the hydro simulation when only the core masses were taken into account.For frequencies greater than 3×10 −6 Hz, the exponential model reproduces the shape of the strain in the simulation best.The three toy models produce an SNR similar to the SNR in the simulation (4.00 × 10 −4 for the exponential, 4.16×10 −4 for the Gaussian model, and 4.17×10 −4 for the linear model).We conclude from this that the models capture the strain of the GWs in the specified frequency range moderately well and can be used as simple toy models for the inspiral of our CE simulation.
Assuming that this approach also holds valid in a CE merger scenario, we applied it to the phase leading up to our CE merger simulation in order to estimate the GW signal produced during its inspiral phase.Disregarding any potential signal from the envelope, we took two point mass particles with the masses of our WD and the core of the RG to compute the reduced mass in Eqs.(15a) and (15b) and shrank the orbit according to Eq. ( 14) at a rate that resembles what we observe in our accurate simulation of the successful CE ejection (τ = 1.47 × 10 6 s in the exponential and Gaussian models, and τ = 2.94 × 10 6 s for the linear model).However, instead of stopping the inspiral at a distance of 3.6 R , The sensitivities of LISA (Robson et al. 2019), DECIGO, and BBO (Yagi & Seto 2011) are shown as dashed orange, green and red lines, respectively.The region shaded in gray represents the frequency range spanning between twice the orbital frequency at the beginning of the simulation and twice the orbital frequency at the moment when the star is disrupted (assuming the orbits to be perfectly circular and Keplerian).The black lines show the characteristic strain produced by our toy models.
we let the orbit shrink to a separation of a final = 3.25 × 10 9 cm, equal to the initial orbit used in our simulation of the CE merger.With these toy models, we obtained an estimate of the characteristic strain of the GWs released during the earlier stages of a CE merger ejection.The orbital evolution of our toy models is shown in Fig. 8.As observed in our previous simulations, the orbital frequency shifts from lower values up to the initial values of the CE merger simulation, while the strength of the GWs increases.The duration of this inspiral phase differs between our toy models, as we stopped them the moment they reached a separation equal to the initial orbital separation in our core merger simulation.
A9, page 8 of 11 The more time the system spends orbiting at a certain orbital frequency, the greater the strain of the GWs associated with that frequency.The exponential model has the highest value for the characteristic strain at high frequencies because it spends more time at orbital separations close to the desired final orbit, and the linear model has the lowest SNR because it reaches its final orbit faster.The signal during this stage crosses the peak sensitivity of LISA, and therefore, most of the SNR for this observatory could actually come from this phase depending on the exact rate of the orbital decay.With the chosen toy models, we indeed obtain SNR values for the inspiral phase of 191.1 for the exponential model, 84.5 for the Gaussian, and 10.5 for the linear model (to compute the SNR in these toy models, we ignored the characteristic strain at frequencies higher than twice the orbital frequency at the final orbital separation because this part of the characteristic strain is just numerical noise introduced by the fast Fourier transform on nonperiodic data).The final part of the inspiral phase begins to enter the frequency range of DECIGO and BBO, and because of the planned increase in sensitivity with respect to LISA, they would be even better able to measure this final stage of the CE merger, obtaining DECIGO SNRs of 6784, 2961, and 299 and BBO obtaining SNRs of 13957, 6092, and 616 for the exponential, Gaussian, and linear models, respectively.

Discussion
We computed the SNRs for all our simulations following Eq.( 8), and assuming a polar distance to the source of 1 kpc.This distance is small and provides a best-case scenario for the detection of these events.The SNR is inversely proportional to the distance between detector and source, that is, if the event were to take place at 10 kpc rather than 1 kpc, the SNR would be ten times lower.
In the case of the merger between two MS stars, we obtain an SNR of 0.062, which suggests that these events are not suitable candidates for LISA sources.In Fig. 2 we show that the characteristic frequency of the GWs of this type of mergers is too far from the peak sensitivity of LISA.We started the simulation only a few days prior to the merger, whereas LISA could theoretically observe such a system for a few years before the disruption takes place, building up more signal and producing a higher SNR, but the frequency of these GWs will still be in the low-sensitivity band.The strength of the produced signal would increase in the case of a more massive merger (GW amplitude ∝ M 5/4 chirp ; Creighton & Anderson 2011), where M chirp is the chirp mass of the binary, However, it would also likely mean a larger tidal disruption radius, leading to lower GW frequencies and a decreased LISA sensitivity.The maximum frequency of GWs released by these systems before disruption can be estimated as twice the orbital frequency at the point where the separation between the stars equals the tidal disruption radius.We conclude that even at higher masses or longer observing times, the mismatch between the frequency of the GWs and the peak sensitivity of LISA makes main-sequence star mergers an unlikely GW source for this observatory.
Figure 4 shows that in our simulation of a successful CE ejection, the final orbit reached is not tight enough to produce GWs in a frequency that LISA can detect even at an optimistic distance of 1 kpc.As indicated by Renzo et al. (2021), the final orbit resulting from a CE ejection has to shrink below one solar radius to be in a frequency range where LISA has a chance to detect them.Unfortunately, we do not observe a sufficient orbital shrinkage in this simulation or in any of the other simulations resulting in successful CE ejection that we analyzed (Sand et al. 2020;Ondratschek et al. 2022;Moreno et al. 2022).Therefore, we do not obtain any detectable GW signal from the dynamical phase of successful CE ejections for the LISA mission.Our final orbital separations are in line with those from Iaconi & de Marco (2019), showing that most simulations of successful CE ejections produce final orbits larger than 1 R .
Different setups might shrink the orbit further.All our setups involved giant stars with loosely bound envelopes that require less energy to be transferred from the binary orbit to eject them.It is possible that CE events in which the giant star has a more strongly bound envelope might result in tighter final orbits.There are some indications that the final orbits reached in our simulations are too large compared to observations (Iaconi & de Marco 2019).Many of the close binaries used as LISA verification targets (Stroeer & Vecchio 2006;Kupfer et al. 2018) will have evolved through a successful CE phase.It is not clear whether their orbits shrink to the current tight orbits after the CE or in the dynamical phase of the CE.This further shrinking stage after the dynamical plunge-in could be due to a self-regulated inspiral, as suggested by Ivanova et al. (2013) or other mechanism draining angular momentum from the binary orbit (e.g., a circumbinary disk) such that the GWs are in the frequency band of LISA.
Figure 7 showed that the frequencies of the GWs in our CEmerger simulation are in a band in which LISA is more sensitive than in the case of the merger between two main-sequence stars.Computing the SNR, we obtain a value of 6.6, even though the total simulated time in this case is only a few thousand seconds.The frequency of the post-merger signal of this simulation is close to the optimal frequency range of decihertz GW observatories such as DECIGO and BBO, and the signal from our merger would produce SNRs of 896 and 1865, much higher than LISA due to the much higher planned sensitivity of the detectors, and because their frequency bands are able to best detect the post-merger signal.In Fig. 9 we showed that a large fraction of the GW signal is likely to be produced during the inspiral phase: As the orbital separation decreases, the frequency of the GWs will move from low frequencies (where LISA is not very sensitive) up to higher frequencies, crossing the peak sensitivity of LISA.The exact rate of the orbital decay strongly impacts the SNR, and at the same time, the dynamical plunge-in is one of the lesser known stages of the CE phase.A detection of such an event could provide information about stellar interiors in RG stars and the drag forces that the companion experiences while spiraling into the envelope.For example, different density profiles of the star lead to differences in the temporal evolution of the drag forces, causing very different waveforms and SNRs.Even if we were to observe the electromagnetic counterpart but do not detect the GW signal, a slow inspiral might be ruled out and a rapid plunge-in through the envelope might be indicated.Therefore, accurate models for this inspiral stage are required for generating waveform templates and estimating the detection rates.In turn, future observations may reveal binary dynamics.
Using semi-analytical models and smoothed particle hydrodynamics simulations, Ginat et al. (2020) computed the waveform and characteristic strain of the inspiral phase of CE mergers.Their simulations covered the initial phase of the CE stage, but their models were stopped at larger orbital separations than the initial separation in our simulation, and a core merger A9, page 9 of 11 was assumed to take place.The authors showed that the GWs released until that point might be observable by LISA.In our 3D models, we followed the final disruption in detail, probing down to shorter orbital separation (higher frequencies), and our findings confirm that statement.The disruption on its own already produces a strong signal that could be detected if the source were close enough to LISA.Although we did not compute a detailed model of the late inspiral phase, our estimates using toy models of the inspiral point in the same direction as those of Ginat et al. (2020).With SNRs ranging between approximately 10 and 191 (depending on the model), the signal from the entire CE event could be strong enough for a detection with LISA.We also computed the SNR for the signals of these toy models with DECIGO and BBO.If their final design achieves the planned sensitivity, they might be able to detect the final phase of the inspiral in CE mergers with a much higher SNR than LISA, ranging between approximately 300 and 7000 for DECIGO and between 600 and 14 000 for BBO, allowing them to be found at much larger distances.
Our findings about the detectability of these events with both LISA and the decihertz observatories are expected to hold for different masses for the core of the giant and also the inspiraling companion.Depending on the exact masses, we might enter a regime in which the CE merger of two white dwarfs leads to runaway thermonuclear burning and to an explosive event (Pakmor et al. 2021), which, in contrast to type Ia supernovae, is expected to contain substantial amounts of H and He in its ejecta (this corresponds to the so-called core degenerate supernova scenario of Soker 2011).However, we did not include nuclear burning here.In the event that one of the stars explodes, the GW signal will diverge from our results from that point on.The explosion could release a strong GW signal (Seitenzahl et al. 2015) and the emission would cease immediately thereafter.If the explosion occurred during or after the disruption of the core, mainly the high-frequency components of the GWs will be affected, which does not alter the SNR significantly, and makes the system an interesting candidate for a multimessenger detection.
Changing the mass of the companion changes the tidal disruption radius of the core, and thus the highest GW frequency reached by the system before the disruption.A higher-mass companion such as a neutron star will disrupt the core at a larger distance and will also produce a stronger signal because the GW amplitude is proportional to M 5/4 chirp , but the maximum frequency is more sensitive to the core radius and mass than to the companion mass.
If the companion were to be a neutron star, the formation of a Thorne-Żytkow object might be observed (Thorne et al. 1977;Hirai & Podsiadlowski 2022) because the strength and frequency of the GW signal during the inspiral and merger would not sensibly differ from the one we computed here.In particular, these objects are predicted to show high-frequency features in the post-merger signal, for whose detection decihertz observatories would be very well suited.

Conclusions
We have computed the GW signal released during the merger between two MS stars, a successful CE ejection, and a CE merger in order to study their detectability with LISA.The simulations were carried out with the 3D MHD code arepo.
The GWs from mergers between main-sequence stars are found to be outside the sensitivity range of LISA at distances equal to or larger than 1 kpc.This result holds true for every merger between main-sequence stars, regardless of their masses.The tidal disruption radii of stars in binary systems can be used to approximate the cutoff frequency of the GWs released by these systems during their mergers.
From our simulation of a successful CE ejection, we find that the final separation between the core of the giant star and the companion is too large to release GWs in the frequency range where LISA could detect them.This result applies to all successful CEs that we simulated, but observations seem to suggest smaller final orbits, which would move the final systems closer to the frequency range of LISA.We also find that during successful CE ejections, the GWs emitted by the envelope contribute to the whole signal mostly in the lowest frequency range.This plays a very small role for the detection of such a system, and thus the signal associated with the envelope can safely be ignored.
Computing the GWs originating from the last stage of a CEmerger in 3D HD simulations of the interaction of the core of a giant star with a WD companion, we find their frequency inside the frequency range of LISA.They might therefore be detectable depending on the distance.We also used toy models to estimate the contribution to the GW signal from the inspiral phase and showed that it further increases the chances of detecting these systems.Such a detection would help to better understand the CE stage.The time evolution of the GWs gives an accurate depiction of how the orbital frequency changes as the companion inspirals through the envelope.GWs could serve as probes to study the interior of giant stars, helping us model the drag forces experienced by the companion and the density profile of the star.The exact conditions required for a successful ejection of the envelope or those leading to a merger are not fully known; detections of these events would place constraints on these cases.These events would probably have a detectable electromagnetic counterpart that could make them a suitable target for multimessenger astronomy.

Fig. 1 .
Fig. 1.Evolution of the main-sequence star merger.Top (a): Strain of the GW signal h over time (+ polarization in the z-direction, h = A z + /D).The amplitude is shown for an observer situated at a distance D. Bottom (b): Orbital separation over time.The dashed black line indicates the last stable orbital separation before disruption at 3.5 × 10 11 cm.

Fig. 2 .
Fig. 2. Characteristic strain of the GW signal released by the merger of two main-sequence stars (in blue) assuming a distance of 1 kpc to the source along the polar axis.The LISA sensitivity from Robson et al. (2019) is shown by the dashed orange line.The region shaded in gray represents the frequency range spanning between twice the orbital frequency at the beginning of the simulation and twice the orbital frequency at the moment the star is disrupted (assuming the orbits to be perfectly circular and Keplerian).

Fig. 3 .
Fig. 3. Time evolution of the orbital separation during the successful CE ejection.In blue we show the data from the arepo simulation.The black lines show the fit to these data with an exponential, a Gaussian, and a linear orbital decrease.

Fig. 4 .
Fig. 4. Characteristic strain of the GW signal released by the successful common-envelope ejection (in blue) assuming a distance of 1 kpc to the source.The LISA sensitivity from Robson et al. (2019) is shown by the dashed orange line.The vertical dashed gray line represents twice the orbital frequency at the beginning of the CE (assuming a perfectly circular and Keplerian motion).The dashed black line represents twice the orbital frequency after the ejection of the envelope.The dotted black lines represent the first four harmonics of the final orbital frequency ( f n = n • f orb,final ).

Fig. 5 .
Fig. 5. Characteristic strain of the successful CE ejection.The characteristic strain of the GWs released by the entire simulation (the core of the RG, its envelope, and the companion) is shown in blue.In green we show the characteristic strain of the GW signal produced by only the core of the RG and the companion.The dashed gray line represents twice the orbital frequency at the beginning of the simulation.The black lines show the characteristic strain produced by toy models (see Sect. 3.4).The LISA sensitivity fromRobson et al. (2019) is shown by the dashed orange line.

Fig. 6 .
Fig. 6.Evolution of the CE merger.Top: Orbital separation over time.The dashed black line indicates the computed tidal disruption radius at 2.4 × 10 9 cm.Middle: Hydrodynamical evolution of the core merger resulting from a CE merger.The density is color-coded.The computation of the gravitational waves starts at the onset of steady mass transfer.The more massive companion shreds mass from the less massive, less dense core of the RG.The separation between the two objects slowly decreases until it reaches the point of tidal disruption at ∼840 s.The companion accretes a fraction of the mass from the former RG core (∼0.11M ), while the remaining mass forms a disk-like structure around it.Bottom: Strain of the GW signal over time released by the CE merger (cross-polarization in the z-direction, h = A z + /D).

Fig. 7 .
Fig. 7. Characteristic strain of the GW signal of the CE merger (in blue) assuming a polar distance of 1 kpc to the source.The sensitivities of LISA(Robson et al. 2019), DECIGO, and BBO(Yagi & Seto 2011) are shown as dashed orange, green, and red lines, respectively.The region shaded in gray represents the frequency range spanning between twice the orbital frequency at the beginning of the simulation and twice the orbital frequency at the moment when the star is disrupted.

Fig. 8 .Fig. 9 .
Fig. 8. Time evolution of the orbital separation of our toy models.The orbit shrinks from the initial orbital separation of the successful CE simulation to the initial orbital separation of the CE merger.

Table 1 .
Summary of the simulations.