# The gravitational-wave detection of exoplanets orbiting white dwarf binaries using LISA

Nicola Tamanini, Camilla Daniëlski

► To cite this version:

Nicola Tamanini, Camilla Daniëlski. The gravitational-wave detection of exoplanets orbiting white dwarf binaries using LISA. Nature Astronomy, 2019, 3 (9), pp.858-866. 10.1038/s41550-019-0807-y. hal-01965353

HAL Id: hal-01965353

https://hal.science/hal-01965353v1

Submitted on 8 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The gravitational-wave detection of exoplanets orbiting white dwarf binaries using LISA

Nicola Tamanini<fnref n="1" /> and Camilla Daniëlski²,³

So far, around 4,000 exoplanets have been discovered orbiting a large variety of stars. Owing to the sensitivity limits of the currently used detection techniques, these planets populate zones restricted either to the solar neighbourhood or towards the galactic bulge. This selection problem prevents us from unveiling the true galactic planetary population and is not set to change for the next two decades. Here, we present a detection method that overcomes this issue and that will allow us to detect massive exoplanets using gravitational-wave astronomy. We show that the Laser Interferometer Space Antenna (LISA) mission can characterize new circumbinary exoplanets orbiting white dwarf binaries everywhere in our Galaxy—a population of exoplanets so far completely unprobed—as well as detecting extragalactic bound exoplanets in the Magellanic Clouds. Such a method is not limited by stellar activity and, in extremely favourable cases, will allow LISA to detect planets down to 50 Earth masses.

In the last 20 years the field of extrasolar planets has witnessed an exceptionally fast development, revealing an incredibly diverse menagerie of planetary companions. These discoveries changed the place that our Solar System occupies in the galactic context and allowed us to develop a deeper understanding of the planetary population around us. For instance, we now know that hot Jupiters are rare, super-Earths are ubiquitous<fnref n="1" /> and that there is a gap in the radius distribution of small planets². Nonetheless, our knowledge is restricted to the solar neighbourhood because the most successful detection techniques, such as radial velocity and transit, can only observe bright stars close to us. Differently, gravitational microlensing can observe farther away towards the galactic bulge, but it does not provide enough targets to develop a robust statistics of the bulge population. Determining if what we see is as a result of a selection bias or not is extremely important, and it cannot be assessed through the usual techniques for at least the next two decades.

## Exoplanets beyond the solar neighbourhood

A breakthrough in the detection of exoplanets beyond the solar neighbourhood is going to be possible only through gravitational wave (GW) astronomy. The Laser Interferometer Space Antenna (LISA) mission³, planned for launch in the early 2030s, will enable us to indirectly probe for the first time the population of giant exoplanets orbiting detached double white dwarf (DWD) binaries, everywhere in the Milky Way and in the nearby Magellanic Clouds. Given that about half of the stellar population resides in multiple stellar systems⁴,⁵, and that ~95% of the stars will become a white dwarf⁶, the LISA circumbinary planets survey will shed light on the final fate of an exoplanet as well as provide a galactic statistic of these objects.

Among the 90+ circumbinary systems currently known in the solar proximity, two tens are P-type, meaning systems with planet(s) orbiting both stars in the binary. Of these, only six systems have a white dwarf as a binary component. Usually, the second companion is a low-mass star⁷ or, for one specific case, a pulsar⁸. Today, no exoplanets have been discovered around DWDs, irrespective of the compactness of the binary.

Owing to the intrinsic faintness of these DWDs, only few tens are known by spectroscopic and variability surveys⁹, but substantial progress in the detection of these sources is expected via GWs. The LISA mission, working in the gravitational-wave low-frequency range between 0.1 mHz and 1 Hz, will detect around $25 \times 10^3$ compact DWDs within and outside the Milky Way⁹,¹⁰, some of which could be perturbed by the presence of a third gravitationally bound stellar companion¹¹ or planetary companion.

The existence of this planetary population in our Galaxy is far from being excluded⁷,¹², yet there is no observational proof of its presence. We show here that LISA represents an important step forward in the science of planetary formation and evolution. In case of positive planetary detections, LISA will provide crucial elements: (1) to understand under which conditions a planet can survive the most critical phases of a close binary evolution¹²,¹³; (2) to set constraints on binary mass-loss and the dynamical aspects that directly follow¹²,¹⁴; and (3) to check whether a second generation of planets exists, that is, planets that form from the stellar material ejected during the binary common-envelope phase(s)¹⁵. Conversely, in case of no detection all over the Milky Way, LISA will enable us to set unbiased constraints on planetary evolution theories, and in particular on the fate of exoplanets bound to a binary that undergoes two common-envelope phases¹².

In what follows we assess the potential of LISA to discover new P-type circumbinary exoplanets through their perturbation on the GW signal emitted by the DWD. This method is fully original and relies on the large DWD population of GW sources to be heard by LISA, which makes it more powerful and interesting than other ideas of direct detection of exoplanets through GWs¹⁶⁻²⁰. Finally, in the spirit of the new era of multi-messenger astronomy, we discuss the possibilities that could open for the field of exoplanets when standard electromagnetic (EM) techniques will work in synergy with GW astronomy.

## Detecting exoplanets with LISA

Galactic DWDs with periods shorter than 1 h emit almost monochromatic GWs in the LISA frequency band. If a third companion

<footnote n="1">¹Max-Planck-Institut für Gravitationsphysik, Albert-Einstein-Institut, Potsdam-Golm, Germany. ²AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, Gif-sur-Yvette, France. ³Institut d'Astrophysique de Paris, CNRS, UMR 7095, Sorbonne Université, Paris, France. e-mail: nicola.tamanini@aei.mpg.de; camilla.danielski@cea.fr</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1 | Geometry of the DWD-CBP three-body system.** a, Source reference frame (primed coordinates). b, Observer reference frame (non-primed coordinates) with the source reference frame rotated.</figcaption></figure>

object orbits the DWD then the centre of mass (CoM) of the DWD is on a Keplerian orbit. Through the Doppler effect, the motion of the DWD CoM induces an observable imprint on the GW wave-form measured by LISA¹¹. See Methods for the dynamical model employed to describe the three-body system (Fig. 1) and details regarding the Fisher matrix parameter estimation with LISA.

The three parameters of the circumbinary planet (CBP), which can be recovered with LISA, are the planetary period *P*, its initial phase *φ*₀ and the parameter

$$K = \left( \frac{2\pi G}{P} \right)^{\frac{1}{3}} \frac{M_p}{(M_b + M_p)^{\frac{2}{3}}} \sin i \quad (1)$$

which depends on the CBP mass *M*<sub>p</sub>, its orbit inclination *i*, and the binary total mass *M*<sub>b</sub>. *G* is the gravitational constant. In Fig. 2 we show how the 1σ relative uncertainties on *K* and *P*, and the 1σ absolute uncertainties on *φ*₀ (in radians), vary as a function of *P*. We consider two values for the GW frequency emitted by a DWD: a representative frequency at *f*₀ = 1 mHz (dashed lines), at which the majority of DWDs detected by LISA are expected⁹, and a higher frequency *f*₀ = 10 mHz (solid lines), where events with higher signal-to-noise ratio (SNR) will be detected. All numbers are linearly scaled with respect to the SNR of the DWD detected by LISA and the mass of the CBP. This implies that to find the precision with which the parameters are measured, the numbers reported in Fig. 2 must be divided by the SNR of each individual DWD event and by the CBP mass, measured in *M*<sub>1</sub> (Jupiter's mass).

Figure 2 shows that the error estimations on *K* and *φ*₀ are better for planetary periods *P* comparable to the LISA nominal lifetime. For shorter planetary periods the errors on *K* and *φ*₀ smoothly decrease with the increasing of *P*, while for planetary periods longer than LISA's nominal duration the precision steeply worsens by a few orders of magnitude. On the other hand, the uncertainty on *P* increases smoothly with the planetary period up to the LISA nominal mission lifetime, after which it rapidly worsens similarly to the behaviour of *K* and *φ*₀. An explanation for the different behaviour of *P*, as well as for the appearance of the peak at 1 year, is provided in the Methods.

Analogously to EM radial velocity techniques, we cannot recover the mass of the companion directly from the three GW parameters *K*, *P* and *φ*₀, implying that we have no means of knowing if the perturbing object is a planet, a star or a different object. Nevertheless, if both *K* and *P* are measured with sufficient precision, then additional EM information on the total DWD mass and on the companion's

<figure><img src="image_3.png" /><figcaption>**Fig. 2 | LISA estimation of planetary parameters.** Relative accuracy on *K* and *P*, and absolute accuracy on *φ*₀, for LISA measurements of planetary orbital parameters. Solid curves denote DWD with *f*₀ = 10 mHz, whereas dashed curves are for *f*₀ = 1 mHz. All results have been scaled with respect to the DWD SNR and the planetary mass measured in Jupiter masses. The peak at 1 year is caused by degeneracies in the GW signal appearing when the planetary orbit is an exact multiple of LISA's orbital period around the Sun.</figcaption></figure>

orbit inclination can help in determining the mass of the companion through equation (1). It is, furthermore, possible to obtain some estimates on the allowed *M*<sub>p</sub> range if the chirp mass of the binary is also measured with LISA, which is usually the case at least for high-frequency DWDs²¹. In fact, by assuming that the symmetric mass ratio *η* of DWDs detected by LISA cannot be lower than a certain value, we can derive upper and lower bounds for the total DWD mass *M*<sub>b</sub>, and from this range of values we can estimate *M*<sub>p</sub>sin(*i*) through equation (1). Since from GWs alone we do not obtain any information on *i*, we can only derive lower limits for *M*<sub>p</sub> by considering sin(*i*) = 1. Values of sin(*i*) different from 1 will yield higher values of *M*<sub>p</sub>. If the lower mass limit estimated in this way is below 13 *M*<sub>1</sub> (the deuterium burning limit), then the companion object is probably an exoplanet that would need to be confirmed by EM follow-ups (see 'Synergy with EM observations and implications'). On the other hand, if the lower bound of this estimate is above 13 *M*<sub>1</sub>, then we can confidently exclude the possibility that the third orbiting object is a planet. Furthermore, if we also set a minimum
<figure><img src="image_4.png" /></figure>

**Fig. 3 | Selection functions of both LISA and of main EM exoplanetary projects.** Selection functions of LISA (gravitational waves, red lines): planetary mass $M_p$ versus distance from the (binary) star. The DWD GW frequency $f_0$ and the LISA SNR are marked in the legend. The typical selection curves of other exoplanet detection experiments are plotted for comparison: Gaia (astrometry, green) line$^{76,77}$ is measured for a $0.5 M_\odot$ M-dwarf at 25 pc; VLT-ESPRESSO (radial velocity, cyan) line$^{52}$ is measured for a velocity amplitude of 10 cm s$^{-1}$ for a $0.8 M_\odot$ star; ARIEL (transiting, orange) line$^{53}$ is based on the Tier 1 sample (that is, separations lower than 1 au); WFIRST (microlensing, brown) line$^{54}$ corresponds to the $3\sigma$ detection; ELT-PCS (imaging, dark blue) line$^{55}$; JWST (imaging, lime) line$^{56}$; and SPHERE (imaging, purple) line$^{55}$. Note that some of the lines overlap. Grey markers correspond to the currently known exoplanets as retrieved from https://exoplanetarchive.ipac.caltech.edu. The dashed horizontal line marks the deuterium burning limit at $M_p = 13 M_\odot$.

allowed value for $\sin(i)$, we can find an upper bound on $M_p$ that would confirm the presence of an exoplanet from the GW signal alone if it does not exceed $13 M_\odot$. In this case, no EM counterpart would be needed to confirm the GW detection of a CBP.

Given the discussion above, we are particularly interested in systems for which both *K* and *P* can be measured with a relative statistical uncertainty better than some accuracy value that we arbitrarily set to 30%. In what follows we will define a 'detection' of a CBP with GWs if our estimated precision on both *K* and *P* is better than this number. Figure 3 shows the region in the mass-separation parameter space where LISA will have the possibility of detecting CBPs, according to the definition of detection above. We see that LISA will be more effective at observing exoplanets with a separation from the DWD between roughly 1 and 2 au (Fig. 3). Depending on the SNR and frequency of the DWD, LISA will be able to detect exoplanets down to ~1 $M_\odot$ and will be efficient in the range between 0.01 and 3 au, roughly.

Although the occurrence rate of CBPs around DWDs is at present unknown, if these planets physically exist, and they orbit even a few per cent of DWDs detectable by LISA, we will observe up to hundreds of new exoplanets, assuming they are sufficiently massive to be heard.

Our results are obtained for binaries composed of two white dwarfs, each of mass $M_\star = 0.23 M_\odot$ (chirp mass $M_c = 0.2 M_\odot$). For DWDs with higher (chirp) masses the perturbation due to a planet will in general be weaker, although the SNR of the binary will be higher for equal distance, symmetric mass ratio $\eta$ and GW frequency. The same reasoning applies to different types of stellar binaries visible with LISA, such as neutron star-neutron star and white dwarf-neutron star binaries, for which sufficiently massive CBPs could be observed as well.

In our analysis we assume for simplicity no eccentricity for both stellar and planetary orbits. Note, however, that the eccentricity of

the planetary orbit constitutes another measurable parameter with LISA, if it differs significantly from zero$^{11}$. For eccentric DWDs, which may form in the presence of a hierarchically bound third companion$^{22,23}$, it might also be possible to measure the individual masses of the stars directly from GW observations$^{24}$, yielding a more accurate and precise estimation of the CBP mass without the need for complementary EM observations.

Moreover, we assume the DWDs to be detached, namely with no accretion. LISA is also expected to detect a few thousands of accreting DWDs$^{25}$, for which, however, it might be difficult to directly measure the chirp mass with GWs and to disentangle an eventual CBP signal from accreting effects. Nevertheless, individual masses can instead be recovered with joint GW-EM observations, for instance with LISA-Gaia for which ~50 such events can be characterized$^{26}$. Furthermore, accreting DWDs are expected to emit specific EM signatures that might facilitate the identification of an EM counterpart$^{27}$, enhancing the multi-messenger potential of these systems.

In our analysis we assume a 4 year nominal LISA mission duration$^3$. However, for a maximal extension of 10 years, we expect the number of DWD detections to roughly double$^9$, and the SNR of each individual DWD to increase as the square root of time$^{21,28}$. This translates into both a higher number of CBP detections with GWs, and a more precise characterization of the systems already observed in 4 years. A 10 year LISA mission will also be more sensitive to CBPs with periods between 4 and 10 years. In fact, we expect the rapid deterioration of the accuracy with which the planetary parameters are recovered present at periods longer than 4 years (Fig. 2) to shift towards longer periods, in particular around 10 years.

Our results are always rescaled with the SNR of the DWD, implying that all effects due to the distance of the GW source from the Earth have been factored out. LISA will detect DWDs at all distances within the Milky Way, and even at intergalactic distances up
to M31¹⁰ (Andromeda galaxy). Note, however, that DWDs at farther distances will generally be observed with a lower SNR, implying that the largest majority of CBPs detected by LISA will nonetheless be located close to us. Although only systems in the solar neighbourhood ($d. \lesssim 3$ kpc) can be repeatedly observed with joint GW and EM observations, CBPs farther away can be detected by GW alone. Given that no EM detection of extragalactic planets bound to their star(s) has been confirmed yet (see refs. ²⁹,³⁰ for detection of the population of unbound planets), LISA will open a new window for the quest of extragalactic exoplanets through GW astronomy.

The other advantage of having produced SNR rescaled results is that we can directly make predictions not only for different DWD systems detected with the same LISA configurations, but also assess the potential for exoplanetary searches of other LISA-like missions. For instance, for DWDs with SNR=100 and $f_0=10$ mHz, LISA will be able to hear CBPs of $M_p \ge 1 M_\odot$. Some of these higher-frequency DWDs will be measured with SNR as high as 600, meaning that, for these specific systems, LISA will be able to detect CBPs down to $50 M_\odot$ (~0.16 $M_\odot$). Similarly, a future more sensitive GW observatory, which can improve LISA SNR values of at least one order of magnitude, will be able to observe CBPs down to a few Earth masses. Note, however, that our Fisher matrix analysis does not take into account the discreteness of LISA frequency measurements. In particular, low-mass CBPs might in reality provide a frequency drift in the GW signal that is comparable to, or lower, than the frequency resolution of LISA. This implies that results obtained for high SNRs can only be considered as an optimistic estimation at best, and that further analyses will be needed to refine these predictions.

## Synergy with EM observations and implications

Observing the EM counterpart of a GW exoplanet detection is necessary to fully characterize a planetary system, in term of stars plus planet.

Recent estimates⁹ suggest that out of the $25 \times 10^3$ DWDs expected to be detected in 4 years by LISA (see Fig. 4 for the location in the Galaxy of such a population³¹), only up to a maximum of 100 will be observed with instruments such as Gaia and the Large Synoptic Survey Telescope (LSST). Large upcoming ground-based facilities might improve these numbers. When a GW event is heard (see 'Detecting exoplanets with LISA'), EM observations of the same system allow us to measure the binary mass. Moreover, if the mass of each binary component is determined, it is possible to directly estimate their stellar radii and ages by using the mass-radius relations⁷ and the age-luminosity law³², respectively (the latter only if the stars are resolved). We remark that the identification of which DWD has a planetary companion, within the LISA sky error box, is left for future studies.

Concerning the planetary characterization, the CBP itself has its typical signature that can be spotted in the EM spectrum by various well-known detection techniques (Fig. 3), which, on average, mainly observe systems relatively close to us. Transit, radial velocity, direct imaging and astrometry methods cover, on average, distances up to $d. \sim 3$ kpc, while gravitational microlensing goes up to $d. \sim 8$ kpc (see Fig. 4 for the location in the Galaxy of the known exoplanets host star). Among these, while the transit and radial velocity methods are affected by stellar activity, direct imaging is limited by the angular resolution of the instrument. Microlensing events are not reproducible and astrometry is severely limited by the time span of the survey. The CBP detection method using GWs proposed here is conceptually similar to the radial velocity method, but has the advantages that it can observe everywhere in the galaxy, is not affected by the activity of the stars and does not need any observational pointing.

After a CBP is heard via GWs, and if we assume that the planetary system complies with the requirements specific for an EM follow-up (for example magnitude and system configuration), joint

GW-EM observations will allow us to resolve the circumbinary system in terms of periods, masses (of the stars+planet) and inclination of the orbits. See Methods for more details on the parameters that can be retrieved in the specific cases of synergy between GWs and radial velocity, transit, astrometry, direct imaging and gravitational microlensing.

As previously mentioned, a 10 year mission can detect and characterize exoplanets on wider orbits, making LISA even more compatible with EM direct imaging techniques, and specifically with large ground-based cameras (for example, the Planetary Camera and Spectrograph for the European Extremely Large Telescope, ELT-PCS, Fig. 3) or the successor of the James Webb Space Telescope (JWST). Imaging of CBPs around DWDs can be used to test the presence of a second generation of exoplanets in the outer regions of a planetary system, and consequently to provide constraints on migration theories. Emission spectra of these objects will furthermore allow us to estimate their temperature and the main molecular component of their atmosphere, making direct connections to chemical element distributions in the atmosphere of white dwarfs. This would also allow to better understand the observed white dwarf pollution effect.⁷ On the other hand, if an existing CBP accretes mass after a common-envelope stage, it becomes brighter, further decreasing the already low planet-to-white dwarfs contrast, meaning that also first-generation, more mature exoplanets can be imaged.

A peculiarity of white dwarfs is their size. Their small radius, $R_{WD} \approx 1 R_\odot$, makes them excellent targets for detecting transiting exoplanets, and in particular Earth-like planets³³. However, Jupiter-like planets, which are a factor ~7 larger than a typical white dwarf, generates a complete eclipse during their transit, fully occulting the stellar light. This facilitates their detection, but prevents atmospheric studies on the exoplanet in question through transmission spectroscopy. For specific configurations, though, atmospheric feature information could be retrieved by using an external background source. If the planet is in close-in orbit and has an exomoon (assuming it is not removed through an ejection resonance while undergoing the planetary migration phase³⁴), then transit spectroscopy will allow to study the atmosphere of the satellite, due to the small size of white dwarfs. For instance, an Europa transit ($R_E \sim 0.24 R_\odot$) in front of a typical white dwarf ($M_* = 0.6 M_\odot$, $R_{WD} = 1.57 R_\odot$) yields a decrease in brightness of approximately ~5%.

## Considerations

CBPs orbiting DWDs are proven to exist in theory. One theoretical example is given by the evolution of the Kepler 1647 system (two solar-mass stars of F- and G-type with a circumbinary gas giant companion³⁵), which could, under a specific scenario, become a white dwarf-white dwarf pair that keeps a bound exoplanet after two common-envelope phases¹². It is important to mention though that after one inner binary common-envelope phase, a CBP can either be ejected by the system⁷,¹⁴ or survive the common-envelope phase(s) while undergoing an orbital expansion/shrinkage, generally coupled with an increased eccentricity¹². The kind of binary evolution, that is, tidal evolution, the energy transfer between the common envelopes and the inner cores, and the mass-loss rate, is directly responsible for the fate of their planetary companion. The handful of currently known post-common-envelope exoplanets, together with the many unknowns related to the complex common-envelope phase³⁶, prevent us from examining the parameters' space expected to be found for CBPs orbiting a DWD. Especially when effects such as planetary migration (see for example ref. ³⁷ and references therein), second-generation formation¹⁵ and the type of common-envelope expansion³⁸ might play a key role. The detection and characterization of CBPs around DWDs is consequently crucial for pinning down the planetary evolution phases while its binary evolves, and for providing observational constraints of the stellar common-envelope phase.
<figure><img src="image_5.png" /></figure>

**Fig. 4** | Location of known planets and expected LISA DWDs. Left: The plot represents both known stars hosting planetary companion(s) (triangles, retrieved from https://exoplanetarchive.ipac.caltech.edu) and the number of DWD detections expected with LISA³¹ (colour bar on the right). Right: Zoom-in on the Solar region. The top legend shows the exoplanet detection method used for the discovery and its colour code. Note that the white dwarf binary population is built assuming spherical symmetry³¹, and that data overlay a face-on black-and-white image of the Milky Way for galactic location reference purposes. We also highlight the known CBPs (P-type, red triangles) and the horizon (500 pc, white circle) for exoplanet detection⁵⁷ of the Gaia nominal mission (5 years). Dashed lines represent the galactocentric coordinates of the Sun. Credit for the background image of the Milky Way in both panels: NASA / JPL-Caltech / R. Hurt (SSC/Caltech).

We note that the orbit of a CBP is stable when the semi-major axis *a* is approximately *P* ≥ 4.5 *P*<sub>b</sub> for circular binaries (CBPs are not expected to exist below this limit<sup>39,40</sup>), which in our specific case corresponds to a minimal planet–binary separation of *a*<sub>min</sub> ~ 0.0007 au and *a*<sub>min</sub> ~ 0.003 au when accounting for a chirp mass of *M*<sub>c</sub> = 0.2 *M*<sub>⊕</sub> (that is, *M*<sub>b</sub>=0.46 *M*<sub>⊕</sub>), and for a binary period *P*<sub>b</sub>=3 min and *P*<sub>b</sub>=30 min, respectively. These distances are well below the minimal separation considered for CBPs in this work, since we always assume *P* ≫ *P*<sub>b</sub> in order for higher-order orbital effects, for example Kozai–Lidov resonances, to be negligible.

White dwarfs in general are common as Sun-like stars and may also provide a source of energy for planets for gigayear (Gyr) durations. Cool white dwarfs (T<sub>eff</sub> < 6,000 K) have an habitable zone that endures for up to 8 Gyr (ref. ³³) located at only few solar radii from the star. This habitable zone regularly moves inwards as the star ages, and hence cools³³. Planets located in the habitable zone of a white dwarf must have surely migrated inwards after the stellar giant phase⁴¹⁻⁴³. The detection of such objects, or any close-in object not necessarily in the habitable zone, would help constraining migration theories of planets around post-common-envelope binaries. For a binary case and with hotter DWDs (both at T<sub>eff</sub> > 6,000 K), the habitable zone scales outwards to account for the energy contribution of both stars and their temperatures.

Conclusion

We have presented an original observational method that employs GWs to detect exoplanets. The conceptual idea is similar to the radial velocity detection technique, but it is notably not affected by stellar activity. Our results show that LISA will allow us to verify the existence of exoplanets orbiting detached white dwarf binaries in our Galaxy, as well as in other nearby galaxies. The discovery of such objects will statistically increase the actual sample of post-main-sequence planets, filling an area of the planetary Hertzsprung–Russell (HR) diagram currently not explored. Such a population will

be unbiased and valid all over the Milky Way. Specifically, LISA will provide observational constraints on both planets that can survive two common-envelope stellar evolution phases and on a possible second-generation (or maybe third) planet population. Depending on the SNR of the observation, LISA will have the potential to detect exoplanets with mass down to 50 *M*<sub>⊕</sub>. On the other hand, in a scenario where LISA will not detect any of these circumbinary planets anywhere in the Milky Way, we will still be able to set strong unbiased constraints on planetary evolution and dynamical theories.

Given that in the next few years the Transiting Exoplanet Survey Satellite (TESS) is expected to observe ~500,000 bright eclipsing binaries⁴⁴, a hefty increase in the number of observed P-type circumbinary planets, mostly around main-sequence binaries, is expected. These discoveries are important for collecting a more robust statistics about such systems in the solar neighbourhood, offering the unique opportunity to improve our understanding of the planets' dynamics in binary systems. This is particularly relevant for close (≤10 au) binary systems, the ancestors of the systems observable by LISA, whose stellar evolution is qualitatively and quantitatively different from wide-orbit binaries. Together with these upcoming evidences, further developments in planet-star interactions and binary evolution theoretical studies (for example refs. ¹²,⁴⁵,⁴⁶) need to be performed in the forthcoming years, in order to thoroughly assess the potential of LISA to detect exoplanets, especially before the European Space Agency (ESA) definitively approves the mission design in the early 2020s³. In fact, although the low-frequency part (<4 mHz) of the LISA sensitivity curve is not expected to deteriorate thanks to the achievements of LISA Pathfinder⁴⁷,⁴⁸ (from which one could also hope for a better low-frequency performance of the instrument), the high-frequency part (>4 mHz) could still evolve depending on the forthcoming technological developments and choices.

The characterization of circumbinary exoplanets detected by LISA, through both GW and EM observations, will deepen our
knowledge of the planet-star interaction during the critical phases of stellar evolution, and will enable us to find the connecting thread that runs from formation through to the very end of an extrasolar planet. GWs are probably not the answer to the ultimate question of life, the Universe and everything, but they might, after all, constitute the key to find Magrathea⁴⁹.

## Methods

**Dynamical model and modified GW phase.** The three-body systems we are interested in are composed of a DWD emitting GWs in the LISA frequency band and a planet orbiting the DWD on an outer orbit. We assume that the separation between the planet and the DWD is much greater than the separation between the two stars forming the DWD. Moreover, the period of the planetary orbit is always assumed to be much longer than the period of the inner DWD orbit. These assumptions imply that in first approximation these three-body systems can be treated as two separate two-body problems. In particular, the internal orbit of the DWD is not perturbed by the CBP, and the planet and the CoM of the DWD are orbiting each other on Keplerian orbits. For the sake of simplicity in our investigation we consider both these orbits as circular. The extension to elliptic orbits should be straightforward and will be left for future more in-depth investigations.

Given these assumptions, in the $(x', y', z')$ reference frame where the direction $\hat{z}'$ is perpendicular to the planetary orbital plane (Fig. 1), the distance vector between the planet and the CoM of the DWD is given by

$$ \mathbf{r}(t) = (R \cos \varphi(t), R \sin \varphi(t), 0) \quad (2) $$

where $t$ is time and $R$ and $\varphi$ are

$$ R^3 = G(M_b + M_p) \left(\frac{P}{2\pi}\right)^2 \quad \text{and} \quad \varphi = \frac{2\pi t}{P} + \varphi_0 \quad (3) $$

with $M_b, M_p, P$ and $\varphi_0$ the total mass of the binary, the mass of the planet, the period and the initial phase of the planetary orbit, respectively. To find the motion in the $(x,y,z)$ reference frame of a general observer whose $\hat{z}$-direction points towards the source, we apply two rotations: a first rotation by $i$ (the inclination angle between $\hat{z}'$ and $\hat{z}$—the line of sight) to bring the $\hat{z}'$ direction to point towards the $\hat{z}$ direction, and a subsequent second rotation by $\sigma$ around the $z'$ axis (which after the first rotation coincides with the $z$ axis), making $\hat{x}'$ and $\hat{y}'$ to point in the same directions of $\hat{x}$ and $\hat{y}$ (Fig. 1). Note, however, that the second rotation around $\hat{z}$ is degenerate with a phase redefinition of the planetary orbit, that is, it is degenerate with $\varphi_0$ and can thus be ignored. This holds as long as one assumes the CBP orbit to be circular and ceases to be true as soon as non-zero eccentricity is considered. The motion in the reference frame of a general observer is thus given by

$$ \mathcal{R}(i)\mathbf{r} = (R \cos \varphi(t), R \cos i \sin \varphi(t), -R \sin i \sin \varphi(t)) \quad (4) $$

where $\mathcal{R}$ is a rotation matrix by the angle $i$.

For what concerns our scope, we are interested in the circular motion of the CoM of the DWD around the common CoM of the three-body system. The distance vector connecting the DWD CoM to the three-body system CoM is

$$ \mathbf{r}_b = \frac{M_p}{M_b + M_p} \mathcal{R}(i) \mathbf{r} \quad (5) $$

We are only interested in the motion along the line of sight from the observer to the system (for any practical purposes the direction pointing from the observer to the CoM of the three-body system and the direction pointing from the observer to the DWD CoM will be assumed to coincide). Since the $z$ axis of the observer is aligned along the line-of-sight direction, the $z$ component of the motion, given by

$$ z_b = -\frac{M_p}{M_b + M_p} R \sin i \sin \varphi(t) \quad (6) $$

is the one we need to consider. The velocity of the DWD CoM along the line of sight is then given by

$$ v_{z,b} = -K \cos \varphi(t) \quad (7) $$

where we defined the parameter

$$ K = \left( \frac{2\pi G}{P} \right)^{\frac{1}{3}} \frac{M_p}{(M_b + M_p)^{\frac{2}{3}}} \sin i \quad (8) $$

In the reference frame comoving with its CoM, the DWD emits almost monochromatic GWs at some specific frequency $f_{GW}$. Since this frequency is

changing on timescales much longer if compared with the observational timescale, an expansion around the initial observed frequency will suffice in describing their time evolution²¹. We can thus assume

$$ f_{GW}(t) = f_0 + f_1 t + \mathcal{O}(t^2) \quad (9) $$

where $f_0$ is the initial observed frequency, $f_1$ is the time derivative of $f_{GW}$ evaluated at the initial time and we neglected second- and higher-order terms. The GW frequency in the observer reference frame changes instead due to the Doppler effect that, as long as the dynamics is non-relativistic, gives

$$ f_{\text{obs}}(t) = \left( 1 + \frac{v_{z,b}(t)}{c} \right) f_{\text{GW}}(t) \quad (10) $$

where $c$ is the speed of light. Finally, the phase at the observer of the GW can be obtained by integrating the frequency $f_{\text{obs}}$:

$$ \Psi_{\text{obs}}(t) = 2\pi \int f_{\text{obs}}(t') dt' + \Psi_0 \quad (11) $$

where $\Psi_0$ is a constant initial phase. This is thus the GW phase measured by LISA from which the parameters characterizing the DWD and the planetary perturbation can be extracted. Note that there should also be an additional correction to equation (11) since the time between the observer and the source is related by $t_{\text{obs}} = t_{\text{src}} + z_b/c$. This is, however, negligible for the systems considered in this work¹¹.

**LISA parameter estimation.** The three arms of LISA constitute a pair of two-arm detectors outputting two linearly independent signals $H_{1,II}(t)$. Assuming that the noise in each independent channel is stationary and Gaussian, these two signals can be expressed in the common amplitude-and-phase form as²⁸

$$ h_{1,II}(t) = \frac{\sqrt{3}}{2} A_{1,II}(t) \cos[\Psi_{\text{obs}}(t) + \Phi_{1,II}^{(p)}(t) + \Phi_{D}(t)] \quad (12) $$

where

$$ A_{1,II}(t) = [A_{+}^2 F_{1,II}^{+2}(t) + A_{\times}^2 F_{1,II}^{\times2}]^{1/2} \quad (13) $$

$$ \Phi_{1,II}(t) = \tan^{-1} \left( \frac{A_{\times} F_{1,II}^{\times}(t)}{A_{+} F_{1,II}^{+}(t)} \right) \quad (14) $$

$$ \Phi_D(t) = \frac{2\pi f_D(t) R_{\text{Earth}}}{c} \sin\theta_S \cos\left(\bar{\varphi}_0 + \frac{2\pi t}{P_{\text{Earth}}} - \varphi_S\right) \quad (15) $$

In these expressions $R_{\text{Earth}}=1$ au and $P_{\text{Earth}}=1$ year are the mean distance from the Sun and the orbital period of the Earth. $F_{1,II}^{\times\times}$ are the antenna pattern functions depending on the source angular position $(\theta_S, \varphi_S)$, the orientation of its orbit $(\theta_L, \varphi_L)$ and the LISA configuration. The quantities $A_{\times,\times}$ are instead constant amplitudes that depend on the physical parameters and orientation of the source. Equations (12)–(15) provide the GW signal measured by LISA irrespectively of the sky location of the source, that is, for arbitrary values of $(\theta_S, \varphi_S)$. The full expressions for all the quantities appearing in these equations can be found in the literature²¹,²⁸,⁵⁰. We basically follow the set up described in ref. ²⁸, adding to the GW signal the perturbation due to the CBP.

To extract information on the parameters characterizing the three-body system under consideration, we employ matched filtering techniques²¹,²⁸. We assume that the GW signal measured by LISA has the time-dependent expression given in equation (12). This is characterized by 11 parameters, namely $\ln(A), \Psi_{i,0}, f_{i0}, f_1, \theta_S, \varphi_S, \theta_L, \varphi_L, K, P, \varphi_0$ that we collectively denote by $\lambda_i$. The SNR of the signal can be computed as

$$ \text{SNR}^2 = \frac{2}{S_n(f_0)} \sum_{a=1,II} \int_0^{T_{\text{obs}}} dt \, h_{\alpha}(t) h_{\alpha}(t) \quad (16) $$

where $T_{\text{obs}}$ is the LISA observational time and $S_n(f_0)$ is the LISA one-sided spectral density noise computed at $f_0$. In all computations we fix $T_{\text{obs}} = 4$ years in agreement with the nominal mission requirements¹. For all the systems we consider here, the SNR as defined in equation (16) is approximately the same as the SNR of an equivalent DWD without the perturbation of the planet. For this reason in our analysis we freely compare this SNR with values reported in the literature for non perturbed DWDs. For GW signals with high SNR, the uncertainties and correlations on the parameters can instead be estimated from the covariance matrix $\Sigma_{ij}$ given by the inverse of the Fisher information matrix

$$ \Sigma_{ij} = \langle \Delta\lambda_i \Delta\lambda_j \rangle = (\Gamma^{-1})_{ij} \quad (17) $$
The standard estimator for the statistical error on the parameter $\lambda_i$ is thus given by $\sqrt{\sum_{i=1}^n \bar{\sigma}_i^2}$. The Fisher information matrix itself can be computed as

$$ \Gamma_{ij} = \frac{2}{S_n(f_0)} \sum_{\alpha=1,\text{II}} \int_0^{T_{\text{obs}}} dt \frac{\partial h_{\alpha}(t)}{\partial \lambda_i} \frac{\partial h_{\alpha}(t)}{\partial \lambda_j} \quad (18) $$

Note that we are simplifying equations (16) and (18) by considering $S_n(f_0)$ to be a constant and thus by taking it out of the time integral. This approximation is justified for almost monochromatic signals whose frequency does not depart considerably from $f_0$. Consequently, in analogy with the works of refs. 21,28, in our analysis we always replace $S_n(f_0)$ with the SNR of the source by using equation (16). This allows us to scale all results with the SNR of each event, without referring to any particular configuration of LISA.

The accuracy with which the parameters associated with the DWD orbit, namely $\Psi_0, f_0, f_1, \theta_S, \phi_S, \theta_L, \phi_L$, can be observed by LISA has already been investigated in several works$^{21,28,51}$. For this reason we focus our analysis on exploring the possibility of measuring the additional parameters coming from the perturbation due to the planet, namely $K, P$ and $\phi_0$. Moreover, we restrict the sampled parameter space by considering only specific values for some of the DWD parameters. In analogy with the examples in refs. 21,28, we set the orbital geometry of the DWD and its sky location by assuming

$$ \Psi_0 = 0; \quad \theta_S = \arccos(0.3); \quad \phi_S = 5; \quad \theta_L = \arccos(-0.2); \quad \phi_L = 4 \quad (19) $$

Different choices of the parameters above do not alter the qualitative conclusions derived in our investigation, although clearly the DWD SNR will change if these parameters change. By considering representative values from the expected population of DWDs observed by LISA, we assume equal mass DWDs with a chirp mass of $M_c = 0.2 M_\odot$. The value of $M_c$, together with the frequency $f_0$, yields $f_1$ as

$$ f_1 = \frac{96}{5} \pi^{8/3} f_0^{11/3} \left( \frac{GM_c}{c^3} \right)^{5/3} \quad (20) $$

The value of the amplitude A is irrelevant for our analysis since we always report results in terms of SNRs. The most favourable orientation for the detectability of the planetary perturbation on the GW signal consists in having $i = \pi/2$. We set $i$ to this value, keeping in mind that different inclinations will always be expected to give worse results. Analogously we set $\phi_0 = \pi/2$. The results might change if other values of $\phi_0$ are chosen, but at least for periods shorter than the LISA mission duration this choice should not affect them. The remaining parameters, namely $f_0, P$ and $K$, are varied in order to explore the parameter space. Since we analyse only systems for which the planetary inclination $i$ and DWD total mass $M_c$ (through $M_c = 0.2 M_\odot$ and $\eta = (M_c/M_\odot)^{5/3} = 1/4$) are fixed, the choice of the planetary mass directly determines the value of $K$. For this reason all results are reported by referring to the planetary mass $M_p$. We provide results for a LISA nominal mission duration of 4 years (ref. 3) and neglect the presence of a possible second planet or tertiary star, which are commonly found around binaries orbiting DWDs$^{52}$.

We, moreover, compared the errors on the DWD parameters with the expected results recovered without the presence of a CBP. We verified that the error estimates for $\Psi_0, f_0$ and $f_1$ deteriorate at maximum by a factor of two for CBPs close to $13 M_p$, while the other DWD parameters remain unaffected. Less massive CBPs will have a lower effect, while more massive objects can even bias the measured value of DWD parameters$^{11}$. Finally, we verified that the LISA accuracy on K and P parameters improves only marginally by fixing the sky location, as in the case where an EM counterpart is spotted.

**Details of Fig. 2.** The precision on the CBP orbital parameters, as estimated with the methodologies described above, is reported in Fig. 2.

First of all we have verified that as long as our assumptions hold, errors on the planetary parameters change linearly as $M_p$ changes. The CBP mass can thus be factored out from the final estimation of the measurement precision on the parameter, similarly to what we have done with the SNR. This is why in Fig. 2 we present results scaled by both the DWD SNR and the CBP mass, and why in Fig. 3 the LISA selection function for SNR=100 and $f_0 = 1$ mHz coincides with the one for SNR=10 and $f_0 = 10$ mHz.

The behaviour of the uncertainty on P in Fig. 2 differs from the behaviour of K and $\phi_0$. An explanation for this can be found by looking at the form of the signal in the Fourier domain (see Fig. 5 of ref. 11). In fact, the spread of the signal over different frequency bins due to the Doppler modulation induced by the CBP directly determines the orbital period of the planet. This spread is easier to measure with respect to other details of the GW signal that are used to determine K and $\phi_0$, together with the other DWD parameters. This implies that the CBP period P will always be well measured as long as it is shorter than the LISA observational time. For longer CBP periods this spread is less effective since only a fraction of the orbit is observed and the signal has no time to be fully spread, which, in turn, implies that P is recovered with less precision. Of course the longer

the CPB period, the lower the spread is and thus the measurement of P is less precise. Note that on top of the frequency spread given by the Doppler modulation due to the CBP, there is also the usual spread given by the orbital motion of LISA around the Sun. It is, however, easy to distinguish the two spreads since we know the orbital period of LISA is exactly 1 year. Similarly, the peak at 1 year on K and $\phi_0$ is due to the degeneracy between the motion of the DWD around the three-body CoM and the motion of LISA around the Sun. We have indeed checked that other smaller peaks appear at multiples of 1 year, confirming that this effect is due to the degeneracy between the two orbital motions.

**Synergy with radial velocity.** The radial velocity technique allows for the determination of the planetary orbital period, P, the semi-amplitude of the radial velocity curve, K, and the eccentricity, e (which we did not take into account in this work), the longitude of the periastron, $\omega$, and the time of periastron passage $T_0$. The drawback of this technique is that the radial velocity signal can be strongly biased by the presence of stellar activity, causing a false-positive detection of a planetary companion.

Note that radial velocity also allows for the determination of the stellar mass. In the case of spectroscopic binaries (and if both stars are bright enough), each resolved stellar spectrum enables the estimation of the effective temperature and surface gravity, and consequently the determination of the stellar mass by using the white dwarf mass–radius relations$^{53,54}$. In the case of eclipsing spectroscopic binaries it is possible to solve for the mass of each member of the binary straight away.

The Very Large Telescope (VLT)/Echelle SPectrograph for Rocky Exoplanets and Stable Spectroscopic Observations (ESPRESSO) is designed to explore a new mass domain, corresponding to rocky planets down to the Earth mass in the habitable zone of solar-type stars, with a radial velocity precision down to the 10 cm s⁻¹ level$^{55}$. Its selection curve is shown in Fig. 3. Such an instrument will be powerful enough to be able to follow up all possible potential exoplanets discovered by LISA, at least within its horizon (V magnitudes as faint as 20 to 21 in dark-sky conditions, https://www.eso.org/sci/facilities/paranal/instruments/ espresso/ESPRESSO_User_Manual_P102.pdf).

For exoplanet detection, radial velocities and GWs are conceptually similar methods, consequently the planetary parameters retrieved are the same. However, GWs are not affected by stellar activity and can observe at larger distances. Radial velocity is thus not expected to improve the CBP characterization already provided by GWs, but it will be extremely useful to measure the stellar masses and to detect the EM counterpart, which would consequently confirm potential GW detections.

**Synergy with transit.** The transit method$^{56}$ alone allows for the determination of the planetary radius, $R_p$, the orbital inclination, $i$, and the semi-major axis, a. Transit surveys, though, are affected by the issue of astrophysical false-positives, which requires complementary observations to constrain the presence of a planet around a single star system. Conversely, a planetary detection in a circumbinary system is unambiguous due to the geometry (Fig. 1) and dynamics of the system itself. When the binary and the planetary orbits are coplanar$^{57,58}$, the transits are only possible on eclipsing binaries, and when there is a misalignment, transits are still possible, but with gaps in the transit sequence and asymmetries in the transit profile$^{51,59}$. Although, this might not be the case for a very short binary period such as as the one we are considering in this work. Furthermore, because the planet is transiting a moving target, both enhanced transit timing$^{9,60}$ and transit duration variations$^{61,62}$ are expected$^{58}$.

Transit measurements, like radial velocities, are contaminated by stellar activity. However, a study on Kepler data showed that white dwarfs are photometrically stable to better than 1% on 1-h to 10-day timescales$^{63}$, a key aspect for this kind of observations.

We plot the selection function of the ESA's Atmospheric Remote-Sensing Infrared Exoplanet Large-survey (ARIEL; launch in 2028) mission$^{64}$, based on its Tier-1 sample. The goal of ARIEL is to investigate the atmospheres of several hundreds of planets through transit and eclipse spectroscopy.

For the transiting hot-Jupiter whose mass and semi-major axis fall inside the LISA selection function (Fig. 3), it will be possible to determine $M_p, i$ and P (together with the eccentricity not considered here).

**Synergy with astrometry.** Notably, astrometry allows to determine all the orbital elements of a detected planet, and yields a direct measurement of $M_p/M_\odot$ from which it is possible to determine the planetary mass if the binary mass is known. Astrometry is sensitive to faint targets, long-period planets and it is not affected by stellar activity. On the other hand, this technique is limited by the time span of the survey (which puts limits on the period of exoplanets that can be detected) and, more importantly, by the distance of the system, which directly affects the precision on the individual measurements and hence the ability to resolve the system itself. Generally, at a fixed distance and stellar mass, the farther and the more massive the planet, the higher the astrometric signature.

No planets have been detected yet through this method, though Gaia, at the end of the nominal mission, will deliver a catalogue of 21,000 (±6,000) high-mass (1–15 $M_\odot$) long-period planets at up to distances of 500 pc (ref. 65) and 500 circumbinary gas giants within 200 pc (ref. 66). Gaia is also expected to detect tens
or hundreds of planets around single white dwarf $M_p > ∼1 M_p$) in long-period orbits. We report in Fig. 3 the typical Gaia selection function⁶⁷ updated with the most recent single-transit astrometric uncertainty model⁶⁸.

In the case of joint GW-astrometric observations, we should be able to determine $M_p$ from the combined measurements of $P$, $i$ and the ratio $M_p/M_b$.

**Synergy with direct imaging.** Direct imaging of extrasolar planets is currently extremely challenging, especially at close projected separations, due to the high planet-to-star contrast. Consequently, the method is currently sensitive to young giant self-luminous planets, orbiting at wide separations around low-luminous stars.

With only one epoch it is possible to measure the planetary angular separation, with more epochs over much longer timescales we can directly retrieve, $i$, $a$ and $P$ (ref. ⁶⁹). Other parameters such as $M_p$, the actual effective temperature and gravity are derived from the observed photometry using age-dependent relationships, meaning that an error on the age of the system translates in a large error on $M_p$ and $R_p$.

The advantage of imaging white dwarfs is that such objects are natural clocks and their age, well constrained, can be derived straightforwardly from their luminosity⁷⁰. Obviously this would imply to be able to resolve both components in the binary system. Furthermore, white dwarfs are ~10⁴ times fainter than their progenitors, and the contrast ratio between DWDs-CBP is lower after each common-envelope phase.

No imaging of any post-common-envelope binaries has been done yet, although giant planets have been found around a couple of these kinds of system with eclipse-time variations⁷¹,⁷². These planets have large separations and can be imaged with upcoming large ground-based telescopes such as the ELT. The ELT will also allow us to perform spectroscopic observations of an unresolved planet and its exomoon(s), (if existing), and resolve the emitted combined light when the moon is hidden by the planet (that is, during the moon eclipse).

We plot in Fig. 3 the selection curves for ELT-PCS⁶⁹ and the Spectro-Polarimetric High-contrast Exoplanet REsearch instrument (SPHERE)⁶⁸ for comparison. While with SPHERE the synergy with LISA is inexistent (at the planetary mass level), the synergy with PCS covers mostly planets within $a=1–2$au, whose $M_p$ varies in the 1–13 (and above) $M_i$ range. If both the DWD total mass and the inclination $i$ are measured through direct imaging, then a joint observation with LISA will yield the planetary mass $M_p$.

One of the advantage of this technique is that it is not affected by the variability of the central star, it is though limited by the distance of the CBP system and by the angular resolution of the instruments.

**Synergy with gravitational microlensing.** Gravitational microlensing alone provides the star-to-planet mass ($M_s/M_p$) ratio and projected planetary separation in units of Einstein ring. However, it also allows to retrieve the distance $d$ of the system, $M_s$, $M_p$ and the semi-major axis $a$ in physical units, when working in synergy with adaptive optics observations and/or variations of the alignment between the foreground-background stars, due to a parallax effect. Microlensing detections are possible even for small planetary masses, and are more suitable for targets located between the centre of the Galaxy and the observer. The downside is that every event is rare and unique, hence not reproducible.

Currently, only one circumbinary system has been discovered using gravitational microlensing⁷³. NASA's Wide Field Infrared Survey Telescope (WFIRST) mission (scheduled for launch in 2028) has both microlensing and direct-imaging capabilities. It is expected that a total of ~1,400 bound (to single stars) exoplanets with mass greater than ~0.1 $M_⊕$, including ~200 with mass $≤3 M_⊕$ (ref. ⁷⁴) will be found through microlensing, and the technique will also vastly increase the number of CBPs detected⁷⁵. We plot in Fig. 3 the selection curve for WFIRST.

Note that there will be strong complementarity with LISA for CBPs with a separation up to ~2au. A joint detection with GWs would be useful to confirm the presence of a CBP by strengthening the lower limit for the planetary mass.

## Data availability

The data that support the plots within this paper and other findings of this study are available from the authors upon reasonable request.

Received: 28 November 2018; Accepted: 7 May 2019;

Published online: 8 July 2019

## References

1. Winn, J. N. in *Handbook of Exoplanets* (eds Deeg, H. J. & Belmonte, J. A.) 1949–1966 (Springer, 2018).

2. Fulton, B. J. et al. The California-Kepler Survey III. A gap in the radius distribution of small planets. *Astron. J.* **154**, 109 (2017).

3. Amaro-Soeane, P. et al. Laser Interferometer Space Antenna. Preprint at https://arxiv.org/abs/1702.00786 (2017).

4. Raghavan, D. et al. A survey of stellar families: multiplicity of solar-type stars. *Astrophys. J. Suppl.* **190**, 1–42 (2010).

5. Duchêne, G. & Kraus, A. Stellar multiplicity. *Annu. Rev. Astron. Astrophys.* **51**, 269–310 (2013).

6. Althaus, L. G., Córsico, A. H., Isern, J. & García-Berro, E. Evolutionary and pulsational properties of white dwarf stars. *Astron. Astrophys. Rev.* **18**, 471–566 (2010).

7. Veras, D. Post-main-sequence planetary system evolution. *R. Soc. Open Sci.* **3**, 150571 (2016).

8. Sigurdsson, S., Richer, H. B., Hansen, B. M., Stairs, I. H. & Thorsett, S. E. A young white dwarf companion to pulsar B1620-26: evidence for early planet formation. *Science* **301**, 193–196 (2003).

9. Korol, V. et al. Prospects for detection of detached double white dwarf binaries with Gaia, LSST and LISA. *Mon. Not. R. Astron. Soc.* **470**, 1894–1910 (2017).

10. Korol, V., Koop, O. & Rossi, E. M. Detectability of double white dwarfs in the local group with LISA. *Astrophys. J.* **866**, L20 (2018).

11. Robson, T., Cornish, N. J., Tamanini, N. & Toonen, S. Detecting hierarchical stellar systems with LISA. *Phys. Rev. D* **98**, 064012 (2018).

12. Kostov, V. B., Moore, K., Tamayo, D., Jayawardhana, R. & Rinehart, S. A. Tatooine's future: the eccentric response of Kepler's circumbinary planets to common-envelope evolution of their host stars. *Astrophys. J.* **832**, 183 (2016).

13. Zuckerman, B., Melis, C., Klein, B., Koester, D. & Jura, M. Ancient planetary systems are orbiting a large fraction of white dwarf stars. *Astrophys. J.* **722**, 725–736 (2010).

14. Veras, D. & Tout, C. A. The great escape – II. Exoplanet ejection from dying multiple-star systems. *Mon. Not. R. Astron. Soc.* **422**, 1648–1664 (2012).

15. Schleicher, D. R. G. & Dreizler, S. Planet formation from the ejecta of common envelopes. *Astron. Astrophys.* **563**, A61 (2014).

16. Ferrari, V., Berti, E., D'Andrea, M. & Ashtekar, A. Gravitational waves emitted by extrasolar planetary systems. *Int. J. Mod. Phys. D* **9**, 495–509 (2000).

17. Berti, E. & Ferrari, V. Excitation of g-modes of solar-type stars by an orbiting companion. *Phys. Rev. D* **63**, 064031 (2001).

18. Ain, A., Kastha, S. & Mitra, S. Stochastic gravitational wave background from exoplanets. *Phys. Rev. D* **91**, 124023 (2015).

19. Cunha, J. V., Silva, F. E. & Lima, J. A. S. Gravitational waves from ultra-short period exoplanets. *Mon. Not. R. Astron. Soc.* **480**, L28 (2018).

20. Wong, K. W. K., Berti, E., Gabella, W. E. & Holley-Bockelmann, K. On the possibility of detecting ultra-short period exoplanets with LISA. *Mon. Not. R. Astron. Soc.* **483**, L33–L36 (2019).

21. Takahashi, R. & Seto, N. Parameter estimation for galactic binaries by LISA. *Astrophys. J.* **575**, 1030–1036 (2002).

22. Thompson, T. A. Accelerating compact object mergers in triple systems with the Kozai resonance: a mechanism for "prompt" Type Ia supernovae, gamma-ray bursts, and other exotica. *Astrophys. J.* **741**, 82 (2011).

23. Seto, N. Highly eccentric Kozai mechanism and gravitational-wave observation for neutron star binaries. *Phys. Rev. Lett.* **111**, 061106 (2013).

24. Valsechi, F., Farr, W. M., Willems, B., Deloye, C. J. & Kalogera, V. Tidally-induced apsidal precession in double white dwarfs: a new mass measurement tool with LISA. *Astrophys. J.* **745**, 137 (2012).

25. Kremer, K., Breivik, K., Larson, S. L. & Kalogera, V. Accreting double white dwarf binaries: implications for LISA. *Astrophys. J.* **846**, 95 (2017).

26. Breivik, K. et al. Characterizing accreting double white dwarf binaries with the Laser Interferometer Space Antenna and Gaia. *Astrophys. J.* **854**, L1 (2018).

27. Nelemans, G., Yungelson, L. R. & Portegies Zwart, S. F. Short-period AM CVn systems as optical, X-ray and gravitational wave sources. *Mon. Not. R. Astron. Soc.* **349**, 181–192 (2004).

28. Cutler, C. Angular resolution of the LISA gravitational wave detector. *Phys. Rev. D* **57**, 7089–7102 (1998).

29. Sumi, T. et al. Unbound or distant planetary mass population detected by gravitational microlensing. *Nature* **473**, 349–352 (2011).

30. Dai, X. & Guerras, E. Probing extragalactic planets using quasar microlensing. *Astrophys. J. Lett.* **853**, L27 (2018).

31. Korol, V., Rossi, E. M. & Barausse, E. A multimessenger study of the Milky Way's stellar disc and bulge with LISA, Gaia, and LSST. *Mon. Not. R. Astron. Soc.* **483**, 5518–5533 (2019).

32. Luyten, W. J. White Dwarfs (Minneapolis, University of Minnesota, 1970).

33. Agol, E. Transit surveys for Earths in the habitable zones of white dwarfs. *Astrophys. J. Lett.* **731**, L31 (2011).

34. Spalding, C., Batygin, K. & Adams, F. C. Resonant removal of exomoons during planetary migration. *Astrophys. J.* **817**, 18 (2016).

35. Kostov, V. B. et al. Kepler-1647b: the largest and longest-period Kepler transiting circumbinary planet. *Astrophys. J.* **827**, 86 (2016).

36. Ivanova, N. et al. Common envelope evolution: where we stand and how we can move forward. *Astron. Astrophys. Rev.* **21**, 59 (2013).

37. Turrini, D., Nelson, R. P. & Barbieri, M. The role of planetary formation and evolution in shaping the composition of exoplanetary atmospheres. *Exp. Astron.* **40**, 501–522 (2015).
38. Alexander, M. E., Chau, W. Y. & Henriksen, R. N. Orbital evolution of a singly condensed, close binary, by mass loss from the primary and by accretion drag on the condensed member. *Astrophys. J.* **204**, 879–888 (1976).

39. Holman, M. J. & Wiegert, P. A. Long-term stability of planets in binary systems. *Astron. J.* **117**, 621–628 (1999).

40. Pilat-Lohinger, E., Funk, B. & Dvorak, R. Stability limits in double stars. A study of inclined planetary orbits. *Astron. Astrophys.* **400**, 1085–1094 (2003).

41. Debes, J. H. & Sigurdsson, S. Are there unstable planetary systems around white dwarfs? *Astrophys. J.* **572**, 556–565 (2002).

42. Livio, M., Pringle, J. E. & Wood, K. Disks and planets around massive white dwarfs. *Astrophys. J. Lett.* **632**, L37–L39 (2005).

43. Faedi, F., West, R. G., Burleigh, M. R., Goad, M. R. & Hebb, L. Detection limits for close eclipsing and transiting substellar and planetary companions to white dwarfs in the WASP survey. *Mon. Not. R. Astron. Soc.* **410**, 899–911 (2011).

44. Quarles, B., Satyal, S., Kostov, V., Kaib, N. & Haghighipour, N. Stability limits of circumbinary planets: is there a pile-up in the Kepler CBPs? *Astrophys. J.* **856**, 150 (2018).

45. Mustill, A. J. et al. Main-sequence progenitor configurations of the NN Ser candidate circumbinary planetary system are dynamically unstable. *Mon. Not. R. Astron. Soc.* **436**, 2515–2521 (2013).

46. Portegies Zwart, S. Planet-mediated precision reconstruction of the evolution of the cataclysmic variable HU Aquarii. *Mon. Not. R. Astron. Soc.* **429**, L45–L49 (2013).

47. Armano, M. et al. Sub-femto-g free fall for space-based gravitational wave observatories: LISA Pathfinder results. *Phys. Rev. Lett.* **116**, 231101 (2016).

48. Armano, M. et al. Beyond the required LISA free-fall performance: new LISA Pathfinder results down to 20µHz. *Phys. Rev. Lett.* **120**, 061101 (2018).

49. Adams, D. *The Hitchhiker's Guide to the Galaxy* (Pan Books, 1979).

50. Cornish, N. J. & Larson, S. L. LISA data analysis: source identification and subtraction. *Phys. Rev. D* **67**, 103001 (2003).

51. Tokovinin, A., Thomas, S., Sterzik, M. & Udry, S. Tertiary companions to close spectroscopic binaries. *Astron. Astrophys.* **450**, 681–693 (2006).

52. Pepe, F. et al. ESPRESSO: the next European exoplanet hunter. *Astron. Nachr.* **335**, 8–20 (2014).

53. Tinetti, G. et al. A chemical survey of exoplanets with ARIEL. *Exp. Astron.* **46**, 135–209 (2018).

54. Penny, M. T. et al. Predictions of the *WFIRST* microlensing survey. I. Bound planet detection rates. *Astrophys. J. Suppl. Ser.* **241**, 3 (2019).

55. Lagrange, A.-M. Direct imaging of exoplanets. *Phil. Trans. R. Soc. Lond. Ser. A* **372**, 20130090 (2014).

56. Beichman, C. A. et al. Imaging young giant planets from ground and space. *Publ. Astron. Soc. Pac.* **122**, 162 (2010).

57. Perryman, M., Hartman, J., Bakos, G. Á. & Lindegren, L. Astrometric exoplanet detection with Gaia. *Astrophys. J.* **797**, 14 (2014).

58. Althaus, L. G., García-Berro, E., Isern, J. & Córsico, A. H. Mass-radius relations for massive white dwarf stars. *Astron. Astrophys.* **441**, 689–694 (2005).

59. Renedo, I. et al. New cooling sequences for old white dwarfs. *Astrophys. J.* **717**, 183–195 (2010).

60. Deeg, H. J. & Alonso, R. in *Handbook of Exoplanets* (eds Deeg, H. J. & Belmonte, J. A.) 633–657 (Springer, 2018).

61. Pierens, A. & Nelson, R. P. Orbital alignment of circumbinary planets that form in misaligned circumbinary discs: the case of Kepler-413b. *Mon. Not. R. Astron. Soc.* **477**, 2547–2559 (2018).

62. Foucart, F. & Lai, D. Assembly of protoplanetary disks and inclinations of circumbinary planets. *Astrophys. J.* **764**, 106 (2013).

63. Martin, D. V. & Triaud, A. H. M. J. Planets transiting non-eclipsing binaries. *Astron. Astrophys.* **570**, A91 (2014).

64. Martin, D. V. Circumbinary planets – II. When transits come and go. *Mon. Not. R. Astron. Soc.* **465**, 3235–3253 (2017).

65. Holman, M. J. & Murray, N. W. The use of transit timing to detect terrestrial-mass extrasolar planets. *Science* **307**, 1288–1291 (2005).

66. Armstrong, D. et al. Placing limits on the transit timing variations of circumbinary exoplanets. *Mon. Not. R. Astron. Soc.* **434**, 3047–3054 (2013).

67. Kostov, V. B. et al. Kepler-413b: a slightly misaligned, Neptune-size transiting circumbinary planet. *Astrophys. J.* **784**, 14 (2014).

68. Liu, H.-G., Wang, Y., Zhang, H. & Zhou, J.-L. Transits of planets with small intervals in circumbinary systems. *Astrophys. J.* **790**, 141 (2014).

69. Hermes, J. J. et al. When flux standards go wild: white dwarfs in the age of Kepler. *Mon. Not. R. Astron. Soc.* **468**, 1946–1952 (2017).

70. Sahlmann, J., Triaud, A. H. M. J. & Martin, D. V. Gaia's potential for the discovery of circumbinary planets. *Mon. Not. R. Astron. Soc.* **447**, 287–297 (2015).

71. Winget, D. E. & Kepler, S. O. Pulsating white dwarf stars and precision asteroseismology. *Annu. Rev. Astron. Astrophys.* **46**, 157–199 (2008).

72. Qian, S.-B. et al. A circumbinary planet in orbit around the short-period white dwarf eclipsing binary RR Cae. *Mon. Not. R. Astron. Soc.* **422**, L24–L27 (2012).

73. Beuermann, K., Dreizler, S. & Hessman, F. V. The quest for companions to post-common envelope binaries. IV. The 2:1 mean-motion resonance of the planets orbiting NN Serpentis. *Astron. Astrophys.* **555**, A133 (2013).

74. Bennett, D. P. et al. The first circumbinary planet found by microlensing: OGLE-2007-BLG-349L(AB)c. *Astron. J.* **152**, 125 (2016).

75. Luhn, J. K., Penny, M. T. & Gaudi, B. S. Caustic structures and detectability of circumbinary planets in microlensing. *Astrophys. J.* **827**, 61 (2016).

76. Casertano, S. et al. Double-blind test program for astrometric planet detection with Gaia. *Astron. Astrophys.* **482**, 699–729 (2008).

77. Sozzetti, A. et al. Astrometric detection of giant planets around nearby M dwarfs: the Gaia potential. *Mon. Not. R. Astron. Soc.* **437**, 497–509 (2014).

## Acknowledgements

We thank E. Berti, A. Buonanno, V. Korol, P.-O. Lagage, C. Miller, A. Petiteau, E. M. Rossi and G. Tinetti for their suggestions and comments. We are particularly thankful to V. Korol for providing the DWD data beyond Fig. 3 of ref. 31. C.D. acknowledges support from the LabEx P2IO, the French ANR contract 05-BLAN-NT09-573739. We acknowledge the use of the Python package mw plot (https://pypi.org/project/mw-plot/). This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.

## Author contributions

N.T. performed the theoretical and numerical studies to produce and analyse the results presented in the paper. C.D. assessed the feasibility within the exoplanetary context, and investigated the synergies with EM observations. Both authors interpreted the results, studied the implications and wrote the paper.

## Competing interests

The authors declare no competing interests.

## Additional information

Reprints and permissions information is available at www.nature.com/reprints.

Correspondence and requests for materials should be addressed to N.T. or C.D.

Peer review information: Nature Astronomy thanks Jonathan Gair and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© The Author(s), under exclusive licence to Springer Nature Limited 2019