© ESO 2022

1 Introduction

Whilst the number of detected exoplanets has increased sharply in recent years, relatively few have been found in the outer regions of planetary systems. This is likely due to the limitations of current detection techniques; radial-velocity (RV) and transit methods are insensitive to wide-separation planets, and direct imaging has yet to achieve the contrasts and sensitivities required to really probe outer-planet populations. However, debris discs provide an alternative, indirect method for exploring system outskirts; discs in these regions can be imaged at high resolution, and their locations and morphologies can bear signatures of dynamical interactions with planets (e.g. Wyatt et al. 1999; Wyatt 2008; Krivov 2010; Matthews et al. 2014). These disc features can be used to infer the presence and properties of otherwise undetectable planets, yielding unique insights into the outer regions of planetary systems.

A debris disc is therefore a planet probe, provided that planet-disc interactions are either ongoing or have previously occurred. Such interactions may be common, for a number of reasons. Firstly, the edges of both the asteroid and Kuiper belts are sculpted by dynamical interactions with Solar System planets (e.g. Moons 1996; Malhotra 2019), so it is reasonable to expect similar processes to occur in extrasolar systems. This is consistent with the emerging picture that many extrasolar debris discs may have edges consistent with planet sculpting (e.g. Faramaz et al. 2021; Marino 2021). Secondly, many extraso-lar systems contain several distinct debris discs at different radial locations, and the broad cavities between them are most simply explained by multi-planet clearing (Faber & Quillen 2007; Su et al. 2013). Thirdly, some debris discs bear features that are directly attributable to interactions with known companions in the same system (e.g. Mouillet et al. 1997; Lagrange et al. 2010). Finally, the fact that we observe extrasolar debris discs at all could imply that planet-disc interactions occur; the observed dust is thought to be released through collisions between larger debris, and a popular explanation is that planetary perturbations excite the random velocities of debris sufficiently for such collisions to be destructive. This is known as ‘planet-stirring’ (Mustill & Wyatt 2009), although alternative stirring mechanisms are also possible. If disc-planet interactions are therefore common, and possibly ubiquitous, then studying a large sample of debris discs should offer insights into the population of unseen planets.

In this paper we, use dynamical arguments to constrain planet properties from debris discs around 178 stars. We consider the planet(s) required to truncate and shape each of these discs and also investigate whether planets are required to stir them; the dynamical arguments we use are outlined in Fig. 1. We aim to explore the expected population of wide-separation exoplanets in order to gain insight into the formation and evolution histories of planetary systems and determine the feasibility of detecting perturbing planets in the near future. Attempts to constrain planets from debris disc features are common in the literature, but most target individual systems (e.g. Mouillet et al. 1997; Chiang et al. 2009; Backman et al. 2009; Read et al. 2018; Booth et al. 2021; Pearce et al. 2021); we consider a large sample of secure debris-disc stars and infer properties of the planet population in a consistent manner. Some of the planets that we predict should soon be detectable with the James Webb Space Telescope (JWST), which will finally provide observational tests of much of the interaction theory discussed over recent decades.

The paper layout is as follows. Section 2 describes the system sample, stellar parameters, planet-detection constraints, and debris-disc parameters that we use. Section 3 considers the planets required to sculpt our discs, and Sect. 4 determines whether planets are also needed to stir them. We discuss our results in Sect. 5, including comparisons of our predicted planet populations to both known planets and those inferred to be forming in protoplanetary discs, as well as identifying systems most favourable for future planet searches. We conclude in Sect. 6. The appendix provides a catalogue of our planet predictions and notes for specific systems.

Fig. 1 Cartoon showing the dynamical arguments used in this paper to constrain unseen planets from debris discs. At least one planet is assumed to have sculpted the inner edge of each disc (brown); this could be a single large planet (red circle) or multiple smaller planets (black circles). If the disc is eccentric, then we ascribe this to perturbations by a planet on an eccentric orbit (red line). The fact that we see dust implies that debris undergoes destructive collisions, which requires the disc to have been stirred; stirring may be performed by unseen planets (although other mechanisms are also possible). Planets can induce other debris morphologies in addition to those shown here, such as warps, resonant structures, and gaps; however, for this paper we only consider the above dynamical arguments.

2 System sample and debris disc data

To infer planet properties from debris discs, we first require information about the discs and their host stars. We also require planet-detection limits for later comparison with our predicted planets. Section 2.1 describes the system sample, stellar properties and planet-detection limits, and Sect. 2.2 the debris-disc data. We also quote 1σ uncertainties on these data; we will later use Gaussian uncertainty propagation to derive 1σ uncertainties on any values calculated from these data.

2.1 System sample, stellar parameters, and planet-detection limits

2.1.1 System sample

Our sample of debris disc stars is drawn from three large, high-contrast, L′-band imaging surveys: the NaCo Imaging Survey for Planets around Young stars (NaCo-ISPY; Launhardt et al. 2020), the Large Binocular Telescope Interferometer Exozodi Exoplanet Common Hunt (LEECH; Stone et al. 2018), and the L′-band Imaging Survey for Exoplanets in the North (LIStEN; Musso Barcucci et al. 2021). LEECH and LIStEN were conducted using LMIRCam on the Large Binocular Telescope Interferometer (LBTI), while ISPY used NaCo on the Very Large Telescope (VLT). The target selection criteria for these surveys are described in their respective papers, but can roughly be summarised as follows: targets were restricted to nearby stars including low-mass star-forming regions (distances <160 pc), exclude extreme spectral types (those later than M5 or earlier than B8), were observable by our telescopes (declinations between –73 and +70°), had good adaptive-optics performance (K-band magnitudes below 10), and were not too young (ages ≳10 Myr), focussing on pre-main-sequence stars up to several 100 Myr but also including main-sequence stars with gigayear-scale ages. For each star in these samples, we performed a spectral energy distribution (SED) fitting (Sect. 2.2.3) and selected only those stars with excess infrared emission that confirmed with high confidence the existence of at least one debris disc; stars with no significant excesses, or possible blending with background objects, were discarded. In total we selected 178 debris disc targets, of which 136 were in the ISPY survey, 33 in LEECH (including eight stars overlapping with ISPY), and 22 in LIStEN (including five overlapping with ISPY). These 178 stars are listed in Table A.1.

2.1.2 Stellar parameters

For each star, its distance d was derived from its Gaia Data Release 2 parallax according to the formalism described by Bailer-Jones et al. (2018). The effective temperature T_eff and luminosity L_* of the star, as well as the fractional luminosity and blackbody temperature T_bb of its associated debris disc(s), were derived by fitting simultaneous stellar (phoenix; Husser et al. 2013) and blackbody models to the observed photometry and spectra (see Sect. 2.2.3). The photometry was obtained from multiple catalogues and publications, including 2MASS, APASS, ffiPPARCOS/Tycho-2, Gaia, AKARI, WISE, IRAS, Spitzer, Herschel, JCMT, and ALMA. In some cases photometry was excluded, for example due to saturation or confusion with other objects. The fitting method used synthetic photometry of grids of models to find best fits with the MultiNest code (Feroz et al. 2009; Buchner et al. 2014; Yelverton et al. 2019). For uncertainties on effective temperature T_eff, we adopted a lower limit of 70 K (even if the formal fitting uncertainty was smaller); this lower limit was derived from the rms scatter of the difference between photometric and spectroscopic effective temperatures.

Stellar metallicities were taken from Gáspár et al. (2016) where available (133 stars); for the remaining stars we adopted [Fe/H] = -0.03 ± 0.19, which is the median and standard deviation of the 133 stars with metallicities. L′ magnitudes were derived by weighted interpolation of the 2MASS (Cutri et al. 2003) JHK and WISE (Cutri et al. 2013) W1–3 fluxes to the central wavelength of the L′ filter (3.80 µm), taking into account saturated flux measurements and upper limits.

To derive stellar ages, we first checked each target for membership of known associations using the banyan Σ tool¹ (Gagné et al. 2018). If the membership probability was >80%, then we assigned the stellar age t_* to be the mean age of the association². In total, we assigned association ages to 78 of our 178 stars. For our remaining 100 field stars, ages were assigned by compiling various literature estimates (particularly from large surveys with diverse age-determination methods). If these various estimates yielded contradictory ages, then we considered possible reasons before assigning an age (for example, the luminosity could be corrupted by unresolved binarity, or some methods may be inappropriate for the spectral type and age, or the literature uses an outdated association membership that appears unlikely given newer photometry, parallaxes, radial velocities and astrometry). If multiple valid age estimates were available, then we adopted an age (and conservative uncertainty range) that accounted for the spread in estimates.

To validate our adopted ages, and to derive stellar masses, we performed Hertzsprung-Russell diagram isochrone fits. We adopted evolutionary models from the MIST project³ (Dotter 2016; Choi et al. 2016; Paxton et al. 2011, 2013, 2015, 2018). Our method consisted of finding combinations of star mass M_*, age and metallicity that reproduced the measured effective temperature and luminosity, using the MIST models (despite having metallicity estimates, we still needed to sample metallicity to predict T_eff and L_* properly). The likelihood of the predicted T_eff and L_* given the data were modelled as a two-dimensional normal distribution, with standard deviations set equal to the 1σ uncertainties on T_eff and L_* (the two were assumed to be uncorrelated, based on inspection of the T_eff and L_* samples produced by the MultiNest fits). We assumed a Chabrier (2003) single-star prior on mass. For field stars, we assumed a 1/t_* prior on age (uniform in log t_*); for stars deemed to belong to associations, we substituted the 1 /t_* prior with a normal one based on the association mean age and its standard deviation. For metallicity, we assumed a normal prior based on the Gáspár et al. (2016) measurements. The posterior distribution of mass, age and metallicity was derived using the emcee sampler⁴. The resulting stellar masses (the medians of the posterior distributions), along with their 1σ uncertainties, were adopted for each star. The corresponding ages and uncertainties (which are naturally large for stars close to the main sequence) were not adopted; however, they were used to verify the consistency of our adopted ages with the isochrone fits⁵. Only for one star (HD 218340), for which we could not find any literature age, did we adopt the isochrone age directly. The corresponding metal-licities were also not used; instead we adopted the Gáspár et al. (2016) measurements, which were typically very similar.

Close binary stars were not excluded a priori, and were only avoided in certain circumstances (for example, where they would hamper coronagraphic observations). However, for each system we checked a posteriori for close binarity (separations < 1 arcsec) using our L′-band images, the 9th Catalogue of Spectroscopic Binary Orbits (SB9; Pourbaix et al. 2004), and the Washington Double Star Catalogue (WDS; Mason et al. 2020). Whilst these resources may be incomplete (in particular for very close companions with separations < 0.2 arcsec), they are the most complete repository to verify whether our calculated luminosity might be affected by previously unresolved sources, which would also affect our mass and age determination (see above). Moreover, physical binarity would affect the system dynamics and consequently our derived planet parameters, although additional data are often required to show physical association (e.g. Pearce et al. 2015). Targets where such close sources could affect our stellar parameters are flagged by asterisks in Table A.1, and some possible effects of binarity on our results are discussed in Sect. 5.3.3.

Our adopted distances, luminosities, effective temperatures, masses and ages are listed in Table A.1, along with their 1σ uncertainties and the references and associations considered in our age determination. For our 178 targets, the distances range from 3.2 to 160 pc (median 42 pc), luminosities from 4.7 × 10^-3 to 390 L_⊙ (median 3.6 L_⊙), stellar temperatures from 3000 to 12 000 K (median 6700 K), masses from 0.22 to 3.7 M_⊙ (median 1.4 M_⊙), and ages from 10 Myr to 6 Gyr (median 200 Myr).

2.1.3 Planet-detection limits

To get upper limits on the masses of undetected planets, we first extract contrast-detection-threshold curves from the L′-band images. For comparability with other surveys, we employ a classical 5σ detection threshold; as a function of angular separation, we determine the mean contrast required for a point-spread-function (PSF) shaped planet to exceed this 5σ threshold by injecting ‘fake planets’ into the images as described in Launhardt et al. (2020). We then determine the physical relation between this contrast-detection threshold and planet mass, given the star's distance, age and L′ magnitude, together with the AMES ‘hot start’ COND (no photospheric dust opacity; Baraffe et al. 2003) and DUSTY (maximal dust opacity; Chabrier et al. 2000) evolutionary model isochrones for the L′ passband. The mass-detection limits from the two models were then stitched together according to their validity range (COND: < 1400 K, DUSTY: > 1700 K) based on the best-fit COND effective temperatures, with a linear interpolation between the two validity ranges. Our resulting L′ contrast limits have median and standard deviation of 6.8 ± 1.3 mag at 0.15 arcsec and 9.6 ± 1.2 mag at 0.5 arcsec, with median and standard deviation mass-detection limits of ${40}_{-10}^{+130}$ M_Jup at 10 au and $8_{-4}^{+15}$ M_Jup at 100 au. These detection limits are not always as constraining as those from other instruments (e.g. SPHERE), but the strength of our data arises from the size and homogeneity of our sample of secure debris disc stars.

We also consider Gaussian uncertainty ranges on our mass-detection curves. For each image, the uncertainty on the azimuthally averaged 5σ contrast curve was taken to be the azimuthal variance of the contrasts. To derive 1σ uncertainties on the mass-detection limits, this contrast variance was combined with uncertainties on other parameters (propagated via Monte Carlo sampling); these uncertainty sources are the L′ magnitude of the star, its distance, and its age (the age in most cases being the dominating uncertainty factor). Throughout the paper we use the mass-detection limits calculated above, together with these associated 1σ uncertainties.

Several of our systems already have low-mass stellar or brown dwarf companions published, and those that could affect our debris disc dynamics are described in Appendix B. In addition, some potentially planet-mass companions are present in our images, but these are undergoing verification and will not be discussed in this paper.

2.2 Debris disc data

Constraining planet parameters from a debris disc requires knowledge of the debris morphology. Our models use the spatial distribution of large parent bodies (e.g. asteroids), which we infer from observations of dust thermal emission; we therefore derived our disc properties from ALMA, Herschel, or SED data. Thermal emission is used because, unlike scattered light, it is expected to represent the underlying parent-body morphology reasonably well; thermal emission traces larger grains that are less affected by non-gravitational forces, and its shallower radial dependence (∝r^−1/2 rather than ∝r⁻² for scattered light, where r is the stellocentric distance) makes it better at constraining the outer edges of broad discs. Hence we do not use scattered light data in our modelling. Whilst this means some additional features are neglected, such as the sweeping morphology of ‘The Moth’ (HD 61005; Hines et al. 2007) in scattered light, considering only thermal emission means that all of our disc morphologies (and resulting planet constraints) will be derived in a consistent manner.

For each system, debris disc parameters are taken from the best-available data source; ALMA data are preferred, then Herschel, and finally unresolved SED data. Only one debris model, derived from one data source, is used per disc. If multiple discs are present, then we only consider the outermost disc. For each disc we aim to predict planets located interior to the disc inner edge; to do this, the most important disc parameters are the locations and shapes of the disc inner and outer edges. We describe these edges using four parameters: the pericentre q_i and apocentre Q_i of an ellipse tracing the disc inner edge, and similarly the pericentre q_o and apocentre Q_o of an ellipse tracing the outer edge (see Fig. 2). These parameters are derived using different models according to whether the system has ALMA, Herschel or SED data available, as described in Sects. 2.2.1–2.2.3.

Fig. 2 ALMA Band 6 (1.3 mm) image of HD 202628 (Faramaz et al. 2019), as an example of the ALMA data used in this paper (Sect. 2.2.1). North is up, and east is left. The debris disc is well resolved and clearly eccentric, with the host star (crosshairs) offset from the geometric centre of the disc. Considering the pericentre, qi, and apocentre, Qi, of the disc inner edge, plus equivalent values qo and Qo for the outer edge, yields the minimum mass of a single planet required to sculpt and stir the disc within the system lifetime (Sects. 3.1 and 4.2). This minimum-mass planet would have the orbit shown by the grey line; alternatively, a higher-mass planet located farther inwards could also sculpt and stir the disc. It is unclear whether the bright point in the north-east part of the ring is a background object or associated with the system (Faramaz et al. 2019).

2.2.1 Systems with ALMA data

32 of our 178 systems have reliable ALMA disc models in the literature. For these systems, we take the disc parameters from literature ALMA fits, which are expected to be the most-reliable disc parameters currently available. These are also the only systems for which we consider asymmetric disc models, and these systems are expected to yield the tightest planet constraints of all systems in our sample. An example ALMA image that we use is shown in Fig. 2.

If several literature ALMA models are available for one disc, then we picked one of them. In choosing, we considered whether the model allows the disc to be asymmetric; some discs have well-constrained eccentricities, and we favour these models over axisymmetric fits because the former potentially allow planet eccentricities to be constrained. We also consider the data and fit quality, and whether system-specific considerations have been included. Systems with disc parameters derived from ALMA data are identified in Table A.1, and the literature sources for individual disc models are given in Appendix B. Some of our ALMA fits have not yet been published but come from the upcoming REsolved ALMA and SMA Observations of Nearby Stars (REASONS) analysis (Sepulveda et al. 2019; Matrà et al., in prep.); their spatial dust model follows Eq. (1) in Matrà et al. (2020), and we refer the reader to that paper for details of the fitting process for these systems.

Since we use ALMA models from a range of literature sources, the fitted disc edges are not defined in the same way across all ALMA systems. Some works, such as Lieman-Sifry et al. (2016), use a broad disc model with sharp edges (similar to our Herschel models in Sect. 2.2.2); in these cases we use their fitted inner and outer edges. Others, such as REASONS, model the dust spatial morphology as Gaussian; in this case we took the disc width as the fitted full width at half maximum and calculated the edges accordingly.

Fig. 3 Herschel PACS 100 µm image of HD 1461, as an example of the Herschel data used in this paper (Sect. 2.2.2). North is up, and east is left. From left to right, the plots show the observation, observation after PSF subtraction, best-fitting disc model, and residuals after model subtraction. Contours show ±2 and ±3 times the root-mean-square intensity. Whilst the disc is resolved, it clearly has a lower resolution than the example ALMA image in Fig. 2. Our Herschel disc models are axisymmetric; the offset between the disc model and the origin reflects the imperfect Herschel pointing rather than a physical offset between the disc and the star.

2.2.2 Systems with Herschel data but no ALMA data

Of the 146 systems without ALMA data, many still have resolved Herschel data from the Photodetecting Array Camera and Spectrometer (PACS), to which disc models can be fitted. The lower Herschel resolution makes these disc parameters less reliable than those in ALMA systems, but they are still favoured over those derived from unresolved SED data alone (Sect. 2.2.3). Whilst our PACS data have shorter wavelengths than ALMA (~ 100 µm for PACS vs. ~ 1 mm for ALMA), both instruments are expected to trace the cold populations of parent planetesimals (e.g. Krivov 2010; Pawellek et al. 2019). The main differences between disc parameters inferred from Herschel and ALMA data are expected to be dominated by differences in resolution, rather than different disc morphologies at different wavelengths; for this reason, we can compare planets inferred from ALMA discs to those from Herschel discs.

Using data from the Herschel Science Archive⁶, we identify every system in our sample with PACS 70 or 100 µm data. If level 2.0 and 2.5 data were available, we chose 2.5 (which combines scans to create a single image). For these systems, a parametric debris disc model was fitted to each individual 70 and 100 µm image. The model is similar to that in Yelverton et al. (2019), and consists of a star plus a single axisymmetric disc with sharp inner and outer edges. The locations of these edges were treated as free parameters, along with the disc orientation, disc flux and star position (to account for non-perfect Herschel pointing). The surface brightness profile was fixed as r^-1.5, and the discs were assumed to be optically and geometrically thin. We note that, unlike in Yelverton et al. (2019), no unresolved inner component was included (discussed later in this section).

The fitting procedure is described in Yelverton et al. (2019), and was conducted using publicly available code⁷. We refer the reader to that paper for a detailed description of the fitting process, but briefly it involved using a Markov chain Monte Carlo method to fit the disc-plus-star model, which had been convolved with a PSF estimated using a calibration star. The disc fitting was only performed if PSF subtraction yielded significant residuals around the target star, which could be indicative of a resolved debris disc. Since each available 70 and 100 µm dataset was considered separately, the analysis yields multiple Herschel models for some systems: one for each 70 µm image, and one for each 100 µm image.

Given the fitting results, we conducted a visual inspection to remove models where the fitting was likely confused by background objects. We also omitted models where the image-minus-model flux within 15arcsec of the image centre had a greater than 5% probability of having a non-Gaussian distribution in an Anderson-Darling test, which could be indicative of a poor fit. The surviving models were deemed reliable. For systems with multiple reliable Herschel models, we selected one fit; we chose 70 µm fits over 100 µm owing to better resolution, and if multiple models still remained (i.e. there were multiple observations at the same wavelength), then we chose the fit to the observation with the longest duration. This yielded 35 systems with reliable Herschel disc models (in addition to those resolved with ALMA), and the edges of these fitted discs, as well as the wavelengths of the PACS images used, are listed in Table A.1. An example Herschel image and disc model that we use is shown in Fig. 3.

The outer edges of the Herschel discs appear resolved and well constrained (bottom plot of Fig. 4), but the majority of the inner edges are unresolved and less reliable (top plot of Fig. 4). This will later lead to large uncertainties on planet parameters derived from Herschel data. The inner-edge determination is further compounded by the modelling not including an unresolved inner component; were separate inner discs also present, then we may underestimate the locations of the inner edges of the outer discs in the Herschel fits. This underestimation would not invalidate our later calculations of the minimum planet masses required to sculpt the inner edges, because the required masses increase with inner edge location (Sect. 3); underestimating the inner edges would mean the planets must be larger than the minimum masses we calculate, so these minimum masses still hold. However, underestimating the inner edges could mean that we overestimate the planet masses required to stir Herschel discs, as discussed in Sect. 5.3.2.

Fig. 4 Fitted debris disc parameters for systems where Herschel PACS data are used. The discs are modelled as axisymmetric with sharp inner and outer edges, as described in Sect. 2.2.2. The outer edges of these discs appear resolved (bottom plot), but the majority of the inner edges are unresolved (top plot).

2.2.3 Systems with only SED data

The majority of our systems (111 out of 178) are neither resolved by ALMA nor Herschel, and for these systems we infer disc properties from SED data alone. These data provide only rough disc locations, with no information on disc widths or degrees of asymmetry. The planet parameters later predicted for individual SED systems may therefore be unreliable; however, we will show that our predicted distributions of these sculpting planets across all SED systems are probably reasonable, since they are similar to those derived from more accurate ALMA and Herschel data. For systems with ALMA or Herschel data, we also use SED fits in the self-stirring calculations in Sect. 4.1.

The SED fitting process is described in Launhardt et al. (2020), and we refer the reader to that paper for details and data sources (see also Yelverton et al. 2019). For each system, the SED is fitted with a stellar model, plus one or two blackbody components representing one or two debris discs. An example is shown in Fig. 5. To estimate the radial location of the disc(s), the temperature of each blackbody component T_BB is first converted to a blackbody radius,

$$r_{BB}=1\text{au}{\left(\frac{T_{BB}}{278\text{K}}\right)}^{-2}{\left(\frac{L_{\ast }}{L_{\odot }}\right)}^{1/2}.$$(1)

However, this blackbody radius is often a poor estimate of the ‘true’ disc location, since it fails to account for grain optical properties and blowout size; in line with Booth et al. (2013) and Pawellek & Krivov (2015), we therefore use a correction factor Γ to better estimate the disc location as r_SED ≡ Γr_BB. The correction factor depends on stellar luminosity as $$\Gamma =A{\left(\frac{L_{\ast }}{L_{\odot }}\right)}^B,$$(2)

where Launhardt et al. (2020) used values for A and B derived from resolved Herschel data (A = 7.0, B = -0.39; Pawellek 2017), with Γ limited to 4 for stars with L_* < 4 L_⊙. However, this trend does not appear to hold for disc radii derived from better-resolved ALMA data (Matrà et al. 2018), and Pawellek et al. (2021) have since argued that a shallower function is better for discs around 1–20 L_⊙ stars: $$\Gamma =\left(2.92\pm 0.50\right)\times {\left(\frac{L_{\ast }}{L_{\odot }}\right)}^{-0.13\pm 0.07}.$$(3)

In Fig. 6, we plot these two Γ relations, and also show the ratios of ALMA- or Herschel-derived radii to blackbody radii for the systems in our sample.

The figure shows that discs around less luminous stars tend to be located exterior to where their blackbody temperatures imply, and the discrepancy is smaller for more luminous stars, in agreement with other studies (e.g. Booth et al. 2013; Pawellek & Krivov 2015; Marshall et al. 2021). Neither Γ correction factor provides a strong fit across our entire sample, but since the form in Eq. (3) (Pawellek et al. 2021) is derived from well-resolved ALMA data, we chose to use this Γ correction over the Herschel-derived factor from Pawellek (2017) (although we again limit Γ to a maximum of 4). It should be noted that, for our sample at least, the ALMA and Herschel data could actually be in agreement, suggesting a Γ steepness in line with Pawellek (2017) but with a smaller offset (e.g. Γ ~ 4(L_*/ L_⊙)^-039); however, we stress that we do not account for the various nuances involved in the estimation of Γ, and that this is an active area of research, so we take Γ from Eq. (3) rather than that roughly calculated from our sample. We also note that Γ is simply a trend used to estimate disc size; factors omitted here, such as disc width, may also be important, so we expect a reasonable scatter in Fig. 6 regardless of the Γ prescription used.

For each system without resolved ALMA or Herschel data, we estimate the location of the (outermost) disc by taking the blackbody radius and scaling it by Γ from Eq. (3) (where we assert Γ ≤ 4 as in Launhardt et al. 2020). This gives a rough location for the discs in our 111 SED systems. Since we have no width information, we are forced to set the inner and outer edges to be equal (i.e. Q_i = q_i = Q_o = q_o = r_SED), and we discuss the effect of this throughout the paper. Disc inclination is also unconstrained for SED systems.

Fig. 5 SED data for one example system, HD 37484 (Sect. 2.2.3). Black circles show photometric data, triangles are upper limits, and the grey shaded region is the Spitzer/InfraRed Spectrograph (IRS) spectrum. The best-fit model (black line) comprises the stellar flux (dotted green line) plus additional emission from two debris discs: a cold disc with blackbody temperature 60 ± 7 K (dashed blue line) and a warm disc with blackbody temperature 130 ± 10 K (dash-dot red line). This paper considers only the outermost (cold) discs.

Fig. 6 Ratios of resolved debris disc radii to blackbody radii for the systems in our sample, as described in Sect. 2.2.3. Red diamonds and yellow squares show systems where resolved disc radii were derived from ALMA and Herschel data, respectively. Each system is plotted once, using only the best available resolved data as described in Sect. 2.2. The lines show Γ correction factors; the solid line is derived from ALMA data (where we have limited Γ to a maximum of 4) and is that used in this paper (Pawellek et al. 2021), and the dashed line is from Herschel data (Pawellek 2017) as used by Launhardt et al. (2020). We note that the vertical axis uses the central radius of the resolved disc, defined as halfway between pericentre of the disc inner edge and the apocentre of the outer edge. The outlier on the left is the M-type Foma-lhaut C (NLTT 54872). The two systems in the bottom right of the plot may be complicated by either an underestimated cold disc location due to an additional warm component (β UMa: HD 95418; Moerchen et al. 2010) or background contamination (γ Oph: HD 161868; Holland et al. 2017).

3 Planet constraints from disc sculpting

Using the star and debris disc parameters from Sect. 2, we now use dynamical arguments to infer the properties of unseen planets that could be interacting with our discs. In this section we consider planets that could be sculpting these discs; later, in Sect. 4, we examine planets required to stir each disc. The discs in our sample are typically located at tens or hundreds of au, and it is generally assumed that planets exist interior to such discs and sculpt their inner edges. In this section we consider the planets required to sculpt the disc inner edges, assuming that the shapes and locations of these edges result from interactions with unseen planets.

For each system, we assume that the sculpting planet(s) orbits just interior to the disc inner edge, in the same plane as the disc, and that the planet orbits have not evolved significantly over the stellar lifetime; this simple model should provide reasonable planet constraints for the great majority of systems, although we note that more complicated interactions are also possible, such as those proposed for Fomalhaut (Faramaz et al. 2015; Pearce et al. 2021), HD 107146 (Pearce & Wyatt 2015; Yelverton & Kennedy 2018; Sefilian et al. 2021), or the Solar System (Gomes et al. 2005). However, these more complicated interactions often produce system-specific features, so to ensure generality across our sample we only consider basic interactions in this paper. We also again emphasise that this section concerns planets sculpting the inner edges of debris discs; for broad discs with radial structure (e.g. HD 15115, HD 92945, HD 107146; MacGregor et al. 2019; Marino et al. 2019, 2018) we are not constraining planets that potentially produce gaps in those broad discs, and care must therefore be taken when comparing our results to others to ensure the same planets are being considered.

If the disc inner edges are truncated by planets, then these planets must be sufficiently massive to sculpt the edges within the lifetime of the star. However, the specific planet masses required depend on the system architecture, history, and debris properties and morphology; many different mass estimation methods exist in the literature, each making different assumptions about the nature and degree of debris sculpting, and each predicting different masses for the planet(s) responsible (e.g. Wisdom 1980; Gladman & Duncan 1990; Faber & Quillen 2007; Mustill & Wyatt 2012; Pearce & Wyatt 2014; Morrison & Malhotra 2015; Nesvold & Kuchner 2015; Shannon et al. 2016; Lazzoni et al. 2018; Regály et al. 2018). In this analysis we use two different sets of assumptions to calculate the minimum mass of the sculpting planet(s), which should provide reasonable upper and lower bounds on the minimum planet masses required in each system. We first use a modified version of the Pearce & Wyatt (2014) method, which assumes that a single planet has ejected 95% of unstable debris from near the disc inner edge (Sect. 3.1). This method is chosen because it is valid for planets on either circular or eccentric orbits, with the latter favoured in some systems with asymmetric discs. We then employ an alternative method based on Faber & Quillen (2007) and Shannon et al. (2016), which assumes that multiple, equal-mass planets on circular orbits have removed 50% of unstable debris (Sect. 3.2). This second method is used because it yields conservative lower limits on planet masses, which are generally an order of magnitude below those predicted by the Pearce & Wyatt (2014) method (since a single, high-mass planet is less efficient at clearing than multiple, lower-mass planets, and it takes longer to remove 95% of unstable debris than 50%). Taken together, these two methods likely bound the minimum planet masses required to sculpt the disc inner edges, with the actual required masses depending on the specific system architecture and degree of clearing.

Throughout this section we refer to an example system, 49 Cet (HD 9672), to demonstrate the concepts used, before applying those methods to all systems in our sample. Figure 7 shows the various arguments used to constrain planets around 49 Cet, as described below. We note that the grain dynamics in this particular system may be complicated by the presence of gas (much less than in a protoplanetary disc, but still potentially enough to couple to dust; Moór et al. 2019). Nonetheless, this system should serve as an adequate example for our purposes.

Fig. 7 Masses and locations of the smallest planets needed to sculpt or stir the outer debris disc of 49 Cet (HD 9672), according to different models. The system hosts two debris discs; the outer disc has its resolved inner edge at 62 ± 4 au (line 1), and the unresolved inner disc is located around 12.2 ± 0.6 au (line 2). Lines 3a and 3b show 5σ detection limits; 3a is our NaCo-ISPY limit, and 3b is the SPHERE H2-band limit from Choquet et al. (2017) (with uncertainties derived as in Sect. 2.1.3). The red circle is the minimum-mass planet predicted by the Pearce & Wyatt (2014) model, which assumes that the outer disc is sculpted by a single planet (Sect. 3.1). Lines 4, 5, and 6 show the constraints from the Hill radius (Eq. (6)), diffusion time (Eq. (7)), and secular time (Eq. (8)), respectively, used to predict the properties of this single planet; if the outer disc is sculpted by one planet, then that planet is predicted to lie in the white region. The black circles denote the minimum-mass planets predicted by the alternative Shannon et al. (2016) model, which assumes that the outer disc is sculpted by multiple, equal-mass planets spanning the gap between the inner and outer discs (Eqs. (12) and (15) in Sect. 3.2). Line 7 is the minimum planet mass required to stir the outer disc (Eq. (23) in Sect. 4.2); we note that this is a separate analysis from the sculpting analysis because we do not necessarily assume that the disc is stirred by the sculpting planet(s). The 1σ uncertainties are shown as error bars around points and dark grey regions around lines. Similar analyses to those shown here are applied to all systems in our sample. We note that grain dynamics could be complicated by the presence of gas in this system (much less so than in a protoplanetary disc but still potentially enough to couple to dust; Moór et al. 2019). The PYTHON program used to produce this plot is publicly available for download8.

3.1 Disc sculpting by a single planet

We first considered a scenario where one planet sculpts the inner edge of each disc. We calculated the minimum mass of the planet required to do this, as well as its maximum semi-major axis and minimum eccentricity. The PYTHON program used for this analysis is publicly available for download⁸.

3.1.1 Disc sculpting by a single planet: Method

A single planet would destabilise non-resonant debris on either side of its orbit, and the larger the planet mass, semi-major axis, and eccentricity, the wider the unstable region. Pearce & Wyatt (2014) show that non-resonant material orbiting within five eccentric Hill radii of planet apocentre would eventually be ejected. Rearranging their Eqs. (9) and (10) yields the mass a single planet would need to sculpt the inner edge of an external disc, as a function of planet semi-major axis a_p and eccentricity e_p: $$M_{\text{p}}\approx 8.38M_{\text{Jup}}\left(\frac{M_{\ast }}{M_{\odot }}\right){\left[\frac{Q_i}{a_{\text{p}}\left(1+e_{\text{p}}\right)}-1\right]}^3\left(3-e_{\text{p}}\right),$$(4)

where Q_i is the apocentre of an ellipse tracing the disc inner edge. The apocentre is used because it equals the semi-major axis of the innermost disc particle if the disc eccentricity arises through a secular interaction with an internal planet (Pearce & Wyatt 2014); if the disc is axisymmetric, then Q_i is just the disc inner edge radius.

Equation (4) shows that the larger the planet eccentricity, the wider the unstable region. However, in addition to removing unstable debris, a coplanar eccentric planet would also drive surviving debris into an eccentric structure, apsidally aligned with the planet orbit (Wyatt et al. 1999; Faramaz et al. 2014; Pearce & Wyatt 2014); this means that a_p and e_p cannot be varied independently in Eq. (4) because an eccentric planet with apocentre close to the disc inner edge would also make the disc eccentric, which may not be compatible with observations. We can use this idea to eliminate e_p in Eq. (4), whilst retaining the dependence on planet eccentricity. If an internal eccentric planet sculpts the disc then, provided the planet eccentricity is not too high, the disc inner and outer edge dimensions are related to the perturbing planet orbit through $$a_{\text{p}}{e}_{\text{p}}\approx \frac{2}{5}\left(Q_i-q_i\right)\approx \frac{2}{5}\left(Q_o-q_o\right),$$(5)

recalling that q_i is the pericentre of an ellipse tracing the disc inner edge, and q_o and Q_o are the pericentre and apocentre of an ellipse tracing the disc outer edge, respectively (Eq. (5) is derived from Eqs. (1) and (5)–(8) in Pearce & Wyatt 2014). We note that the disc inner edge would be more elliptical than the outer edge, and for a broad, moderately eccentric disc, Eq. (5) can be used to assess whether that eccentricity could be driven by a single eccentric planet. Combining Eqs. (4) and (5) yields the general planet mass (as a function of semi-major axis) required to sculpt an external disc, where that disc may or may not be eccentric: $$M_{\text{p}}\approx 8.38M_{\text{Jup}}\left(\frac{M_{\ast }}{M_{\odot }}\right){\left[\frac{Q_i}{a_{\text{p}}+0.4\left(Q_i-q_i\right)}-1\right]}^3\times \left(3-0.4\frac{Q_i-q_i}{a_{\text{p}}}\right).$$(6)

This is line 4 in Fig. 7. The equation holds for both circular planets, and those with eccentricities up to moderate values⁹.

For the resolved discs in our sample, the vast majority do not show significant asymmetries. The deprojected edges of these discs do not deviate significantly from circles, and they have usually been fitted with axisymmetric models. This does not necessarily mean that these discs are axisymmetric; many could be asymmetric, but we lack the resolution and sensitivity necessary to detect this. Nonetheless, for this analysis any discs without significant asymmetries are treated as axisymmetric. In this case a circular planet could sculpt the disc inner edge; terms relating to planet eccentricity in Eqs. (5) and (6) disappear (i.e. Q_i = q_i and Q_o = q_o). For almost all resolved discs in this analysis, the eccentricity of the sculpting planet is therefore set to zero. Likewise, for unresolved discs we have no information on widths or asymmetries, so for these discs we also assume axisymmetry and non-eccentric sculpting planets (using Q_i = q_i = Q_o = q_o = r_SED in the above equations). However, there are several resolved discs that do display significant eccentricities, for which we can also constrain the eccentricity of the minimum-mass sculpting planet.

Equations (5) and (6) yield the sculpting planet mass and eccentricity as functions of its semi-major axis. However, this equation alone does not yield a lower bound on the planet mass required to sculpt a disc, which tends to zero as the planet apoc-entre approaches the disc inner edge. We can, however, combine this equation with a timescale argument to calculate the minimum planet mass required; if a planet sculpts a disc, then it must be massive enough to clear debris within the stellar lifetime. For an eccentric planet there are two timescales involved; the diffusion timescale sets how quickly debris close to the planet is removed via scattering, and the secular timescale determines how quickly more distant debris particles are driven onto eccentric, planet-crossing orbits. Pearce & Wyatt (2014) show that the time taken for a planet to clear 95% of unstable debris is 10 times the longer of the secular and diffusion timescales at the outer edge of the unstable region. Rearranging their equation 18 yields the planet mass required to clear the unstable region within the stellar lifetime t_* if the process is set by scattering: $$M_p\ge 0.331M_{\text{Jup}}\left(\frac{a_{\text{p}}}{\text{au}}\right){\left(\frac{Q_i}{\text{au}}\right)}^{-1/4}{\left(\frac{t_{\ast }}{\text{Myr}}\right)}^{-1/2}{\left(\frac{M_{\ast }}{M_{\odot }}\right)}^{3/4}.$$(7)

This is line 5 in Fig. 7. Similarly, if the clearing is limited by secular evolution, then rearranging Eq. (17) in Pearce & Wyatt (2014) yields the minimum planet mass required to clear the unstable region as $$M_p\ge 0.0419M_{\text{Jup}}{\left(\frac{a_{\text{p}}}{\text{au}}\right)}^{-1}{\left(\frac{Q_i}{\text{au}}\right)}^{5/2}{\left(\frac{t_{\ast }}{\text{Myr}}\right)}^{-1}\times {\left(\frac{M_{\ast }}{M_{\odot }}\right)}^{1/2}{\left[b_{3/2}^{\left(1\right)}\left(\frac{a_{\text{p}}}{Q_i}\right)\right]}^{-1},$$(8)

where $$b_s^{\left(j\right)}\left(\alpha \right)\equiv \frac{1}{\pi }{\displaystyle {\int }_0^{2\pi }\frac{\cos \left(j\psi \right)}{{\left(1-2\alpha \cos \psi +{\alpha }^2\right)}^s}\text{d}\psi }$$(9)

is a Laplace coefficient (Murray & Dermott 1999). Equation (8) is line 6 in Fig. 7. This secular timescale can be important for highly eccentric planets, but is less important for low-eccentricity planets because these do not significantly increase debris eccentricities through secular interactions. Furthermore, for planetary-mass objects the secular timescale is shorter than the diffusion timescale unless the planet semi-major axis is much smaller than the disc inner edge location (Eqs. (17) and (18) in Pearce & Wyatt 2014), and in this case the planet is unlikely to sculpt the disc unless the former has very high eccentricity. For the systems considered in this analysis we will find that the required eccentricities of the sculpting planets are low even for asymmetric discs (Sect. 3.1.2), and so the clearing times are set by the diffusion timescale rather than the secular timescale for all of our systems.

Equations (7) and (8) put lower bounds on the sculpting planet mass as functions of planet semi-major axis, and Eq. (6) relates planet mass and semi-major axis. We can therefore find the minimum planet mass required to sculpt a disc by equating Eq. (6) and the larger of Eqs. (7) and (8) (for our sample, the larger is always Eq. (7)). This is done numerically. This also yields the semi-major axis of the minimum-mass perturbing planet, which is the maximum semi-major axis that a perturbing planet could have. These minimum mass and maximum semi-major axis predictions are shown for 49 Cet (HD 9672) as the red circle in Fig. 7. Finally, Eq. (5) gives the eccentricity of the minimum-mass planet (in cases where significant disc asymmetry is observed), which is the minimum possible eccentricity of a perturbing planet; in general, since we model the planet peri-centre as being internal to the pericentre of the disc inner edge, Eq. (5) also shows that the eccentricity of a sculpting planet must be at least $$e_p\gtrsim \frac{4e_{\text{i}}}{5-e_{\text{i}}},$$(10)

where e_i is the eccentricity of an ellipse tracing the disc inner edge.

3.1.2 Disc sculpting by a single planet: Results

We apply the above method to every system in our sample, to constrain the minimum masses of planets required to sculpt the inner edges of our discs (assuming that one planet is responsible in each case), as well as the maximum semi-major axes and minimum eccentricities of those planets. The results for each system are listed in Table A.2. The predicted minimum masses are shown on the top plot of Fig. 8, with a median of 0.4 M_Jup and first and third quartiles of 0.2 and 0.8 M_Jup, respectively. These values were calculated from our whole sample, but the distribution is similar if we instead only consider more reliable constraints from systems where the discs are resolved; for systems with ALMA- or Herschel-resolved discs, the minimum sculpting planet masses have a median of 0.3 M_Jup and first and third quartiles of 0.1 and 0.9 M_Jup, respectively. This suggests that the ‘typical’ individual planet required to sculpt a debris disc is around half a Jovian mass; such bodies are unlikely to be detectable in the outer regions of systems with current instruments, but could be in the near future (as discussed below).

The left plot of Fig. 9 shows the minimum masses and maximum semi-major axes of the sculpting planets. The two are correlated because the more distant discs require more distant planets to sculpt them, and these sculpting planets would need higher masses because the clearance timescale increases with stellocentric distance (line 5 in Fig. 7). For each system, comparing the predicted semi-major axis of the minimum-mass sculpting planet and the location of the disc inner edge yields a median ratio of 81% (first and third quartiles: 78 and 84%, respectively), which is comparable to Neptune and the Kuiper Belt (although care is required with this comparison due to the complicated history of our Solar System; Levison et al. 2008). We note that we have only considered minimum-mass planets here; higher-mass planets located farther inwards could also sculpt the discs to the same degree (line 4 in Fig. 7). For the few discs with significant asymmetries in thermal emission, the sculpting planets typically require minimum eccentricities of around 0.1. These eccentricities are listed in Table A.2.

Figure 10 shows the detectability of our predicted sculpting planets, given our survey sensitivities. The figure shows the minimum sculpting-planet masses, divided by our L′-band AMES-COND/DUSTY mass-detection limits at the apocentres of those planets (if we have multiple contrast curves, for example if the target is in more than one of our surveys, then we use the most constraining detection limit at that location). Any planet lying above the horizontal line marking unity could therefore be detectable. The left plot shows our predictions if each disc is sculpted by a single planet (the right plot, for multi-planet sculpting, will be discussed in Sect. 3.2.2). We see that, to 1σ, all predicted planets on the left plot are consistent with being below the detection limits; this means that our discs could each be sculpted by a single, unseen planet. However, some of these planets would be at the limits of detectability; for 23 systems, the predicted planet masses are greater than 10% of our detection limits (to 1σ), including eight systems with reliable ALMA or Herschel disc data. These planets would be detectable following a factor of 10 improvement in mass-detection limits. The planets in a further 23 systems could be above 10% of our detection limits too, but their 1σ uncertainties mean they could also be below 10% of these limits. The conclusion is that, if a single planet sculpts each disc, then many such planets could be observable in the near future; we discuss planet-detection possibilities further in Sect. 5.2. However, we also note that the sculpting planets could be less massive than those predicted here if sculpting is not performed by a single planet, but instead by multiple planets; that scenario is considered in the following section.

Fig. 8 Minimum masses of the planet(s) predicted to exist in our debris disc systems, for three different scenarios: a single planet sculpts the inner edge of the disc (top plot), multiple planets sculpt the inner edge of the disc (middle), or a single planet stirs the disc (bottom). The solid bars show predicted planet parameters in all of our systems, and the hatched bars show only systems where the disc parameters come from resolved images (ALMA or Herschel); the results for resolved and unresolved discs are similar for sculpting planets (top two plots) but differ for stirring planets because the results for unresolved discs are conservative lower bounds (Sect. 4.2.1). For the multi-planet scenario (middle plot), we note that the plot shows the mass of each individual planet.

3.2 Disc sculpting by multiple planets

The previous section considered a scenario where one planet sculpts each debris disc, and identified the minimum-mass planet required in each case. However, multiple planets can clear debris more efficiently than a single planet could; if disc sculpting is instead performed by multiple planets, then these planets could be less massive than those predicted in Sect. 3.1.

In this section, we aim to identify the smallest planets that could possibly sculpt each disc, by considering a somewhat-idealised scenario where multiple, equal-mass planets remove unstable debris. To this end we also relax our definition of a ‘sculpted’ disc; here we consider the planets required to remove just 50% of unstable material, rather than 95% as in Sect. 3.1. This lowers the required planet masses still further, although the mass reduction is predominantly caused by the use of multiple planets as argued later in this section. The planet masses predicted in this section are the very minimum possible¹⁰, and it is likely that the actual sculpting planets have masses more similar to those calculated in Sect. 3.1. Nonetheless, calculating such extreme lower bounds is useful; it demonstrates that, even if a single planet predicted to sculpt a disc can be ruled out by observations, the disc could still be sculpted by multiple, lower-mass planets that evade detection.

3.2.1 Disc sculpting by multiple planets: Method

We explore the scenario of Faber & Quillen (2007) and Shannon et al. (2016), where multiple, equal-mass planets on circular orbits span a region interior to the disc, and clear that region of debris (the black circles in Fig. 7 show such planets for 49 Cet as an example). The planet orbits are assumed to be spaced such that each adjacent pair of orbits are separated by K mutual Hill radii, where K is a constant and the mutual Hill radius is $$R_{\text{H}}=\frac{a_1+a_2}{2}{\left(\frac{M_1+M_2}{3M_{\ast }}\right)}^{1/3};$$(11)

here the two planets have masses M₁ and M₂, and semi-major axes a₁ and a₂. If each planet has the same mass M_p,n, and the innermost and outermost planets have orbital radii r₁ and r₂, respectively, then Eq. (11) can be used to show that the number of planets between r₁ and r₂ is $$n_{\text{p}}=1+\frac{\log \left(r_2/r_1\right)}{\log \left[\frac{1+\eta \left(M_{\text{p},\text{n}}\right)}{1-\eta \left(M_{\text{p},\text{n}}\right)}\right]},$$(12)

where $$\eta \left(M_{\text{p,n}}\right)\equiv 0.0430K{\left(\frac{M_{\text{p,n}}}{M_{\text{Jup}}}\right)}^{1/3}{\left(\frac{M_{\ast }}{M_{\odot }}\right)}^{-1/3}.$$(13)

Equations (12) and (13) are equivalent to Eq. (5) in Shannon et al. (2016), only ours are written in terms of Jupiter masses, and with K (the number of mutual Hill radii between each pair of adjacent planet orbits) explicitly stated. We note that the logarithms in Eq. (12) can have any matching base.

Equations (12) and (13) yield the number of planets of a given mass that can fit between radii r₁ and r₂, but we still need to establish how massive those planets would have to be to clear that region of debris. As in Sect. 3.1, the planets must be massive enough to clear unstable debris within the stellar lifetime. Shannon et al. (2016) ran numerical simulations of this setup, and found that the time taken to clear 50% of debris depends on planet mass, star mass, outer radius r₂ and orbit spacing constant K: $$t_{\text{clear}}\approx \tau \left(K\right){\left(\frac{r_2}{\text{au}}\right)}^{3/2}{\left(\frac{M_{\text{p,n}}}{M_{\oplus }}\right)}^{-1}{\left(\frac{M_{\ast }}{M_{\odot }}\right)}^{1/2},$$(14)

where τ(K) ≈ 4 Myr if K = 20, or 2 Myr if K = 16. For this analysis we assume that each adjacent planet pair is separated by K = 16 mutual Hill radii, for several reasons; first, Eq. (12) fails if K³ ≥ 12 M_*/M_p,n, and we found that some of our systems encounter this problem if K = 20, but not if K = 16. Second, we are interested in finding the smallest planets that could possibly sculpt each disc, and taking K = 16 rather than 20 reduces the planet masses required by a factor of 2¹¹. We therefore estimate the minimum mass of a planet required to clear debris in a multi-planet system as $$M_{\text{p,n}}\ge 6.29\times {10}^{-3}{M}_{\text{Jup}}{\left(\frac{t_{\ast }}{\text{Myr}}\right)}^{-1}{\left(\frac{r_2}{\text{au}}\right)}^{3/2}{\left(\frac{M_{\ast }}{M_{\odot }}\right)}^{1/2}.$$(15)

This is equivalent to Eq. (4) in Shannon et al. (2016), only ours is written in terms of Jupiter masses, and we use K = 16 rather than 20.

Equation (15) shows that the minimum planet mass required to sculpt a disc depends on the orbital radius of the outermost planet, which Shannon et al. (2016) argue is roughly the inner edge of the disc. Whilst such a planet would also clear debris beyond its orbit, and would actually be located interior to the disc edge as in Sect. 3.1, for this analysis we follow Shannon et al. (2016) and assume that the outermost planet resides at the inner edge of our disc. The effect of this is that clearing may occur slightly faster in reality than assumed here, but the discrepancy is expected to be small. In the few cases where the disc is asymmetric, we approximate the outermost planet orbit to be circular with radius equal to the pericentre of the disc inner edge because this provides a lower limit for the required planet mass; whilst the assumption of circular planet orbits may not hold for asymmetric discs, the low degrees of asymmetry in our discs make this a reasonable approximation. We therefore use r₂ = q_i.

Shannon et al. (2016) designed this model for two-disc systems, where planets clear a gap between a warm and a cold disc. They therefore argue that r1 in Eq. (12) can be set to the location of the inner disc (r₁ = r_warm), allowing the calculation of the number of planets in the gap. However, this number should be treated with caution because it depends very strongly on model assumptions (equal-mass planets extending all the way in to the warm disc). Owing to this, plus the fact that the majority of our systems host only one detected disc, we chose not to calculate the total number of planets expected from the Shannon et al. (2016) multi-planet model. Despite this, the minimum masses of these planets are expected to be robust; the clearing timescale is set only by the outermost planet, so we can still calculate the masses required if multiple, equal-mass planets separated by K = 16 mutual Hill radii sculpt the disc, without knowing the actual number of such planets (provided it is greater than one). We calculate these masses for each system in the following section (we note that, for illustrative purposes only, we did calculate the number of planets between the inner and outer disc of 49 Cet in Fig. 7, using r_warm = 12.2 ± 0.6 au).

Fig. 9 Masses and locations of the smallest planets required to sculpt our debris disc inner edges, assuming the planets lie interior to those discs (Sect. 3). Left plot: minimum planet masses required if each disc is sculpted by a single planet (Pearce & Wyatt 2014). More massive planets located farther inwards could also sculpt the discs to the same degree. Right plot: minimum individual planet masses required if each disc is instead sculpted by multiple, equal-mass planets (Shannon et al. 2016), where the outermost planet location is fixed at the disc inner edge. The single-planet model (left plot) requires larger planet masses to sculpt discs than the multi-planet model does (right plot), mainly because multiple planets are more efficient at removing debris. Symbols denote the disc data sources, where the best-available source has been used for each system; systems with ALMA data (red diamonds) are expected to yield the most reliable planet constraints, followed by those with Herschel data (yellow squares), then finally those with only unresolved SED data (black circles). The larger uncertainties on the multi-planet masses are driven by uncertainties on the disc inner edge locations; the multi-planet model has a stronger dependence on this than the single-planet model does.

Fig. 10 Detectability of the predicted disc-sculpting planets (Sect. 3). Horizontal axes show the minimum planet masses required to sculpt the discs if debris clearing is performed by one planet (left plot) or by multiple, equal-mass planets (right plot). The vertical axis shows sculpting planet masses divided by our L′-band detection limits at the apocentres of those planets (assuming the best-case scenario, where the systems are face-on to the observer). The dashed lines show unity. Symbols have the same meanings as in previous figures. If the discs are each sculpted by one planet (left plot), then some of these planets will lie at or close to the current detection limits, and a factor of 10 improvement in mass detection limits could yield many more predicted planets. Conversely, if instead each disc is sculpted by multiple planets (right plot), then these planets could be significantly smaller and many might not be observable even with significant improvements in detection limits.

3.2.2 Disc sculpting by multiple planets: Results

We now apply the above method to every system in our sample, to constrain the minimum masses of planets required to sculpt the inner edges of discs in the most conservative scenario (sculpting by multiple, equal-mass planets on circular orbits). The results for each system are listed in Table A.2. The minimum masses of the individual planets in the multi-planet scenario are shown on the middle plot of Fig. 8, with a median of 0.01 M_Jup (3 M_⊕) and first and third quartiles of 2 × 10^-3 M_jup (0.8 M_⊕) and 0.05 M_jup (16 M_⊕), respectively. As in the single-planet case, the distribution is similar if instead we only consider the more reliable constraints from systems where the discs are resolved; for those systems the minimum masses have a median of 9 × 10^-3 M_Jup (3 M_⊕) with first and third quartiles of 8 × 10^-4 M_jup (0.3 M_⊕) and 0.05 M_jup (17 M_⊕), respectively.

The planets required to sculpt each disc in the multi-planet scenario are on average 40 times less massive than a single planet would have to be to sculpt the same disc (first and third quartiles: 20 and 80 times less massive). This difference is predominantly caused by multiple planets being more efficient at clearing debris than a single planet is; the requirement that only 50% of unstable debris is removed in the multi-planet model (rather than 95% in the single-planet model) cannot account for the difference, since Yabushita (1980) show that the number of particles surviving repeated close encounters with a single planet decreases over time as t^-1 or t^-2, and propagating this through Eq. (7) shows that the required mass of a single planet would only decrease by a factor of 2–3 due to this effect. Similarly, Eq. (14) shows that the multi-planet masses required for sculpting depend strongly on the inter-planetary spacing; a modest reduction in the number of mutual Hill radii between planets from K = 20 to 16 lowers the required planet masses by a factor of 2, implying that the presence of multiple planets strongly influences the interaction timescales. Nonetheless, the multi-planet model represents a somewhat-idealised scenario, so the minimum planet masses found in this section are conservative lower bounds. In reality, it is likely that the smallest planets capable of sculpting a 'typical’ cold debris disc lie between the super-Earths predicted in this section and the Jupiters predicted by the single-planet model (Sect. 3.1); the former may remain undetectable in the outer regions of systems for the near future, but the latter could be observable soon.

The multi-planet model assumes that the outermost planet is located at the inner edge of the disc, and we show the predicted semi-major axes of these planets on the right plot of Fig. 9; as in the single-planet case, the planet masses required to sculpt a disc increase with the disc inner edge radius. Given our L′-band detection limits at the inner edges, the right plot of Fig. 10 shows the detectability of these planets; the minimum planet masses required in the multi-planet-sculpting scenario all lie below our detection limits, so all of our discs could be sculpted by multiple unseen planets. In fact, for the majority of systems such planets could lie very far below current detection limits; if multiple planets sculpt these discs, then only one system (HD 126062) requires the outermost planet mass to be greater than 10% of our detection limit (to 1σ), and since no ALMA or Herschel data are available for this system these specific constraints are not particularly reliable. A further three systems also have minimum outer-planet masses compatible with 10% of the detection limits (to 1σ), although their uncertainties mean they could also be below 10% of these limits. Of these, just one (CPD-72 2713) has resolved disc data. Therefore, in a worst-case scenario, even a tenfold improvement in mass detection limits would not necessarily allow the observation of significantly more sculpting planets. However, it must again be recalled that the multi-planet model produces a conservative lower bound on planet mass, and the (less idealised) single planet model suggests that a tenfold improvement in mass detection limits could yield a substantial number of planets.

To summarise, in Sect. 3 we considered the planets required to sculpt the inner edges of our debris discs. We showed that, if single planets sculpt each disc, then many of those planets could be detectable with upcoming instrumentation. Conversely, if each disc is instead sculpted by multiple planets, then it is possible for those planets to evade detection even following a tenfold improvement in mass detection limits. Future planet searches should bear this in mind; even if observations can rule out an individual planet required to sculpt a disc, that does not mean planet sculpting is not occurring, since that disc could still be sculpted by multiple, lower-mass planets lying below the detection limits. However, this scenario is a somewhat-idealised case, where all planets have equal, minimal masses and are optimally spaced to facilitate clearing. Even if multiple planets are present, then (by analogy with Jupiter and Neptune in the Solar System) we suspect that clearing would often still be dominated by a single planet, and those planets would be more akin to those predicted in Sect. 3.1 (Neptune- to Jupiter-mass planets, as opposed to the super-Earths predicted in this section).

4 Planet constraints from debris stirring

The previous section discussed the planets required to sculpt the inner edges of debris discs. We now use a different physical argument to constrain unseen planets: the requirement that the discs have been somehow stirred, such that collisions between large bodies are sufficiently common and violent to produce the dust that we see. We first show that disc self-gravity alone is insufficient to bring about these collisions in a number of systems, unless the disc masses were unfeasibly high (Sect. 4.1). This could imply that planets perform the stirring instead, and in Sect. 4.2 we show that unseen planets would indeed be capable of stirring our discs.

4.1 Self-stirring by debris discs

For a debris disc to be self-stirred, its mass must be sufficient to drive debris onto crossing orbits and induce destructive collisions at its outer edge within the system lifetime (Kenyon & Bromley 2001, 2008, 2010; Kennedy & Wyatt 2010; Krivov & Booth 2018). In this model the largest bodies in the disc, those with radii s_max, perturb and induce destructive collisions between smaller bodies. However, the size of these largest bodies is unclear. We conduct a new analysis, combining the results of Krivov & Booth (2018) and Krivov & Wyatt (2021), to find the minimum-possible s_max required to stir a disc; in turn, this provides a self-consistent lower bound on the minimum disc mass required for self-stirring. The python program used for this analysis is publicly available for download⁸.

4.1.1 Self-stirring by debris discs: Method

Rearranging Eq. (28) of Krivov & Booth (2018) yields a lower limit on the disc mass required to initiate self-stirring within the system lifetime: $$M_{\text{disc,stir}}\gtrsim 1.62\times {10}^{-3}{M}_{\oplus }{\left(\frac{t_{\ast }}{\text{Myr}}\right)}^{-1}{\left(\frac{M_{\ast }}{M_{\odot }}\right)}^{-1/2}\times \left(\frac{a_2-a_1}{\text{au}}\right){\left(\frac{a_1+a_2}{\text{au}}\right)}^{5/2}\times {\gamma }^{-1}{\left(\frac{{\rho }_{\text{s}}}{1\text{g}{\text{cm}}^{-3}}\right)}^{-1}{\left[\frac{{\upsilon }_{\text{frag}}\left(s_{\text{weak}}\right)}{m{s}^{-1}}\right]}^4{\left(\frac{s_{\max }}{\text{km}}\right)}^{-3}.$$(16)

Here a₁ and a₂ are the innermost and outermost semi-major axes of the debris disc particles, respectively, γ parametrises the eccentricity of the stirring bodies (1 ≤ γ ≤ 2), ρ_s is the planetes-imal material density, and v_frag the collision speed required to fragment the weakest colliding body (which has radius s_weak). Equation (16) strongly depends on v_frag and s_max; Krivov & Booth (2018) used 30 m s^-1 and 200 km, respectively, as examples, but the actual values are unclear. We now reconsider these values, in order to place self-consistent lower bounds on the minimum disc masses required for self-stirring.

To proceed, we must first consider the mass of a debris disc. Defining n(s) as the number of particles of radius s in the range s → s + ds, the debris disc size distribution is often modelled as a single power law n(s) ∝ s^−q, with q ≈ 3.5 expected from destructive collisions (Dohnanyi 1969). However, Krivov & Wyatt (2021) argue that debris is unlikely to follow a single power law all the way from dust up to the largest planetesi-mals, and that a three-slope model is more realistic; in this case there exist transitional sizes s_mm and s_km, where the size distribution slope changes due to different physics affecting debris. For grains smaller than s_mm ~ 1 mm, radiation effects are important and the size distribution is defined to be n(s) ∝ s−^q. Larger bodies are unaffected by radiation forces, so debris larger than s_mm but smaller than the largest colliding body s_km have their size distribution set by destructive collisions, with n(s) ∝ s^−qmed. Finally, bodies larger than s_km are not yet involved in this collisional cascade; it takes time to stir large bodies such that they undergo destructive collisions, so bodies larger than (time-dependent) s_km are primordial and not yet colliding. Hence the size distribution between the largest fragmenting body s_km and the largest body s_max goes as $n\left(s\right)\propto s^{-q_{\text{big}}}$. As in Krivov & Wyatt (2021), we fixed q = 3.5, q_med = 3.7, and q_big = 2.8 ± 0.1, noting that these values describe indices and should not be confused with the disc peri-centre locations q_i and q_o. Integrating over the size distribution yields the total disc mass in the three-slope model as $$M_{\text{disc}}=M_{\text{dust}}\frac{4-q}{4-q_{\text{big}}}{\left(\frac{s_{\text{km}}}{s_{\text{mm}}}\right)}^{4-q_{\text{med}}}{\left(\frac{s_{\max }}{s_{\text{km}}}\right)}^{4-q_{\text{big}}},$$(17)

where M_dust is the mass in grains smaller than s_mm. This can be derived from observations: $$M_{\text{dust}}\approx \frac{F_v{d}^2}{\kappa \left(\lambda \right)B_v\left(T_{\text{BB}}\right)},$$(18)

where F_v is the disc specific flux per unit frequency range at observation wavelength λ, κ(λ) ≡ 1.7 cm² g^-1 (850 µm/λ) is the assumed dust opacity (noting that this involves various factors and is subject to uncertainty; e.g. Krivov & Wyatt 2021), and B_v(T_BB) is the blackbody emission intensity at v for dust with temperature T_BB, where T_BB is found from the SED disc fits (Sect. 2.2.3).

The total disc mass (Eq. (17)) therefore depends on the unknown sizes of both the largest body (s_max) and the largest colliding body (s_km). Krivov & Wyatt (2021) show that s_km is set by s_max; combining their Eqs. (17) and (19) yields $${\left(\frac{s_{\text{km}}}{\text{km}}\right)}^{q_{\text{big}}-q_{\text{med}}}=2.66\times {10}^{-3}{\left(\frac{t_{\ast }}{\text{Myr}}\right)}^{-0.59}\times {\left[\frac{{M^{\prime }}_{\text{disc}}\left(s_{\max }\right)}{M_{\oplus }}\right]}^{-0.59}{\left(\frac{M_{\ast }}{M_{\odot }}\right)}^{-0.79}{\left(\frac{a_1+a_2}{\text{au}}\right)}^{2.56}$$(19)

where ${M^{\prime }}_{\text{disc}}\left(s_{\max }\right)$ is the mass that would arise from Eq. (7) if s_km were set to 1 km. So if we knew the size of the largest body s_max, we could uniquely calculate the corresponding disc mass from Eq. (17).

We now reconsider the disc mass required for self-stirring (Eq. (16)), which depends strongly on s_max and the fragmentation speed of the weakest colliding body, v_frag(s_weak). These are unknown, but by assuming the aforementioned size distribution from Krivov & Wyatt (2021), we can calculate v_frag(s_weak) as a function of s_max and therefore express Eq. (16) as a function of a single unknown, s_max. Equation (1) in Krivov et al. (2018) gives the material strength $Q_{\text{D}}^{\ast }$ of a colliding particle as $$Q_{\text{D}}^{\ast }\left(s,v_{\text{col}}\right)=9.13\text{J}{\text{kg}}^{-1}{\left(\frac{v_{\text{col}}}{\text{m}{\text{s}}^{-1}}\right)}^{1/2}\times \left[{\left(\frac{s}{m}\right)}^{-0.36}+{\left(\frac{s}{\text{km}}\right)}^{1.4}\right],$$(20)

which depends on particle size s and collision speed v_col, approximating the material as basalt. This is a v-shaped function of particle size, where the weakest-possible colliding particle would have radius 110 m. Since the largest colliding particle has radius s_km, the radius of the weakest colliding particle s_weak will therefore be the smaller of s_km and 110 m. We can calculate the fragmentation speed of this particle; since destructive collisions would occur between equal-sized colliders at the onset of the collisional cascade, we can take $v_{\text{frag}}\left(s\right)=\sqrt{8Q_{\text{D}}^{\ast }\left[s,v_{\text{frag}}\left(s\right)\right]}$ from Krivov et al. (2005), which yields $$v_{\text{frag}}\left(s_{\text{weak}}\right)=17.5\text{m}{\text{s}}^{-1}{\left[{\left(\frac{s_{\text{weak}}}{\text{m}}\right)}^{-0.36}+{\left(\frac{s_{\text{weak}}}{\text{km}}\right)}^{1.4}\right]}^{2/3}.$$(21)

So for a given s_max, we can determine the size of the largest colliding particle s_km, which then yields the size of the weakest colliding particle s_weak, which in turn yields the fragmentation speed v_frag(s_weak) of that particle¹².

Combining the above arguments, for a given s_max we can uniquely solve the disc mass (Eq. (17)) and determine whether that mass is sufficient to stir the disc (Eq. (16)). Furthermore, inspection of Eqs. (16)–(21) show that the disc mass M_disc monotonically increases with s_max, while the mass required for stirring, M_discstir, monotonically decreases with s_max. This means that there will be a single critical value of s_max where the disc mass is just sufficient for stirring; for any higher s_max the disc can be self-stirred, and for lower s_max it cannot. In turn, this critical s_max corresponds to the minimum mass that a disc requires to self-stir.

Using Eqs. (16)–(21), we can therefore compute the minimum masses that our discs would require to be self-stirred. We do this by equating Eqs. (16) and (17), and numerically solving the result for the critical s_max required for self stirring (recalling that s_km and v_frag(s_weak) are also functions of s_max). We note that we are considering the most optimistic case here, where the largest planetesimals are only just massive enough to stir the disc; even so, we will show that some of these required masses could still be unfeasibly large. Such discs would need to be stirred by other mechanisms, such as by unseen planets (Sect. 4.2).

For axisymmetric discs we set a₁ and a₂ in Eq. (16) to the inner and outer disc edges, respectively. If the disc is eccentric then we assume that the eccentricity is caused by a secular interaction with an eccentric planet; in this case the innermost and outermost debris semi-major axes will equal the apocentre of the disc inner edge and the pericentre of the disc outer edge, respectively (Eqs. (5) and (7) in Pearce & Wyatt 2014). We therefore generally set a₁ = Q_i and a₂ = q_o in Eq. (16). Also, as in Musso Barcucci et al. (2021), we set γ ≡ 1.5 and use ρ_s = 1 g cm^-3 (Krivov & Booth 2018). Finally, since the required disc mass is proportional to the disc width (a₂ - a₁ in Eq. (16)), and no width information is available from disc SEDs, we omit systems with only SED data from this analysis.

4.1.2 Self-stirring by debris discs: Results

Across our resolved discs, we find that the median s_max required for stirring is 14 km (first and third quartiles: 12 and 18 km). This means that, for a typical disc to be self-stirred, it must contain bodies larger than 14 km; for these values of s_max, the largest colliding bodies would have median size s_km = 55 m (first and third quartiles: 31 and 94 m), and the median fragmentation speeds of the weakest colliding bodies v_frag(s_weak) = 7.0 ms^-1 (first and third quartiles: 6.6 and 7.8 ms^-1). We note that these minimum s_max and v_frag(s_weak) values are much smaller than the values of 200 km and 30 ms^-1 respectively used by Krivov & Booth (2018) and that their values are not necessarily consistent with each other; for example, if the largest planetesi-mal in the disc of 49 Cet (HD 9672) really had s_max = 200 km, then the fragmentation velocity of the weakest colliding body v_frag(s_weak) would be just 7 m s^-1. This difference in u_frag(s_weak) is very important due to the $v_{\text{frag}}^4$ term in Eq. (16). In general, we therefore argue that Eq. (28) in Krivov & Booth (2018) (Eq. (16) here) should not be used with ‘standard’ values of s_max = 200 km and v_frag(s_weak) = 30 m s^-1 to calculate the minimum disc mass required to self-stir, but that these values should instead be calculated on a system-specific basis using a method akin to that in this paper.

If our resolved discs are self-stirred, then their minimum possible masses have a median of 110 M_⊕ (first and third quartiles: 12 and 680 M_⊕). These minimum masses are shown in Fig. 11 and listed in Table A.2. Of our 67 resolved discs, 26 would need masses greater than 100 M_⊕ to self-stir (to 1σ), including six requiring masses greater than 1000 M_⊕. Table 1 lists the systems where the disc masses required for self-stirring are compatible with being 1000 M_⊕ or greater (to 1σ), However, Krivov & Wyatt (2021) argue that debris disc masses greater than 100-1000 M_⊕ may be unfeasible, since these could exceed the total mass in solids available in their parent protoplanetary discs. This implies that some, and potentially many, of our debris discs cannot be self-stirred because this would require unreal-istically high disc masses. In general, Krivov & Wyatt (2021) argue that the resolution to the ‘debris disc mass problem’ is that the largest planetesimals in some debris discs are smaller than previously assumed (s_max ~ 1 km, rather than ~ 100 km), but in this case these discs would be incapable of self-stirring. To summarise, even if we assume the optimal case where the largest planetesimals are just large enough to stir the discs, it is likely that a reasonable fraction of our discs cannot be self-stirred without exceeding their allowed mass budget. If correct, these discs would need to be stirred by other mechanisms. Alternatively, it may be possible to self-stir these discs without requiring unreasonable masses, if we assume a very different initial size distribution (e.g. Najita et al. 2022).

Figure 11 also shows that discs requiring larger masses to self stir also need larger planets to sculpt their inner edges (Sect. 3). This is expected because both the sculpting planet mass and self-stirring disc mass increase with increasing disc size, and decrease with increasing system age. For the discs requiring potentially unfeasible masses to self-stir, this means that large planets are typically expected in those systems anyway, whatever the stirring mechanism. An obvious hypothesis is therefore that planets stir the discs that cannot be self-stirred, and these could be the same planets that sculpt the disc inner edges (although not necessarily). In the following section we calculate the planets required to stir the debris discs, and show that known or observationally allowed companions could stir all discs in our sample.

Fig. 11 Self-stirring constraints in systems with resolved debris discs. The vertical axis shows the minimum masses of discs required for those discs to be self-stirred (Sect. 4.1). Horizontal axes show the minimum planet masses required to sculpt the inner edges of those discs (Sect. 3) if debris clearing is performed by one planet (left plot) or by multiple, equal-mass planets (right plot). Symbols have the same meanings as in previous figures. Background colours denote debris disc masses that may be unfeasibly large (red), borderline (yellow), or allowed (green), according to the quantity of solids that may have been available at the protoplanetary disc stage (Krivov & Wyatt 2021). It may not be possible for discs in the red and yellow regions to be self-stirred, in which case other stirring mechanisms are required. Since these discs also typically require large planets to sculpt their inner edges, it is likely that planet-stirring operates in these systems, as shown in Sect. 4.2. The requirement for large planets to both stir and sculpt these discs makes systems in the yellow and red regions prime targets for future planet searches, as discussed in Sect. 5.2.

4.2 Stirring the debris discs with planets

The previous section examined the disc masses required for self-stirring, and showed that some would be unfeasibly high. There are two main alternatives to explain these discs; either they were ‘pre-stirred’ (i.e. debris orbits were already crossing at the end of the protoplanetary disc phase; Wyatt 2008), or unseen planets are responsible for stirring (Mustill & Wyatt 2009). Here, we explore the latter scenario, where each disc is stirred by a single planet. Secular perturbations by this planet would cause debris eccentricities and orientations to evolve, eventually causing debris orbits to cross and initiating the collisional cascade. Since the secular interaction does not change the debris and planet semi-major axes, if the primordial debris has very low eccentricities then the perturbing planet would have to be eccentric to stir a disc. For this analysis we assume that primordial debris is unexcited, and therefore that any stirring is driven by a planet with non-zero eccentricity. Whilst a single planet could be responsible for both stirring and sculpting the disc, we treat these two analyses separately; the stirring constraints derived in this section are independent from the sculpting constraints from Sect. 3.

4.2.1 Stirring the debris discs with planets: Method

We assume that the stirring planet is located interior to the inner edge of the disc. The stirring timescale decreases with planet mass M_p,stir, and increases with the outer radius of the disc; a lower bound on stirring planet mass is therefore the planet mass that is just sufficient to stir the disc outer edge within the stellar age. Mustill & Wyatt (2009) give timescales for an eccentric planet to stir debris at a given radius; rewriting their Eq. (15) in our units yields the time for a planet to stir the outer edge of an external disc as $$t_{\text{p,stir}}=5.07\times {10}^{-5}\text{Myr}{\left(\frac{a_2}{\text{au}}\right)}^{9/2}{\left(\frac{a_{\text{p}}}{\text{au}}\right)}^{-3}\times \frac{{\left(1-e_{\text{p}}^2\right)}^{3/2}}{e_{\text{p}}}{\left(\frac{M_{\text{p,stir}}}{M_{\text{Jup}}}\right)}^{-1}{\left(\frac{M_{\ast }}{M_{\odot }}\right)}^{1/2},$$(22)

where a₂ is the semi-major axis of the outermost disc particle.

If the planet stirs the disc then the timescale from Eq. (22) must be shorter than the stellar age. Also, since we assume that the planet is interior to the inner edge of the disc, we have the additional constraint that a_p(1 + e_p) < Q_i, where Q_i is the apoc-entre of the disc inner edge (for axisymmetric discs this is just the inner edge location). Combining these two constraints with Eq. (22) yields the minimum mass of an eccentric planet required to stir an external disc as $$M_{p,stir}\gtrsim 5.07\times {10}^{-5}{M}_{\text{Jup}}{\left(\frac{q_o}{\text{au}}\right)}^{9/2}{\left(\frac{Q_{\text{i}}}{\text{au}}\right)}^{-3}\times {\left(\frac{t_{\ast }}{\text{Myr}}\right)}^{-1}{\left(\frac{M_{\ast }}{M_{\odot }}\right)}^{1/2}\frac{{\left(1+e_{\text{p}}\right)}^3{\left(1-e_{\text{p}}^2\right)}^{3/2}}{e_{\text{p}}}.$$(23)

To make Eq. (23) general for both axisymmetric and eccentric discs, we set a₂ = q_o, the pericentre of the disc outer edge, because if disc eccentricity is driven by an eccentric planet, then this is the semi-major axis of the outermost debris particle (Eq. (7) in Pearce & Wyatt 2014). We note that the minimum-mass stirring planet would have apocentre at the disc inner edge; planets located farther inwards would require larger masses to stir the disc.

Equation (23) depends on planet eccentricity; the higher the eccentricity, the lower the planet mass required for stirring. However, this eccentricity is unknown. If the planet also sculpts the disc inner edge then a lower bound on planet eccentricity is given by Eq. (10), but this is only a lower limit because a planet with higher eccentricity located farther inwards could also sculpt the disc to the same degree (Eq. (5)). Also, since the great majority of our discs are either unresolved or axisym-metric, we cannot bound the eccentricity of a stirring planet in most systems. In order to estimate the minimum planet mass required to stir the disc, we must therefore assume a maximum planet eccentricity. Denoting the eccentricity term in Eq. (23) as $f\left(e_p\right)\equiv {\left(1+e_{\text{p}}\right)}^3{\left(1-e_{\text{p}}^2\right)}^{3/2}/e_{\text{p}}$, we see that for small planet eccentricities $f\left(e_{\text{p}}\right)\propto e_{\text{p}}^{-1}$, whilst for larger eccentricities the scaling is roughly ƒ(e_p) ∝ (1 − e_p). The transition between these two regimes occurs around e_p ≈ 0.3. For the planet-stirring analysis we therefore assume that the eccentricity of the stirring planet is less than 0.3 because for smaller eccentricities the minimum planet mass is very sensitive to the exact choice of eccentricity, whilst for larger eccentricities the dependence is less severe. A maximum eccentricity of 0.3 is also consistent with the planet eccentricities required to sculpt the inner edges of the discs (Sect. 3.1). We therefore assume a maximum planet eccentricity of 0.3 for this analysis; for comparison, if we had assumed a maximum eccentricity of 0.1 then the minimum planet masses required for stirring would be doubled. We note that we expect the great majority of planets to actually have eccentricities smaller than 0.3, so assuming 0.3 as an upper limit will set a lower bound on the stirring masses required.

Assuming that the stirring planets have eccentricities smaller than 0.3, we use Eq. (23) with e_p = 0.3 to find the minimum planet masses required to stir the discs. For resolved discs we take Q_i and q_o from the observed disc edges, and for discs with only SED data we calculate conservative lower limits on the stirring planet masses by assuming that those discs have zero width, are located at the corrected blackbody radii, and are axisym-metric (i.e. Q_i = q_o = r_SED; non-zero widths would increase the required stirring planet masses). Using this method we plot the planet-stirring constraints for 49 Cet (HD 9672) in Fig. 7 as an example (line 7); if a single planet sculpts this disc then that planet could also stir it, but if multiple planets instead sculpt the disc, then those planets may be too small to perform the stirring themselves.

4.2.2 Stirring the debris discs with planets: Results

We now apply the above method to all of our systems, and list the resulting minimum planet masses for stirring in Table A.2. In Fig. 12, we compare these minimum stirring masses to our L′-band observational upper limits on planet masses, and also to the disc masses required for self-stirring from Sect. 4.1.

The top plot of Fig. 12 shows the minimum planet masses required to stir the discs, divided by our L′-band planet detection limits at the semi-major axes of those planets. The plot shows that all of our discs could be stirred by unseen planets (i.e. the minimum planet masses required are all consistent with being below the detection limits, to 1σ), Almost all of the required planet masses are orders of magnitude too small to be detectable with current instrumentation, with only one system (β Leo: HD 102647) needing a planet mass greater than 10% of our L′-band detection limit (to 1σ) to stir the disc (this system is discussed below). For a further 17 systems, the stirrers could have masses above 10% of our detection limits too, but their large uncertainties mean they could also be below 10% of the detection limits (to 1σ), This leaves 160 of our 178 systems where the minimum planet masses required for stirring are well below current detection limits; therefore, even order-of-magnitude improvements in mass-detection limits may not be sufficient to confirm or exclude planetary stirrers, unless those planets also sculpt the discs (see below).

The bottom plot of Fig. 12 compares the minimum planet masses required for planet-stirring to the minimum disc masses required for self-stirring. Generally, discs requiring larger masses to self-stir also need larger planets to planet-stir, since both mass requirements increase with disc size and decrease with stellar age. However, the discs that require unfeasibly high masses to self-stir could alternatively be stirred by moderate-mass planets; for example, the outer disc of 49 Cet (HD 9672) requires a disc mass of at least 2300 ± 200 M_⊕ (7.3 ± 0.7 M_Jup) to self-stir, but could alternatively be stirred by a 1.2 ± 0.5 M_Jup planet with eccentricity ≤ 0.3. Table 1 shows that all systems potentially requiring disc masses above 1000 M_⊕ to self-stir could instead be stirred by planets with masses of 1 MJup or less (to 1σ), Since all minimum-mass stirring planets are consistent with being below the detection limits (Fig. 12, top plot), it appears that whilst not all debris discs can be self-stirred (unless they have very high masses, argued against by Krivov & Wyatt 2021), all discs (in our sample at least) could be planet-or companion-stirred. This could point towards general conclusions about the relative prevalence of self-stirring and planet-or companion-stirring across all debris discs. However, we urge caution on this point; many of our discs could still be self-stirred with masses below 100-1000 M_⊕.

The bottom plot of Fig. 8 shows the distribution of minimum planet masses required to stir our discs. Whilst these masses appear considerably smaller than those required to sculpt the discs, the stirring planet mass distribution is biased towards unre-alistically small values because SED discs were modelled as having zero width. It is still useful to provide these estimates as conservative lower bounds for the stirring planet masses in specific systems without resolved disc data, but more realistic estimates of the ‘typical’ planet masses required for stirring arise by considering only discs resolved by ALMA or Herschel (although we note that these also may not be fully representative, since ALMA and Herschel observations are biased towards the largest, brightest discs). For these systems, the median planet mass required for stirring is 0.05 M_Jup (20 M_⊕), with first and third quartiles of 0.02 M_Jup (5 M_⊕) and 0.6 M_Jup, respectively.

Figure 13 compares the planet masses required to stir the discs to those required to sculpt the discs, for the two scenarios where either one or multiple planets perform the sculpting. The left plot shows that, for all but two of our systems, if a single planet sculpts the disc then that planet could also stir the disc (i.e. the minimum mass required by a single planet to sculpt the disc is greater than that required to stir the disc, or the two are consistent to σ), The two exceptions (both Herschel discs) are β Tri (HD 13161) and β Leo (HD 102647). However, β Tri is a close binary with a circumbinary disc, and Kennedy et al. (2012) show that the binary would be capable of stirring the disc. For β Leo, Churcher et al. (2011b) argue that multiple belts exist, with one or two inner belts plus an outer disc spanning 15–70 au (compared to the 7–79 au that we use in our single-disc model); it is therefore likely that we underestimate the outer disc inner edge, and taking this to be 15 rather than 7 au reduces the required stirring planet mass by an order of magnitude, below that needed to sculpt the disc (the general effect of underestimating the disc inner edges in our Herschel models is discussed in Sect. 5.3.2). In conclusion, if the inner edges of our discs are each sculpted by a single planet, then these planets could also fully stir the discs (except for β Tri, where stirring could instead be performed by the binary).

Conversely, the right plot of Fig. 13 shows that, if multiple planets sculpt each disc, then for eight systems the minimum-mass sculpting planets would be too small to also stir the discs (to These systems are 49 Cet (HD 9672), β Tri (HD 13161), ψ⁵ Aur (HD 48682), HD 92945, β Leo (HD 102647), HD 107146, HD 131835, and HD 138813. If the inner edges of these systems are carved by multiple planets, then those planets must either be larger than the minimum masses derived in Sect. 3.2 (in which case they could also stir the discs), or the stirring is performed by some other process(es). Some of these discs could be self-stirred without requiring disc masses greater than 1000 M_⊕, so for these systems a mixture of self-stirring and multi-planet sculpting could be occurring. Alternatively, the planet(s) responsible for stirring the discs could be different to those sculpting them; larger planets located farther in could stir the discs, whilst smaller ones at the disc edges could sculpt them (discussed further in Sect. 5.1.2). Another possibility is that, since some of these discs contain gaps that could host planets, planets located in the gaps could be responsible for stirring, whilst other planets sculpt the inner edges. However, two of the above discs (49 Cet and HD 138813) have no observed gaps, and would require masses of more than 1000 M_⊕ to be self-stirred, which are possibly unfeasibly high; for these systems it seems likely that the sculpting planets are larger than the minimum masses derived assuming multi-planet sculpting, and that these planets also stir the discs, or that other interior planets stir the discs. We note that both of these discs also have gas detected (Hughes et al. 2017; Lieman-Sifry et al. 2016), so dust-gas interactions could also be occurring in them (e.g. Pearce et al. 2020). In summary, if one planet sculpts each of our discs, then the discs could also all be stirred by that planet or a known companion; conversely, if the discs are sculpted by multiple planets, then for some systems the minimum planet masses required for sculpting are smaller than those needed for stirring.

Fig. 12 Constraints on the planet masses required to stir debris discs. Horizontal axes show the minimum planet masses required for stirring if those planets have eccentricities smaller than 0.3, as described in Sect. 4.2 (higher eccentricities would require lower masses). The vertical axis definitions are as follows. Top plot: minimum stirring-planet masses, divided by our L′-band planet detection limits at the apocentres of those planets. The minimum stirring-planet masses are all consistent (to 1σ) with lying below the detection limits (i.e. below the dashed line, which shows unity), so all of the discs may potentially be stirred by unseen planets. The stirring planet predictions for systems where discs are unresolved (black circles) are conservative lower limits since these discs are modelled as having zero width. Bottom plot: minimum masses of the resolved discs required for those discs to be self-stirred, from Fig. 11. Many systems would require disc masses greater than 1001000 M⊕ to self-stir, which may be unfeasibly high (Krivov & Wyatt 2021); however, the top plot shows that these systems could still be stirred by unseen planets. This suggests that self-stirring is not the dominant stirring mechanism in at least some debris discs. No unresolved discs are shown on the bottom plot because they were omitted from the self-stirring calculation owing to a lack of width information. Also, we note that a few systems on the bottom plot are missing from the top plot because their predicted planet locations are not covered by our observations. The uncertainties on stirring planet masses are large owing to the very strong dependence of Eq. (23) on disc edge locations.

Fig. 13 Comparison of the planet masses required to stir and sculpt debris discs. The vertical axis shows the minimum planet masses required to stir each disc (Sect. 4.2). Horizontal axes show the minimum planet masses required to sculpt the inner edges of those discs (Sect. 3) if sculpting is performed by one planet (left plot) or by multiple, equal-mass planets (right plot). Symbols are the same as in previous figures, and the dashed lines denote 1:1 relations. If one planet sculpts each disc inner edge (left plot), then for almost all systems that planet could also stir the disc (either the sculpting planet mass is greater than that required for stirring, or the two are consistent to 1σ), The only exception is β Tri (HD 13161), which is a binary system where the binary could instead stir the disc (for the other labelled system, β Leo (HD 102647), we likely underestimate the disc inner edge and therefore overestimate the required stirring planet mass). Conversely, if multiple planets sculpt each disc (right plot), then for several systems the minimum-mass sculpting planets would be too small to stir those discs (to 1σ), some, but not all, of those discs could feasibly be self-stirred, so either the planets in those systems are more massive than required by the multi-planet sculpting model, or disc sculpting and stirring are performed by different entities. We note that the stirring planet masses derived from SEDs of unresolved discs (black circles) are conservative lower limits since these discs were assumed to have zero width for this analysis; in this case, the planet masses required for sculpting and stirring have similar or identical dependences on system parameters, resulting in the tight relationships between the two. As on previous plots, the large uncertainties are due to the strong dependences of Eqs. (15) and (23) on disc edge locations.

5 Discussion

We used sculpting and stirring arguments to predict planet parameters in 178 debris-disc systems. In this section we compare our predictions to known planets and those inferred from protoplanetary discs (Sect. 5.1), consider the systems best-suited for future planet searches (Sect. 5.2), and discuss several caveats that must be considered when interpreting our results.

5.1 Comparisons to known planet populations

Here we compare our planet predictions to both the detected planet population, and to planets inferred to be forming in protoplanetary discs.

5.1.1 Predicted planets versus detected and forming planet populations

Our debris disc analysis predicts that planets of at least super-Earth- to Jovian-masses commonly reside between 10–100au. We show our planet predictions from sculpting and stirring arguments in Fig. 14, compared to the known planetary population. Our predicted planets occupy a different parameter space to known exoplanets, in a region that is difficult to probe with current instrumentation and consequently hosts few known objects. However, our analysis predicts that this parameter space should harbour a rich family of planets, some of which should be detectable with next-generation instruments (Sect. 5.2). Whilst dissimilar to known exoplanets, our predicted planets are more similar to those in the outer Solar System; Uranus and Neptune in their current positions would be capable of sculpting and stirring about half of our discs, and many of these systems may well contain Uranus or Neptune analogues. For the remaining systems, we typically require either ice giants at 2–4 times the distances of Uranus and Neptune or Jovian planets at 5–10 times the distance of Jupiter; such bodies have no known Solar System analogues.

Whilst dissimilar to the known exoplanets, our predicted planets are similar to those inferred to be forming in protoplan-etary discs; the crosses on the right plot of Fig. 14 show the young planets inferred to be carving gaps in ALMA-resolved protoplanetary discs (Lodato et al. 2019). Both planetary and non-planetary origins have been suggested for such gaps (e.g. Okuzumi et al. 2016; Dullemond et al. 2018; Richert et al. 2018; Zhang et al. 2018; Pinte et al. 2020), but if these gaps are created by forming planets then we may be witnessing planets being born that will later join the population responsible for sculpting and stirring known debris discs. Of course, this picture could be complicated by planetary migration both in and beyond the gas disc phase, as predicted from population synthesis models (e.g. Alibert et al. 2013; Pfyffer et al. 2015; Mordasini 2018). However, it is interesting that the planets potentially forming in protoplan-etary disc gaps, if they remained in situ, overlap very strongly with those that we infer to exist from debris disc observations. If these inferred forming planets are destined to migrate to inwards of 10 au as suggested (e.g. Lodato et al. 2019; Wang et al. 2021), then either a different population of planets is required to later interact with debris discs, or these forming planets must quickly sculpt and stir the debris discs before migrating inwards (see Sect. 5.1.2), or debris discs are neither sculpted nor stirred by planets. Alternatively, since the forming planets already have the masses and locations of planets inferred from debris discs, this could be evidence that the forming planets do not migrate as far as previously thought.

5.1.2 Predictions for systems with known planets

Planets have already been detected in some of our systems. Here we compare these to our inferred planets, to assess whether the known planets could be those interacting with our discs.

Of our 178 debris disc systems, 15 host at least one planet residing interior to the disc inner edge. These planets were detected using a variety of techniques, and specific details for individual systems are given in Appendix B. In Fig. 15, we compare the masses and semi-major axes of these detected planets to those required to sculpt the associated discs; we consider single-planet sculpting here (Sect. 3.1) because it predicts the required planet mass as a function of its semi-major axis (rather than assuming that the outermost planet resides exactly at the disc inner edge, as in the multi-planet model). If several planets have been detected in a system then we only consider the outermost planet, and we estimate the masses of RV planets by assuming that they share the same inclination as the associated debris disc.

Figure 15 shows that, whilst many known planets would have masses sufficient to sculpt their associated discs were the planets located farther out in those systems, these detected planets are mostly located too far inwards to be sculpting the discs in situ; the mass required by a sculpting planet sharply increases as its semi-major axis decreases (line 4 in Fig. 7). The possible exceptions are 51 Eri (HD 29391) and HD 69830, where the outermost planets may be capable of sculpting the disc inner edges in situ; for the remaining systems, if disc inner edges are sculpted by planets then these planets probably reside exterior to the known ones. However, for systems above the horizontal line in Fig. 15, it is also possible that the detected planets formed in the outer regions, sculpted and stirred the discs, and have since migrated inwards to their current locations. This picture is consistent with the expectation that young planets inferred to reside in gaps in ALMA-resolved protoplanetary discs will later migrate inwards of 10 au (see Sect. 5.1.1), although in this case some planets would have to be considerably more massive than ourpredictions to rapidly sculpt or stir debris discs before migrating.

Whilst most of the known planets are unlikely to be sculpting discs in situ, a larger fraction of them could be stirring the discs. The circles in Fig. 15 mark systems where the known planet would, in their present location, have sufficient mass to stir the disc within the system age according to Eq. (22) (assuming the planet eccentricities are smaller than 0.3, as in Sect. 4.2). This is possible because the planet masses required for stirring typically have a much shallower dependence on semi-major axis than the masses required for sculpting (compare lines 4 and 7 in Fig. 7), since it is often easier to drive debris orbit crossings at the disc outer edge than it is to eject all material from the disc inner edge. Some detected planets could therefore still be interacting with the discs by stirring them, even if they are located too far inwards to sculpt them; this emphasises the general point that planets responsible for stirring discs need not be located at the disc inner edges, but could reside much closer to the star.

Fig. 14 Planets predicted from our debris analyses, compared to known planets and those inferred to be forming in protoplanetary discs (Sect. 5.1.1). Small black circles are known exoplanets from Exoplanet.eu, and large circles are Solar System planets. The left plot shows the planetary populations inferred from our debris analyses; green diamonds and blue crosses are the planets required to sculpt the discs if sculpting is performed by one or multiple planets, respectively (Sects. 3.1 and 3.2). Red circles are planets required to stir the discs (only systems with resolved discs are shown; Sect. 4.2). The plot shows that debris disc morphologies imply planetary populations that are poorly explored with current instruments, so the hypotheses that planets sculpt and stir debris discs are consistent with the lack of planet detections in these systems. Right plot: Comparison of our planets inferred from debris discs to those inferred to be forming in protoplanetary discs (black crosses; Lodato et al. 2019). To make the plot clearer, we have replaced our individual planet predictions with rough areas representing the three populations. Our planetary population inferred from debris discs well matches the population inferred to be forming in protoplanetary discs; if the forming planets remain in situ, then they could later join the population responsible for sculpting and stirring known debris discs. Alternatively, it is possible that the forming planets rapidly sculpt debris discs before migrating inwards. Error bars are omitted for clarity but can be seen in our other figures.

5.2 Potential for future observations

Many of our predicted planets lie just below current detection limits, so could be observable in the near future. We now discuss the detectability of our predicted planets, and identify some of the more promising targets for future direct-imaging observations.

To detect a planet via direct imaging, it must be bright enough to lie above the contrast curve at its projected separation from the star. This naturally favours massive planets, and there is also a bias towards younger planets because these are more luminous (Baraffe et al. 2003; Phillips et al. 2020). For our predicted planets these two considerations are linked because younger systems require higher-mass planets to sculpt and stir the discs within the system lifetime; our systems requiring the most-massive planets are therefore also the youngest, so these should be strongly favoured for future observations. For a planet to be detectable it must also have a wide angular separation from its host star at the moment of observation; this not only favours planets on wide orbits, but also those with low sky-plane inclinations (because these would be well separated from their stars in projection for greater fractions of their orbits). Since we assume that our predicted planets are coplanar with the debris discs, we therefore also favour systems with low disc-to-sky inclinations as targets for future observations.

Figure 16 shows the minimum masses, maximum projected separations, ages and inclinations of the planets that we predict to be sculpting our discs (assuming that one planet is responsible in each case; Sect. 3.1). We only plot systems where the disc has been resolved by ALMA or Herschel, since these provide the most-reliable planet constraints. The figure shows that we predict planets of up to several Jupiter masses at arcsec-ond separations, that our most-massive predicted planets are also young, and that a reasonable number are expected on low-inclination orbits that maximise the chances of observations. These compare favourably with the observational capabilities of the upcoming JWST; 0.1−1 M_Jup planets at arcsecond separations around young stars such as those in the β Pic and TW Hya moving groups would have detection thresholds of 50% when observed for 1 hour with the Mid-Infrared Instrument (MIRI; Beichman et al. 2010; Carter et al. 2021). Other upcoming instruments should also be able probe our planet parameter space; these include the Enhanced Resolution Imager and Spec-trograph (ERIS; Davies et al. 2018) at the VLT, due to start operations in 2022. Considered the successor of NaCo, it will provide enhanced adaptive optics and improved sensitivity in the J to M bands. In addition, the grating vector Apodizing Phase Plate (gvAPP; Boehle et al. 2018) should enable high-contrast imaging observations, and should complement JWST by probing small angular separations, where some of our most-massive planets are predicted to reside. Further ahead, the Mid-infrared E-ELT Imager and Spectrograph (METIS; Brandl et al. 2018) should detect planets down to several Earth masses at small angular separations (Quanz et al. 2015), thanks to its 39 m mirror that increases angular resolution by a factor of ~4–5 compared to current telescopes. The expected background-limited magnitude for a one-hour observation in L′ with METIS is 21.2 mag (Bowens et al. 2021), which is several orders of magnitude deeper than current NaCo limits (~ 16−17 mag, Launhardt et al. 2020). Taken together, many of our predicted planets should be observable by upcoming instrumentation.

Table 2 shows all ALMA and Herschel systems where if the discs are each sculpted by one planet, then that planet would be more massive than 0.5 M_Jup (to Many of these systems would be good targets for future planet searches. In particular, we can be reasonably confident that large planets are interacting with the discs of 49Cet (HD 9672), HD 121617, HD 131835, and HD 138813; these all require large sculpting planets (>1 M_Jup if one planet, ≳0.1 M_Jup if multiple planets), are unlikely to have disc masses sufficient for self-stirring, and would need ≳ 0.1 M_Jup planets for planet-stirring. These systems are also all young, and the discs of HD 121617 and HD 138813 have low-to-moderate inclinations, making these good candidates for future observations. We note, however, that these systems all have gas detected, so the grain dynamics could also be more complicated than assumed in this paper (e.g. Pearce etal. 2020). Another promising target is HR 8799 (HD 218396), which is young, has a low inclination and requires a well-separated, ≳0.1 M_Jup planet to sculpt its disc (in addition to the four known planets; see Appendix B.62).

Figure 16 shows our predicted minimum-mass planets if the discs are each sculpted by a single planet. If we instead consider our more conservative, multi-planet model (Sect. 3.2) then the sculpting planets could be less massive by an order of magnitude or more, whilst their projected maximum separations would be similar to those plotted. Some of these would still be in the parameter space probed by upcoming instruments, but some would not be. Non-detections of the above planets would therefore not necessarily rule out planet sculpting in some cases, but could potentially differentiate between the scenario where one planet dominates the truncation from that where sculpting is performed by multiple planets in tightly packed arrangements (or scenarios where sculpting was performed by a massive planet that has since migrated inwards). Nonetheless, the smaller planets predicted by the multi-planet sculpting model are conservative lower bounds on the required planet masses for disc sculpting, and in many cases would be incapable of stirring the discs (Sect. 4.2); in practice, we would therefore expect upcoming instruments such as JWST to detect planets in many of our systems.

Fig. 15 Comparison of detected and predicted planets in our debris disc systems, as described in Sect. 5.1.2. For systems with a detected planet, the vertical axis shows the mass of one detected planet divided by the mass required by a single planet to sculpt the disc inner edge (Sect. 3.1). We calculate the masses of RV planets by assuming that they are copla-nar with the discs. The horizontal axis shows the semimajor axis of the detected planet, divided by the semimajor axis that the minimum-mass sculpting planet would have. For systems with multiple known planets, only the outermost detected planet is shown. Dashed lines denote unity. With the possible exceptions of 51 Eri (HD 29391) and HD 69830, the detected planets are unlikely to be massive enough to sculpt the disc inner edges in situ; whilst some are larger than the minimum masses required for sculpting, these detected planets lie far interior to the discs and so would need significantly greater masses to sculpt the discs in situ (line 4 on Fig. 7). This means that either additional undetected planets also reside in these systems, or the known planets formed in the outer regions and rapidly sculpted the discs before migrating inwards. Conversely, some detected planets are massive enough to have stirred the discs in situ, despite their smaller semimajor axes; systems where detected planets could stir the discs are circled. Arrows show lower limits on detected planet masses for RV planets with no disc inclination data, and the remaining symbols were defined on previous figures.

5.3 Caveats

Finally, we discuss several caveats that must be considered when interpreting our results: the potential effects of non-zero disc mass (Sect. 5.3.1), Herschel's spatial resolution (Sect. 5.3.2) and close binaries (Sect. 5.3.3) on our planet predictions.

Fig. 16 Identification of our predicted planets with the best chances of future detection (Sect. 5.2). Points show the minimum masses and maximum projected separations of the planets predicted to be sculpting the inner edges of our discs if one planet is responsible in each system (Sect. 3.1). Only systems where the disc is resolved by ALMA (diamonds) or Herschel (squares) are shown. Left plot: systems coloured by age. Right plot: systems coloured by disc inclination to the sky. The best targets have predicted planets with high masses, which are well separated from their stars, are young (so the planets are brightest), and have low-inclination discs (so a coplanar planet would be separated from its star in projection over a large fraction of its orbit). JWST will be capable of imaging 0.1−1 MJup planets at approximately arcsecond separations (Beichman et al. 2010; Carter et al. 2021), and this should be complemented at small angular separations by upcoming instruments such as VLT/ERIS (Davies et al. 2018), so many of these predicted planets could be observable in the near future. The plots show the single planets required to sculpt each disc, but if multiple planets were instead responsible (Sect. 3.2), then those planets could be at least an order of magnitude less massive than those shown here. The planets required to stir each disc are typically smaller than those shown (Sect. 4.2). We note that some systems are missing from the right plot because they lack inclination information.

5.3.1 Omission of disc mass in truncation and planet-stirring models

Whilst disc mass is explicitly included in our self-stirring model (Sect. 4.1), it is omitted from our truncation and planet-stirring models (Sects. 3 and 4.2). These models therefore implicitly assume that the mass of the planet(s) is significantly greater than that of the disc. We assume this because debris disc masses are fundamentally unknown (e.g. Krivov & Wyatt 2021), and because the lower limits on perturbing-planet mass arise if the disc has zero mass (i.e. the planet need not be massive enough to overcome disc gravity). In reality, planets interacting with a massive debris disc would migrate as they scattered particles (Ida et al. 2000; Kirsh et al. 2009; Pearce & Wyatt 2014, 2015), and disc self-gravity could allow debris to resist deformation and stirring (Löhne et al., in prep.). The result would be that truncating planets could clear a significantly wider region (and hence be located much farther inwards) than our massless-disc models imply (Friebe et al. 2022), or that the planets required to sculpt and stir discs could be more massive than our predictions.

5.3.2 Limitations of Herschel versus ALMA data

The disc widths derived from Herschel PACS are systematically larger than those from ALMA, as expected from the lower resolution of the former. This, combined with our neglecting any additional inner emission in the Herschel fitting (Sect. 2.2.2), implies that the inner edges of Herschel discs are likely underestimated. This would lead to the masses and locations of sculpting planets also being underestimated, as appears to be the case from Fig. 9, where the planet predictions for Herschel systems are offset from those for ALMA systems. Similarly, underestimating the disc inner edges would both increase the predicted planet masses needed for planet-stirring, and (for wide discs) decrease the disc masses required for self-stirring. Both of these effects appear to occur, as evidenced by the distributions of ALMA and Herschel points on the bottom plot of Fig. 12. Whilst observational biases are also present (the largest, brightest discs are likely to be those selected for ALMA observations, and these would be expected to need large masses if they were to self-stir), the reader should be aware of the above caveats when considering the dynamical predictions for Herschel systems.

5.3.3 Effects of close binaries on predictions

Several of our systems are possible close binaries, which could impact our dynamical predictions in several ways. Our measured disc parameters should be unaffected by binarity because dust emission is determined by temperature, which is set by the stellar luminosity; since our luminosities correspond to the combined star luminosities in binary systems, the dust masses and locations that we use should not depend on binarity. However, the effective central mass is affected by close binarity, which in turn affects our dynamical calculations; for such systems, both the discs and our predicted planets are assumed to be circumbinary, and so the star masses M_* used in our equations should be the combined binary mass. Close binarity could also affect our stellar mass estimates, as well as our ages for field stars, which would in turn affect our derived planet masses. Finally, potentially the largest effect is that some binaries could sculpt and/or stir their debris discs, removing the need for planets in those systems altogether. Whilst it is unlikely that our potential binaries could sculpt their debris discs, it is possible for some to have stirred them (e.g. β Tri; HD 13161).

6 Conclusions

We analysed 178 debris-disc systems in the ISPY, LEECH, and LIStEN planet searches, using dynamical arguments to constrain the properties of unseen planets from disc morphologies. We show that a typical’ cold debris disc requires a Neptune- to Saturn-mass perturbing planet located at 10–100 au, with many needing Jupiter-mass perturbers. All of our discs could be both sculpted and stirred by either known companions or unseen planets, which is consistent with current non-detections of planets in many of these systems. Our predicted stirring planets are often consistent with those required to sculpt the discs, suggesting that in many systems a single planet may be responsible for both sculpting and stirring debris. Conversely, we show that many discs cannot be self-stirred without requiring unreasonable disc masses. Our predicted planet populations are quite different from those currently detected, both in the Solar System and beyond, but they are very similar to those inferred to be forming in pro-toplanetary discs; this could imply that these are the same planet population if the forming planets do not undergo significant migration. For systems with known planets, few could be sculpting their associated discs in situ, although many could have done so in the past and since migrated inwards. Conversely, many of the known planets could be stirring their associated discs in situ. Finally, we provide a catalogue of planet predictions for our 178 systems, show that JWST should be able to find many of these planets, and identify the most promising targets for future planet searches.

Acknowledgements

We thank Luca Matrà for sharing unpublished REASONS modelling results, Anne-Lise Maire and Élodie Choquet for their HR 8799 and 49 Cet planet-detection limits respectively, Virginie Faramaz for the HD 202628 ALMA image, and Johan Olofsson for helpful comments. We also thank André Müller and Ingo Stilz, whose previous efforts made this work possible, and Nicole Pawellek and Sebastián Marino for useful discussions. We are grateful to the anonymous referee, whose comments and suggestions improved the paper. TDP, MB and AVK are supported by the Deutsche Forschungsgemeinschaft (DFG) grants Kr 2164/13-2, Kr 2164/14-2, and Kr 2164/15-2. GMK is supported by the Royal Society as a Royal Society University Research Fellow. GC thanks the Swiss National Science Foundation for financial support under grant number 200021_169131. TH acknowledges support from the European Research Council under the Horizon 2020 Framework Program via the ERC Advanced Grant Origins 83 24 28. ECM acknowledges support from the Swiss National Science Foundation (SNSF); this work has been carried out in part within the framework of the National Centre for Competence in Research PlanetS supported by SNSF. We made use of the Data & Analysis Center for Exoplanets (DACE) at the University of Geneva (CH), a platform of the Swiss National Centre of Competence in Research (NCCR) PlanetS. The DACE platform is available at https://dace.unige.ch.

Appendix A System data

Table A.1 lists the parameters that we use for each system, and Table A.2 gives the planet and disc constraints that we calculate.

Appendix B Notes on specific systems

This section gives additional notes for some systems. We mainly discuss dynamical considerations that may affect our system-specific planet predictions, such as the presence of known planets or companions.

B.1 CPD-72 2713

The ALMA data are from Moór et al. (2020), and the disc model from REASONS (Matrà et al., in prep.). Moór et al. (2020) argue that the disc can be self-stirred, using Herschel data and the Krivov & Booth (2018) method with S_max = 200 km and v_frag = 30 ms⁻¹. However, our analysis with ALMA data and self-consistent S_max and v_frag values shows that it cannot be, unless the disc mass is unfeasibly high.

B.2 GJ 581

There are Neptunes and super-Earths very close to the star (< 0.1 au: Bonfils et al. 2005; Udry et al. 2007; Mayor et al. 2009; Robertson et al. 2014). The outermost planet is at 0.074 ± 0.001 au with M_p sm(i_p) = 5.7 ± 0.4 M_⊕; Fig. 15 shows that this planet cannot be sculpting or stirring the disc in situ.

B.3 GJ 649

An eccentric planet with M_p sm(i_p) = 0.33 ± 0.03 M_jup is detected at 1.14 ± 0.04 au (Johnson et al. 2010), although its eccentricity is uncertain (Wittenmyer et al. 2013). Figure 15 shows that the planet could be stirring the disc, and is massive enough to have sculpted it (but not in situ).

B.4 HD 1461

There are probably two super-Earths inside 0.1 au (Rivera et al. 2010; Diaz et al. 2016). The outermost planet is at 0.112 ± 0.004 au with M_p sm(i_p) = 5.6 ± 0.7 M_⊕; Fig. 15 shows that this planet cannot be sculpting or stirring the disc in situ.

B.5 HD 3003 (β³ Tuc)

This is a multiple system with up to six components (Tokovinin 2018).

B.6 HD 9672 (49 Cet)

The ALMA data are from Hughes et al. (2017), and the disc model from REASONS (Matrà et al., in prep.). There is a large CO gas mass coincident with the cold disc, which could complicate grain dynamics (Zuckerman et al. 1995; Moór et al. 2019). There is tentative evidence of a warp or spiral arms in the gas disc (Hughes et al. 2017). Various planet constraints are shown for this specific system in Fig. 7.

B.7 HD 10647 (q¹ Eri)

The ALMA data are from Lovell et al. (2021) and show that the disc is asymmetric with a sharp inner edge. Lovell et al. (2021) model this as a symmetric disc with a possible resonant clump, but they also show that a non-resonant, eccentric planet with a_pe_p = 2 – 4 au could instead drive the observed asymmetries. If the latter scenario applies, then these planet constraints agree with ours. A planet lies at 2.0 ± 0.2 au with M_p sin(i_p) = 0.9 ± 0.2 M_jup and eccentricity 0.15 ± 0.08 (Butler et al. 2006; Marmier et al. 2013), but Fig. 15 shows that this inner planet cannot be sculpting or stirring the disc in situ.

B.8 HD 10700 (τ Cet)

The ALMA data are from MacGregor et al. (2016) and show an axisymmetric disc. At least four super-Earths lie within 1.5 au (Tuomi et al. 2013; Feng et al. 2017), the outermost at 1.33 ± 0.02 au with M_p sm(i_p) = 4 ± 1 M_⊕ (Feng et al. 2017); Fig. 15 shows that this planet cannot be sculpting or stirring the disc in situ. There are possibly more planets present (Dietrich & Apai 2021), and a 1 – 2 M_jup planet is predicted at 3 – 20 au from astrometry (Kervella et al. 2019); this planet could sculpt and stir the disc in situ. Lawler et al. (2014) show that the planetary system would be stable if an additional planet with mass below that of Neptune exists outside 5 au, which is consistent with the planet(s) we infer.

B.9 HD 13161 (β Tri)

The disc is circumbinary (binary separation 0.3 au, eccentricity 0.4; Pourbaix 2000). Kennedy et al. (2012) show that the binary could stir the disc.

B.10 HD 15115

The ALMA data are from MacGregor et al. (2019) and show no significant asymmetries, although the disc has a possible eccentricity of 0.06 in the H-band (Sai et al. 2015). ALMA finds a 14 ± 6 au gap at 59 ± 5 au, which MacGregor et al. (2019) predict is carved by an embedded 0.2 M_Jup planet at 58 au. A migrating planet of comparable mass could also produce the observed gap (Friebe et al. 2022). If a separate planet sculpts the disc inner edge, then we predict that planet has mass ≥ 0.85 ± 0.09 M_Jup and semi-major axis ≤ 34 ± 4 au (our inner planet is larger than the gap-clearing planet because we assume a higher degree of debris clearing). Alternatively, if the gap is created through a secular interaction with a precessing planet near the disc inner edge (e.g. Pearce & Wyatt 2015; Yelverton & Kennedy 2018; Sefilian et al. 2021) then that planet would need a mass of at least 0.05 M_Jup (10 M_⊕) to carve the gap (i.e. for 10 secular times to have elapsed within the system age). A large asymmetry is seen in scattered light (Kalas et al. 2007; Schneider et al. 2014), which does not necessarily imply an asymmetric parent ring (Mazoyer et al. 2014). There is a possible disc interaction with an M star (Kalas et al. 2007).

B.11 HD 17094 (µ Cet)

The disc is circumbinary, where the binary has a separation of 3 au (Trilling et al. 2007). It is possibly a quadruple system (Richichi et al. 2000). The binary is unlikely to be sculpting the disc at ~ 100 au.

B.12 HD 17925 (EPEri)

The central mass is possibly a close binary (two K-types at 0.07au separation; Rodriguez et al. 2015), which would make the 6 - 25 au disc circumbinary. If so, the binary could stir the disc, but probably not sculpt it.

B.13 HD 21997

The ALMA data and disc model are from REASONS (Matrà et al., in prep.). The disc is axisymmetric (Kóspál et al. 2013; Moór et al. 2013). There is a very high CO gas mass, which could affect grain dynamics. It is questionable whether the disc can be self-stirred, as its mass would have to be ${580}_{-80}^{+60}{\text{M}}_{\oplus }$ (also see Moór et al. 2015 and Krivov & Booth 2018).

B.14 HD 22049 (ϵ Eri)

The ALMA data are from Booth et al. (2017), and we use their Gaussian model. The disc is probably axisymmetric in SMA and Spitzer images (MacGregor et al. 2015; Backman et al. 2009), although there may be low-significance clumps. A ${0.66}_{-0.09}^{+0.12}{\text{M}}_{\text{Jup}}$ planet lies much closer in at 3.52 ± 0.04 au, which may be significantly inclined relative to the disc (Llop-Sayson et al. 2021). Figure 15 shows this planet could be stirring the disc, and is massive enough to have sculpted it (but not in situ). Booth et al. (2017) predict a single sculpting planet at > 45 au with mass ≤ 1.3 M_Jup, consistent with our prediction of a ≥ 0.19 ± 0.04 M_Jup planet at $\le {52.7}_{-1.0}^{+0.9}\text{au}$. They also show that a planet at 48 au could be in 3:2 and 2:1 resonance with the disc inner and outer edges, respectively, which is also consistent with our predicted planet.

B.15 HD 27290 (γ Dor)

Our Herschel model fits a single broad disc spanning 14 ± 7 to 280 ± 40 au, but modelling by Broekhoven-Fiene et al. (2013) shows that the system could be fitted by a warm ring at several au plus a cool disc spanning 55 to 400 au, or two narrow rings at 70 and 190 au, or any configuration between.

B.16 HD 29391 (51 Eri)

A directly imaged object with a mass of 2 – 12 M_Jup (depending on whether hot- or cold-start models are used) exists with a semi-major axis of ${11}_{-1}^{+4}\text{au}$, and possibly a non-zero eccentricity (Macintosh et al. 2015; De Rosa et al. 2020). Figure 15 shows that this object could be both sculpting and stirring the disc in situ.

B.17 HD 32297

The ALMA data are from MacGregor et al. (2018), and the disc model from REASONS (Matrà et al., in prep.). There is possibly some unresolved, non-axisymmetric thermal structure. A swept-back structure is visible in scattered light, similar to ‘The Moth’ (HD 61005; Schneider et al. 2014). There is also a large CO gas mass (Moór et al. 2019), which could affect grain dynamics.

B.18 HD 35850 (AF Lep)

This is a close binary, with separation 0.021 au (Rodriguez & Zuckerman 2012); the binary is unlikely to be sculpting the disc at ~ 50 au.

B.19 HD 38206

The ALMA data are from Booth et al. (2021) and show a possibly asymmetric disc (eccentricity ${0.25}_{-0.09}^{+0.10}$). If a single planet both stirs and sculpts the disc, then Booth et al. (2021) predict it to have mass ${0.7}_{-0.3}^{+0.5}{\text{M}}_{\text{Jup}}$, semi-major axis ${76}_{-13}^{+12}\text{au}$ and eccentricity ${0.34}_{-0.13}^{+0.20}$; our predicted planet is at a similar location but is three times more massive because Booth et al. (2021) use the shorter multi-planet sculpting timescale from Shannon et al. (2016) rather than the single-planet timescale from Pearce & Wyatt (2014).

B.20 HD 48682 (ψ⁵ Aur)

There are potential brightness asymmetries in the Herschel images (Hengst et al. 2020).

B.21 HD 50554

A large, eccentric companion is present with M_p sin(i_p) = 5 ± 3 M_Jup, semi-major axis 2 ± 1 au, and eccentricity ~ 0.5 (Fischer et al. 2002; Perrier et al. 2003). Figure 15 shows that this planet could be stirring the disc, and is massive enough to have sculpted it (but not in situ).

B.22 HD 53143

A potential asymmetry is visible in scattered light (Schneider et al. 2014, although also see Kalas et al. 2006). Possible clumps may be visible in micron grains, and the disc is highly asymmetric in unpublished ALMA data (Stark et al. 2020).

B.23 HD 61005 ('The Moth')

The ALMA data are from MacGregor et al. (2018), and the disc model from REASONS (Matrà et al., in prep.). The disc is very asymmetric in scattered light (Schneider et al. 2014; Olofsson et al. 2016). No significant asymmetries are seen by ALMA, but there are hints of unresolved asymmetries. Esposito et al. (2016) show that the structure could be caused by an inclined, eccentric planet and predict a planet with eccentricity ~ 0.2 and semi-major axis < 50 au, with many possible masses. Our predicted sculpting planet with semi-major axis ${41.0}_{-0.5}^{+0.6}\text{au}$ and mass ${0.70}_{-0.04}^{+0.05}{\text{M}}_{\text{Jup}}$ is consistent with this.

B.24 HD 69830

The disc is close to the star. Three Neptunes exist inside 1 au, with M_p sin(i_p) = 10 – 20 M_⊕ (0.03 – 0.06 M_Jup; Lovis et al. 2006). The outermost planet has semi-major axis 0.6 ± 0.3 au and M_p sin(i_p) = 18 ± 9 M_Jup (we assume 50 per cent uncertainties); Fig. 15 shows that this planet could be both sculpting and stirring the disc in situ. However, the system dynamics may be complicated; the planets may even have migrated through the disc (Payne et al. 2009).

B.25 HD 92945

The ALMA data are from Marino et al. (2019), and the disc model from REASONS (Matrà et al., in prep.). There is a gap at 73 au, and small asymmetries are seen (possibly due to a background galaxy). If the gap is carved by an in situ planet, then Marino et al. (2019) argue that that planet has mass < 0.6 M_Jup (although a migrating planet could be smaller; Friebe et al. 2022). Alternatively, if the gap is created through a secular interaction with planets near the disc inner edge (e.g. Pearce & Wyatt 2015; Yelverton & Kennedy 2018; Sefilian et al. 2021), then the outermost planet would have a semi-major axis of 40 – 45 au and mass of at least 10 M_⊕ (0.03 M_Jup). This is consistent with our prediction of a planet sculpting the disc inner edge with mass ≥ 0.32 ± 0.06 M_Jup and semi-major axis $\le {45}_{-4}^{+1}\text{au}$.

B.26 HD 102647 (β Leo)

There are possibly separate hot, warm and cold debris components present (Churcher et al. 2011b; Defrère et al. 2021). Churcher et al. (2011b) argue for either three discs at 2, 9, and 30 – 70 au or two discs at 2 and 15 – 70 au. Our Herschel model assumes only a single broad disc (7 – 79 au), so we may underestimate the outer disc inner edge. A companion lies at 440 au (Rodriguez et al. 2015).

B.27 HD 107146

The ALMA data are from Marino et al. (2018), and the disc model from REASONS (Matrà et al., in prep.). The disc contains a 39 au wide gap at 76 au (Ricci et al. 2015; Marino et al. 2018). Marino et al. (2018) argue that the gap could be carved by multiple 10 M_⊕ (0.03 M_Jup) planets in situ, and Eq. 23 shows that such planets could stir the disc outer edge even with small eccentricities. The gap could also be carved by a single, migrating planet (Friebe et al. 2022). Alternatively, if the gap is created through a secular interaction with a precessing planet near the disc inner edge (e.g. Pearce & Wyatt 2015; Yelverton & Kennedy 2018; Sefilian et al. 2021) then that planet would need a mass of at least 10 M_⊕ (0.04 M_jup) to carve the gap (i.e. for 10 secular times to have elapsed within the system age). If there were several such planets then they could also sculpt the disc inner edge; alternatively, this could be done by a single planet with mass $\ge {0.44}_{-0.1}^{+0.07}{\text{M}}_{\text{Jup}}$ and semi-major axis $\le {41.6}_{-0.7}^{+1.0}\text{au}$.

B.28 HD 109085 (η Crv)

The ALMA data are from Marino et al. (2017b), and the disc model from REASONS (Matrà et al., in prep.). Marino et al. (2017b) predict a sculpting planet at 75 – 100 au with a mass of 3 – 30 M_⊕ (0.01 - 0.1 M_jup) and eccentricity < 0.08. We predict a $\ge {0.34}_{-0.1}^{+0.05}{\text{M}}_{\text{Jup}}$ sculpting planet at $\le {104}_{-2}^{+3}\text{au}$; our planet is larger because Marino et al. (2017b) use the shorter multi-planet sculpting timescale from Shannon et al. (2016) rather than the single-planet timescale from Pearce & Wyatt (2014).

B.29 HD 110058

The ALMA data are from Lieman-Sifry et al. (2016). The disc inner edge is unresolved but is less than 30 au, with best fit 10 au, so we use 10 au with 100% uncertainties. In scattered light, the disc beyond 40 au shows a strong warp, implying that the inner disc is misaligned with the outer disc (Kasper et al. 2015); if caused by a planet, then that planet should be aligned with either the inner or outer disc. Our predicted ≥ 0 5 ± 0 4 Mjup sculpting planet with semi-major axis ≤ 8 ± 8 au would yield a secular timescale equal to the system age at ~ 40 au, so this planet, if aligned with the inner disc, could be driving the warp.

B.30 HD 111520

The ALMA data are from Lieman-Sifry et al. (2016), and the disc model from REASONS (Matrà et al., in prep.). Draper et al. (2016) find a strong brightness asymmetry in scattered light, but argue against a high disc eccentricity.

B.31 HD 113337

One planet with M_p sin(i_p) = 3 1 ± 0 2 M_jup lies at 1 03 ± 0 02 au (Borgniet et al. 2014), and a candidate planet with M_p sin(i_p) = 7.2 ± 0.5 M_jup may lie at 4.8 ± 0.2 au (Borgniet et al. 2019). We take the candidate planet as the outermost planet; Fig. 15 shows that this planet could be stirring the disc, and is massive enough to have sculpted it (but not in situ). There is a probable misalignment between the planet(s) and the disc (Xuan et al. 2020), and an M-star companion is also present at 4000 au (Reid et al. 2007); the misalignment could be driven by the planet(s) and stellar companion (Xuan et al. 2020).

B.32 HD 114082

Wahhaj et al. (2016) predict a planet at ~ 25 au with mass ≤ 1 M_jup; this is similar to our predicted sculpting planet of mass ≥ 1.1 ± 0.2 Mjup and semi-major axis 22 ± 4au.

B.33 HD 115617 (61 Vir)

The ALMA and JCMT data are from Marino et al. (2017a). Three super-Earths or Neptunes with minimum masses 0.02 – 0.07 M_jup lie within 0.5 au (Vogt et al. 2010), although the outermost is not confirmed by HARPS (Wyatt et al. 2012). We consider the outermost (unconfirmed) planet at 0.476 ± 0.001 au, with M_p sin(i_p) = 23 ± 3 M_⊕ (Feng et al. 2017); Fig. 15 shows that this planet cannot be sculpting or stirring the disc in situ. Marino et al. (2017a) predict a stirring planet at 10 – 20 au with mass > 10 M_⊕ (0.03 M_jup), if its eccentricity is less than 0.1. This is similar to our stirring planet prediction (mass $\ge {0.02}_{-0.03}^{+0.04}{\text{M}}_{\text{Jup}}$), although the planet would have to be slightly more massive (≥ 0.06 ± 0.02 M_jup) to also sculpt the disc inner edge.

B.34 HD 120326

A Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE) image shows a ring at ~ 50 au and a possible second structure at 130 au (Bonnefoy et al. 2017); this scattered-light ring is slightly larger than the 31 ± 6 au SED-derived radius that we use.

B.35 HD 121191

The ALMA data are from Moór et al. (2017) and Kral et al. (2020), and the disc model from REASONS (Matrà et al., in prep.).

B.36 HD 121617

The ALMA data are from Moór et al. (2017), and the disc model from REASONS (Matrà et al., in prep.).

B.37 HD 123889

A companion lies at projected separation 200 au (Janson et al. 2013), whilst the disc is at 34 ± 8 au; it is possible that the two are interacting.

B.38 HD 127762 (γ Boo)

The disc is circumbinary, where the binary separation is 2 au (Yelverton et al. 2019). The binary is unlikely to be sculpting the disc at 130 ± 40 au.

B.39 HD 128311 (HNBoo)

At least two Jovian planets lie within 2 au, possibly in 2:1 resonance (Butler et al. 2003; Vogt et al. 2005). We consider the outermost at 1.8 ± 0.1 au with M_p sm(i_p) = 3.2 ± 0.3 M_Jup. Figure 15 shows that this planet could be stirring the disc, and is massive enough to sculpt it (but not in situ); this agrees with Moro-Martín et al. (2010) who show that the detected planets are unlikely to eject planetesimals outside 4 au (compared to the disc at 140 ± 70 au).

B.40 HD 131511 (DE Boo)

The disc is circumbinary, where the binary comprises 0.79 and 0.45 M_⊙ stars with a semi-major axis of 0.2 au (Jancart et al. 2005). The binary is aligned with the disc (Kennedy 2015), but is unlikely to be sculpting it (the disc is at ~ 100 au).

B.41 HD 131835

The ALMA data and disc model are from REASONS (Matrà et al., in prep.). In scattered light, the disc is wide and possibly comprises three concentric rings, with gaps at 46 – 57 au and 71 – 85 au (Feldt et al. 2017). Our sculpting planet predicted from ALMA data would be at ≤ 30 ± 2 au, which could be coincident with a tentative asymmetric inner ring seen in scattered light.

B.42 HD 138813

The ALMA data and disc model are from REASONS (Matrà et al., in prep.).

B.43 HD 139006 (α CrB)

The disc is circumbinary, where the binary has semi-major axis 0.2 au and eccentricity 0.37 (Tomkin & Popper 1986). Kennedy et al. (2012) show that the binary could stir the disc. An inner belt may also be resolved (Moerchen et al. 2010), but is not factored in to our Herschel single-disc model; our cold disc inner edge may therefore be underestimated.

B.44 HD 139664 (g Lup)

The ALMA data and disc model are from REASONS (Matrà et al., in prep.). The disc is possibly asymmetric in scattered light (Schneider et al. 2014, although also see Kalas et al. 2006).

B.45 HD 143894 (π Ser)

JCMT data may show a very large (560 au) structure, but could be confused by the background (Holland et al. 2017). This large size is comparable to our SED radius of 240 ± 70 au.

B.46 HD 146181

The ALMA data are from Lieman-Sifry et al. (2016). Their median ALMA model gives the inner edge as ${73}_{-19}^{+14}\text{au}$ with width < 50 au, whilst their best-fitting model gives the inner edge as 83 au with width 20 au. We therefore take a width of 20 au, and assume uncertainties of 50%.

B.47 HD 146897

The ALMA data are from Lieman-Sifry et al. (2016). Their median ALMA model gives the inner edge as ${60}_{-13}^{+11}\text{au}$ with width < 50 au, whilst their best-fitting model gives the inner edge as 71 au with width 20 au. We therefore take a width of 20 au, and assume uncertainties of 50%.

B.48 HD 156623

The ALMA data are from Lieman-Sifry et al. (2016). The inner edge is unresolved but < 20 au with best fit 10 au, so we use 10 au with 100% uncertainties.

B.49 HD 161868 (γ Oph)

The ALMA data and disc model are from REASONS (Matrà et al., in prep.).

B.50 HD 172555

A giant impact may have recently occurred in this system (Johnson et al. 2012; Schneiderman et al. 2021), so our dynamical models assuming a low-excitation debris disc may not apply.

B.51 HD 181296 (η Tel)

A brown dwarf is present beyond the disc, at 200 au (Lowrance et al. 2000); this could be interacting with the disc, which is at several tens of au.

B.52 HD 181327

The ALMA data and disc model are from REASONS (Matrà et al., in prep.). There is possibly some unresolved exterior emission from another ring or a broad population (Marino et al. 2016). The disc is asymmetric in scattered light (Stark et al. 2014). Nesvold & Kuchner (2015) predict a sculpting planet at 55 – 82 au, similar to our predicted sculpting planet with semi-major axis ≤ 53.9 ± 0.4 au and mass ≥ 1.6 ± 0.1 M_Jup.

B.53 HD 191089

The ALMA data and disc model are from REASONS (Matrà et al., in prep.). From scattered light images, Ren et al. (2019) give an upper limit on the ring eccentricity as 0.3. Churcher et al. (2011a) predict a ≥ 2.6 M_⊕ planet at 28 au stirs the disc, similar to our prediction of a ≥ 5 ± 2 M_⊕ (0.017 ± 0.006 M_Jup) stirring planet. If the disc is sculpted by a single planet, then Ren et al. (2019) predict a < 5 M_Jup planet at > 20 au, consistent with our prediction of a $\ge {0.89}_{-0.06}^{+0.07}{\text{M}}_{\text{Jup}}$ planet at ≤ 28 ± 1 au.

B.54 HD 193571 (k¹ Sgr)

The disc is circumbinary, where the binary is an A0 star with a 0.3 – 0.4 M_⊙ companion at ~ 11 au (Musso Barcucci et al. 2019). Our models show that this companion cannot sculpt the disc at 120 ± 30 au, but could stir it.

B.55 HD 196544 (ι Del)

The disc is circumbinary, where the binary has separation 0.2 au (Rodriguez & Zuckerman 2012). The binary is unlikely to be sculpting the disc at 100 ± 30 au.

B.56 HD 197481 (AU Mic)

The ALMA data are from Daley et al. (2019), and the disc model from REASONS (Matrà et al., in prep.). Two transiting Neptunes lie in the plane of the disc, inside 0.1 au (Plavchan et al. 2020; Martioli et al. 2021). The outermost planet has mass > 2 M_⊕ and is located at 0.110 ± 0.002 au (Martioli et al. 2021); this mass is a lower bound so we assume a 50 per cent uncertainty. Figure 15 shows that this planet cannot be stirring or sculpting the disc.

B.57 HD 202628

The ALMA data are from Faramaz et al. (2019) and show the disc to be narrow and eccentric (e = 0.09). Since Faramaz et al. (2019) use offset circles (rather than ellipses) to model the disc edges, for most of our calculations we approximate the edges as ellipses with eccentricity 0.09 and semi-major axes equal to the offset circle radii. However, this approximation fails when calculating the self-stirring disc mass from Eq. 16, so for that equation we take the semi-major axes of the innermost and outermost disc particles to be the radii of the inner and outer offset circles, respectively. Our predicted sculpting planet (mass ≥ 0.33 ± 0.06 M_Jup, semi-major axis ≤ 116 ± 2 au, eccentricity ≥ 0.09 ± 0.01) agrees with various literature predictions; Faramaz et al. (2019) predict mass 0.8 M_Jup, semi-major axis 110 au and eccentricity 0.1, Thilliez & Maddison (2016) predict mass 3 M_Jup, semi-major axis 100 au and eccentricity 0•2, and Nesvold & Kuchner (2015) predict mass < 15 M_Jup and semi-major axis 86 – 160 au.

B.58 HD 206860 (HN Peg)

A companion T dwarf (0.02 M_⊙, 20 M_Jup) lies at 800 au (Luhman et al. 2007). This is unlikely to be interacting with the disc at 32 ± 8 au.

B.59 HD 207129

The ALMA data and disc model are from REASONS (Matrà et al., in prep.). There are some hints at asymmetry in Herschel images (Löhne et al. 2012), but we use an axisymmetric ALMA model. If a single planet sculpts the disc then we predict that planet to have mass $\ge {0.21}_{-0.07}^{+0.03}{\text{M}}_{\text{Jup}}$ and semi-major axis ≤ 109 ± 3 au, which agrees with the prediction of Thilliez & Maddison (2016) that the planet is at 100 au (they also expect the planet to have eccentricity 0.06 based on the possible disc asymmetry). A companion lies at 900 au (Rodriguez et al. 2015), which could potentially interact with the disc at 130 – 170 au.

B.60 HD 213398 (β PsA)

A stellar companion lies at 1000 au (Rodriguez et al. 2015), which is unlikely to be interacting with the disc at 60 ± 30 au.

B.61 HD 216435 (τ¹ Gru)

A planet with M_p sin(i_p) = 1.3 ± 0.1 M_jup lies at 2.6 ± 0.2 au (Jones et al. 2003); Fig. 15 shows that this planet could be stirring the disc, and is massive enough to sculpt it (but not in situ).

B.62 HD 218396 (HR 8799)

The ALMA data are from Faramaz et al. (2021). We use their model 3, with a smoothly rising inner edge and smoothly falling outer edge, as expected from a disc containing a high eccentricity component (e.g. a scattered disc). Emission therefore extends both inwards and especially outwards of the radii we use, so our sculpting and stirring planet masses may be underestimates. A scattered population was also inferred by Geiler et al. (2019). Four giant planets are detected (Marois et al. 2008, 2010), with semi-major axes of 15 – 70 au and masses 6 – 7 M_jup (Wang et al. 2018). We consider the outermost planet at 70.8 ± 0.2 au with mass 5.8 ± 0.5 M_jup; Fig. 15 shows that this planet could be stirring the disc, and is massive enough to sculpt it (but not in situ). This is consistent with Booth et al. (2016), Read et al. (2018) and Geiler et al. (2019), who also argue that the outermost known planet cannot sculpt the disc in situ. Read et al. (2018) predict a fifth planet, showing that a 0.1 M_jup planet at 138 au could reproduce the disc morphology (although other planets are also possible). Our predicted sculpting planets are larger than this because our models assume that planets have cleared most unstable debris; however, unless the observed planets have since migrated, the presence of a significant scattered disc implies that debris clearing is ongoing in this system, in which case the planets could be less massive than our predictions. Planet detection limits were taken from Maire et al. (2015).

B.63 HD 223340

The target is component C in a quadruple system involving HD 223352 (which is component A). The target is separated from AB by 3100 au in projection (Phillips 2011).

B.64 HD 223352 (δ Scl)

The target is component A in a quadruple system involving HD 223340 (which is component C). The projected separation of AB is 164 au (Phillips 2011), so B could affect the disc around A (at 26 ± 8 au). B is two stars (De Rosa et al. 2011).

B.65 NLTT 54872 (Fomalhaut C)

The ALMA data are from Cronin-Coltsmann et al. (2021) and show the disc to be narrow and possibly eccentric $\left(e={0.04}_{-0.02}^{+0.03}\right)$. Since Cronin-Coltsmann et al. (2021) use an offset symmetrical torus (rather than ellipses) to model the disc edges, for most of our calculations we approximate the edges as ellipses with eccentricity 0.04 and semi-major axes equal to the radii of the offset torus edges. However, this approximation fails when calculating the self-stirring disc mass from Eq. 16, so for that equation we take the semi-major axes of the innermost and outermost disc particles to be the radii of the offset torus edges, respectively. Cronin-Coltsmann et al. (2021) argue that, if the disc edges are bound by the 2:1 and 3:2 resonances with a planet, then that planet would lie at 17 – 20 au. This is consistent with our prediction that a single sculpting planet would have semi-major axis $\le {21.6}_{-0.6}^{+0.5}\text{au}$ and mass ≥ 0.048 ± 0.003 M_jup (15.3 ± 0.9 M_⊕). The target is part of the Fomalhaut triple system.

References

All Tables

All Figures

Fig. 1 Cartoon showing the dynamical arguments used in this paper to constrain unseen planets from debris discs. At least one planet is assumed to have sculpted the inner edge of each disc (brown); this could be a single large planet (red circle) or multiple smaller planets (black circles). If the disc is eccentric, then we ascribe this to perturbations by a planet on an eccentric orbit (red line). The fact that we see dust implies that debris undergoes destructive collisions, which requires the disc to have been stirred; stirring may be performed by unseen planets (although other mechanisms are also possible). Planets can induce other debris morphologies in addition to those shown here, such as warps, resonant structures, and gaps; however, for this paper we only consider the above dynamical arguments.

In the text

Fig. 2 ALMA Band 6 (1.3 mm) image of HD 202628 (Faramaz et al. 2019), as an example of the ALMA data used in this paper (Sect. 2.2.1). North is up, and east is left. The debris disc is well resolved and clearly eccentric, with the host star (crosshairs) offset from the geometric centre of the disc. Considering the pericentre, qi, and apocentre, Qi, of the disc inner edge, plus equivalent values qo and Qo for the outer edge, yields the minimum mass of a single planet required to sculpt and stir the disc within the system lifetime (Sects. 3.1 and 4.2). This minimum-mass planet would have the orbit shown by the grey line; alternatively, a higher-mass planet located farther inwards could also sculpt and stir the disc. It is unclear whether the bright point in the north-east part of the ring is a background object or associated with the system (Faramaz et al. 2019).

In the text

Fig. 3 Herschel PACS 100 µm image of HD 1461, as an example of the Herschel data used in this paper (Sect. 2.2.2). North is up, and east is left. From left to right, the plots show the observation, observation after PSF subtraction, best-fitting disc model, and residuals after model subtraction. Contours show ±2 and ±3 times the root-mean-square intensity. Whilst the disc is resolved, it clearly has a lower resolution than the example ALMA image in Fig. 2. Our Herschel disc models are axisymmetric; the offset between the disc model and the origin reflects the imperfect Herschel pointing rather than a physical offset between the disc and the star.

In the text

	Fig. 4 Fitted debris disc parameters for systems where Herschel PACS data are used. The discs are modelled as axisymmetric with sharp inner and outer edges, as described in Sect. 2.2.2. The outer edges of these discs appear resolved (bottom plot), but the majority of the inner edges are unresolved (top plot).
In the text

Fig. 5 SED data for one example system, HD 37484 (Sect. 2.2.3). Black circles show photometric data, triangles are upper limits, and the grey shaded region is the Spitzer/InfraRed Spectrograph (IRS) spectrum. The best-fit model (black line) comprises the stellar flux (dotted green line) plus additional emission from two debris discs: a cold disc with blackbody temperature 60 ± 7 K (dashed blue line) and a warm disc with blackbody temperature 130 ± 10 K (dash-dot red line). This paper considers only the outermost (cold) discs.

In the text

Fig. 6 Ratios of resolved debris disc radii to blackbody radii for the systems in our sample, as described in Sect. 2.2.3. Red diamonds and yellow squares show systems where resolved disc radii were derived from ALMA and Herschel data, respectively. Each system is plotted once, using only the best available resolved data as described in Sect. 2.2. The lines show Γ correction factors; the solid line is derived from ALMA data (where we have limited Γ to a maximum of 4) and is that used in this paper (Pawellek et al. 2021), and the dashed line is from Herschel data (Pawellek 2017) as used by Launhardt et al. (2020). We note that the vertical axis uses the central radius of the resolved disc, defined as halfway between pericentre of the disc inner edge and the apocentre of the outer edge. The outlier on the left is the M-type Foma-lhaut C (NLTT 54872). The two systems in the bottom right of the plot may be complicated by either an underestimated cold disc location due to an additional warm component (β UMa: HD 95418; Moerchen et al. 2010) or background contamination (γ Oph: HD 161868; Holland et al. 2017).

In the text

Fig. 7 Masses and locations of the smallest planets needed to sculpt or stir the outer debris disc of 49 Cet (HD 9672), according to different models. The system hosts two debris discs; the outer disc has its resolved inner edge at 62 ± 4 au (line 1), and the unresolved inner disc is located around 12.2 ± 0.6 au (line 2). Lines 3a and 3b show 5σ detection limits; 3a is our NaCo-ISPY limit, and 3b is the SPHERE H2-band limit from Choquet et al. (2017) (with uncertainties derived as in Sect. 2.1.3). The red circle is the minimum-mass planet predicted by the Pearce & Wyatt (2014) model, which assumes that the outer disc is sculpted by a single planet (Sect. 3.1). Lines 4, 5, and 6 show the constraints from the Hill radius (Eq. (6)), diffusion time (Eq. (7)), and secular time (Eq. (8)), respectively, used to predict the properties of this single planet; if the outer disc is sculpted by one planet, then that planet is predicted to lie in the white region. The black circles denote the minimum-mass planets predicted by the alternative Shannon et al. (2016) model, which assumes that the outer disc is sculpted by multiple, equal-mass planets spanning the gap between the inner and outer discs (Eqs. (12) and (15) in Sect. 3.2). Line 7 is the minimum planet mass required to stir the outer disc (Eq. (23) in Sect. 4.2); we note that this is a separate analysis from the sculpting analysis because we do not necessarily assume that the disc is stirred by the sculpting planet(s). The 1σ uncertainties are shown as error bars around points and dark grey regions around lines. Similar analyses to those shown here are applied to all systems in our sample. We note that grain dynamics could be complicated by the presence of gas in this system (much less so than in a protoplanetary disc but still potentially enough to couple to dust; Moór et al. 2019). The PYTHON program used to produce this plot is publicly available for download8.

In the text

Fig. 8 Minimum masses of the planet(s) predicted to exist in our debris disc systems, for three different scenarios: a single planet sculpts the inner edge of the disc (top plot), multiple planets sculpt the inner edge of the disc (middle), or a single planet stirs the disc (bottom). The solid bars show predicted planet parameters in all of our systems, and the hatched bars show only systems where the disc parameters come from resolved images (ALMA or Herschel); the results for resolved and unresolved discs are similar for sculpting planets (top two plots) but differ for stirring planets because the results for unresolved discs are conservative lower bounds (Sect. 4.2.1). For the multi-planet scenario (middle plot), we note that the plot shows the mass of each individual planet.

In the text

Fig. 9 Masses and locations of the smallest planets required to sculpt our debris disc inner edges, assuming the planets lie interior to those discs (Sect. 3). Left plot: minimum planet masses required if each disc is sculpted by a single planet (Pearce & Wyatt 2014). More massive planets located farther inwards could also sculpt the discs to the same degree. Right plot: minimum individual planet masses required if each disc is instead sculpted by multiple, equal-mass planets (Shannon et al. 2016), where the outermost planet location is fixed at the disc inner edge. The single-planet model (left plot) requires larger planet masses to sculpt discs than the multi-planet model does (right plot), mainly because multiple planets are more efficient at removing debris. Symbols denote the disc data sources, where the best-available source has been used for each system; systems with ALMA data (red diamonds) are expected to yield the most reliable planet constraints, followed by those with Herschel data (yellow squares), then finally those with only unresolved SED data (black circles). The larger uncertainties on the multi-planet masses are driven by uncertainties on the disc inner edge locations; the multi-planet model has a stronger dependence on this than the single-planet model does.

In the text

Fig. 10 Detectability of the predicted disc-sculpting planets (Sect. 3). Horizontal axes show the minimum planet masses required to sculpt the discs if debris clearing is performed by one planet (left plot) or by multiple, equal-mass planets (right plot). The vertical axis shows sculpting planet masses divided by our L′-band detection limits at the apocentres of those planets (assuming the best-case scenario, where the systems are face-on to the observer). The dashed lines show unity. Symbols have the same meanings as in previous figures. If the discs are each sculpted by one planet (left plot), then some of these planets will lie at or close to the current detection limits, and a factor of 10 improvement in mass detection limits could yield many more predicted planets. Conversely, if instead each disc is sculpted by multiple planets (right plot), then these planets could be significantly smaller and many might not be observable even with significant improvements in detection limits.

In the text

Fig. 11 Self-stirring constraints in systems with resolved debris discs. The vertical axis shows the minimum masses of discs required for those discs to be self-stirred (Sect. 4.1). Horizontal axes show the minimum planet masses required to sculpt the inner edges of those discs (Sect. 3) if debris clearing is performed by one planet (left plot) or by multiple, equal-mass planets (right plot). Symbols have the same meanings as in previous figures. Background colours denote debris disc masses that may be unfeasibly large (red), borderline (yellow), or allowed (green), according to the quantity of solids that may have been available at the protoplanetary disc stage (Krivov & Wyatt 2021). It may not be possible for discs in the red and yellow regions to be self-stirred, in which case other stirring mechanisms are required. Since these discs also typically require large planets to sculpt their inner edges, it is likely that planet-stirring operates in these systems, as shown in Sect. 4.2. The requirement for large planets to both stir and sculpt these discs makes systems in the yellow and red regions prime targets for future planet searches, as discussed in Sect. 5.2.

In the text

Fig. 12 Constraints on the planet masses required to stir debris discs. Horizontal axes show the minimum planet masses required for stirring if those planets have eccentricities smaller than 0.3, as described in Sect. 4.2 (higher eccentricities would require lower masses). The vertical axis definitions are as follows. Top plot: minimum stirring-planet masses, divided by our L′-band planet detection limits at the apocentres of those planets. The minimum stirring-planet masses are all consistent (to 1σ) with lying below the detection limits (i.e. below the dashed line, which shows unity), so all of the discs may potentially be stirred by unseen planets. The stirring planet predictions for systems where discs are unresolved (black circles) are conservative lower limits since these discs are modelled as having zero width. Bottom plot: minimum masses of the resolved discs required for those discs to be self-stirred, from Fig. 11. Many systems would require disc masses greater than 1001000 M⊕ to self-stir, which may be unfeasibly high (Krivov & Wyatt 2021); however, the top plot shows that these systems could still be stirred by unseen planets. This suggests that self-stirring is not the dominant stirring mechanism in at least some debris discs. No unresolved discs are shown on the bottom plot because they were omitted from the self-stirring calculation owing to a lack of width information. Also, we note that a few systems on the bottom plot are missing from the top plot because their predicted planet locations are not covered by our observations. The uncertainties on stirring planet masses are large owing to the very strong dependence of Eq. (23) on disc edge locations.

In the text

Fig. 13 Comparison of the planet masses required to stir and sculpt debris discs. The vertical axis shows the minimum planet masses required to stir each disc (Sect. 4.2). Horizontal axes show the minimum planet masses required to sculpt the inner edges of those discs (Sect. 3) if sculpting is performed by one planet (left plot) or by multiple, equal-mass planets (right plot). Symbols are the same as in previous figures, and the dashed lines denote 1:1 relations. If one planet sculpts each disc inner edge (left plot), then for almost all systems that planet could also stir the disc (either the sculpting planet mass is greater than that required for stirring, or the two are consistent to 1σ), The only exception is β Tri (HD 13161), which is a binary system where the binary could instead stir the disc (for the other labelled system, β Leo (HD 102647), we likely underestimate the disc inner edge and therefore overestimate the required stirring planet mass). Conversely, if multiple planets sculpt each disc (right plot), then for several systems the minimum-mass sculpting planets would be too small to stir those discs (to 1σ), some, but not all, of those discs could feasibly be self-stirred, so either the planets in those systems are more massive than required by the multi-planet sculpting model, or disc sculpting and stirring are performed by different entities. We note that the stirring planet masses derived from SEDs of unresolved discs (black circles) are conservative lower limits since these discs were assumed to have zero width for this analysis; in this case, the planet masses required for sculpting and stirring have similar or identical dependences on system parameters, resulting in the tight relationships between the two. As on previous plots, the large uncertainties are due to the strong dependences of Eqs. (15) and (23) on disc edge locations.

In the text

Fig. 14 Planets predicted from our debris analyses, compared to known planets and those inferred to be forming in protoplanetary discs (Sect. 5.1.1). Small black circles are known exoplanets from Exoplanet.eu, and large circles are Solar System planets. The left plot shows the planetary populations inferred from our debris analyses; green diamonds and blue crosses are the planets required to sculpt the discs if sculpting is performed by one or multiple planets, respectively (Sects. 3.1 and 3.2). Red circles are planets required to stir the discs (only systems with resolved discs are shown; Sect. 4.2). The plot shows that debris disc morphologies imply planetary populations that are poorly explored with current instruments, so the hypotheses that planets sculpt and stir debris discs are consistent with the lack of planet detections in these systems. Right plot: Comparison of our planets inferred from debris discs to those inferred to be forming in protoplanetary discs (black crosses; Lodato et al. 2019). To make the plot clearer, we have replaced our individual planet predictions with rough areas representing the three populations. Our planetary population inferred from debris discs well matches the population inferred to be forming in protoplanetary discs; if the forming planets remain in situ, then they could later join the population responsible for sculpting and stirring known debris discs. Alternatively, it is possible that the forming planets rapidly sculpt debris discs before migrating inwards. Error bars are omitted for clarity but can be seen in our other figures.

In the text

Fig. 15 Comparison of detected and predicted planets in our debris disc systems, as described in Sect. 5.1.2. For systems with a detected planet, the vertical axis shows the mass of one detected planet divided by the mass required by a single planet to sculpt the disc inner edge (Sect. 3.1). We calculate the masses of RV planets by assuming that they are copla-nar with the discs. The horizontal axis shows the semimajor axis of the detected planet, divided by the semimajor axis that the minimum-mass sculpting planet would have. For systems with multiple known planets, only the outermost detected planet is shown. Dashed lines denote unity. With the possible exceptions of 51 Eri (HD 29391) and HD 69830, the detected planets are unlikely to be massive enough to sculpt the disc inner edges in situ; whilst some are larger than the minimum masses required for sculpting, these detected planets lie far interior to the discs and so would need significantly greater masses to sculpt the discs in situ (line 4 on Fig. 7). This means that either additional undetected planets also reside in these systems, or the known planets formed in the outer regions and rapidly sculpted the discs before migrating inwards. Conversely, some detected planets are massive enough to have stirred the discs in situ, despite their smaller semimajor axes; systems where detected planets could stir the discs are circled. Arrows show lower limits on detected planet masses for RV planets with no disc inclination data, and the remaining symbols were defined on previous figures.

In the text

Fig. 16 Identification of our predicted planets with the best chances of future detection (Sect. 5.2). Points show the minimum masses and maximum projected separations of the planets predicted to be sculpting the inner edges of our discs if one planet is responsible in each system (Sect. 3.1). Only systems where the disc is resolved by ALMA (diamonds) or Herschel (squares) are shown. Left plot: systems coloured by age. Right plot: systems coloured by disc inclination to the sky. The best targets have predicted planets with high masses, which are well separated from their stars, are young (so the planets are brightest), and have low-inclination discs (so a coplanar planet would be separated from its star in projection over a large fraction of its orbit). JWST will be capable of imaging 0.1−1 MJup planets at approximately arcsecond separations (Beichman et al. 2010; Carter et al. 2021), and this should be complemented at small angular separations by upcoming instruments such as VLT/ERIS (Davies et al. 2018), so many of these predicted planets could be observable in the near future. The plots show the single planets required to sculpt each disc, but if multiple planets were instead responsible (Sect. 3.2), then those planets could be at least an order of magnitude less massive than those shown here. The planets required to stir each disc are typically smaller than those shown (Sect. 4.2). We note that some systems are missing from the right plot because they lack inclination information.

In the text