The quest for Magrathea planets. II. Orbital stability of exoplanets formed around double white dwarfs

TL;DR

This study investigates the orbital stability of second-generation giant planets formed around double white dwarfs (DWDs) in circumbinary, P-type configurations, comparing non-resonant and resonant architectures over a few million years using a hybrid symplectic N-body framework tailored to circumbinary dynamics. By sampling 2–5 planet systems around three DWD binaries and tracking catastrophic events with metrics like Normalized Angular Momentum Deficit (NAMD), orbital spacing, and center-of-mass shifts, the authors map how multiplicity and resonance influence long-term stability. They find that two-planet systems are most robust, while three to five-planet architectures frequently undergo instability, ejecting planets and sometimes disrupting the system; resonant chains can help in some low-multiplicity cases but typically break for higher multiplicities. The work also assesses detectability with the LISA mission, providing a formula to estimate gravitational-wave frequency shifts induced by planets and highlighting that many single-planet configurations around DWDs could be viable LISA sources, though multi-planet detections remain challenging. Overall, the results inform expectations for the survival of Magrathea planets and guide GW-based exoplanet searches in the LISA era, while acknowledging model limitations such as neglected relativistic effects and binary orbital evolution.

Abstract

Planetary formation might occur at different stages of the stellar evolution of compact binaries. In recent years, the formation of second-generation planets has been tested in circumbinary discs formed by the ejection of stellar material from double white dwarf (DWD) binaries. In these environments, planets ranging from sub-Neptunian to giant masses can form and migrate to within 1 au of the central binary. Nevertheless, studies on the orbital stability of such systems have yet to be undertaken. In this work, we use N-body simulations to study the stability of multi-planet systems formed around compact DWDs in both resonant and non-resonant configurations over timescales of a few million years. We track the occurrence of catastrophic events and employ a variety of metrics, e.g., orbital spacing, centre-of-mass variations and Normalized Angular Momentum Deficit, to explore the outcomes of their evolution. Furthermore, we evaluate the potential for detecting these systems in their final configurations with the Laser Interferometer Space Antenna (LISA) mission by measuring the overall gravitational-wave frequency shift amplitude induced by their planets. Our results show that planets orbiting DWDs can be stable over the studied timescales. While planetary systems starting with two-planets are more likely to survive unaltered, planetary systems with three, four or five planets, experience catastrophic events that cause them to lose some of their original planets, ending up hosting only two surviving planets in the majority of cases. This increases the number of two-planet systems by 122% with respect to their initial abundance and creates a single-planet population amounting to 7% of the totality of systems. The majority of these single-planet systems are potential candidates for LISA. Concerning multi-planet systems, a handful of systems could be detected.

Figures (16)

• Figure 1: Population A: schematic representation of the percentage of simulated planetary systems (1000 simulations in total for DWD$_2$ and DWD$^*_4$) based on two specific properties: whether they experience or not at least one catastrophic event (ejection/scattering/collision) and whether they are dynamically cold or hot. The number in each box is related to the percentage of systems which have the combination of parameters specified on the x an y axes. The 14.1% of systems displayed in upper right corner include the 0.2% of systems that are disrupted (see text for details).
• Figure 2: Population A: $S_m^f/S_m^s$ as a function of $S_s^f/S_s^s$ (top left panel) with a zoom-in of the region where both ratios are equal or very close to unity (top right panel) and $\Delta_{\rm CoM}$ as a function of $S_s^f/S_s^s$ (bottom left panel) with a zoom-in of the region where $S_s^f/S_s^s$ is equal or very close to unity and $\Delta_{\rm CoM}$ is equal or very close to zero (bottom right panel). DWD$_2$ and DWD$^*_4$ systems are displayed as circles and triangles, respectively. Disrupted systems are not displayed in this figure as it was not possible to compute their NAMD and $\Delta_{\rm CoM}$.
• Figure 3: DWD$^*_4$ systems (Population A): initial and final $a$ (empty green circles, filled points, respectively). Systems IDs #331-350 (on the vertical axis, top panel), System IDs #441-450 (on the vertical axis, bottom panel). Each point represents a planet and its size is proportional to the planetary mass; the color-map gives a measure of the degree of stability. Note that in some systems the green circles overlap with the filled points.
• Figure 4: Population B: schematic representation of the percentage of simulated Population B planetary systems (500 simulations) based on two specific properties: whether they experience or not at least one catastrophic event (ejection/scattering/collision) and whether they are dynamically cold or hot. The number in each box is related to the percentage of systems which have the combination of parameters specified on the axes. The 47% of systems displayed in upper right corner include the 9% of systems that are disrupted (see text for details).
• Figure 5: Population B: $S_m^f/S_m^s$ as a function of $S_s^f/S_s^s$ (top left panel) with a zoom-in of the region where both ratios are equal or very close to unity (top right panel) and $\Delta_{\rm CoM}$ as a function of $S_s^f/S_s^s$ (bottom left panel) with a zoom-in of the region where $S_s^f/S_s^s$ is equal or very close to unity and $\Delta_{\rm CoM}$ is equal or very close to zero (bottom right panel). Disrupted systems are not displayed in this figure as it was not possible to compute their NAMD and $\Delta_{\rm CoM}$.