Statistics of Magrathea exoplanets beyond the Main sequence. Simulating the long-term evolution of circumbinary giant planets with TRES.

The focus of this work is to quantify the survival rate of circumbinary planets (CBPs) in the context of the binary host evolution, from the zero age main sequence (ZAMS) to one Hubble time (i.e. the age of the Universe).

Multiple studies have shown that stars in binary and multiple systems are ubiquitous throughout the Milky Way (Duchêne & Kraus 2013; Tokovinin 2021), but very little is known about exoplanets around evolved binary systems.

Among the ∼45 confirmed circumbinary substellar objects (SSOs) known to date, seven planets orbit one WD star and eight planets have one host that has completed a first giant branch phase.

The evidence of a larger number of planetary systems in evolved binaries than single WDs has been argued not to be a coincidence. A work by Kostov et al. (2016) found that CBPs have more chances of survival when orbiting compact binary systems than single stars of similar mass.

Our work is set within the broader context of study for the development of the planetary detection science case of the LISA mission. This development includes the LISA detection prospects of SSOs orbiting DWDs through Bayesian analysis. Indeed, LISA will be able to detect giant planets through gravitational waves in the entire Milky Way (Tamanini & Danielski 2019, Danielski et al. 2019) and also in the Large Magellanic Cloud (Danielski & Tamanini 2020).

The powerhouse of our simulations is the publicly available TRES package (Toonen et al. 2016). This software simulates hierarchical triple star systems, with a thorough analytical treatment including secular orbital evolution, and various stellar interactions (e.g. tides, CE evolution, stellar winds, and supernovae and associated natal kicks). The single stellar evolution is modelled by SeBa code (Toonen & Nelemans 2013).

We further developed TRES and SeBa so that we can simulate SSOs as well.

• Criterion of Holman & Wiegert (1999) for CBP (P-type orbits) stability
• Time-dependent internal structure constants: k₂ (apsidal motion constant) and gyration radius, for both SSOs and stars from Claret (2019)
• Atmospheric photoevaporation via energy-limited mass loss due to XUV radiation
• SSO mass-radius relation from Chen & Kipping (2017) and initial spin velocity from Bryan et al. (2018)

We simulated two populations of CBPs (Pop. A / B) with 10500 systems each, differing only for the planetary priors used.

The systems are required to be secularly stable at ZAMS. We placed an upper limit of CPU time for computation to avoid cluttering the simulation of the entire population with anomalous systems. 

The simulation of those systems in which any of the components either merge, or become unbound, or initiate a phase of stable mass transfer, are automatically stopped. The same happens if the triple became dynamically unstable. Systems undergoing common-envelope phases were not stopped. We set the maximum simulation time to 13.5 Gyr (one Hubble time).

Inner binary stars

• Semimajor axis: Log₁₀-Uniform distribution for a_in ∈ (15 R_☉; 2200 R_☉ ) (e.g. Duchêne & Kraus, 2013)
• Primary mass (M₁): Kroupa’s mass function n(M) ∼ M^−2.35 , for M₁ ∈ (0.95, 10) M_☉
• Mass ratio (secondary): Uniform distribution for q = M₂/M₁ ∈ (0, 1] (e.g. Moe & Di Stefano 2017; Duchêne & Kraus 2013) with a lower limit of 0.95 M_☉ on M₂
• Eccentricity: Thermal distribution f(e) = 2e , for e_in ∈ (0, 0.95) (Heggie 1975)

CBPs (Danielski et al. 2019; Katz et al. 2022)

• Semimajor axis: Log₁₀-Uniform / Uniform (Pop A / B) distribution for a_out ∈ (0.17 au; 200 au) 
• Mass: Uniform distribution for M_P ∈ (0.2, 16) M_Jup
• Eccentricity: Beta / Uniform (Pop A / B) distribution, for e_out ∈ (0, 0.95) (Bowler et al. 2020)
• Inclination (orbit): cosine-Uniform distribution

The results of the two simulated populations show that between 23% - 32% of the CBPs survive the entire evolution and become Magrathea planets. 

They reside on rather wide orbits around tight, mostly circularised, inner binaries. CBP eccentricity is not a decisive parameter by itself to determine the Magratheas survival. They have 10% more prograde than retrograde orbits. 

The Merged category includes all systems where the inner binary stars merge into one single object, via merger or collision. They are around 1/3 of our total sample and involve heavy progenitors. 

The Destabilised planets most characterising feature is the very high eccentricity, on average higher than all the other categories. They live on the brink of instability as the binaries are wide and the CBPs weaky hierarchical. These systems became unstable for either dynamical instabilities, orbit-overlappings or SN-kicks. 

Planet survival seems linked to the stellar mass-loss kick, which is weaker for compact binaries than for single stars. In all categories, excluding the Destabilised, we see an expansion of the planetary final semimajor axis of a factor x4, due to the adiabatic stellar winds. 

Also, our results show no planetary-mass selection effect for any of the system categories, indicating that the mass does not play a role in the survival of CBPs (in the gas giant range). 

Dynamical instabilities can be a source of free-floating planets, unbinding the SSO from the host stars. Following Sutherland & Fabrycky (2016) we estimate that around 2% of our total simulated CBPs become free-floating. Direct collisions are significantly more rare, but not impossible; we find destabilisation possible even at WD stages, thus setting the prerequisites for WD pollution. 

Veras et al. (2011) showed the break-down of the adiabatic stellar wind assumption for very wide-orbit circumstellar planets or Oort-cloud analogues. Similarly, the common-envelope (CE) phase mass loss can trigger orbit disruption, on timescales of 10^4-5 yr (600-3000 au). 

Given our maximum initial orbital separation of 200 au for the planets, our simulated systems are typically not disrupted neither by CE mass loss nor by stellar winds. 

The DWD mergers happening within one Hubble time are particularly interesting as, when their period falls between 3 to 60 min, they can be detected by LISA (Korol et al. 2017). LISA has the potential to detect the giant planets (M_P > 0.2 M_Jup) whose periods are shorter than the LISA observing lifetime T_obs ≈ 4-10 yr.

We have no planet in the Merged-DWDs category with period < 4 yr: the shortest is 13.8 yr (4.6 au). 
Similarly, for the Magrathea planets the shortest period is 43 yr (9.6 au) in Pop. A, and 20 yr (6 au) in Pop. B.

This lack of planets with LISA-compatible periods might result from statistical bias due to the limited number of systems simulated. Destabilisation and dynamical chaos, for instance in planet-planet scattering, can be responsible for migrating the planet inwards. 
Furthermore, the presence of a second-generation dust disk (post-CE) can also be responsible for inward migration of the surviving planets (Perets 2010), similarly to what is theorised during planetary formation.

Main takeaways of our analysis:

• 23% to 32% of all giant CBPs survive for one Hubble time to become Magrathea planets (respectively Pop. A and B) and they tend to reside on wide orbits around light hosts.
• A large fraction (∼ 33%) of the systems end up with a merger in the inner binary.
• The CBP mass in the simulated giant planet range, does not appear to be correlated with its host system category, and does not play a role in the survival.
• Photoevaporation has a negligible impact on the majority of our giant planets, due to the adiabatic expansion of their orbits during every stellar wind phase. 
• CBPs prefer prograde orbits during their evolution, especially the Magratheas.
• Unstable systems are short-lived, they can potentially create small percentage of free-floating planets, and set the conditions for WD pollution.
• The evolution of single giant planet systems do not produce any CBP in the LISA sensitivity range. 

All results and full discussion can be found in Columba et al., 2023, A&A - DOI: https://doi.org/10.1051/0004-6361/202345843

Contact me at:  gabriele.columba@phd.unipd.it

Bowler, B. P., Blunt, S. C., & Nielsen, E. L. 2020, ApJ, 159, 63
Bryan, M. L., Benneke, B., Knutson, H. A., Batygin, K., & Bowler, B. P. 2018,Nature Astronomy, 2, 138
Chen, J. & Kipping, D. 2017, ApJ, 834, 17
Claret, A. 2019, A&A, 628, A29
Claret, A. & Torres, G. 2019, ApJ, 876, 134
Danielski, C., Korol, V., Tamanini, N., & Rossi, E. M. 2019, A&A, 632, A113
Danielski, C. & Tamanini, N. 2020, International Journal of Modern Physics D, 29, 2043007
Duchêne, G. & Kraus, A. 2013, ARA&A, 51, 269
Heggie, D. C. 1975, MNRAS, 173, 729

Holman, M. J. & Wiegert, P. A. 1999, AJ, 117, 621
Katz, M. L., Danielski, C., Karnesis, N., et al. 2022, Monthly Notices of the Royal Astronomical Society, 517, 697
Korol, V., Rossi, E. M., Groot, P. J., et al. 2017, Monthly Notices of the Royal Astronomical Society, 470, 1894
Kostov, V. B., Moore, K., Tamayo, D., Jayawardhana, R., & Rinehart, S. A. 2016, ApJ, 832, 183
Moe, M. & Di Stefano, R. 2017, ApJS, 230, 15
Perets, H. B. 2010, arXiv e-prints, arXiv:1001.0581
Sutherland, A. P. & Fabrycky, D. C. 2016, ApJ, 818, 6
Tamanini, N. & Danielski, C. 2019, Nature Astronomy, 3, 858
Tokovinin, A. 2021, Universe, 7, 352
Toonen, S., Hamers, A., & Portegies Zwart, S. 2016, Computational Astrophysics and Cosmology, 3, 6
Toonen, S. & Nelemans, G. 2013, A&A, 557, A87
Toonen, S., Nelemans, G., & Portegies Zwart, S. 2012, A&A, 546, A70

Veras, D., Wyatt, M. C., Mustill, A. J., Bonsor, A., & Eldridge, J. J. 2011, MNRAS, 417, 2104