© The Authors 2024

1. Introduction

Binary systems formed by a neutron star (NS) and a white dwarf (WD) are one of the most abundant double compact object systems in the Galaxy (Nelemans et al. 2001), second only to WD-WD systems. If the two compact objects are born with a short orbital separation, the release of energy in the form of gravitational waves (GWs) can make the binary system merge on timescales shorter than the age of the Universe. This process can be accelerated if the binary has undergone a common-envelope (CE) phase, where the drag in the envelope drives the orbital shrinking and may even lead to a prompt merger. When the orbit has decreased to the point where the WD fills its Roche lobe, a mass transfer episode commences. The stability of this mass transfer depends primarily on the mass of the WD. For sufficiently low-mass WDs, the mass transfer will be stable, and the system can be observed as an ultra-compact X-ray binary (Nelson et al. 1986). If the WD mass is above a critical value of M_WD = 0.2 M_⊙ (Bobrick et al. 2017), the mass transfer will be unstable and the system will merge on a dynamical timescale, leading to a variety of electromagnetic transients (Metzger 2012; Fernández et al. 2019; Zenati et al. 2019, 2020; Bobrick et al. 2022; Kaltenborn et al. 2023) and a burst of GWs, similar to the ones released by binary black hole (BH) mergers (Abbott et al. 2016) and binary NS mergers (Abbott et al. 2017) that have been detected by the LIGO¹-Virgo-KAGRA² (LVK) Collaboration (Aasi et al. 2015; Acernese et al. 2014; Akutsu et al. 2021), albeit in a different frequency range.

The GWs released by these systems prior to the merger make these binaries a likely source (van Zeist et al. 2023; Abdusalam et al. 2020) for the LISA³ mission (Robson et al. 2019). The GWs released by the merger itself will have higher frequency and will be better captured by proposed deci-hertz observatories such as space-based DECIGO⁴ (Seto et al. 2001) and BBO⁵ (Harry et al. 2006) or the Moon-based LGWA⁶ (Harms et al. 2021). The post-merger signal which could originate as accretion onto the NS excites oscillation modes of the NS (Schutz 2008), would be in the frequency range of the LVK network (Aasi et al. 2015; Acernese et al. 2014; Akutsu et al. 2021).

The electromagnetic signal released during these mergers has been proposed to explain several fast-evolving transients such as AT2018kzr (Gillanders et al. 2020) or AT2018cow (Metzger 2022), as well as long gamma-ray bursts (GRBs) such as GRB 211211A (Yang et al. 2022; Zhong et al. 2023) or GRB 060614 (King et al. 2007). NS-WD mergers are also predicted to result in magnetars, and were suggested to cause some fast radio bursts (FRBs) such as FRB 190824 and FRB 180916.J0158+65 (Zhong & Dai 2020).

In this work, we aim to study such a merger by the means of magnetohydrodynamic (MHD) simulations. Previous studies have looked into several different aspects of these mergers. Fernández et al. (2019) and Zenati et al. (2019) used two-dimensional (2D) hydrodynamic simulations to study the accretion disk resulting from the disruption of the WD and the nuclear reactions that would take place on it. Paschalidis et al. (2011) performed the first fully general relativistic simulations of mergers between NSs and scaled-down “pseudo-WDs” and studied whether the remnant would collapse into a BH. Zenati et al. (2020) and Bobrick et al. (2022) ran three-dimensional (3D) simulations of NS-WD mergers using smooth particle hydrodynamics (SPH) starting from realistic initial conditions. The enormous range of scales required for this simulations, spanning over ten orders of magnitude in density, forced previous works to revert to two-dimensional models or a relatively low resolution when simulating these events. Furthermore, none of the existing hydro simulations have yet included magnetic fields, which could be decisive in shaping the merger and post-merger dynamics (Zhong & Dai 2020). In this paper, we study the dynamics of the disruption and present the first high-resolution, 3D MHD simulation of a NS-WD merger.

2. Methods

For our simulations, we use the AREPO code (Springel 2010; Pakmor et al. 2011, 2016; Weinberger et al. 2020). AREPO uses the finite volume method to solve the 3D, fully compressible equations of ideal MHD with gravity:

$$\begin{array}{cc}& \frac{\partial \rho }{\partial t}+\mathrm{\nabla }·(\rho \mathbf{u})=0,\hfill \end{array}$$(1)

$$\begin{array}{cc}& \frac{\partial (\rho \mathbf{u})}{\partial t}+\mathrm{\nabla }·(\rho \mathbf{u}\otimes \mathbf{u}+p_{\mathrm{tot}}\mathbb{I}-\mathbf{B}\otimes \mathbf{B})=\rho \mathbf{g},\hfill \end{array}$$(2)

$$\begin{array}{cc}& \frac{\partial E}{\partial t}+\mathrm{\nabla }·[(E+p_{\mathrm{tot}})\mathbf{u}-\mathbf{B}(\mathbf{B}·\mathbf{u})]=\rho (\mathbf{g}·\mathbf{u}),\hfill \end{array}$$(3)

$$\begin{array}{cc}& \frac{\partial \mathbf{B}}{\partial t}+\mathrm{\nabla }·(\mathbf{u}\otimes \mathbf{B}-\mathbf{B}\otimes \mathbf{u})=\mathbf{0},\hfill \end{array}$$(4)

where ρ denotes the density, u the velocity, B the magnetic field⁷, p_tot = p + |B|²/2 the total pressure, p the gas pressure, E the total energy density per unit volume, and g the gravitational acceleration. Equations (1)–(4) are discretized on an unstructured, moving, Voronoi mesh as described in Pakmor et al. (2016). The boundary conditions are periodic for all variables except the gravitational potential, which we compute by assuming that the density outside of the box domain is 0. Hyperbolic fluxes are computed at cell interfaces using the five-wave HLLD solver of Miyoshi & Kusano (2005), and the strength of magnetic monopoles is kept under control using the Powell scheme (Powell et al. 1999). The grid is generated at each timestep from a set of mesh-generating points that move with the local velocity of the fluid with small corrections to keep the mesh regular, leading to a nearly Lagrangian scheme. The resolution of the simulation can be changed locally by refining or de-refining the mesh following a set of criteria. The main criterion enforces the same mass for all cells within a factor of two of a chosen target mass and a maximum volume ratio between neighboring cells of a factor of five. On top of these criteria, we employed additional refinement criteria which are described in detail in Sects. 2.1 and 2.2. Self-gravity is included with a tree-based algorithm for the gas cells, and via exact force computation for the NS. We performed these simulations using Newtonian gravity, as the processes and dynamics that we aim to resolve are sufficiently distant from the NS. The equation of state (EoS) used in our simulations is the Helmholtz EoS (Timmes et al. 2000). We assume that the simulation takes place in a regime where the gas is completely transparent to neutrinos, and they leave the system without interactions. Therefore, neutrino transport simplifies to a cooling process. We take into account thermal neutrino cooling by implementing⁸ the fits from Itoh et al. (1996). The energy lost to GW is negligible and not taken into account, as it is not dynamically relevant for the part of the merger that we simulate. We computed the GWs in the same way as Morán-Fraile et al. (2023) by calculating the approximate quadrupole radiation from Newtonian gravity. We did not use a nuclear network as nuclear processes are predicted to play a minor role in the dynamical evolution of the disruption (Zenati et al. 2019).

2.1. Stellar setup

We took a 1.4 M_⊙ NS as a primary star in our binary system, and for the companion we took a 1 M_⊙ WD composed of an equal-by-mass mixture of carbon and oxygen with a central density of ρ_c = 3.4  ×  10⁷g cm⁻³ and an initial radius of 5.3  ×  10⁸ cm (defined as the radius containing 99.9% of the stellar mass). The WD model was generated with a uniform initial temperature of 5  ×  10⁷ K. We then obtained the density profile by assuming hydrostatic equilibrium and mapped it into 3D using a healpix-based algorithm following Pakmor et al. (2012) and Ohlmann et al. (2017). We used the Helmholtz EoS to compute the values of the internal energy for the AREPO cells given the temperature. We added an initial dipole magnetic field with an polar surface value of 10⁵ G to the WD. The magnetic dipole field is set up along the z-axis, and we imposed a cutoff radius in the center of the star with a size equal to two times the radius of the smallest AREPO cell to avoid an arbitrarily large value for the magnetic field in the center of the WD.

This WD is relaxed in isolation to ensure hydrostatic equilibrium, as deviations from equilibrium could appear due to discretization errors when mapping the one-dimensional profiles to the unstructured 3D mesh of AREPO. For the relaxation we place the star in a box of 10¹⁰ cm side length with a uniform background grid of density ρ = 10⁻⁵ g cm⁻³. The total relaxation time is chosen to be ten times the dynamical timescale of the star, corresponding to 11 s. During the first half of this time we apply a damping force to reduce the spurious velocities that could have arisen during the mapping. During the second half of the relaxation, the star is left to evolve on its own (with no damping force being applied), and we check that the density, pressure, and internal energy profiles remain stable.

Due to the difference in the dynamical timescales of a NS and a WD, it is computationally unachievable to resolve the hydrodynamic evolution of both the NS and the WD in one simulation. We therefore model the NS as a gravitation-only point particle with a mass of 1.4 M_⊙. To limit the pressure gradient next to our point mass NS, we used a softened gravitational potential. We found that low-resolution inside the softening length produces issues with energy conservation from the beginning of the simulation. At the same time, increasing the resolution close to the NS becomes prohibitively expensive for small values of the softening length. Experimentally, we increased the resolution while decreasing the size of the softening length, and we found a good compromise between computational expense and energy conservation by using a softening length of 8.4  ×  10⁷ cm and enforce a minimum of 40 cells per dimension within the softening length.

The gravitational potential inside the volume enclosed by the softening length deviates from a Newtonian potential. Therefore, we consider this region unresolved, and we report the values obtained by the simulation within this volume solely as a reference as they may deviate from realistic ones. As a point mass, the only characteristic for a NS that we assume is its mass and we do not model any other physical property. Thus, we cannot give it a magnetic field with the current implementation of MHD in AREPO. We could set a dipole configuration in the background around it, but as the mass contained in the background is so small compared with the one of the WD, as soon as accretion onto the NS starts, the initial magnetic field around the NS will become dynamically irrelevant. Therefore, we did not establish an initial magnetic field around the NS.

2.2. Binary setup

The target mass resolution of our simulation is m_cell = 5  ×  10⁻⁷ M_⊙, which results in ∼2  ×  10⁶ AREPO cells before further refinement. Due to the extreme pressure and density gradients in the neighborhood of the NS, we further increased the resolution in a disk-like region around the NS where a large fraction of the gas ends up after the disruption of the WD. We define this region as a cylinder centered on the NS with a radius of 6  ×  10⁸ cm (approximately half of the tidal disruption radius of the WD) and a height of 5  ×  10⁷ cm. In this region, we imposed the same refinement and derefinement conditions as within the softening length of the NS, which in turn means ensuring a minimum of ∼6 million particles in this region. We allowed for further refinement if the mass of the cells inside this region exceeds our target m_cell. With this setup, the typical total number of cells in the simulation is ∼8.1 million.

At the beginning of the simulation, the stars were placed in a Keplerian circular orbit at a distance of 2.27  ×  10⁹cm, which corresponds to 1.5 times the separation at Roche-lobe overflow (RLOF, Eggleton 1983). The WD is initially set to solid-body rotation in co-rotation with the NS at the point of RLOF. The orbit is then artificially shrunk on a timescale longer than the dynamical timescale of the WD so that potential tidal effects can still take place. This is done using the scheme described by Pakmor et al. (2021) that removes angular momentum from the binary system. We stopped shrinking the orbit at the point of RLOF (corresponding to an orbital separation of 1.52  ×  10⁹cm), when a steady mass transfer sets in. At this point, the system has a total angular momentum of L = 8.3  ×  10⁵⁰g cm² s⁻¹. From here on, the system was left to evolve on its own and the WD disruption takes place on a few dynamical timescales. As the outer layers of the WD cannot be fully resolved, there is some numerical mass loss during the orbit-shrinking stage. Part of this mass will be accreted by the NS before the onset of the steady accretion stream. Furthermore, as explained in Bobrick et al. (2017), the tidal force exerted by the NS deforms the WD making it not completely spherical and its gravitational potential not perfectly Keplerian. This results in our initial orbit having a small eccentricity (∼0.001), which also contributes to a small amount of mass transfer before the end of the orbit shrinking process. The amount of mass accreted into the softening length at the time that we stop the orbit shrinking is 10⁻⁴ M_⊙. This initial mass transfer is dynamically irrelevant, but during this stage the amplification of the magnetic fields around the NS already begins. Table 1 shows the initial parameters and important timescales of the numerical setup for easy reference.

3. Results

The hydrodynamical evolution of the merger is summarized in Fig. 1⁹, and the temperature is shown in Fig. 2. After the end of the artificial orbital shrinkage (t = 0 s), a steady mass transfer developed. Moreover, during this stage, a small fraction of mass from the WD left the binary through the Lagrange point, L₂. The transfer of mass and angular momentum makes the orbit shrink further, and after four complete orbits (t = 80 s, Fig. 1b) the WD becomes closer to the NS than its tidal disruption radius and is torn apart at an orbital separation of 1.14  ×  10⁹ cm. By this point, 0.005 M_⊙ of gas has been unbound. This takes place mainly through the L₂ point, with a smaller fraction of the gas being also ejected along the polar axis of the NS. The disruption unbinds another 0.01 M_⊙ of material. The disrupted stellar material initially forms a spiral arm and then starts to be accreted onto the NS (Fig. 1d).

Fig. 1. Hydrodynamical evolution of merger. The panels show density slices in the orbital plane from the beginning of the simulation (panel a) until we stop it at t = 242 s (panel j). The region with a softened gravitational potential is shown in a blue dashed line.

The infalling material intersects itself while accreting, developing shear instabilities and shocks that heat the gas and transport angular momentum outwards. In this transient process, the orbit of the disrupted material circularises and a disk is formed due to the tidal interaction of the turbulent, infalling material with the NS (Figs. 1e,f). The temperature in the shocks reaches values up to T = 9 × 10⁸ K (Fig. 2), but the hottest regions of the shocks remain at rather low densities of < 10⁶ g cm⁻³. Therefore, only minor nuclear burning is expected to take place in these shocks. The regions close to the NS heat up rapidly, reaching temperatures of 10¹⁰ K inside the softening length. In this region, nuclear reactions are essentially unavoidable. By t ≃ 160 s, the disrupted material has formed a disk-like structure around the NS. As the NS does not have a surface and its gravitational potential is softened, the accreted mass accumulates inside the volume of the softening length. In our simulation, any “acretion” is not onto the NS, but into this volume. During the 150 s following the disruption of the WD (t < 230 s), 0.056 M_⊙ of material is accreted into the gravitationally softened volume at an average rate of 3.7  ×  10⁻⁴ M_⊙ s⁻¹. The density inside this region keeps increasing as material is accreted. At the time of the WD disruption, the maximum density is ρ = 2  ×  10⁷ g cm⁻³, and at the end of the simulation it has reached a value of ρ = 1.4 × 10⁹ g cm⁻³. As we did not model nor resolve the physical processes around and at the NS, we observe an artificial accumulation of material inside the softened volume. Part of this material should be accreted onto the NS, potentially through an accretion disk. During several stages of the simulation the combination of temperatures and densities inside the softened volume are suitable for a strong nuclear burning, which would also impact the density and temperature profiles in the neighborhood of the softening length.

Fig. 2. Slice through orbital plane showing the temperature at selected times.

The tidal interaction of the NS with the nonuniform accretion disk generates shocks, producing spiral density waves (Figs. 1i,j) that transfer angular momentum (Fig. 3a). The heat dissipation produced by the shocks speeds up the transfer of angular momentum, increasing the accretion rate (Fig. 3c). This culminates in a rapid accretion event at t = 242 s, where we observe a peak accretion rate of 3  ×  10⁻³ M_⊙ s⁻¹. At this point, 0.077 M_⊙ of gas has been accreted into the gravitationally softened volume. Inside the gravitationally softened volume, the gas rotates as a solid body and transitions to a quasi-Keplerian rotation in the accretion disk. The distribution of angular momentum and mass during the simulation can be seen in Figs. 3a,b. Most of the mass and angular momentum after the disruption of the WD is concentrated at a distance from the NS equal to the tidal disruption radius of the WD. It is initially slowly pushed out of this region by the turbulent motions and the shocks. When the spiral waves are developed they transport momentum outwards very efficiently, as can be seen by the rapid change in the angular momentum distribution between t = 180 s and t = 242 s.

Fig. 3. Evolution of the disrupted material around the NS. Panel a: Z component of the angular momentum with respect to the distance to the NS for selected times. Panel b: distribution of gas mass with respect to the distance to the NS for selected times. The distance at which the WD is disrupted is shown by the black dashed line. Panel c: rate in M⊙ s−1 at which mass is accreted onto the gravitationally softened volume.

At t = 242 s we stop the simulation as we cannot resolve the accretion and energy errors start accumulating in the simulation. Considering every AREPO cell with a total energy larger than zero as unbound material, by the time we end the simulation, 0.94 M_⊙ of gas remain gravitationally bounded when only taking into account potential and kinetic energy. This number decreases to 0.93 M_⊙ when also including magnetic and internal energy in the total energy.

3.1. Magnetic field amplification

The evolution of the magnetic energy in the simulation is shown in Fig. 4. As explained in Sect. 2.1, the WD is given an initial magnetic field seed with a value of 10⁵ G at the equator. As soon as mass transfer begins, shear within the accretion stream amplifies the magnetic field around the NS, while compression of the magnetic field lines amplifies it in the shocked regions. The amplification process is complex and involves a combination of several mechanisms including magneto-rotational instability (MRI; discussed in more detail in Appendix A); Kelvin-Helmholtz instabilities and dynamo processes, similar to what is observed in 3D simulations of common-envelope evolution (Ondratschek et al. 2022); main-sequence star mergers (Schneider et al. 2019); binary white dwarf mergers (Zhu et al. 2015); or binary neutron star mergers (Obergaulinger et al. 2010; Kiuchi et al. 2015; Aguilera-Miret et al. 2020). This amplification process begins as soon as we start the simulation, and the NS begins to accrete the pseudo-vacuum surrounding it as well as part of the mass lost by the WD during this stage. By the time the artificial orbital shrinking phase finishes, the magnetic energy (E_B) of the system has increased from an initial value of 2  ×  10⁴¹ erg up to a value of 7.0  ×  10⁴⁵ erg, with the fields reaching values of up to 5.6  ×  10¹¹ G inside the gravitationally softened volume around the NS.

Fig. 4. Evolution of magnetic field throughout the simulation. Top: slices through orbital plane, centered on the NS, showing the absolute magnitude of the magnetic field for selected snapshots in color. Bottom: evolution of total magnetic energy over time (blue, dashed line) as well as the energy of the different components of the magnetic field in spherical coordinates (orange, green, and red for the radial, azimuthal, and polar components, respectively).

This initial amplification of the magnetic energy due to the NS accreting the pseudovacuum is unavoidable. To ensure that the amplification is not boosted by our artificial orbit shrinkage, we have performed two more simulations. In both of them we have placed the stars at a distance of 1.1 times the distance of RLOF. In the first one, we let the stars evolve freely (keeping constant orbital separation), and in the second one we shrank the orbit further down to the distance of RLOF before letting it evolve freely. At this initial distance of 1.1 times the separation at RLOF, the accretion of mass lost from the WD by numerical effects is the main source of magnetic amplification rather than the accretion of the pseudo-vacuum. In both simulations, the total magnetic energy is amplified to values of ∼1  ×  10⁴⁶ erg on a timescale of 25 s, matching the magnetic energy generated in our fiducial simulation at this separation, but generating it on a shorter timescale. The simulation left to evolve freely at a separation equal to 1.1 of the separation at RLOF remains at a constant orbital separation, and after this initial fast amplification the growth of the magnetic energy slows down and plateaus at around 2  ×  10⁴⁶ erg. In contrast, in the simulation where we shrank the orbit down to the point of RLOF, a steady mass transfer set in and the magnetic fields continued to be amplified with the evolution of the total magnetic energy matching the one of our fiducial simulation. Hence, we conclude that this initial phase of magnetic amplification previous to the start of our simulation is self-consistent and independent of the artificial orbital shrinking.

The amplification of the magnetic field takes place dominantly in a toroidal direction around the NS, following the accretion stream, as evidenced by the ϕ component of the magnetic energy in Fig. 4. When a steady mass transfer sets in, the amplification of the magnetic field speeds up as the NS begins to accrete the outermost layers of the WD. At this initial stage, the e-folding growth time of the magnetic energy peaks at a value of ∼20 s. As the orbit shrinks and the mass transfer rate increases, the magnetic energy keeps steadily growing, with an abrupt increase shortly after the disruption of the WD due to the turbulent motions. At this point, the total magnetic energy has increased to 6.7  ×  10⁴⁷ erg, reaching maximum local values of up to 5  ×  10¹² G inside the gravitationally softened region around the NS. During the period of rapid accretion onto the NS driven by the spiral waves (t ∼ 240 s), there is a second phase of fast amplification of the magnetic fields with E_B raising up to 5  ×  10⁴⁸ erg and reaching maximum values of 2.3  ×  10¹³ G next to the NS.

3.2. Jet-launching

The magnetic fields generated throughout the entire event power two polar outflows launched in opposite directions (Fig. 5). After the onset of mass transfer, regions with a magnetic-to-fluid pressure ratio (p_b/p_f) larger than one develop perpendicular to the orbital plane (third column in Fig. 5), where the magnetic pressure exceeds the dynamic pressure. In these regions, the Maxwell stresses dominate the dynamics, driving the collimated ejecta (note the white-to-blue shades in the upper panels of Fig. 6). As a result, since the beginning of the simulation, a small fraction of the accreting mass is ejected along the polar axis of the NS following the direction of the steepest magnetic pressure gradient. This outflow becomes faster and more collimated as the energy associated with the radial component of the magnetic field (E_Br in Fig. 4) increases. At t  ∼  140 s (Figs. 5e–h) the outflows become distinguishable as highly magnetized, high-entropy gas ejected perpendicular to the orbital plane with velocities close to 1% of the speed of light, c. By t  ∼  170 s (Figs. 5i–l), a bipolar jet structure has developed with characteristic velocities of ∼0.05c, reaching maximum velocities of 0.08c. In spite of the large entropy of the jets, the (relativistic) specific enthalpy, h = 1 + ϵ/c² + p_f/(ρc²) (where ϵ is the internal energy per unit mass), is very close to unity, indicating the ejection of a cold and highly magnetized jet. The computed values of h provide a proxy for the asymptotic Lorentz factor of the fluid, Γ_∞  ∼  1.005 − 1.01 (Γ_∞  ∼  h, according to Bernoulli’s law along streamlines). Likely, this is a lower bound, since after the formation of a low-density polar funnel, matter further ejected through the jet will be accelerated to higher terminal speeds. The strength of the jets varies in time and evolves unevenly in the two directions. During the late period of accretion dominated by the spiral waves, the outflow above the orbital plane weakens, while the outflow below the orbital plane greatly increases in strength with typical velocities v ∼ 0.08c and a maximum velocity of v ∼ 0.1c at t = 242 s (Figs. 5q–t). The outflows are initially launched at a distance of ∼4.4  ×  10⁸ cm from the NS, which roughly corresponds to five times its softening length. In the ejecta, p_b/p_f  ∼  10, and locally p_b/p_f reaching values of ∼1000 close to the gravitationally softened volume around the NS. The energy and luminosity of the jets is shown in the upper panel of Fig. 7. There, we measure a non-monotonically growing jet energy along with a growing luminosity, with imprinted variability timescales down to ∼25 s (Fig. 7 bottom panel, showing the ratio jet-energy-to-jet-luminosity, which provides a idea of the timescale of changes in the jet luminosity). Likely, this timescale is an upper bound of the smallest variability timescales induced by the jet interacting with the surrounding, stratified medium, dictated by the time it takes Alfvén or sound waves to cross the jet transverse radius (∼1–2 s; see, e.g., Aloy et al. 2002; Matsumoto et al. 2017 for analogous jets propagating in stellar environments). Capturing the jet–environment interaction requires a very large numerical resolution on the jet, which we are lacking here, as outside of our disk-like region surrounding the NS the resolution criteria fall back to target a mass per cell of m_cell = 5  ×  10⁻⁷ M_⊙. Also, in our case, it depends on the operative definition of matter belonging to the jet. While the impact of small variations of the radial velocity threshold adopted to define the jet on its energy is negligible, the choice of the threshold may induce moderate changes for the luminosity variability timescale. Shorter timescales may be induced by the dynamics of the accretion/ejection of matter in the vicinity of the (here unresolved) NS, but likely at later times (in analogy to what is found for jets induced in collapsars Aloy et al. 2000; Matzner 2003; Morsony et al. 2010; Gottlieb et al. 2020). The jets eject material at an average rate of ∼3 × 10⁻⁴ M_⊙ s⁻¹ and have opening half-angles oscillating between ∼22° and 25° for the top jet and between 18° and 22° for the bottom jet. As the accretion rate increases, the region of gas with the highest p_b/p_f draws closer to the gravitationally softened volume around the NS, and by t  ∼  242 s the material is launched from the outer regions of the softened volume.

Fig. 5. Slices through midplane perpendicular to the orbital plane of the binary, centered on the NS before and after the formation of the polar outflows. The first column shows the radial velocity of the gas as a fraction of the speed of light, (vr/c). The second column shows the entropy of the gas. The third column shows the magnetic-to-fluid pressure ratio, pb/pf. The fourth column shows the absolute magnitude of the magnetic field in the gas.

Fig. 6. Slices through midplane perpendicular to the orbital plane of the binary, centered on the NS at two times during the jet evolution. The top row shows the rϕ component of the Maxwell stress tensor in the gas normalized by the centrifugal component of the Reynolds stress tensor. The bottom row shows the specific enthalpy of the gas.

Fig. 7. Time evolution of energy and luminosity of the jets. Top: energy in the jets is shown by the plain blue line and the jets luminosity is shown in a blue dashed line. Bottom: energy divided by luminosity over time.

As gravity inside the softening length around the NS diverges from a realistic gravitational potential, we must be careful when interpreting processes that have an origin close to the softened volume. Also, nuclear reactions can massively change the energy budget of the central region. As the jet is initially launched at a reasonably large distance from the NS, we consider its formation a realistic physical effect rather than an artifact of the numerical model. However, as the region of acceleration draws closer to the NS, we must consider the jet-launching at later stages unresolved. Moreover, we used a nonrelativistic version of AREPO, and therefore gas moving at velocities close to 0.1c is no longer accurately treated, contributing to our decision to stop the simulation during the fast accretion event that takes place at t  =  242 s.

3.3. Gravitational waves

Figure 8 shows the GWs released during the simulation as seen from the top of the orbital plane, which is the viewing angle that results in the loudest signal, providing a best-case scenario for their observation. The cross (h_×) and plus (h₊) polarizations of the GWs have been computed using the approximate quadrupole radiation in the same way as in Morán-Fraile et al. (2023), and the characteristic strain is defined as $h_c(f)=2f\stackrel{\sim }{h}(f),$ where $\stackrel{\sim }{h}$ is the Fourier transform of the time-domain amplitude of the GW:

$$\begin{array}{c}\hfill h(f)=\frac{1}{R}\sqrt{|h_{+}{(f)|}^2+{|h_{\times }(f)|}^2},\end{array}$$(5)

Fig. 8. Gravitational waves produced in simulation. Top: strain of GW signal h times the distance D to the source, observed from the z-direction (polar axis), over time. Bottom: characteristic strain of GW signal assuming a distance of 1Mpc to the source along the polar axis. The sensitivities of LISA (Robson et al. 2019); DECIGO and BBO (Yagi & Seto 2011); LIGO, VIRGO, and KAGRA (Abbott et al. 2020); and LGWA (Dupletsa et al. 2023) are shown as dashed orange, green, red, purple, brown, and yellow lines, respectively. The gray shaded area denotes the frequencies below two times our initial orbital frequency where signal coming from the inspiral phase of the binary is expected to be present.

with

$$\begin{array}{c}\hfill h(f)={\displaystyle {\int }_{-\infty }^{\infty }}{e}^{-2\pi ift}h(t)\mathrm{d}t.\end{array}$$(6)

During the accretion of the outer parts of the WD, the amplitude of the waves oscillates in the time domain (Fig. 8a) mainly due to the eccentricity of the binary orbit at this stage, but it does not significantly rise or fall prior to the merger, while the frequency remains approximately constant. The emission reduces abruptly shortly after the tidal disruption of the WD. Only higher frequency components due to the lack of perfect axial symmetry of the accretion disk remain in the signal (Fig. 8a).

The characteristic strain (Fig. 8b) shows a strong peak at f = 0.1 Hz, which roughly corresponds to twice the orbital frequency of the binary at t = 0 s (f_orb = 4.8  ×  10⁻² Hz). As expected from binary star dynamics, the frequency of the GWs released corresponds to twice the orbital frequency of the system. This can be used to explain some of the main features seen in the characteristic strain. Because we only simulate the merger from an initial orbital frequency of ∼0.05 Hz, we miss the GWs released during the previous stages of the evolution of the binary, which would emit GWs mainly at frequencies below the peak. NS-WD binaries are expected to release GWs below this frequency during their inspiral stage (van Zeist et al. 2023).

In our simulation, the orbital frequency of the binary does not change substantially before the tidal disruption of the WD, and the amplitude of the GWs decreases sharply after the disruption, explaining the narrow peak that we observe in the characteristic strain (Fig. 8b). The higher frequency components of the GW signal produced primarily after the disruption have a much lower amplitude than the GWs released prior to the merger. Between f = 0.2 Hz and f = 0.5 Hz, they show a mainly constant characteristic strain ∼50 times smaller than at the peak. For higher frequencies, the strain drops down another order of magnitude, and at frequencies higher than 5 Hz we become limited by the sampling rate of our simulation. Frequencies on the order of some hertz and higher may arise from material close to and inside the softened volume, and they should therefore be considered unresolved.

The characteristic frequencies of the signal released during and after the merger make this event an ideal target for deci-hertz observatories such as BBO and DECIGO, which, for a source at a distance of 1Mpc, would produce a signal-to-noise ratio (S/N) of approximately 115 and 50, respectively. The S/N is defined as

$$\begin{array}{c}\hfill S/N=\sqrt{4{\int }_0^{\infty }\frac{{|h(f)|}^2}{S_n(f)}\mathrm{d}f},\end{array}$$

where S_n is the noise power spectrum density of the corresponding detector given by Yagi & Seto (2011). For the moon-based LGWA, the signal is not loud enough at 1 Mpc, and the merger would have to take place significantly closer to be detected.

4. Discussion and summary

We performed the first 3D-MHD simulations of the merger between a 1 M_⊙ COWD and a 1.4 M_⊙ NS. We find that the disruption of the WD takes place a few dynamical timescales after reaching the separation of RLOF. The mass transfer from the WD to the NS and the ensuing disruption of the WD amplify any pre-existing magnetic fields by several orders of magnitude, producing magnetic fields up to 10¹³ G.

The strong magnetic fields generated during the event launch two mildly relativistic jets perpendicular to the orbital plane with velocities on the order of 0.05 − 0.1c. These values should be taken as order-of-magnitude estimates because several aspects of the final outcome of these mergers are expected to depend on the physics that are not resolved in our simulations.

As we model the NS as a point mass particle with a softened potential, not only gravity inside the softened volume deviates from a realistic gravitational potential, but we also do not provide any physical model for accretion onto the NS, which would potentially be limited by the Eddington accretion luminosity, dramatically changing the conditions in the neighborhood of the NS. As the conditions around the NS allow for nuclear burning, we have run a simulation in which we couple the MHD solver to a 13 isotope nuclear network following the same methods as in Gronow et al. (2021) and Pakmor et al. (2021). In this simulation, the production of ⁵⁶Ni is observed both outside and inside the softening length. While the first effect is resolved and physical, the burning inside the softening length is unresolved but may still be physical. While we must exercise caution when interpreting these results, which will be studied in more detail in follow up publications, they do provide support for the idea that mergers between NS and WD result in iron rich transients such as AT2018kzr (Gillanders et al. 2020).

As reported in Sect. 5 the jet is initially launched sufficiently far from the NS, so we consider its origin physical, but as accretion into the softened potential increases, the jet becomes faster and its acceleration region draws closer to the NS. A shorter softening length will lead to larger flow accelerations close to the NS, potentially generating stronger magnetic fields, which may alter the properties of the jet. Moreover, because the velocities of the ejected material reach up to 0.1c, our Newtonian treatment becomes less accurate and a relativistic solver is needed at some point (e.g., Lioutas et al. 2022).

The merger of our two compact stellar remnants releases GWs in the deci-Hertz frequency range. We observe that GWs at the moment of WD disruption are released mainly at a frequency of 0.1 Hz, and the post merger releases higher frequency GWs, making these mergers an ideal target for proposed deci-Hertz observatories such as DECIGO, BBO, and LGWA. We only report GWs up to frequencies of ∼5 Hz as higher frequency components will mainly originate from the regions of the simulation with characteristic frequencies smaller than that which correspond to the softened volume. To obtain accurate GW signals of higher frequencies, accretion on to the NS must be modeled. With a hydrodynamically resolved NS, one could excite oscillation modes that would result in GWs with frequencies up to kilohertzs, which would be in the range of the LVK network (Schutz 2008).

The long-term evolution of the merger remnant is difficult to predict. At the point we ended the simulation, there was 0.94 M_⊙ of gravitationally bound gas. If this material were to be accreted onto the NS, it could lead to an accretion-induced collapse into a BH as we approach the maximum mass for an NS depending on the EoS. The jets have steadily ejected material at a rate of 3  ×  10⁻⁴ M_⊙ s⁻¹. We cannot predict the long-term evolution of this rate, as the launching region enters the unresolved volume at the end of the simulation. Additionally, the amount of mass ejected by the jet is likely dependant on unresolved or non-modeled physics. In particular, accurately modeling the accretion onto the NS can dramatically change the outcome. The accretion would release a great amount of the gravitational energy budget that is contained on the disk that and could contribute to the heating or ejection of the material as well as the release of neutrinos.

Nuclear reactions also inject energy into the gas, reducing the amount of gravitationally bound material. Depending on the amount of mass that these non-modeled effects can unbind, the NS may survive or collapse into a BH.

Despite the limitations of our current simulations, we show that a bipolar jet is self-consistently launched. The duration of the jetted ejection is difficult to estimate from the period of time computed so far. A rough extrapolation of the mass accretion rate onto the unresolved central compact object implies that the bound mass around the unresolved NS may be swallowed on scales of few thousand seconds. We find that the variability of the jet-luminosity happens down to timescales of ∼25 s. This is probably an upper bound induced by the moderate numerical resolution of the jets combined with the effective criteria employed to define the matter in the jet. The power at (significantly) smaller variability timescales would be induced by the changes close to the neutron star, which we cannot resolve.

The jet may be sustained for a sizeable fraction of the accretion time. The transient produced by this merger may share a number of qualitative similarities with other jetted ejections resulting from, for example, the tidal disruption of main-sequence stars by black holes (e.g., Mimica et al. 2015), with some key differences. Firstly, due to the compactness of the WD (compared to a 1 M_⊙ star), the disruption takes place very close to the NS and, as a result, a large fraction of the disrupted material is still gravitationally bound. In analogy to jetted tidal disruption events (TDEs), the accretion of most of the WD mass onto the central object on scales of thousands of seconds will likely produce a very luminous transient. Secondly, the jets launched during this event may asymptotically reach velocities of the order of the escape velocity at the jet basis.

Assuming the NS does not collapse into a BH despite the accretion, we expect asymptotic jet velocities close to those observed in objects such as galactic microquasars (≳0.9c, e.g., Mirabel & Rodríguez 1998), with typical Lorentz factors of two to ten. As a result, the spectrum of this event could be characterized by an initial soft X-ray or UV flash that will evolve into progressively lower frequencies. This would make this event resemble some fast blue optical transients (FBOTs; see, e.g., Drout et al. 2014; Ho et al. 2023; Prentice et al. 2018; Coppejans et al. 2020). These transients have also been observed in the outskirts of galaxies, suggesting an origin related to old stellar populations, for which a NS-WD merger would be a good match. Depending on the properties of the interstellar medium surrounding the NS-WD merger, a bright, weeks-to-months-long radio signal may be observed (see Mimica et al. 2015). The luminosity of this transient could be larger than that of typical kilonovae (∼ 10⁴¹ erg s⁻¹; Metzger 2017), not only because of the larger mass and velocity of the ejecta, but also because the simulated event does not appear to harbor optimal conditions for the production of r-process yields. The presence of lanthanides in kilonovae obscures these events due to their large contribution to the opacity. If, however, the NS eventually does collapse into a BH, significantly higher (ultrarelativistic) jet speeds could be produced, resulting in some long or possibly ultra-long GRB-like transient (Thöne et al. 2011; Levan et al. 2014).

We obtain an average neutrino flux of 2  ×  10⁴⁸ erg s⁻¹, which would put this transient one order of magnitude below the fluxes predicted for supernovae. However, we caution that in our simulations the neutrino luminosity is dominated by the unresolved conditions next to the NS.

The magnetically driven jet observed in our simulations is promising when it comes to explaining some recently discovered astrophysical transients. Therefore, a better numerical resolution of the material in the jet is desirable, and this can easily be achieved by changing the refinement criteria applied in our AREPO simulations. Future modeling work on mergers between WDs and NSs includes the addition of nuclear reactions to our simulations (for models including nuclear reactions, see Bobrick et al. 2022), paying special attention to the treatment of nuclear burning inside the unresolved regions, and a more accurate modeling of the physics around the NS, ideally including accretion onto the compact object.

Acknowledgments

The authors thank Valeriya Korol, Eva Laplace and Maryam Modjaz for their valuable input. J.M-F. is fellow of the International Max Planck Research School for Astronomy and Cosmic Physics at Heidelberg (IMPRS-HD) and acknowledges financial support from IMPRS-HD. This work was supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 759253 and 945806, the Klaus Tschira Foundation, and the High Performance and Cloud Computing Group at the Zentrum für Datenverarbeitung of the University of Tübingen, the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grant no INST 37/935-1 FUGG. M.A.A. acknowledges the support through the grant PID2021-127495NB-I00 funded by MCIN/AEI/10.13039/501100011033 and by the European Union, the Astrophysics and High Energy Physics program of the Generalitat Valenciana ASFAE/2022/026 funded by MCIN and the European Union NextGenerationEU (PRTR-C17.I1), and the grant CIPROM/2022/13 from the Generalitat Valenciana. Georgios Lioutas acknowledges support by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 759253 and by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (ERC Advanced Grant KILONOVA No. 885281). Georgios Lioutas also acknowledges support by the Deutsche Forschungsgemeinschaft (DFG, German ResearchFoundation) – Project-ID 138713538 – SFB 881 (“The Milky Way System”, subproject A10) and by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – MA 4248/3-1.

References

Appendix A: MRI

In order to analyze how well we resolved the MRI amplification in our simulation, we computed the fastest growing mode following Rembiasz et al. (2016) and compare it to the size of the grid cells in Fig. A.1. We find that in the regions closest to the NS, where most of the magnetic field amplification takes place (see Fig. 4), the fastest growing mode of the MRI is well resolved.

Fig. A.1. Number of cells per wavelength λMRI of the fastest growing MRI mode. Top: Slices through the orbital plane, centered on the NS. Bottom: Slices perpendicular to the orbital plane, centered on the NS.

All Tables

All Figures

	Fig. 1. Hydrodynamical evolution of merger. The panels show density slices in the orbital plane from the beginning of the simulation (panel a) until we stop it at t = 242 s (panel j). The region with a softened gravitational potential is shown in a blue dashed line.
In the text

	Fig. 2. Slice through orbital plane showing the temperature at selected times.
In the text

Fig. 3. Evolution of the disrupted material around the NS. Panel a: Z component of the angular momentum with respect to the distance to the NS for selected times. Panel b: distribution of gas mass with respect to the distance to the NS for selected times. The distance at which the WD is disrupted is shown by the black dashed line. Panel c: rate in M⊙ s−1 at which mass is accreted onto the gravitationally softened volume.

In the text

Fig. 4. Evolution of magnetic field throughout the simulation. Top: slices through orbital plane, centered on the NS, showing the absolute magnitude of the magnetic field for selected snapshots in color. Bottom: evolution of total magnetic energy over time (blue, dashed line) as well as the energy of the different components of the magnetic field in spherical coordinates (orange, green, and red for the radial, azimuthal, and polar components, respectively).

In the text

Fig. 5. Slices through midplane perpendicular to the orbital plane of the binary, centered on the NS before and after the formation of the polar outflows. The first column shows the radial velocity of the gas as a fraction of the speed of light, (vr/c). The second column shows the entropy of the gas. The third column shows the magnetic-to-fluid pressure ratio, pb/pf. The fourth column shows the absolute magnitude of the magnetic field in the gas.

In the text

	Fig. 6. Slices through midplane perpendicular to the orbital plane of the binary, centered on the NS at two times during the jet evolution. The top row shows the rϕ component of the Maxwell stress tensor in the gas normalized by the centrifugal component of the Reynolds stress tensor. The bottom row shows the specific enthalpy of the gas.
In the text

	Fig. 7. Time evolution of energy and luminosity of the jets. Top: energy in the jets is shown by the plain blue line and the jets luminosity is shown in a blue dashed line. Bottom: energy divided by luminosity over time.
In the text

Fig. 8. Gravitational waves produced in simulation. Top: strain of GW signal h times the distance D to the source, observed from the z-direction (polar axis), over time. Bottom: characteristic strain of GW signal assuming a distance of 1Mpc to the source along the polar axis. The sensitivities of LISA (Robson et al. 2019); DECIGO and BBO (Yagi & Seto 2011); LIGO, VIRGO, and KAGRA (Abbott et al. 2020); and LGWA (Dupletsa et al. 2023) are shown as dashed orange, green, red, purple, brown, and yellow lines, respectively. The gray shaded area denotes the frequencies below two times our initial orbital frequency where signal coming from the inspiral phase of the binary is expected to be present.

In the text

	Fig. A.1. Number of cells per wavelength λMRI of the fastest growing MRI mode. Top: Slices through the orbital plane, centered on the NS. Bottom: Slices perpendicular to the orbital plane, centered on the NS.
In the text