Circumbinary exoplanets and brown dwarfs with the Laser Interferometer Space Antenna

We explore the prospects for the detection of giant circumbinary exoplanets and brown dwarfs (BDs) orbiting Galactic double white dwarfs binaries (DWDs) with the Laser Interferometer Space Antenna (LISA). By assuming an occurrence rate of 50%, motivated by white dwarf pollution observations, we built a Galactic synthetic population of P-type giant exoplanets and BDs orbiting DWDs. We carried this out by injecting different sub-stellar populations, with various mass and orbital separation characteristics, into the DWD population used in the LISA mission proposal. We then performed a Fisher matrix analysis to measure how many of these three-body systems show a periodic Doppler-shifted gravitational wave perturbation detectable by LISA. We report the number of circumbinary planets (CBPs) and (BDs) that can be detected by LISA for various combinations of mass and semi-major axis distributions. We identify pessimistic and optimistic scenarios corresponding, respectively, to 3 and 83 (14 and 2218) detections of CBPs (BDs), observed during the length of the nominal LISA mission. These detections are distributed all over the Galaxy following the underlying DWD distribution, and they are biased towards DWDs with higher LISA signal-to-noise ratio and shorter orbital period. Finally, we show that if LISA were to be extended for four more years, the number of systems detected will be more than doubled in both the optimistic and pessimistic scenarios. Our results present promising prospects for the detection of post-main sequence exoplanets and BDs, showing that gravitational waves can prove the existence of these populations over the totality of the Milky Way. Detections by LISA will deepen our knowledge on the life of exoplanets subsequent to the most extreme evolution phases of their hosts, clarifying whether new phases of planetary formation take place later in the life of the stars.


Introduction
In an epoch in which the field of exoplanets is moving at a fast pace and groundbreaking discoveries are made, very little is known about the ultimate fate of planetary systems.In the Milky Way more than ∼97% of the stars will turn into a white dwarf (WD), meaning that the vast majority of the more than known 3000 planet-hosting stars will end their life as WDs.Can their planets survive stellar evolution?Theoretical models indicate that a planet can endure the host-star evolution if it avoids engulfment or evaporation throughout the red giant or/and the asymptotic giant branch phases (e.g.Livio & Soker 1984;Duncan & Lissauer 1998;Nelemans & Tauris 1998), where survival itself depends, among various parameters, on the initial semi-major axis and planetary mass (Villaver & Livio 2007).For what remains of the planetary system the complex longterm orbital evolution, consequent to stellar evolution, may yield to planet ejections and/or collisions (e.g.Debes & Sigurdsson 2002;Veras et al. 2011;Veras 2016;Mustill et al. 2018).Besides, if migration or scattering occurs towards the proximity of the Roche limit, strong tidal forces can further crush the planetary cores (Farihi et al. 2018), like in the case of the planetesimal found shattering around WD 1145+017 (Vanderburg et al. 2015).Such a fragmentation process consequently enables the formation of a debris disc, made of metal-rich planetary material, which could in turn accrete onto the WD, polluting its atmosphere (e.g.Jura et al. 2009;Farihi et al. 2010;Farihi 2016;Veras 2016;Brown et al. 2017;Smallwood et al. 2018).A&A 632, A113 (2019) metal-rich material accreting onto these stars must be present.There are several sources of WD pollution proposed in the literature.This WD pollution could originate from planetary material (i.e. from circumstellar debris discs as previously explained), moons via planet-planet scattering (Payne et al. 2016(Payne et al. , 2017)), or comets (Caiazzo & Heyl 2017).It could also be from perturbations created by eccentric high-mass planets, which drive substantial asteroids or minor bodies to the innermost orbital region around the star (which in some cases is within the stellar Roche limit), thereby yielding to tidal fragmentation (e.g.Frewen & Hansen 2014;Chen et al. 2019;Wilson et al. 2019).The last source of WD pollution is currently preferred in the community.Pollution of WDs in wide binaries may also be caused by Kozai-Lidov instabilities, which can cause the orbit of objects such as planets, to intersect the tidal radius of the WD, causing their distruction (Hamers & Portegies Zwart 2016;Petrovich & Muñoz 2017).
Overall WD pollution studies support the evidence of dynamically active planetary systems orbiting WDs.Nonetheless, because of the intrinsic low luminosity of these stars, no planets have been detected yet around single WDs, but an intact planetesimal has been observed inside a debris disc belonging to a WD (Manser et al. 2019).

Generations of circumbinary post-common envelope exoplanets
Contrary to the single star case, P-type exoplanets (Dvorak 1986) have been detected orbiting binary stars in which the higher mass component has already grown to be a WD (i.e. the mass of its progenitor is M * 10 M ); the second component is usually a low-mass star that will become a giant later in its life (i.e.NN Ser, HU Aqr, RR Cae, UZ For, and DP Leo; Beuermann et al. 2010Beuermann et al. , 2011;;Qian et al. 2011Qian et al. , 2010Qian et al. , 2012;;Potter et al. 2011).These discoveries prove that planets can survive at least one common envelope (CE) phase, i.e. a shared stellar atmosphere phase typical of close binary stars, which happens when one of the binary components becomes a giant (see Sect. 3.1 for more detail on this phase).Surviving planets (in this case first generation planets) are usually called post-main sequence exoplanets or post-CE exoplanets, and they are more likely to survive around evolving close binary stars than around evolving single stars (Kostov et al. 2016).Only a small amount is known about these planets, but they are extremely interesting as they provide a link between planetary formation and fate, as well as constraints on tidal, binary mass loss, and radiative process (Veras 2016).Detection and study of these bodies can also provide us with further information about planetary formation processes.There is the interesting hypothesis that some of these known post-CE planets belong to a "new generation", i.e. they have formed after the first CE phase (e.g.Zorotovic & Schreiber 2013;Völschow et al. 2014).A study by Kashi & Soker (2011) has shown that because of angular momentum conservation and further interaction with the binary system, 1-10% of the ejected envelope does not reach the escape velocity.This material remains bound to the binary system, falls back on it, flattens, and forms a circumbinary disc, which could provide the necessary environment for the formation of a second generation of massive exoplanets (Perets 2010;Völschow et al. 2014;Schleicher & Dreizler 2014).On the other hand, as already mentioned for the single star case, some sub-stellar bodies (e.g.first generation exoplanets, asteroids, and comets), whose orbit is small and/or eccentric enough, could be tidally disrupted during the CE phase, creating a circumbinary disc of rocky debris, out of which new terrestrial exoplanets can grow (Farihi et al. 2017).In both cases photoheating from the binary, photoionisation, radiation pressure, and differences in the magnetic field, would likely be responsible for influencing the discs in different ways, causing second generation planets to differ from first generation planets (Perets 2010;Schleicher & Dreizler 2014;Veras 2016).
Another possibility is the existence of a hybrid generation: first generation planets that survive the first CE phase and may have been subject to mass loss throughout the whole process.Either way, the resulting planet/planetesimal could now accrete on the disc material, producing more massive planets on higher eccentricity (Armitage & Hansen 1999;Perets 2010).The outcome would be a planet with a first generation inner core and second generation outer layers.In this case the formation of a giant planet could be faster than for first generation giant planets.
The same hypothesis is similarly applicable if the binary overgoes a second CE phase, i.e. the low-mass star overflows its Roche lobe and shares its atmosphere with the WD companion.After this stage we might have either a third generation of exoplanets forming around a double white dwarf (DWD) system, a hybrid generation, or previous surviving generations.To date no exoplanets are known orbiting a DWD (Tamanini & Danielski 2019), and the only circumbinary exoplanet known orbiting a system with two post-main sequence stars (i.e. a WD and a millisecond pulsar) is the giant PSR B1620-26AB b; this planet is also the first circumbinary exoplanet confirmed (Sigurdsson 1993;Thorsett et al. 1993).Because the planetary system PSR B1620-26AB is the result of a stellar encounter in the Milky Way plane (Sigurdsson et al. 2003), it is not directly representative of a standard (i.e.isolated) binary planetary system evolution.
Possibly because of an observational bias, all the post-CE planets discovered until now are giant planets with masses M ≥ 2.3 M J and semi-major axes a ≥ 2.8 au.The most successful technique used for their detection is eclipse timing variation (ETV), which is sensitive to wider planetary orbits and hence requires a long observational baseline to precisely time the eclipses.Also, ETV typically suffers from a lack of cross validation and errors that are not uncommon, for example, the lack of accurate timing in the instrumentation used or procedures used to place the recorded times onto a uniform timescale corrected for light travel time (Marsh 2018).Any small inaccuracy or analysis imprecision could lead to uncertainties in the validity of a planet in the system, with the planets potentially being the wrong interpretation of the Applegate mechanism (Applegate 1992).
Recently Tamanini & Danielski 2019 showed the possibility to detect Magrathea-like (Adams 1979) planets, i.e. circumbinary exoplanets orbiting DWDs by using the Laser Interferometer Space Antenna (LISA; Amaro-Seoane et al. 2017) to measure the characteristic periodic modulation in the gravitational wave (GW) signal produced by the DWD.Compared to the classic detection methods, the GW approach has the advantages that this method is (i) able to find exoplanets all over the Milky Way and in other close-by galaxies; (ii) not limited by the magnitude of the WDs, but on the parameters from which the GW depends (see Sect. 3.3); and (iii) not affected by stellar activity, which is an issue reported in electromagnetic (EM) observations.

Brown dwarfs
Brown dwarfs (BDs) are by definition bodies that are not massive enough to fuse hydrogen in their interior stably, but are mas-The focus of this work is to follow up on Tamanini & Danielski (2019) and to quantitatively estimate the LISA detection rates of circumbinary exoplanets as well as circumbinary BDs. Brown dwarfs have masses ranging between the stellar and planetary domain; nevertheless, while the difference with stars is well defined, the separation with planets is still an open subject of discussion.The different nature of these objects could be either based on their intrinsic physical properties (Chabrier & Baraffe 2000) or on their formation mechanism (Whitworth et al. 2007).Furthermore, more recently Hatzes & Rauer (2015) analysed the density versus mass relationship for objects with mass ∼0.01 M J < M < 0.08 M , and identified three distinct regions that are separated by a change in slope in such a relation (at M = 0.3 M J and M = 60 M J ). Above M = 60 M J , but lower than M = 0.08 M , the BDs domain, and below that limit (but above M = 0.3 M J ) the giant planets domain.
Because of this ongoing discussion we hence decided to not limit our analysis to the mass domain reported in Tamanini & Danielski (2019), but to account for a larger mass range, up to the stellar limit.Consequently throughout this manuscript, for simplicity we define a sub-stellar object (SSO) to be a celestial body with mass less than 0.08 M (the hydrogen burning limit, which includes the upper uncertainty by Whitworth 2018).This category is divided between CBPs and BDs.As in Tamanini & Danielski (2019) we define the former as objects with mass M ≤ 13 M J (the deuterium burning limit) and the latter as those with mass 13 M J < M < 0.08 M .For simplification only the mass and no spectroscopic and/or formation mechanism classification are accounted for in this work.
The outline of this manuscript is as follows: in Sect. 3 we present the characteristics of populations used in the investigation, and we summarise the GW detection method discussed in Tamanini & Danielski (2019).In Sect. 4 we report CBPs and BDs detection rates, with their error analyses, for both the LISA nominal mission and for a possible extension of four more years.We discuss the implications of our results in Sect. 5 and we conclude in Sect.6.

Method
To reach the scope of this study we worked throughout two different stages.First, we constructed a population of Galactic detached DWDs with circumbinary exoplanets/BDs.To do so we injected a simulated population of SSOs into a synthetic population of DWDs (Korol et al. 2017).Such a DWD population was specifically designed to study the LISA detectability of these binary WDs, and it was employed in the LISA mission proposal (Amaro-Seoane et al. 2017).Second, we used the method described in Tamanini & Danielski (2019) to measure how many SSOs LISA will be able to detect.
In Sect.3.1 we summarise the most important features of the DWDs population.In Sect.3.2 we provide details about the SSOs population injection process.In Sect.3.3 we summarise the method used for the LISA GW detection of a circumbinary SSOs.

LISA DWD population
Our method relies on the binary population model which Toonen et al. (2012) obtained using binary population synthesis code SEBA, originally developed by Portegies Zwart & Verbunt (1996) and later adapted for DWDs by Nelemans et al. (2001b) and Toonen et al. (2012).The progenitor population is initialised by randomly sampling initial distributions of binary properties with a Monte Carlo technique.Specifically, the mass of the primary star is drawn from the Kroupa initial mass function in the range between 0.95 and 10 M (Kroupa et al. 1993).The mass of the secondary star is derived from a uniform mass ratio distribution between 0 and 1 (Duchêne & Kraus 2013).A log-flat distribution and a thermal distribution are adopted for the initial binary orbital separations and binary eccentricities, respectively (Abt 1983;Heggie 1975;Duchêne & Kraus 2013).The initial binary fraction is fixed to 0.5 value.The SEBA code evolves binaries until both stars turn into WDs and beyond up to the present time.More details and discussion on the sensitivity of the binary population synthesis outcome to the aforementioned assumptions are given in Toonen et al. (2012Toonen et al. ( , 2017)).The adopted model has also been recently tested against observations of both single WDs and WDs in binary systems (including DWDs) in the solar neighbourhood by Toonen et al. (2017).In particular, the adopted model currently better represents the space density of DWDs derived from a spectroscopically selected sample of Maoz et al. (2018).
One of the most impacting assumptions in DWD population synthesis is the prescription for the CE evolution (e.g.Toonen et al. 2017).As mentioned in Sect. 1, CE is a short phase of the binary evolution in which the more massive star of the pair expands and engulfs its companion (Paczynski 1976;Webbink 1984).During the CE phase the binary orbital energy and angular momentum can be transferred to the envelope because of the dynamical friction that the companion star experiences when moving through the envelope.Typically, this process A113, page 3 of 15 A&A 632, A113 (2019) Fig. 1.Geometry of the outer three-body system (DWD+planet/BD) and inner compact two-body system (DWD).The quantities û and û denote the directions perpendicular to the outer and inner orbital planes, respectively.The acronyms LoS and CoM stand instead for line of sight and centre of mass (of the whole three-body system).
is implemented in the binary population synthesis either by parametrising the conservation equation for energy (through the α parameter) or that for angular momentum (through the γ parameter) (see Ivanova et al. 2013, for a review).In particular, the γ-prescription was introduced with the aim to reconstruct the evolution path of observed DWDs by Nelemans et al. (2000); Nelemans & Tout (2005).In the model adopted for this study, γα, we allowed both parametrisations; the γ-prescription was applied unless the binary contains a compact object or the CE is triggered by a tidal instability.It has also been shown that γα model describes observations better than the model in which only α-prescription is employed (Toonen et al. 2012).Future optical surveys such as the Large Synoptic Survey Telescope (LSST; LSST Science Collaboration 2009) will provide large samples of new DWDs that will help to further constrain CE evolution for these systems (Korol et al. 2017).
Next, we distributed DWDs in a Milky Way-like galaxy according to a star formation history.We adopted a simplified Galactic potential composed of an exponential stellar disc and a spherical central bulge.Similarly to Ruiter et al. (2009); Lamberts et al. (2019), we found that the contribution of the stellar halo to the total amount of detectable GW sources is at most of a few percent.Thus, it is not included in this study.We populated the disc according to the star formation rate (SFR) from Boissier & Prantzos (1999) and assumed the current age of the Galaxy to be 13.5 Gyr.To model the bulge of the Milky Way we doubled the SFR in the inner 3 kpc as in Nelemans et al. (2001a).The detailed description of the Galactic model is presented in Korol et al. (2019).Finally, we assigned binary inclination angle i b , drawn from a uniform distribution in cos i b .Thus, each DWD in the catalogue is characterised by seven parameters: m 1 , m 2 , P b , i b , the Galactic latitude l and longitude b, and the distance from Sun d (see Fig. 1).
To obtain a sub-sample of DWDs detectable by LISA we employed the Mock LISA Data Challenge (MLDC) pipeline, designed for the analysis of a large number of GW sources simultaneously present in the data (e.g.Littenberg et al. 2013;Cornish & Robson 2017).This is realised throughout an iterative process that is based on a median smoothing of the power spectrum of the input population to compute the overall noise level (instrument plus confusion from the input population).The resolved sources (i.e.those with S /N > 7) are extracted from the data until the convergence.We adopted the LISA noise curves and orbits according to the latest mission design, the nominal mission duration of four years and the extended mission duration of eight years (Amaro-Seoane et al. 2017).
We find approximately 26 × 10 3 and 40 × 10 3 detached DWDs with S /N > 7 for the nominal four years and extended eight years of the LISA mission duration, respectively.Figure 2 illustrates the distribution of detected DWD in our mock Galaxy showing that GW detections can map both disc and bulge at all latitudes.We represent the mean signal-to-noise ratio (S/N) per bin in colour.
We note that in this work we focus on detached DWD binaries only.In principle, other Galactic binaries composed of compact objects (such as WD -neutron star and double neutron stars) and accreting systems could also host a CBP/BD.However, these are significantly less abundant in the Milky Way (e.g.Nelemans et al. 2001a), and thus would not affect much our estimates.In addition, GW signal of accreting systems contains an imprint of the mass-transfer process, which could affect the detection of circumbinary companions.We leave these investigations for future work.

Exoplanet and brown dwarf injection
Since the WD pollution effect supports evidence of dynamically active planetary systems around single WDs (Sect.2.1) and since no data are available for the binary WD case, we set the WD pollution upper limit occurrence rate (i.e.50% Koester et al. 2014) to be the occurrence rate (O.R.) of the synthetic population of SSOs orbiting DWD.We neglected the presence of an external third star and we assumed that pollution derives from asteroidal or moon material, rather than cometary material.We also rejected exceptions such as the capture of a free-floating planet at thousands of astronomical units, and we assumed that each DWD can harbour only one SSO; we briefly discuss the implications of considering multiple circumbinary objects in Sect. 5.
For the following we note that co-evolution of the binary plus SSO was neglected, and that the SSO population was injected into already formed WD-WD systems in which the stability criterion (P 4.5 P b ) of Holman & Wiegert (1999) was always satisfied.
In accordance with the pollution O.R. employed in this investigation, we set the SSO maximum distance (a) to be the approximate maximum limit for pollution to occur.Given that the maximum distance at which those asteroids reside around DWDs is completely unconstrained (Veras et al. 2019), we assumed 200 au to be a reasonable distance at which the SSO could still perturb asteroids which lie outwards or inwards towards the binary.
We set a uniform SSO inclination in cos i (cf.Fig. 1) and uniform initial phase φ 0 between 0 and 2π.Given that the planet distribution function is unknown and that no compelling physical motivation for a specific model at wide separations exists for these systems, we tested a combination of various semi-major axis a and SSO mass M distributions, commonly presented in the literature, to measure the number of possible detection of both CBPs and BDs.More specifically we defined the semi-major axis distributions as follows: (A) uniform distribution U a (0.1-200 au); and (B) log 10 uniform distribution log U a (0.1-200 au); and (C) log-normal distribution f (x) = A e ln(x)−µ/2σ 2 /(x2πσ), where A is the amplitude, µ the mean of the log-normal, and σ the square root of the variance.Meyer et al. (2018) give more details and specific values of the parameters; and (D) power-law distribution a −0.61 (0.1-200 au) (Galicher et al. 2016).

LISA detection of a third sub-stellar object
To model the perturbation induced by the SSO on the GW signal emitted by the DWDs, we followed the procedure presented in Tamanini & Danielski (2019).Figure 1 shows the geometry of the three-body system under consideration.The motion of the DWD around the centre of mass of the three-body system modulates the GW frequency through the well-known Doppler effect.The resulting frequency observed by LISA is written as where v is the line-of-sight velocity of the DWD with respect to the common centre of mass, while f GW is the GW frequency in the reference frame at rest with respect to the DWD centre of mass.Since the DWDs observed by LISA do not merge before a time much larger than the observational lifetime of the mission, we can effectively model the emitted frequency with a Taylor expansion around a constant value and only keep the first order term where f 0 is the frequency when LISA starts taking data and f 1 is its first derivative evaluated at the same time.The line-of-sight velocity of the DWDs is instead given by where we defined the parameters and the orbital phase both derived assuming an SSO circular orbit.In the expressions above P is the SSO orbital period, M is the SSO mass, M b is the DWD total mass, ϕ 0 is the outer orbital initial phase, and i is the SSO orbital inclination (cf.Fig. 1).The phase of the waveform observed by LISA is then given by where Ψ 0 is a constant initial phase.The main contribution of the Doppler frequency modulation Eq. ( 1) consists in a periodical shift of the GW frequency towards higher and lower values around f 0 .This effect is qualitatively depicted in Fig. 3 in which the Doppler modulation has been extremely exaggerated with respect to the perturbation induced by a SSO on a DWDs.In the real case the modulation timescale, of the order of roughly years, is much longer than the period of the GW produced by the binary, roughly minutes, implying that the effect would not be visible by eye.For each DWDs in our mock catalogue we can thus build a waveform depending on 11 parameters: 8 parameters associated with the DWD, namely ln(A), Ψ 0 , f 0 , f 1 , θ S , φ S , θ L , φ L , and 3 parameters associated with the SSO orbit, namely K, P, ϕ 0 .In this case θ S , φ S , θ L , φ L are the two sky localisation angles and the two angles defining the orbital geometry of the DWDs, respectively (directly related to the inclination i b and polarisation angle ψ b ; see e.g.Cornish & Larson 2003).
To simulate the response of LISA and perform a parameter estimation of the GW waveform, we followed Tamanini & Danielski (2019) again.The full expressions for the two linearly independent signals observed by LISA h I,II (t), including the LISA antenna pattern functions and effects due to its orbital motion, can be found in Cutler (1998); Takahashi & Seto (2002); Cornish & Larson (2003).For the sake of simplicity we are not reporting those expressions in this work.The S/N of each event is computed as the following: where T obs is LISA observational time period and S n ( f 0 ) is the one-sided spectral density noise of LISA computed at f 0 .Parameter estimation is performed by employing a Fisher information approach, where we define the Fisher matrix as Marginalised 1σ errors for each waveform parameter are thus estimated from the square root of the diagonal elements of the covariance matrix, the inverse of the Fisher matrix.

Results
We focus first on the properties of the detected population of SSOs (Sect.4.1), showing also how the numbers improve for an extended eight-year LISA mission (Sect.4.2).We then discuss the recovered accuracy on the waveform parameters in Sect.4.3.

LISA detection of SSOs
As in Tamanini & Danielski (2019) we assume that a SSO (either a CBP or a BD) is detected if both K and P parameters are measured with a relative accuracy better than 30%.For every injected SSO population, defined by a combination of semi-major axis a and mass M distributions (see Sect. 3.2), we counted the number of SSOs whose GW perturbation can be detected by LISA.We report in Table 1 the total number and percentage of circumbinary exoplanets and BDs detected during the nominal LISA mission length.For both CBPs and BDs we identified optimistic, pessimistic, and intermediate scenarios.While the first and second represent the cases in which the highest and lowest numbers of CBPs (or BDs) are detected, the last scenario represents the case with the median number of detections, rounded by excess.Among the available combinations, the B1 scenario, i.e. that whose injected SSO population follows a logarithmic a distribution log U a , and uniform M distribution U M (see Sect. 3.2), is the optimistic case for both CBPs and BDs with 83 and 2218 detections, respectively.The intermediate scenario is represented by C1 (log-normal a ; U M ), and B2 (log U a ; M −1.31 ), for CBPs and BDs with 18 and 316 detections, respectively.The pessimistic scenario is represented by A1 (U a ;U M ) and A2 (U a ; M −1.31 ) with 3 and 14 detections for CBPs and BDs, respectively.We plot in Fig. 4 the location in the Milky Way of the detections for the three CBPs scenarios together with a zoomin on the solar neighbourhood for the optimistic scenario.From Fig. 4 it is easy to understand that LISA will be able to observe CBPs and BDs orbiting DWDs everywhere in the Galaxy.
Furthermore, for the six scenarios selected above, Figs.A.1 and A.2 (currently appearing after the references) show the distribution of detected CBPs and BDs, respectively, over the CBP/BD separation from the DWD (a), the mass of the CBP/BD (M), the CBP/BD orbital inclination (i), the parameter K, the CBP/BD period (P), the DWD period (P b ), the DWD S/N, the distance from the Earth (d), the DWD chirp mass (M c ), and the total DWD magnitude measured in the Gaia G band (G DWD ).To highlight possible observational biases, in Figs.A.1 and A.2 we also show the underlying distribution of injected CBPs/BDs in grey.

Detection rates for an extended LISA mission
We repeated our analysis for an eight-year LISA mission, corresponding to a possible realistic extension beyond the nominal four-year mission; this can also approximately be considered as ten years of mission operations, the maximal envisaged extended duration, with duty cycle of 80% similar to the LISA Pathfinder (Armano et al. 2016).We used the catalogue of 40 × 10 3 DWD detected over the eight years of mission presented in Sect.3.1 injecting SSOs according to the optimistic and pessimistic scenarios only.The total detections of CBPs and BDs, together with the percentage over the total number of DWDs detected by LISA, are reported in Table 2.In the optimistic scenario (B1) we find a total of 215 (4684) detected CBPs (BDs), corresponding to the 0.822% (17.913%) of the total population of detected DWDs, and to an improvement of the 259% (211%) over the detections of the nominal four-year mission.The numbers for the pessimistic scenarios, (A1) for CBPs and (A2) for BDs, are instead 8 (43) detected CBPs (BDs), corresponding to 0.02% (0.107%) of the total DWD population, and to an improvement of the 267% (307%) over the 4 yr detections.
In the hypothesis of eight years of observations, LISA will be able to detect SSOs with longer period P and consequently larger separation a.These are the SSO orbital parameters that present a significant improvement with respect to the four-year case, i.e. for which a larger range of measured values is recovered, instead of only a larger number of detections within the same parameter interval.We plot for comparison in Fig. 5 the distributions (injected and recovered) of these two quantities Each plot shows the location of the binary WD system with a planetary companion (red) and BD (green) detection through GWs.In each panel we also plot the known detected exoplanets's host-star (see legend for colour scheme; data from https:// exoplanetarchive.ipac.caltech.edu).We note that data overlay a face-on black and white image of the Milky Way for Galactic location reference purposes.
for both time frames of eight and four years.In general the longer the LISA observational period, the longer the SSO period and separation that will be recovered.This can be easily visualised in Fig. 5 where the eight-year bulk of detected CBPs (BDs), presents periods up to ∼12 (∼30) yr, compared to only ∼6 (∼10) yr over a four-year mission.A similar trend is observed for the separation a, as of course this is directly related to the period.
A113, page 7 of 15 (B1) -8 yrs Fig. 6.Detections by LISA of circumbinary exoplanets (red) and BDs (blue) in the optimistic (B1) scenario for 4 yr (left panel) and 8 yr (right panel) of observations.The SSOs mass (M, or M sin(i) for those planet whose inclination is not known yet), as a function of the SSO-to-binary separation (a) is shown, together with values of known exoplanets (data were taken from https://exoplanetarchive.ipac.caltech.edu).
The horizontal dotted line corresponds to the limits in mass we explored: 13 and 0.08 M J .Similarly, Fig. 6 shows the mass M of the detected SSOs (B1 scenario), as a function of the semi-major axis a.Both detections obtained in a four-and eight-year time frame are shown next to each other for comparison.We note that during eight-year observations LISA will generally be able to identify a larger number of lighter exoplanets below 2 M J .This is because a longer baseline would allow us to disentangle the gravitational pull of the small exoplanet from the gravitational waveform, and consequently consent to measure K and P with a relative precision better than 30%.Similarly, the SSO range of detectable separations a are roughly doubled in an eight-year mission.Such an increase in the parameter space enables LISA to be more compatible with imaging surveys, but also with the bulk of radial velocity surveys, for which a good overlap is already visible during the nominal mission.During both the four-year and eight-year surveys there is no real comparison with the bulk of the transit population, but this is barely a feature of the constructed SSO population, which we limited at 0.1 au.Synergies are possible between 0.1 au and 1 au, however.

Error distributions of third-body parameters
In this subsection we look at the distributions of 1σ errors for the parameters K and P, i.e. the third-body parameters which are interesting from an observational perspective.We report first the best and average error with which these parameters are recovered for detectable systems, i.e. for systems which already have relative errors on both K and P estimated to be below 30%.Again we focus on the optimistic, median, and pessimistic scenarios, as selected above.We only analyse data collected with a nominal four-year LISA mission.
In the CBP optimistic (B1), median (C1), and pessimistic (A1) scenarios, we obtain an average relative error on K of the 14.9, 14.9, and 21.0%, respectively.The best recovered K values are instead measured with relative errors of 0.86, 2.1, and 10.0%, respectively, for the same three models above.The average relative errors on P are instead 4.4, 10.9, and 16.5% for the three models (B1), (C1), and (A1), respectively.The best recovered P values are measured with an estimated relative error of 0.040, 0.23, and 9.6%, again, respectively, for the same models mentioned above.
For BDs we obtain instead an average relative error on K of 12.0, 12.3, and 11.8%, for the optimistic (B1), median (B2), and pessimistic (A2) scenarios, respectively.The same parameter is recovered with a best relative error of 0.10, 0.16, and 3.4%, respectively, in the same three scenarios.The average relative errors on P are instead 3.3, 3.8, and 5.3%, while the best recovered relative errors are 0.0049, 0.013, and 0.28%, again for the scenarios (B1), (B2), and (A2), respectively.

Discussion
The results presented above show that during the four years of its nominal mission, LISA will be able to detect from a few to a few tens of CBPs down to a few Jupiter masses and up to a few astronomical units in separation.Analogously we find that LISA will likely detect from several to few thousands BDs in roughly the same semi-major axis range.These observations will be of fundamental importance for the field of exoplanetary science.As shown in Fig. 4 in the optimistic scenario LISA detections will be distributed all over the Milky Way, but even in a pessimistic scenario we would be able to detect at least some exoplanet far outside the solar neighbourhood.In our study we only considered a Galactic population of DWDs, but we stress that LISA will be able to observe DWDs even in nearby galaxies (e.g. in the Magellanic Clouds and M31; Korol et al. 2018) and consequently, if conditions are optimal (e.g.high values of S/N, f 0 , M, ...), it could also detect extragalactic bound CBPs/BDs (Tamanini & Danielski 2019), possibly leading to the discovery of the first bound extragalactic SSO.Meanwhile, expanding the exoplanetary census beyond the local Galactic environment with GW observations, will help integrate the information collected (and that will be collected) with current (and future) EM surveys, and it will provide a more robust and unbiased statistic on the life of giant exoplanets.If this population is not detected, given the mass-separation parameter space accessible to LISA, we can confidently say that SSO do not survive a second CE phase and are either destroyed or ejected from the system.But whether or not the population exists, beyond the pure survival rates we will set constraints on the dynamical evolution of the tertiary body consequent to the CE phases and the binary mass ejections.A more robust statistic will also allow us to have a better understanding on the existence and nature of planetary generations, by testing the dynamical stability timescale of the systems and identifying if any correlation between the orbital properties of the systems is present.Inevitably, if the range of parameters detected is large, for instance if exoplanets are both found orbiting short (within the maximum radius of the CE of the progenitors), and wide orbits (where giant planets usually form), depending on the binary cooling time we could gather information on both formation and migration processes.The same reasoning applies to the BDs, for which these further studies would help address their difference from planets.
Our results also suggest that an extended LISA mission, up to eight years, will yield a larger parameter space than the one spanned by the nominal four-year mission, and a more robust statistic.The number of detected CBPs and BDs will more than double, implying an incremental trend which grows more than linearly.This is mainly because a longer observational window allows us to unlock the detections of SSOs with longer period, as clearly shown in Fig. 5.To give a numerical example we note that, in the B1 scenario, over four years LISA will detect 0.32% (8.48%) of DWDs with a CBPs (BDs), while over eight years it will detect 0.82% (17.91%) of them (cf.Tables 1 and 2).This clearly shows that a higher percentage of the underlining SSO population will be detected with a longer observational time period, providing another scientific case for an extension of the LISA mission beyond its nominal four-year lifetime.If the maximal envisaged duration of the mission is considered, namely ten years, then the results within our optimistic scenario suggest that LISA should be able to detect more than ∼250 new CBPs and more than ∼6000 new BDs.
For our analysis we tested detection rates of single SSOs (1 M ⊕ < M < 0.08 M ) orbiting DWDs with the future LISA mission.We assumed circular orbits that satisfy the stability criteria by Holman & Wiegert (1999) and no mass transfer among the two WDs.Galactic DWDs represent a stellar population older than 100 Myr, ergo they are not expected to follow the spiral structure of the Milky Way (see Fig. 4).For this reason, in our fiducial simulation we neglected the spiral structure itself and we distributed binaries in a smooth exponential disc potential with a prominent central bulge.Even when adopting a high resolution numerical simulation for the mass distribution, the contrast between the spiral arms and the background disc in GWs is too low to be detected (Wilhelm et al., in prep.).We also note that the space density of DWDs in the solar neighbourhood is three orders of magnitude lower than that of main sequence stars (e.g.Hollands et al. 2018).This translates into a low detection statistics when comparing GW detection with currently used EM methods for the detection of exoplanets (see Fig. 4, top right panel).However, because the GW signal scales as 1/d instead of 1/d 2 , which is typical of EM observations, exoplanets can be detected farther away than will ever be possible at optical wavelengths (out to the far side of the Milky Way and satellite galaxies e.g. the Magellanic Clouds; Korol et al. 2018Korol et al. , 2019)).
The core composition of WDs could be a relevant element with respect to the O.R. of second generation exoplanets.Given the enrichment with heavy elements of the envelope of CO-progenitors, occurring at the end of the asymptotic giant branch, planets should be more frequent around compact CO-core DWDs (than around DWDs with a He-core WD primary; Zorotovic & Schreiber 2013).The CO WDs should have much higher metallicities than He WDs, which is a characteristic that in a protoplanetary disc would promote the formation of a greater number of high-mass giant planets (Johnson & Li 2012).In our mock population 38% of DWDs are He-He type, 27% are He-CO, 33% are CO-CO, and 2% are CO-ONe (see Toonen et al. 2012 for more details and formalism).We note however that these percentages depend on the adopted stellar evolution tracks in the binary population synthesis and can differ significantly from one model to another.We also note that the absolute magnitudes of DWDs in our mock population are derived from the WD cooling curves of pure hydrogen atmosphere model of A113, page 9 of 15 A&A 632, A113 (2019) Table 3.Among the detected SSOs in all 4 yr scenarios we report the systems with the least massive CBP detected (x, A2), with the highest S /N (y, B1), and with the CBP with longest period (z, C1).Notes.In this table the DWD S/N is denoted as S/N DWD rather than S/N as in the text.Holberg & Bergeron (2006).Thus, by construction all binaries in the simulation are composed of DA1 WDs.

S /N
As an example of LISA capabilities, in Table 3 we report the parameters of the system with the least massive planet detected in this analysis, the system with the highest S/N of the DWDs, and the system with the longest planetary period, for which the sky localisation error boxes measure 1.77 2 , 0.14 2 , and 12.6 2 , respectively.As suggested by Table 3 and confirmed by Fig. A.1, binaries with P b < 10 min are optimal for detecting circumbinary companions.This result was expected because low orbital period DWDs emit high frequencies GWs; it is thus easier to discern the Doppler perturbation in the GW waveform produced by the circumbinary object (Tamanini & Danielski 2019).Moreover, the higher the S/N, the easier it is to detect the same perturbations, meaning that detections of both CBPs and BDs are biased towards DWDs with high S/N and low orbital period (within the global DWD population that LISA will observe).This can be quickly confirmed by Figs.A.1 and A.2 simply by comparing the recovered versus the injected distributions of both S/N and P b .Furthermore high S/N necessarily corresponds to high frequency DWDs for two reasons: the GW amplitude scales as f 2/3 and LISA is more sensitive at f ∼ 10 −2 Hz, where we find DWDs with shortest periods (a few minutes; Amaro-Seoane et al. 2017).Because of this, high frequency and high S/N sources can be detected anywhere in the Galaxy (as shown in Fig. 2), with a peak at ∼8.5 kpc due to the high density of DWD in the bulge (Fig. A.2).
With reference to short period binaries, the time τ these bodies will take before colliding can be approximated by meaning that, for a P b ≤ 10 min and a typical chirp mass M c = 0.2 M (Korol et al. 2017;Tamanini & Danielski 2019), the colliding time is τ 1.558 Myr.Even considering the maximum chirp mass of a DWD detectable by LISA, say 1 M , and its minimum orbital period, say 3 min, the DWD will not merge before 3300 yr.
Our detections present some events far in the tails of the observed distributions.These events are usually associated with a combination of high DWD S/N, high DWD GW frequency, and high SSO mass, which correspond to a stronger perturbation in the GW signal and which are thus easier to detect.Consequently it is not surprising that they can be detected even for unusual values of the SSO parameters.The CPB period distribution in Fig. A.1 for example shows few events with periods around six years, while the bulk of the distribution is set on periods shorter than four years.A numerical example in the C1 case is given by the system with the longest detected planetary period (8.82 yr) which report a S /N = 182, a GW frequency f 0 = 12.53 mHz, M = 11.11M J , and a = 3.67 au.This is even more evident for BDs, for which in the D1 scenario we detect a system with a BD whose orbital period is 19.5 yr, even if the global detections are set at P < 12 yr (see Fig. A.2, where an outlier with P ∼ 17 yr is also appearing in the (B1) scenario).This event is characterised by S /N = 169 f 0 = 12.37 mHz, M = 78.86M J , and a =7.2 au, which shows that such outliers can only be detected for systems with high S/N, high frequency, and high SSO mass.All this shows that, with ideal conditions, LISA could detect CBPs with periods up to P ∼ 10 yr and BDs with periods up to P ∼ 20 yr, with only four-year observations.The bulk of detections however are expected at P < 4 yr for CBPs and at P < 11 yr for BDs (see Figs. A.1 and A.2); only rare events appear at higher orbital periods.Moreover, as noted by Tamanini & Danielski (2019), the Fisher matrix approach adopted in this work might not be reliable for events with extremely high S/N and further more detailed data analysis techniques should be used to determine the real detectability of such systems.
In this work we accounted for only one circumbinary SSO for DWD, but observations show that multiple SSOs can orbit evolved binaries, i.e.NN Ser (b,c), UZ For (b,c), or HU Aqr, which is a system that hosts one giant planet and two BDs (Goździewski et al. 2015).Consequently, since multiple circumbinary objects could co-exist (see Veras & Gänsicke 2015 for a single WD case), our results report lower limits of detections in all the possible scenarios of mass and planet-to-binary separation distributions.We note however that additional SSOs, or even a low-mass star, orbiting the same DWD would complicate the GW signal detected by LISA because of the simultaneous Doppler perturbations of different circumbinary objects.This might worsen the precision with which the SSO orbital parameters are recovered, possibly leading to some detections being missed.Future analyses, which lie outside the scope of the present work, will be needed to explore the detectability of CBPs and BDs in systems with multiple orbiting SSOs or with a third star composing hierarchical triples with the DWD.Alongside with it, also studies on the dynamical stability of multi-planets, similarly to e.g.Mustill et al. (2014); Veras & Gänsicke (2015); Kostov et al. (2016); Mustill et al. (2018), but specific for SSO orbiting DWDs are needed to understand the dynamical grounds of these objects.This might need to take into account also the possibility of co-existing generations of planets, aspect that would necessarily make the analysis computationally expensive.
We mentioned in Sect.3.2 that no specific O.R. for planets orbiting a binary WD is available, therefore we used the atmospheric pollution frequency (25-50%) for single stars, robustly measured by Zuckerman et al. (2010) and Koester et al. (2014).Recently Wilson et al. (2019), using Spitzer and Hubble data, estimated the pollution rate in WDs in wide binaries to be 67 +10 −15 %, consistent within 2σ to the single WDs value measured in the same work (45 ± 4%) and to the rate value applied in our study.These rates are also consistent to the O.R. of planets transiting single WDs measured by van Sluijs & Van Eylen (2018), using K2 data.For Jupiter-size planets and planetary periods between 10.12 and 40 days, the authors calculated a detection probability of 53.3 ± 3.0% within a 68% confidence interval.
Given that the length of the survey was only 40 days, there is no information for larger periods and extrapolation at wider orbits would be highly inaccurate given the lack of physical constraints.We consequently decided the 50% upper hand limit reported by Koester et al. (2014) to be our reference value.The total number of detections in a scenario with 25% O.R. (lower hand limit) can be directly estimated from our results in Table 1 since all numbers scale linearly and can thus just be divided by 2. Specifically, in the optimistic scenario (B1) CBPs and BDs detections would be 42 (0.16) and 1109 (4.24%), respectively.In the CBP pessimistic scenario (A1) instead we would only get one detection in four years of observations, while in the pessimistic BD scenario (A2) we would count seven detections.The same reasoning applies to the eight-year results, for which numbers in Table 2 should just be halved.
Concerning BDs detections, we notice that in the optimistic scenario (B1) they amount to ∼27 times the number of CBPs (versus the 2.3 times factor for the pessimistic A2 case), representing the 8.48% of the binaries DWD population in the Milky Way.Such a result was expected given that, assuming the same mass of the binary, a more massive object would produce a larger motion of centre of mass of the three-body system (cf.Eq. ( 4)), and hence a larger shift in frequency, which is easier to detect.The mass distributions in Figs.A.1 and A.2 show this dependency very clearly.The residuals of the injected versus detected population of BDs, normalised to the injected population, in the B1 case, i.e. that presenting a more robust statistics, goes from 91% for BD masses between 15 and 20 M J to 73% for BD masses between 75 and 80 M J .On average thus SSOs with larger masses have a higher probability of being detected by LISA, as expected.The total BDs (i.e. over the mass range M > 13 M J ) normalised residuals (80%) are indeed smaller that the total CBPs (M ≤ 13 M J ) normalised residuals (96%), again for the optimistic (B1) scenario.
Besides, GWs do not allow for a direct measurement of the mass of the SSO.The mass can be estimated only once both K and P are known, and only after we assume a value (or a range of values) for both the binary mass ratio, and the SSO orbital inclination i, in analogy with radial velocity measurements (Tamanini & Danielski 2019).These considerations imply that without EM counterpart data it will be difficult to discern a CBP from a BD for masses around 13 M J , especially if the GW measurement is not precise enough; the needed level of precision depends on the specific case.Only an independent EM estimation of the binary total mass, the SSO orbital inclination, and the SSO radius will enable us to unambiguously characterise the nature of the GW detected SSO (see Sect. 5.2).
In this investigation we injected SSOs with masses up to 0.08 M .This was justified by the fact that the separation between the nature of planets and BDs is still uncertain.By applying the same reasoning the WD pollution, whose O.R. we used, could be also driven by low-mass BDs, i.e. very massive exoplanets.Because of this we took into account the largest possible mass range to cover both populations.However, had we assumed that the O.R. was only valid for planetary masses (M ≤ 13 M J ), i.e. by abruptly excluding the hypothesis that pollution could be also caused by objects able to at least burn deuterium, the CBPs detection rates would have been higher.On the other hand we would have not detected BDs, as none of them would have been injected.

The LISA duty cycle and system identification
Our study was based on the nominal LISA mission lifetime, i.e. four years of uninterrupted observations.However, during ∼30% of this time, LISA will not be acquiring scientific data because of expected maintenance operations (duty cycle).Nevertheless, even though the total effective observational period will be below four years, a periodic stop of scientific operation should not negatively impact the detection capability of LISA, at least for long-living GW sources if additional data analysis tools are employed (Baghi et al. 2019).Our results should thus not be affected by the duty cycle of LISA, albeit we note that a future dedicated investigation is required to address this aspect fully.
We expect DWDs to be very numerous in the Milky Way.Population synthesis studies predict that ∼10 6 DWDs have periods within the LISA frequency band, (e.g.Korol et al. 2019).Only 1% of these objects will have S /N > 7 and be individually resolved by LISA, while the overlapping signals of the remaining DWDs will sum up to form a Galactic noise background.Using the same mock DWD population as in this work (Fig. 2), Littenberg & Cornish (2019) show that the DWD confusion background is mostly confined between ∼0.4 and ∼4 mHz, meaning that at frequencies f > 4 mHz DWDs can be individually identified.The typical LISA error on sky localisation for DWDs is <10 deg 2 (Korol et al. 2019), although for DWDs with a detected SSOs, which we recall are biased towards higher frequency and higher S/N, the error is much lower.For example in our optimistic scenario (B1) the mean sky location accuracy of the DWDs with a detected CBP (BD) is 0.29 (2.83) deg 2 .The 74.5% (37.4%) of these systems are above f ∼ 4 mHz.This implies that LISA DWDs with a detected CBP/BD have higher chances to be spotted by EM telescopes.

What does the future look like?
Mainly because they are intrinsically faint and physically small, DWDs are difficult targets for optical telescopes.Typically, spectra of DWDs are virtually identical to those of single WDs, while their eclipses are very short (e.g.Rebassa-Mansergas et al. 2019).This drastically limits the observed volume, with the most distant detached DWD detected around ∼2.4 kpc (Burdge et al. 2019).Nonetheless, the number of DWDs detected with EM techniques are expected to increase substantially with the upcoming future all-sky and wide optical surveys, which also cover low Galactic latitudes, for example BlackGem (Bloemen et al. 2015), GOTO (Gravitational-wave Optical Transient Observer; Steeghs 2017), Gaia (Gaia Collaboration 2016), and LSST.There can be different detection strategies to identify the EM counterparts.For example, high cadence (several observations per night) photometric surveys such as ZTF (Zwicky Transient Facility; Bellm 2014) can be used to search for variable sources with an orbital period provided by LISA.Surveys with longer cadence (one observation in a few days) such as PTF (Palomar Transient Factory;Law et al. 2009) or LSST, can also be used for finding DWD EM counterparts.However, in the latter case it will be important to account for the orbital period derivative (also provided by LISA if the system is chirping) caused by GW radiation; see example of retrieving J153932.16+502738.8 in PTF archival data in Burdge et al. (2019).The work by Korol et al. (2017) in particular showed that at least 100 DWDs are expected to have GW counterparts.These predictions give us optimistic prospects for observing exoplanets and BDs around DWDs (but also other A113, page 11 of 15 A&A 632, A113 (2019) stellar compact object binaries), which would require monitoring DWDs for a few years.
Concerning a detection of a SSO signal with EM techniques we refer to the Methods section of Tamanini & Danielski (2019) for specific discussion on the EM synergies with GWs.We stress though that upcoming data from Gaia will highly increase the sample of giant planets in long-period orbits around binaries with FGK-dwarf primaries, located within 200 pc from the Sun (Sahlmann et al. 2015).Gaia will also detect tens or hundreds of planets (M > ∼1 M J ) around single WD that, combined with those that the PLATO 2.0 mission is also expected to find (Rauer et al. 2014), will increase this population statistic (should they exist in the favourable region of parameter space).This will allow us to begin placing limits on the masses of planets that can survive stellar evolution.Furthermore, Gaia, LSST and WFIRST will help detect free-floaters i.e. planets not bound to any star(s), which may help constrain the fraction of ejected planets due to mass loss (Veras et al. 2013;Veras 2016).These observations, combined with a continuous development of long-term dynamical evolution models of planetary systems, will help to acquire a more focused picture on the surviving life of exoplanets.Similarly, new studies on formation of second/third generation bodies orbiting post-CE binaries, as well as accretion studies of first and second generation objects, are important to address both the rates and orbital characteristics of the population investigated in this work and to understand whether the presence of a "surviving generation(s)" inhibits or promotes the formation process of new planets.

Possibly planet detection by LISA around ZTFJ1539+5027
The recently discovered ZTF J153932.16+502738.8, which has an orbital period of ∼7 min at a distance d = 2.34 kpc, is a great example of the potential of the multi-messenger observations with the aforementioned surveys together with LISA (Burdge et al. 2019).According to our detection definition in Sect. 3 we note that a planet with mass M = 1 M J orbiting J153932.16+502738.8 at a separation of 1 au, would not be detected by LISA.Using the measured orbital parameters of J153932.16+502738.8 (Burdge et al. 2019;Littenberg & Cornish 2019), setting Ψ 0 = 0 (initial DWD orbital phase), marginalising over ψ b (the polarisation angle unconstrained by EM observations), and assuming the planetary orbit to have ϕ 0 = π/2 and i = π/2 (most favourable orientation), LISA would measure its parameters as K = 31.4± 39.4 m s −1 (relative accuracy: 126%), P = 1.1043 ± 0.0857 yr (7.8%) and ϕ 0 = π/2 ± 1.08 rad (69%), where the sky location has been fixed to the real one measured for J153932.16+502738.8.We see that, although the planetary period is well constrained, K (and ϕ 0 ) are unconstrained, we would thus not be able to estimate the mass of the planet.We note however that an accurate measure of the planetary period could be extremely useful if taken in combination with other EM observations, which for J153932.16+502738.8 can be easily planned.This highlights the multi-messenger potential of LISA in terms of exoplanetary observations.The situation changes for more massive planets.If we consider a planet with 13 M J at the same separation of 1 au and repeat the analysis above, we find that LISA would be able to measure K = 405.1 ± 38.4 m s −1 (9.5%), P = 1.0967 ± 0.0067 yr (0.61%), and ϕ 0 = π/2 ± 0.0853 rad (5.4%).In this case the planet would be easily detected by LISA with accurate measurements of its orbital parameters.This also implies that any BD orbiting at the same separation from J153932.16+502738.8 would be detected by LISA and its parameters would be measured with even higher precision.
If instead J153932.16+502738.8 would appear at a distance of 9 kpc (implying S /N 37), well beyond the Galactic bulk where the majority of DWDs are expected to be detected by LISA, we would only be able to measure the same 13 M J planet with a relative accuracy of σ K /K = 96% and σ P /P = 10%, which again shows that, even if we will not access any information on the mass of the planet, we would still be able to measure its orbital period quite well.However in this case EM complementary observations will be impossible to obtain with the current instrumentation.
We finally compute the probability that an SSO in our optimistic (B1) population has the same f 0 , f 1 , and η of J153932.16+502738.8 (practically within 10% of these values).Among all 13 086 SSOs present in our population, only 14 (0.11%) have these characteristics, 12 of which are BDs and 2 CBPs.Of the 12 BDs, 5 are detectable by LISA while the 2 CBPs are not detectable.If we project these numbers on J153932.16+502738.8, and we recall that we are assuming a 50% O.R., we find that this system has a 17.9% probability of harbouring BDs detectable by LISA.The same reasoning cannot be performed for CBPs for which we can only conclude that the probability that J153932.16+502738.8 has a circumbinary planet detectable by LISA is very small.

Summary
In this work we quantitatively estimated the detection rates, by the LISA mission led by ESA, of circumbinary SSOs (i.e.planets and BDs) orbiting Galactic detached DWDs.To do so we injected a simulated population of SSOs into a synthetic population of already formed DWDs, with an O.R. of 50% i.e. the observed frequency of polluted WDs (Koester et al. 2014).We then applied the method presented by Tamanini & Danielski (2019) to probe how many systems we can identify to have a SSO perturbing the DWD GW signal because of its gravitational pull.Given that currently no theoretical and observational constraints are present to define such a specific population, we tested various combination of semi-major axis and mass distributions for estimating the number of detections over the course of the LISA nominal mission.Our analysis identified an optimistic and pessimistic scenario for which we counted a total of 63 (2218), and 3 (14) detections of circumbinary exoplanets (BDs) in the Milky Way, respectively.These numbers corresponds to 0.317% (8.482%), and 0.011% (0.054%) of the total DWDs visible by LISA, and these have more than doubled in a time frame of eight-year continuous LISA observations, corresponding to a realistic extended mission.In such a case the range of recovered planetary periods (and semi-major axis) would double for planets and increase almost threefold for BDs.The SSO detections that we found are also biased towards high frequency and high S/N binaries, as expected from basic considerations.The advantages of using the GW method for detection of CBPs and BDs comes from the fact that GWs are not affected by dust extinction and can be measured from all over the Milky Way and the Local Group.In constrast to EM techniques, this method is most efficient in the most dense regions of the Milky Way like the central bulge.A full investigation of a realistic observational strategy, including EM complementary observations, will be performed in future studies.

Fig. 2 .
Fig. 2. Signal-to-noise map of DWDs detected by LISA (4 yr) in the galactocentric Cartesian coordinate system.The colour represents the mean S/N per bin.The red triangle identifies the position of the LISA detector in our simulation.

Fig. 3 .
Fig.3.Qualitative example of a DWD waveform with (blue) and without (orange dashed) the presence of a third body.The Doppler modulation is extremely exaggerated for visualisation purposes.

Table 1 .Fig. 4 .
Fig. 4. Optimistic (top left, B1) with its zoom-in on the solar region (top right, heliocentric coordinates), intermediate (bottom left, C1), and pessimistic (bottom right, A1) scenarios.Each plot shows the location of the binary WD system with a planetary companion (red) and BD (green) detection through GWs.In each panel we also plot the known detected exoplanets's host-star (see legend for colour scheme; data from https:// exoplanetarchive.ipac.caltech.edu).We note that data overlay a face-on black and white image of the Milky Way for Galactic location reference purposes.

Table 2 .
Number of detections (over a total of 40 251 DWDs detected with 8 yr LISA mission) for the best and worst scenarios for both CBPs (B1 and A1) and BDs (B1 and A2).