{'fraction_non_alphanumeric': 0.02652295299556044, 'fraction_numerical': 0.012082932857339008, 'mean_word_length': 4.449432597580746, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Listening to the gravitational wave sound of circumbinary exoplanets

To date more than 3500 exoplanets have been discovered orbiting a large variety of stars. Due 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 new detection method that overcomes this issue and that will allow us to detect gas giant exoplanets using gravitational wave astronomy. We show that the Laser Interferometer Space Antenna (LISA) mission can characterise hundreds of 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 super-Earths down to 10 Earth masses.


Exoplanets beyond the solar neighbourhood
In the last twenty 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 1 , and that there is a gap in the radius distribution of small planets 2 . 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 it is given by a selection bias or not, it is extremely important, and it cannot be assessed through the usual techniques for at least the next two decades.
This breakthrough is going to be possible only through gravitational wave (GW) astronomy. The Laser Interferometer Space Antenna (LISA) mission 3 , planned for launch in the early 2030s, will enable us to indirectly probe for the first time the population of gas giants orbiting detached double white dwarfs (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 4,5 , and that ∼ 95% of the stars will become a white dwarf 6 , 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 6 systems have a white dwarf as binary component. Usually the second companion is a low-mass star 7 or, for one specific case, a pulsar 8 . Today, no exoplanets have been discovered around double white dwarfs, irrespectively of the compactness of the binary.
Due to the intrinsic faintness of these DWDs, only few tens are known by spectroscopic and variability surveys 9 , 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 25x10 3 compact DWDs within and outside the Milky Way 9, 10 , some of which could be perturbed by the presence of a third gravitationally bound stellar companion 11 or planetary companion.
The existence of this planetary population in our Galaxy is far from being excluded 7,12 , 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 (i) to understand under which conditions  Relative accuracy on K and P, and absolute accuracy on ϕ 0 , for a LISA measurements of planetary orbital parameters. Solid curves denote DWD with f 0 = 10 mHz, while dashed curves are for f 0 = 1 mHz. All results have been scaled with respect to the DWD SNR and the planetary mass measured in Jupiter masses. The peak at one 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. a planet can survive the most critical phases of a close binary evolution 12,13 ; (ii) to set constraints on binary mass-loss and the dynamical aspects which directly follow 12,14 ; (iii) to check whether a second-generation of planets exists, i.e., planets that form from the stellar material ejected during the binary common envelope phase(s) 15 . 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 12 .
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 GWs sources to be heard by LISA, which makes it more powerful and interesting than former ideas of direct detection of exoplanets through GWs [16][17][18][19][20] . 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 one hour emit almost monochromatic GWs in the LISA frequency band. If a third companion 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 waveform measured by LISA 11 . We refer to the Methods section for the dynamical model employed to describe the three-body system, 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 ϕ 0 and the parameter which depends on the CBP mass M p , its orbit inclination i, and the binary total mass M b . In Fig. 1 we show how the 1σ relative 2/15 uncertainties on K, P, and the 1σ absolute uncertainties on ϕ 0 (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 0 = 1 mHz (dashed lines), at which the majority of DWDs detected by LISA are expected 9 , and a higher frequency f 0 = 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 in order to find the precision with which the parameters are measured, the numbers reported in Fig. 1 must be divided by the SNR of each individual DWD event and by the CBP mass, measured in M J (Jupiter's mass). Fig. 1 shows that the error estimations on K and ϕ 0 are better for planetary periods P comparable to the LISA nominal lifetime. For shorter planetary periods the errors on K and ϕ 0 smoothly decrease with the increasing of P, while for planetary periods longer than LISA's nominal duration the precision steeply worsen by few order of magnitude. On the other hand, the uncertainty on P increases smoothly with the planetary period till the LISA nominal mission lifetime, after which it rapidly worsen similarly to the behaviour of K and ϕ 0 . An explanation for the different behaviour of P, as well as for the appearance of the peak at one year, is provided in the Methods section.
Analogously to EM radial velocity techniques, we cannot recover the mass of the companion directly from the three GW parameters K, P, ϕ 0 , implying that we have no mean to know 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 orbit inclination, can help determining its mass through Eq. (1). It is furthermore possible to obtain some estimates on the allowed M p range if also the chirp mass of the binary is measured with LISA, which is usually the case at least for high frequency DWDs 21 . 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 b , and from this range of values we can estimate M p sin(i) through Eq. (1). Since from GWs alone we do not obtain any information on i, we can only derive lower limits for M p by considering sin(i) = 1. Values of sin(i) different from 1 will yield higher values of M p . If the lower mass limit estimated in this way is below 13 M J (the deuterium burning limit), then the companion object is probably an exoplanet which would need to be confirmed by EM follow-ups (cf. Sec. 3). On the other hand, if the lower bound of this estimate is above 13 M J , then we can confidently exclude the possibility that the third orbiting object is a planet. Furthermore, if we also set a minimum allowed value for sin(i), we can find an upper bound on M p which would confirm the presence of an exoplanet from the GW signal alone if it does not exceed 13 M J . 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 which 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. Fig. 2 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. From the figure we see that LISA will be more effective at observing exoplanets with a separation from the DWD between roughly 1 and 2 AU. Depending on the SNR and frequency of the DWD, LISA will be able to detect exoplanets down to ∼1 M J and will be efficient in the range between 0.01 and 3 AU, roughly.
Although the occurrence rate of CBPs around DWDs is presently unknown, if these planets physically exist, and they orbit even a few percent 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 by two WDs, each of mass M * = 0.23M (chirp mass M c = 0.2M ). 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 η and GW frequency. The same reasoning applies to different types of stellar binaries visible with LISA, such as neutron star-neutron star and WD-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 29,30 , it might also be possible to measure the individual masses of the stars directly from GW observations 31 , yielding a more accurate and precise estimation of the CBP mass without the need of complementary EM observations.
We moreover assume the DWDs to be detached, namely with no accretion. LISA is expected to detect also few thousands accreting DWDs 32 , 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 characterised 33 . Moreover accreting DWDs are expected to emit specific EM signatures which might facilitates the identification of an EM counterpart 34 , enhancing the multi-messenger potential of these systems.
In our analysis we assume a 4 years 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,35 . This translates into both a higher number of CBPs detections with GWs, and a more precise characterisation of the systems already observed in 4 years. A 10 years LISA mission will also be more sensitive to CBPs with periods between 4 and 10 years. In fact we expect the rapid deterioration present at periods longer than 4 years of the accuracy with which the planetary parameters are recovered (cf. Fig. 1), to shift toward 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 10 (Andromeda galaxy). Although only systems in the solar neighbourhood (d * 3 kpc) can be repeatedly observed with joint GW and EM observations, farther CBPs can be detected by GW alone. Given that no EM detection of extragalactic planets bound to their star(s) has been confirmed yet (see 36,37 for detection of population of unbound planets), LISA will open a new window for the quest of extra-galactic 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 ≥ 1 M J . Some of these higher frequency DWDs will be measured with SNR ∼ 1000, meaning that, for these specific systems, LISA will be able to detect CBPs down to 30 M ⊕ (∼ 0.1 M J ). 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 few Earth masses. These CBPs might as well be located in the habitable zone (HZ) of the DWDs (a ∼ 0.014 AU, to account for the double luminosity) 38 .

Synergy with EM observations and implications
Observing the EM counterpart is important to fully characterise the planetary system. When a planet is first heard by GWs (cf Sec. 2), we can precisely determine M p only if the total DWD mass is measured with EM observations. Moreover when the 4/15 mass of each component in the binary is determined, then it is possible to estimate their radius using mass-radius relations 7 . Recent estimates 9 suggest that out of the 25 × 10 3 DWDs expected to be detected in 4 years by LISA, only up to a maximum of 100 will be observed with instruments such as Gaia and LSST, although upcoming large ground-based facilities might well improve these numbers. We verified also that once the EM counterpart has been spotted, the LISA accuracy on K and P parameters improves only marginally by fixing the sky location.
On the other hand, the CBP itself has its typical signature that can be spotted in the EM by various well-known detection techniques (cf. Fig. 2), which on average mainly observe systems relatively close to us. Transit, radial velocity (RV), direct imaging, and astrometry methods cover distances up to d * ∼ 3 kpc, while gravitational microlensing goes up to d * ∼ 8 kpc. Among these, while the transit and RV 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 timespan of the survey. The CBP detection method by GWs proposed here is conceptually similar to RV, but has the advantages that it can observe everywhere in the galaxy, it is not affected by the activity of the stars, and it is does not need any observational pointing. If we assume that the planetary system complies the requirements specific for the EM detection (e.g. magnitude and system configuration), the synergy with GWs allows us to resolve the circumbinary system in term of periods, masses (of the stars + planet), and inclination of the orbits. We refer to the Methods for more details on the parameters that can be retrieved in the specific synergy cases.
We note that the orbit of a CBP is stable when the semi-major axis a is approximately a 4.5 the binary period P b 39, 40 , which in our specific case corresponds to a minimal separation of a min ∼ 0.0009 AU and a min ∼ 0.004 AU for the binary period P b = 30 min and P b = 1 hour, respectively. These distances are well below the minimal separation considered for CBPs in this work, since we always assume P P b in order for higher order orbital effects, e.g. Kozai-Lidov resonances, to be negligible.
WDs are common as Sun-like stars and may also provide a source of energy for planets for gigayear (Gyr) durations. Cool white dwarfs (T eff <6000 K), have an HZ that endures for up to 8 Gyr 38 located at only few solar radii from the star. This HZ regularly moves inwards with the star aging, and hence cooling 38 . Planets in the HZ of a WD must have consequently migrated inwards after the stellar giant phase [41][42][43] . The detection of such objects, or any close-in object not necessarily in the HZ, would help constraining migration theories of planets around post-common envelope binaries. For a binary case and with hotter DWDs (both at T eff >6000 K), the HZ scales outwards to account for the energy contribution of both stars and their temperatures.
It has been shown that, after one inner binary common envelope (CE) phase, CBPs can either be ejected by the system 7, 14 , or survive the CE phase(s) while undergoing an orbital expansion/shrinkage, generally coupled with an increased eccentricity 12 . The kind of binary evolution, i.e., the tidal evolution, the energy transfer between the CE and the inner cores, and the mass-loss rate, is directly responsible for the fate of their planetary companion. Consequently, the detection and characterisation of CBPs around DWDs is crucial for enabling a comparison with known existing circumbinary systems, and for pinning down the planetary evolution phases, valid all over the Galaxy, while their binary evolves. An example is the circumbinary system Kepler 1647 44 (gas giant orbiting two solar-mass stars of F and G type), which could, under a specific theoretical scenario, become a WD-WD pair that keeps a bound exoplanet after two CE phases 12 .
As previously mentioned, a 10 years mission lifetime can detect and characterise exoplanets on wider orbits, making LISA even more compatible with EM direct imaging techniques, and specifically with large ground-based cameras (e.g., ELT-PCS, cf. Fig. 2) or the successor of JWST. Imaging of CBPs around DWDs can be used to test the presence of a second-generation exoplanets in the outer regions of planetary system, and consequently to provide constraints on migration theories. Emission spectra of these objects will furthermore allow us to estimate its temperature and the main molecular component of its atmosphere, making direct connections to chemical element distributions in the WDs atmosphere. This would also permit to better understand the observed WD pollution effect 7 . On another hand, if an existing CBP accretes mass after a CE stage, it becomes brighter, further decreasing the already low planet-to-WDs contrast, meaning that also first-generation, more mature exoplanets, can be imaged.
Another peculiarity of these stars is their size. Their small radius, R WD ≈ 1 R ⊕ , makes them excellent targets for detecting transiting exoplanets, and in particular Earth-like planets 38 . However, Jupiter-like planets, which are a factor ∼ 7 larger than a typical WD, 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 have an exomoon (assuming it is not removed through an evection resonance while undergoing the planetary migration phase 45 ), then transit spectroscopy will allow to study the atmosphere of the satellite, due to the small size of WDs. For instance, an Europa transit (R E ∼ 0.24 R ⊕ ) in front of a typical WD (M * = 0.6 M , R WD = 1.57 R ⊕ ) yields a decrease in brightness of approximately ∼ 5%.

5/15 4 Conclusion
We have presented an original observational method which employs gravitational waves 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 dwarfs 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 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 super-Earth exoplanets with mass down to 30 M ⊕ . 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 46 , 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, a further development in the planet-star interaction and binary evolution theoretical studies (e.g. 12,47,48 ) needs to be performed in the forthcoming years, in order to thoroughly assess the potential of LISA to detect exoplanets, especially before the definitive ESA approval of the mission design in the early 2020s 3 .
The characterisation of LISA circumbinary exoplanets, by both GW and EM observations, will deepen our knowledge on planet-star interaction during the critical phases of stellar evolution, and will enable us to find the connecting thread that runs from formation 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 49 .

Dynamical model and modified GW phase
The three-body systems we are interested in are composed by 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ẑ is perpendicular to the planetary orbital plane (cf. Fig. 3), the distance vector between the planet and the CoM of the DWD is given by where R and ϕ are with M b , M p , P and ϕ 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ẑ-direction points toward the source, we apply two rotations: a first rotation by i (the inclination angle betweenẑ andẑ -the line of sight) to bring theẑ direction to point toward theẑ direction, and a subsequent second rotation by α around the z axis (which after the first rotation coincides with the z axis) makingx andŷ to point in the same directions ofx andŷ (see Fig. 3). Note however that the second rotation aroundẑ is degenerate with a phase redefinition of the planetary orbit, i.e. it is degenerate with ϕ 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 where R is a rotation matrix by the angle i.
For what concerns our scopes, 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 We are only interested in the motion along the line of sight from the observer to the system 1 . Since the z-axis of the observer is aligned along the line of sight direction, the z-component of the motion, given by is the one we need to consider. The velocity of the DWD CoM along the line of sight is then given by where we defined the parameter In the reference frame comoving with its center of mass, the DWD emits almost monochromatic GWs at some specific frequency f GW . Since this frequency is changing on time scales much longer if compared to the observational time scale, an expansion around the initial observed frequency will suffice in describing their time evolution 21 . We can thus assume 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 which, as long as the dynamics is non-relativistic, gives Finally the phase at the observer of the GW can be obtained integrating the frequency f obs : where Ψ 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.

LISA parameter estimation
The three arms of LISA constitute a pair of two-arm detectors outputting two linearly independent signals h I,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 35 where In these expressions R Earth = 1 AU and P Earth = 1 year are the mean distance from the Sun and the orbital period of the Earth. F +,× I,II are the antenna pattern functions depending on the source angular position (θ S , φ S ), the orientation of its orbit (θ L , φ L ) and the LISA configuration. The quantities A +,× are instead constant amplitudes which depend on the physical parameters and orientation of the source. Eqs. (12)-(15) provide the GW signal measured by LISA irrespectively of the sky location of the source, i.e. for arbitrary values of (θ S , φ S ). The full expressions for all the quantities appearing in these equations can be found in the literature 21,35,50 . We basically follow the set up described in 35 , adding to the GW signal the perturbation due to the CBP.
In order to extract information on the parameters characterizing the three-body system under consideration, we employ matched filtering techniques 21,35 . We assume that the GW signal measured by LISA has the time dependent expression given in Eq. (12). This is characterized by 11 parameters, namely ln(A), Ψ 0 , f 0 , f 1 , θ S , φ S , θ L , φ L , K, P, ϕ 0 which we collectively denote by λ i . The SNR of the signal can be computed as where T obs is LISA observational time and S n ( f 0 ) is LISA one-sided spectral density noise computed at f 0 . In all computations we fix T obs = 4 years in agreement with the nominal mission requirements 3 . For all the systems we consider here, the SNR as defined in Eq. (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 Σ i j , given by the inverse of the Fisher information matrix The standard estimator for the statistical error on the parameter λ i is thus given by √ Σ ii . The Fisher information matrix itself can be computed as

8/15
Note that we are simplifying Eqs. (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 21,35 , in our analysis we always replace S n ( f 0 ) with the SNR of the source by using Eq. (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 ln(A), Ψ 0 , f 0 , f 1 , θ S , φ S , θ L , φ L , can be observed by LISA, has already been investigated in several works 21,35,50 . 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, ϕ 0 . We moreover restrict the sampled parameter space by considering only specific values for some of the DWD parameters. In analogy with the examples in 21, 35 , we set the orbital geometry of the DWD and its sky location by assuming 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 9 , we assume equal mass DWDs with a chirp mass of M c = 0.2 M . The value of M c , together with the frequency f 0 , yields f 1 as The value of the amplitude A is irrelevant for our analysis since we always report results in terms of SNRs. The most favorable orientation for the detectability of the planetary perturbation on the GW signal consists in having i = π/2. We set i to this value, keeping in mind that different inclinations will always be expected to give worse results. Analogously we set ϕ 0 = π/2.  Fig. 1 The precision on the CBP orbital parameters, as estimated with the methodologies exposed above, is reported in Fig. 1.

Details of
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. 1 we present results scaled by both the DWD SNR and the CBP mass, and why in Fig. 2 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. 1 differs from the behaviour of K and ϕ 0 . An explanation for this can be found by looking at the form of the signal in the Fourier domain (cf. Fig. 5 of 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 which are used to determine K and φ 0 , together with the other DWD parameters. This implies that the CBP period P will be always 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 is the spread and thus the less precise the measurement of P. 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 one year. Similarly, the peak at one year on K and ϕ 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 one year, confirming that this effect is due to the degeneracy between the two orbital motions.

Synergy with radial velocity
The radial velocity (RV) 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, ω, and the time of periastron passage T 0 . The drawback of such technique is that the RV signal can be strongly biased by the presence of stellar activity, causing a false positive detection of a planetary companion.
Note that RV 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 WD mass-radius relations 52,53 . In the case of eclipsing spectroscopic binaries it is possible to solve for the mass of each member of the binary straight away.
The 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 RV precision down to 10 cm s −1 level 24 . Its selection curve is shown in Fig. 2. 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) 2 .
For exoplanet detection RVs 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. RV 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 detect the EM counterpart, which would consequently confirm potential GW detections.

Synergy with transit
The transit method 54 alone allows for the determination of the planetary radius, R p , the orbital inclination, i, 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. 3) and dynamics of the system itself. When the binary and the planetary orbits are coplanar 55,56 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 57,58 . Furthermore, because the planet is transiting a moving target, both enhanced transit timing 59, 60 and transit duration variations 61, 62 are expected 58 .
Transit measurements, like RVs, are contaminated by stellar activity. However, a study on Kepler data showed that WDs are photometrically stable to better than 1% on 1-hr to 10-d timescales 63 , key aspect for this kind of observations.
We plot the selection function of the ESA-Atmospheric Remote-Sensing Infrared Exoplanet Large-survey (ARIEL; launch in 2028) mission 25 , based on its Tier-1 sample. ARIEL goal is to investigate the atmospheres of several hundreds planets through transit and eclipse spectroscopy.
For the transiting hot-Jupiter whose mass and semi-major axis fall inside the LISA selection function (Fig. 2) 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 b 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 timespan of the survey (which put limits on the period of exoplanets that can be detected) and, more importantly, by the distance of the system, which directly affect 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 21000 (±6000) high-mass (1 -15 M J ) long-period planets at up to distances of 500 pc 64 and 500 circumbinary gas giants within 200 pc 65 . Gaia is also expected to detect tens or hundreds of planets around single WD (M p >∼ 1 M J ) in long period orbits. We report in Fig. 2 the typical Gaia selection function 22 updated with a most recent single-transit astrometric uncertainty model 23 .
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 giants 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 27 . 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 WDs is that such objects are natural clocks and their age, well constrained, can be derived straightforward from their luminosity 66 . Obviously this would imply to be able to resolve both components in the binary system. Furthermore, WDs are ∼ 10 4 fainter than their progenitors, and the contrast ratio between DWDs-CBP is lower after each CE phase.
No imaging of any post-common envelope binaries has been done yet, though gas giants have been found around a couple of these kinds of systems with eclipse-time variations 67,68 . 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)-(whether existing), and resolve the emitted combined light when the moon is hidden by the planet (i.e., during the moon eclipse).
We plot in Fig. 2 the selection curve for ELT-PCS 27 and SPHERE 27 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 J 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 star-to-planet mass (M * /M p ) ratio and projected planetary separation in units of Einstein ring. However, it allows to retrieve also the distance d * of the system, M * , 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 doable 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 69 . The NASA-WFIRST mission (launch foreseen in 2028) has both microlensing and direct-imaging capabilities. Through microlensing it is expected to find a total of ∼ 1400 bound (to single stars) exoplanets with mass greater than ∼ 0.1 M ⊕ , including ∼200 with mass 3 M ⊕ 26 , and will vastly increase the number of CBPs detected 70 . We plot in Fig. 2 the selection curve for WFIRST. Note that there will be strong complementarity with LISA for circumbinary planets with a separation up to ∼ 2 AU. 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 statement
Data sharing not applicable to this article as no datasets were generated or analysed during the current study.   Figure 3. Geometry of the DWD-CBP three-body system. Left panel: source reference frame (primed coordinates). Right panel: observer reference frame (non-primed coordinates) with the source reference frame rotated.