# Insights into the planetary dynamics of HD 206893 with ALMA

S. Marino, A. Zurlo, V. Faramaz, J. Milli, Th Henning, G. M. Kennedy, L.

Matrà, S. Pérez, P. Delorme, L. A. Cieza, et al.

► To cite this version:

S. Marino, A. Zurlo, V. Faramaz, J. Milli, Th Henning, et al. Insights into the planetary dynamics of HD 206893 with ALMA. Monthly Notices of the Royal Astronomical Society, 2020, 498, pp.1319-1334. 10.1093/mnras/sta2386. insu-03705173

HAL Id: insu-03705173

https://insu.hal.science/insu-03705173

Submitted on 27 Jun 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Insights into the planetary dynamics of HD 206893 with ALMA

S. Marino<sup>1,★</sup> A. Zurlo<sup>2</sup> V. Faramaz, <sup>3</sup> J. Milli, <sup>4</sup> Th. Henning, <sup>1</sup> G. M. Kennedy<sup>5,6</sup> L. Matrà<sup>7</sup>

S. Pérez<sup>8</sup> P. Delorme, <sup>4</sup> L. A. Cieza<sup>2</sup> and A. M. Hughes<sup>9</sup>

Accepted 2020 August 6. Received 2020 July 24; in original form 2020 June 12

## ABSTRACT

Radial substructure in the form of rings and gaps has been shown to be ubiquitous among protoplanetary discs. This could be the case in exo-Kuiper belts as well, and evidence for this is emerging. In this paper, we present ALMA observations of the debris/planetesimal disc surrounding HD 206893, a system that also hosts two massive companions at 2 and 11 au. Our observations reveal a disc extending from 30 to 180 au, split by a 27 au wide gap centred at 74 au, and no dust surrounding the reddened brown dwarf (BD) at 11 au. The gap width suggests the presence of a $0.9M_{\text{Jup}}$ planet at 74 au, which would be the third companion in this system. Using previous astrometry of the BD, combined with our derived disc orientation as a prior, we were able to better constrain its orbit finding it is likely eccentric ($0.14^{+0.05}_{-0.04}$). For the innermost companion, we used radial velocity, proper motion anomaly, and stability considerations to show its mass and semimajor axis are likely in the ranges 4–100$M_{\text{Jup}}$ and 1.4–4.5 au. These three companions will interact on secular time-scales and perturb the orbits of planetesimals, stirring the disc and potentially truncating it to its current extent via secular resonances. Finally, the presence of a gap in this system adds to the growing evidence that gaps could be common in wide exo-Kuiper belts. Out of six wide debris discs observed with ALMA with enough resolution, four to five show radial substructure in the form of gaps.

**Key words:** methods: numerical-techniques: interferometric-planets and satellites: dynamical evolution and stability-circumstellar matter-stars: individual: HD 206893-submillimetre: planetary systems.

## 1 INTRODUCTION

The study of exoplanetary systems has been revolutionized in the last decade with the discovery of thousands of exoplanets and several hundreds of debris discs (analogous to the Kuiper belt), evidenced by short-lived dust that is being replenished via collisions among an underlying population of planetesimals (see reviews by Wyatt 2008; Hughes, Duchêne & Matthews 2018). Some of these systems are known to host both exoplanets and *exo-Kuiper* belts, allowing for a more detailed characterization of their architecture, dynamics, and formation since they provide complementary information (e.g. Moro-Martín et al. 2007, 2010).

As the number of known systems hosting both planets and exo-Kuiper belts grew, studies have tried to find correlations between the two. Some have provided tentative evidence of a possible higher occurrence rate of debris discs (indicative of more mass in the form of planetesimals) in systems hosting low-mass planets detected through radial velocities (RVs; Wyatt et al. 2012; Marshall et al.

2014; Moro-Martín et al. 2015), typically located within 1 au and with discs at tens of au (e.g. Kennedy et al. 2015b; Marino et al. 2017b), but this trend has been recently shown to be not significant (Yelverton, Kennedy & Su 2020). On the other hand, there seems to be an anticorrelation between the presence of massive close-in planets (or stellar metallicity) and detectable debris discs (Greaves et al. 2004; Moro-Martín et al. 2007). More recently, Meshkat et al. (2017) also showed that systems with bright debris discs seem to be more likely to have planets at least a few times more massive than Jupiter at separations of 10–1000 au, where planets and debris generating planetesimals could interact. The origin for these tentative correlations is still unclear and it is likely that many factors during the planet formation process and subsequent dynamical evolution contribute to these.

One way to improve our understanding is to look in detail how planets and debris discs interact. Thanks to ALMA, it has been possible to image tens of debris discs at millimetre wavelengths, typically tracing mm-sized grains unaffected by radiation (Burns, Lamy & Soter 1979) or gas drag forces (e.g. Marino et al. 2020, and references therein), thus tracing the spatial distribution of the parent km-sized planetesimals. ALMA images have revealed at

<footnote n="1"><sup>1</sup>Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="2"><sup>2</sup>Núcleo de Astronomia, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejercito 441, Santiago, Chile</footnote>

<footnote n="3"><sup>3</sup>Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA</footnote>

<footnote n="4"><sup>4</sup>IPAG, Univ. Grenoble Alpes, CNRS, F-38000 Grenoble, France</footnote>

<footnote n="5"><sup>5</sup>Department of Physics, University of Warwick, Gibbet Hill Road, Coventry CV4 7AL, UK</footnote>

<footnote n="6"><sup>6</sup>Centre for Exoplanets and Habitability, University of Warwick, Gibbet Hill Road, Coventry CV4 7AL, UK</footnote>

<footnote n="7"><sup>7</sup>School of Physics, National University of Ireland Galway, University Road, Galway, Ireland</footnote>

<footnote n="8"><sup>8</sup>Departamento de Física, Universidad de Santiago de Chile, Av. Ecuador 3493, Estación Central, Santiago, Chile</footnote>

<footnote n="9"><sup>9</sup>Astronomy Department and Van Vleck Observatory, Wesleyan University, Middletown, CT 06459, USA</footnote>

<footnote>* E-mail: sebastian.marino.estay@gmail.com</footnote>
unprecedented detail asymmetric structures (e.g. β Pic, Fomalhaut, and HD 202628; Dent et al. 2014; MacGregor et al. 2017; Faramaz et al. 2019), annular gaps (HD 107146, HD 92945, and HD 15115, Ricci et al. 2015; Marino et al. 2018b, 2019; MacGregor et al. 2019), and vertical substructure (e.g. β Pic; Matrà et al. 2019), suggesting the presence of as yet unseen low-mass planets.

While most systems with exo-Kuiper belts do not have known planetary mass companions, in a few of these it has been possible to directly image one, thus enabling the study of planet–disc interactions in more detail. There are well-known examples such as β Pic with a massive planet possibly warping the disc (Mouillet et al. 1997; Lagrange et al. 2012, 2019; Matrà et al. 2019); HR 8799 with four giant planets creating a scattered disc and possibly replenishing its warm dust closer in (e.g. Marois et al. 2010; Booth et al. 2016; Zurlo et al. 2016; Read et al. 2018; Wilner et al. 2018; Geiler et al. 2019; Faramaz et al., in preparation); HD 95086's axisymmetric disc implying a low eccentricity of its $4M_{Jup}$ planet (Rameau et al. 2016; Su et al. 2017); and Fomalhaut having a narrow and eccentric planetesimal belt (Kalas, Graham & Clampin 2005; Acke et al. 2012; Boley et al. 2012; MacGregor et al. 2017), implying that its candidate companion on an eccentric orbit has a low mass (~Earth or super-Earth) and is not sculpting the belt (Quillen 2006; Kalas et al. 2008; Chiang et al. 2009; Beust et al. 2014; Faramaz et al. 2015), or is not a compact object but rather the dusty aftermath of a recent planetesimal collision (Gaspar & Rieke 2020). Some exo-Kuiper belt host systems even have companions in the brown dwarf (BD) or low stellar mass regime, suggesting that their likely formation through gravitational instability (Boss 1997, 2003, 2011; Vorobyov 2013) is compatible with the formation of massive Kuiper belt analogues, e.g. HR 2562 (Konopacky et al. 2016), HD 193571 (Musso Barcucci et al. 2019), HD 92536 (Launhardt et al. 2020), and HD 206893 (Milli et al. 2017). The last one is the subject of this paper. For even more massive companions, Yelverton et al. (2019) found a significant lower detection rate of debris discs around binaries, with no discs detected in binaries with separations between 25 and 135 au (comparable to typical debris disc radii; Matrà et al. 2018b). This is likely due to dynamical perturbation inhibiting planetesimal formation or clearing any debris disc formed near those separations.

Located at 40.8 pc (Gaia Collaboration 2018), the F5V star HD 206893 is known to host a companion, HD 206893 B, at a separation of ~11 au (Milli et al. 2017) and a debris disc (Moór et al. 2006; Chen et al. 2014) that was marginally resolved with Herschel (Milli et al. 2017). Given the estimated age of this system of 50–700 Myr, the companion mass is probably in the range 12–50$M_{Jup}$ (Delorme et al. 2017), and thus it is likely a BD. Astrometric follow-up of this companion using VLT/NACO and SPHERE placed some constraints on the period (or semimajor axis) and orientation of the orbit, and set an upper limit to the eccentricity of ~0.5 (Grandjean et al. 2019). The same study using Hipparcos (van Leeuwen 2007) and Gaia Data Release 2 (DR2) data (Gaia Collaboration 2018) revealed a significant proper motion anomaly (PMa) in a direction that cannot be explained by the BD, suggesting the presence of an additional companion closer in (also confirmed in Kervella et al. 2019). This additional companion would also be responsible for an observed RV drift (or stellar acceleration), indicating that this inner companion (HD 206893 C) must have at least a mass of ~15$M_{Jup}$ and orbit with a semimajor axis between 1.4 and 2.6 au (Grandjean et al. 2019). These two companions make this system an ideal target to image its disc with ALMA, in order to both constrain the dynamics of this system and look for additional companions that could shape the distribution of planetesimals. This paper is structured as follows. In Section 2, we present our new ALMA observations

that show evidence of a gap. We then fit these observations using a parametric disc model in Section 3 to constrain the disc orientation and radial structure. Using these constraints, particularly the disc orientation, in Section 4 we improve the previous orbital constraints of HD 206893 B by assuming it is co-planar with the disc. In Section 5, we discuss and summarize the different constraints on companions in this system and potential origins of the gap. Finally, in Section 6 we summarize our results and conclusions.

## 2 OBSERVATIONS

We observed HD 206893 with ALMA in band 7 (average wavelength and frequency of 0.88 mm and 342 GHz, respectively) as part of the cycle 5 project 2017.1.00828.S (PI: A. Zurlo). Observations were taken with both the Atacama Compact Array (ACA) and the main 12-m array, in order to recover the large-scale structure up to sizes of 20 arcsec and the small-scale structures down to 0.3 arcsec, respectively. Details about these observations are summarized in Table 1. The correlator was set up using four spectral windows to study primarily the dust continuum emission in the system. Three of these were centred at 348.4, 336.5, and 334.6 GHz, with a bandwidth of 2 GHz and a channel width of 15.6 MHz. The fourth spectral window was centred at 348.4 GHz and had a bandwidth of 1.875 GHz, with a narrower channel width of 0.488 MHz (0.42 km s⁻<fnref n="1" /> and effective bandwidth of 1.1 km s⁻<fnref n="1" />)<fnref n="1" /> to search for a serendipitous CO J = 3 – 2 detection. The data were calibrated using CASA 5.4 (McMullin et al. 2007) and the standard calibration routines provided by ALMA. Additionally, we applied a phase centre shift to the ACA data that were not correctly centred on HD 206893 at the corresponding epoch according to its Gaia DR2 astrometry. To complement these band 7 observations (both continuum and CO emission), we also retrieved archival band 6 observations (1.3 mm, 222 GHz) that we calibrated using the standard ALMA routines for CASA. These observations are described in Nederlander et al. (2020). Below we present the analysis of continuum and CO observations.

### 2.1 Dust continuum

Continuum images at wavelengths of 0.88 and 1.3 mm are obtained using the TCLEAN task in CASA with Briggs weighting and a robust parameter of 2.0 (to maximize sensitivity). These clean images are presented in Fig. 1. The images are smoothed using a Gaussian tapering<fnref n="2" /> of 0.9 arcsec at 0.88 mm and 0.4 arcsec at 1.3 mm, which leads to a loss of resolution, but allows for an increase of the signal-to-noise ratio (S/N) per beam. The beam size in the tapered images is 1.02 arcsec × 0.91 arcsec at 0.88 mm and 0.92 arcsec × 0.80 arcsec at 1.3 mm. Disc emission is detected at both wavelengths within 4 arcsec (160 au) of the star, distributed over a wide range of radii. At the centre of the images, the star is significantly detected (more clearly seen in non-tapered images; see Section 2.1.2). The stellar flux is consistent with Rayleigh–Jeans extrapolations of its flux at shorter wavelengths (i.e. 30 µJy at 0.88 mm and 13 µJy at 1.3 mm), and thus we attribute it to photospheric emission. We estimate a total integrated flux of 2.68 ± 0.36 and 1.05 ± 0.12 mJy at 0.88 and 1.3 mm, respectively (including 10 per cent absolute calibration uncertainties). These fluxes are computed by integrating all emission within an ellipse of semimajor axis of 5 arcsec (~200 au) and oriented as the disc on the sky (PA = 61° and i = 40°; see Section 3.1). In

<footnote n="1">¹https://help.almascience.org/index.php?/Knowledgebase/Article/View/29</footnote>

<footnote n="2">²This is done in the UV space through the TCLEAN argument uvtaper.</footnote>
<table><caption>**Table 1.** Summary of band 7 (12 m and ACA) and band 6 (12 m) observations. The image rms and beam size correspond to Briggs weighting with a robust parameter of 2.0.</caption><thead><tr><th>Observation</th><th>Dates</th><th>t<sub>sci</sub><br>(h)</th><th>Image rms<br>(μJy)</th><th>Beam size (PA)</th><th>Minimum and maximum<br>baselines (m)<br>(5th and 95th percentiles)</th></tr></thead><tbody><tr><td>Band 7: 12 m</td><td>2018 September 26, 28, and 30; 2019 April 20</td><td>4.3</td><td>9.1</td><td>0.31 arcsec × 0.25 arcsec (−82°)</td><td>55 and 733</td></tr><tr><td>Band 7: ACA</td><td>2017 October 23</td><td>7.9</td><td>110</td><td>5.0 arcsec × 2.7 arcsec (88°)</td><td>9 and 45</td></tr><tr><td></td><td>2018 March 10, April 5 and 7, May 8, 14–18, and 26</td><td></td><td></td><td></td><td></td></tr><tr><td>Band 6: 12 m</td><td>2018 June 27, August 30, September 10 and 17</td><td>4.3</td><td>4.9</td><td>0.70 arcsec × 0.57 arcsec (67°)</td><td>45 and 670</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Continuum clean images at 0.88 mm (12 m + ACA, left-hand panel) and 1.3 mm (right-hand panel) of HD 206893 obtained using Briggs weighting and a robust parameter of 2. Additionally, we applied an ultraviolet (UV) tapering of 0.9 arcsec to the band 7 data and 0.4 arcsec to the band 6 data. The images are also corrected by the primary beam; hence, the noise increases towards the edges. The contours represent three and five times the image rms (17 and 5.3 μJy per beam at the centre of the band 7 and 6 images, respectively). The stellar position is marked with a white plus symbol near the centre of the image (based on Gaia DR2) and the beams are represented by white ellipses in the bottom-left corners (1.02 arcsec × 0.91 arcsec and 0.92 arcsec × 0.80 arcsec, respectively).</figcaption></figure>

both band 6 and band 7 maps, there is evidence of extended emission arising near the star suggesting that the planetesimal disc is wide. The detailed radial structure is analysed in Section 2.1.1.

### 2.1.1 Dust radial structure

In order to study the radial structure of the disc, we azimuthally average the deprojected emission (as in Marino et al. 2016) using the best-fitting disc position angle and inclination (see Section 3). Both band 6 and band 7 profiles are shown in Fig. 2. For this process, we use the band 6 clean image without tapering and the band 7 image with a 0.4 arcsec tapering. This choice results in a similar beam at both wavelengths (0.57 arcsec × 0.50 arcsec at 0.88 mm and 0.70 arcsec × 0.57 arcsec at 1.3 mm) and is a good compromise between spatial resolution and S/N (see Fig. A1 in Appendix A). Based on these profiles, the disc emission is detected from 30 to 180 au, with a peak near 110 au and a local minimum at roughly 70 au. This minimum hints at the presence of a gap at a similar radial distance

compared to HD 107146, HD 92945, and HD15115 (72 ± 3, 73 ± 3, and 59 ± 5 au, respectively; Marino et al. 2018b, 2019; MacGregor et al. 2019). The gap seems to be deeper at 0.88 mm, but this is likely due to the lower resolution at 1.3 mm, which does not resolve well this minimum. Interior to the gap, the disc intensity peaks at around 40 au, but its exact inner edge is uncertain. Note that the emission interior to 20 au is simply consistent with photospheric emission from the star convolved with the beam. Moreover, the disc is not expected to extend interior to 15 au since it would be truncated by HD 206893 B’s chaotic zone if it is on a circular orbit and has a mass of $12M_{Jup}$. Where exactly it should be truncated is uncertain since B could be more massive (up to $50M_{Jup}$) and/or on an eccentric orbit (<0.5). Moreover, secular interactions between the two inner companions could have depleted the disc at regions between 20 and 40 au via secular resonances (see discussion in Section 5.2.2).

Since some debris discs are known to be asymmetric (e.g. Dent et al. 2014; MacGregor et al. 2017; Faramaz et al. 2019; Marino et al. 2019), and expected to be so when interacting with planets on eccentric orbits (e.g. Pearce & Wyatt 2014; Regály et al. 2018),
<figure><img src="image_4.png" /></figure>

**Figure 2.** Azimuthally averaged surface brightness profile computed by deprojecting the emission according to our best-fitting model (Section 3.1) and azimuthally averaging the emission (top). The shaded regions correspond to $1\sigma$ uncertainties. Note that the shaded regions are representative of the uncertainty over a resolution element, i.e. 23 and 29 au for band 7 and band 6, respectively. The vertical dashed line represents the semimajor axis of B and the grey region its chaotic zone if on a circular orbit. The bottom panel shows the difference between the north-east (NE) and south-west (SW) halves of the disc in significance levels (i.e. the difference is divided by the local uncertainty).

we perform two tests in order to search for asymmetries. We first compare the integrated flux of one-half of the disc against the other while varying the angle of the axis that divides the two halves. We find that when comparing the NE and SW halves (divided by the disc minor axis), this difference is maximized and marginally significant ($2\sigma$) at both wavelengths. Averaging both wavelengths, we find that the NE side is $30 \pm 13$ per cent brighter. This difference suggests that the disc could be either asymmetric or instead the measured flux is contaminated by a background submillimetre galaxy (SMG) as found in other ALMA observations of debris discs (e.g. Marino et al. 2017b; Zapata, Ho & Rodríguez 2018). In order to check the radial location where this difference is strongest, we compute azimuthally averaged radial profiles of the two disc halves and subtract them (bottom panel of Fig. 2). We find that the NE side is overall brighter at almost all radii, especially inside the gap and at 140 au, although these differences are not larger than $3\sigma$. They do reveal nevertheless that the flux difference does not arise from a single compact SMG (typically smaller than 0.5 arcsec or 20 au at HD 206893's distance; Simpson et al. 2015; Lindroos et al. 2016; Fujimoto et al. 2017), but rather from a much broader region. In Section 5.2, we will discuss what could be the origin of these asymmetries and how they could be connected with the formation of the gap and the orbit of HD 206893 B.

In order to further constrain the radial structure, we use the PYTHON module FRANKENSTEIN (FRANK; Jennings et al. 2020) that uses Gaussian process to reconstruct the intensity radial profile of a disc. We first subtract the stellar emission from the visibilities by simply subtracting a constant value of 30 µJy at 0.88 mm and 13 µJy at 1.3 mm.<fnref n="3" /> Then, using the derived disc orientation from Section 3.1, we deproject the visibilities and use FRANK to reconstruct the radial

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

**Figure 3.** Disc surface brightness profile computed using FRANK and deprojecting the visibilities according to our best-fitting model (Section 3.1). The shaded regions correspond to $1\sigma$ uncertainties. The vertical dashed line represents the semimajor axis of B and the grey region its chaotic zone if on a circular orbit.

profile at 0.88 and 1.3 mm. Fig. 3 presents the derived profiles using a maximum radius of 400 au and hyperparameters $\alpha = 1.01$ and $w_{\text{smooth}} = 10^{-2}$. These profiles confirm the presence of a gap centred around 70 au, which appears slightly deeper than that in the profiles derived in the image space due to the higher resolving power of FRANK. These profiles also show a clearer inner cavity of $\sim 20$ au mainly because of the subtraction of the stellar emission. A caveat in the derivation of these profiles is that no primary beam correction is taken into account, although its effect is smaller than 10 per cent within 150 au in both band 7 and band 6 data.

## 2.1.2 BD dust emission upper limit

An additional goal of these observations was to constrain the amount of dust surrounding HD 206893 B. Delorme et al. (2017) presented new photometric and spectroscopic measurements with SPHERE, and pointed out that the BD companion is a peculiar object. Neither empirical models combined with absorption by forsterite nor synthetic dusty spectra can describe its very red colour. Since the debris disc is outside the orbit of the companion and it has a very low optical depth, it cannot explain the reddening and extinction of the BD. A possible explanation is that the red colour of this object is produced by the extinction from a circumplanetary disc. At near-infrared (NIR) wavelengths, this small disc is obviously unresolved since diffraction-limited observations lead to 2.0 au resolution at 40 pc in the case of SPHERE/IRDIS (K band). However, it could be detected in thermal emission at longer wavelength.

We do not detect any emission arising from a point source towards the NE at the BD separation of $\sim 0.25$ arcsec (see Fig. 4), where the BD is expected to be (Grandjean et al. 2019). Nevertheless, we can place $3\sigma$ upper limits at 0.88 and 1.3 mm of 27 and 15 µJy, respectively. Using the same radiative transfer tools as in Pérez et al. (2019b), we estimate a dust mass upper limit of $2 \times 10^{-4} M_{\oplus}$ ($2 \times 10^{-2} M_{\text{moon}}$), which is equivalent to a dust-to-planet ratio of $5 \times 10^{-8}$ for a disc with a size $\gtrsim 0.1$ au. A much smaller disc would be optically thick and could hide a higher mass remaining undetected by ALMA. However, this is an unlikely scenario because such a dense and compact dusty disc would quickly become depleted in dust due to collisional evolution and loss processes (e.g. radiation pressure).

<footnote n="3">³The Fourier transform of a point source at the phase centre is a positive and real constant with a value equal to its total flux.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Zoomed clean image at 0.88 mm using Briggs weighting and a robust parameter of 2. The contours represent three and four times the rms (9 µJy per beam). The stellar position is marked with a black plus symbol near the centre of the image (based on Gaia DR2) and the BD position is marked with a black multiplication symbol. The beam is represented by a white ellipse in the bottom-left corner (0.30 arcsec × 0.24 arcsec).</figcaption></figure>

This non-detection of dust surrounding HD 206893 B does not rule out that B could have a massive circumstellar disc in the form of satellites, but that is very depleted in dust. Even if those satellites were capable of producing dust through collisions, for HD 206893 B that dust would hardly be detectable according to theoretical predictions for the collisional evolution of satellite swarms (Kennedy & Wyatt 2011). In such swarms, dust is continually produced via collisions of km-sized satellites (analogous to circumstellar debris discs) and lost due to radiation pressure or Poynting–Robertson drag, in which case it is accreted on the companion. Given HD 206893’s age (>10 Myr) and distance (40.8 pc), the models by Kennedy & Wyatt (2011) predict fluxes below 1 µJy at 1 mm (see their fig. 4 for the case of a system at 10 pc), and thus we cannot rule out that HD 206893 B hosts a massive satellite swarm.

Another possibility to explain the reddening is that it is caused by accretion of dust from the circumstellar debris disc. While the planetesimal belt resides beyond the orbit of the BD, it is still possible that dust is being accreted by the BD, potentially explaining its reddened spectrum. Small dust grains are expected to migrate in through the Poynting–Robertson drag interior to debris discs (e.g. Burns et al. 1979; Wyatt et al. 2005; van Lieshout et al. 2014; Kennedy et al. 2015a), and interact with intermediate planets. Such interactions can lead to trapping in resonance, ejections, and accretion on to the planets (Shannon, Mustill & Wyatt 2015; Bonsor et al. 2018). Using equations (21) and (22) in Bonsor et al. (2018),<fnref n="4" /> we can estimate the rate at which the BD could be accreting dust.

We find that the BD could be accreting µm-sized dust at a rate of ~4 kg s⁻¹. A simple calculation can show that this level of dust accretion cannot lead to the needed reddening, by comparing the total cross-sectional area accreted over the age of this system versus the surface area of the BD. Assuming the accreted dust is all made of grains of 1 µm in radius, has a density of 3 g cm⁻³, and has been accreted for 700 Myr, we find the equivalent dust layer would only cover 0.001 per cent of the surface of the BD, and thus insufficient to produce reddening or extinction at the necessary levels (A<sub>K</sub> ~ 0.5; Delorme et al. 2017). Note that to produce reddening at NIR larger grains might be necessary and hence the total cross-section smaller. Moreover, we neglect here any mixing in the atmosphere of the BD that would tend to remove even more of this dust below the photosphere, by settling. Therefore, we conclude that if dust in the atmosphere is indeed responsible for its reddening, it cannot be supplied in sufficient quantity by the circumstellar disc.

These findings suggest that the physical conditions within the atmosphere of HD 206893 B allow to lift dust (or to prevent it to condense) inside and above its photosphere in quantities that are at least in the highest ranges that were considered as possible by substellar atmosphere models. Indeed, Delorme et al. (2017) showed that the spectra of HD 206893 B could not be correctly fitted without external extinction, unless the dust settling parameters of the models were manually tuned to increase dust opacity within the atmosphere beyond the range that fits other L dwarfs.

## 2.2 CO gas emission

We search for any CO line emission in both band 7 (J = 3 – 2) and band 6 (J = 2 – 1) data. This search is done by first subtracting the continuum emission from the visibilities and imaging the data with TCLEAN. This produces data cubes with channel widths of 0.43 km s⁻¹ (effective bandwidth of 1.1 km s⁻¹ due to Hanning smoothing) in band 7 and 1.3 km s⁻¹ (effective resolution of 2.0 km s⁻¹ due to Hanning smoothing and averaging every two channels) in band 6. The band 7 and 6 cubes have a noise level of 0.7 and 0.3 mJy per beam per channel, respectively, and do not have any significant emission over a single beam in a single channel. In order to search for low-level emission spread over multiple beams and channels, we applied a Keplerian mask as in Matrà et al. (2015, 2017) and Marino et al. (2016). This procedure corrects for the Doppler shift in each individual pixel due to Keplerian motion and centres the emission within a few channels, assuming a disc orientation (derived in Section 3) and sense of rotation – we test both possible directions. After integrating the emission within a deprojected radius of 150 au, we still do not detect any significant emission in the recovered spectra with rms levels of 20 mJy for J = 3 – 2 and 2 mJy for J = 2 – 1. Given this noise level and the fact that real CO emission in Keplerian rotation would appear as a single peak with a width equal to the effective bandwidth, we estimate 3σ upper limits for the line fluxes of 66 mJy km s⁻¹ for J = 3 – 2 and 12 mJy km s⁻¹ for J = 2 – 1.

These CO flux upper limits can be translated to CO masses assuming optically thin emission. We cannot simply assume local thermodynamic equilibrium (LTE) since at these low densities the collisional excitation of rotational levels can be very low and thus non-LTE effects must be taken into account (Matrà et al. 2015). Instead, we use the tool developed by Matrà et al. (2018a) to derive CO mass upper limits for a wide range of collisional partner densities (spanning from the radiation-dominated regime to LTE) and a range of kinetic temperatures from 20 to 200 K. We find that our upper

<footnote n="4">⁴There is a typo in table 1 in Bonsor et al. (2018), where K<sub>ej</sub> should be 5.14 × 10⁸ (Amy Bonsor, private communication).</footnote>
limit on CO $J = 2-1$ emission is the most constraining, with an upper limit of $2.4 \times 10^{-6}$ M$_{\odot}$.

Since CO gas is expected to be released in collisions if planetesimals are rich in CO (e.g. Zuckerman & Song 2012; Dent et al. 2014; Marino et al. 2016, 2020; Kral et al. 2017, 2019; Matrà et al. 2017; Moór et al. 2017), we can use this CO upper limit to place an upper limit in the CO + CO$_2$ ice mass fraction of planetesimals. In steady state, CO gas molecules are photodissociated by interstellar UV photons at a rate equal to the rate at which they are released from planetesimals. The latter is expected to be roughly equal to the product of the mass-loss rate of small dust and the mass fraction of CO + CO$_2$ in solids (as long as CO and CO$_2$ molecules escape before solids are ground down to µm-sized grains). Using equation (2) in Matrà et al. (2017) together with the fractional luminosity of the system (~3 × 10$^{-4}$), stellar mass (1.3 M$_{\odot}$), stellar luminosity (2.9 L$_{\odot}$), disc radius and width (approximately 70 and 100 au), and an expected photodissociation rate of 120 yr (Visser, van Dishoeck & Black 2009), we find an upper limit of 9 per cent for the fractional mass of CO + CO$_2$ in planetesimals. This limit is lower than that in some CO-rich comets, but it is still consistent with the wide distribution of abundances of Solar system comets (Mumma & Charnley 2011). Therefore, the CO non-detection is consistent with HD 206893's exo-Kuiper belt having a volatile composition similar to that of comets.

## 3 PARAMETRIC AXISYMMETRIC DISC MODEL

In order to quantify the location and width of the gap, we fit a parametric disc model with a gap as in Marino et al. (2018b, 2019), combining radiative transfer simulations (RADMC-3D; Dullemond et al. 2017) and a Markov chain Monte Carlo (MCMC) fitting procedure in the visibility space. The surface density of dust is defined by a two-power-law distribution, with an inner edge at $r_{min}$, surface density slopes (power-law index) $\gamma_1$ for $r < r_c$ and $\gamma_2$ for $r > r_c$, and a gap centred at $r_g$ (with $r_{min} < r_g < r_c$) and with a Gaussian profile characterized by a full width at half-maximum (FWHM) $w_g$. In contrast to Marino et al. (2018b, 2019), here we fixed the gap depth to 1 (i.e. the surface density is zero at $r_g$). The surface density is normalized such that the total dust mass is $M_d$. Note that this dust mass only includes the mass in grains smaller than 1 cm, which is set by the assumed dust opacities (see Marino et al. 2019). Since the disc is seen close to face-on and the total S/N is not very high, there is negligible information about the disc vertical structure in these data. We thus choose to fix the vertical extent of the disc to 5 per cent of the radius and use a uniform vertical mass distribution to simplify the model and speed up our simulations. The dust mass together with the assumed opacities set the disc overall brightness at 0.88 mm, while the brightness at 1.3 mm is set by the spectral index $\alpha$, which we leave as a free parameter that is uniform across the disc. The star is modelled with a template spectrum corresponding to an effective temperature of 6500 K and a radius of 1.3 R$_{\odot}$, which leads to a stellar flux of 30 and 13 µJy at 0.88 and 1.3 mm, respectively. In addition, we also fit the disc inclination $i$ and position angle PA, and nuisance parameters such as phase centre offsets, and the position, size (modelled as a 2D Gaussian), and flux of an SMG found 11.5 arcsec towards the SW of the star at both frequencies.

While in previous work we left the depth of the gap as a free parameter, here we opt to leave it fixed after a few tests where we found that if the depth is not fixed, the location and width of the gap are not well constrained. Note that this does not mean that we are

<table><caption>**Table 2.** Best-fitting parameters of the ALMA data using our parametric model. The quoted values correspond to the median, with uncertainties based on the 16th and 84th percentiles of the marginalized distributions.</caption><thead><tr><th>Parameter</th><th>Best-fitting value</th><th>Description</th></tr></thead><tbody><tr><td><i>M</i><sub>d</sub> (<i>M</i><sub>⊕</sub>)</td><td>0.031 ± 0.002</td><td>Total dust mass</td></tr><tr><td><i>r</i><sub>min</sub> (au)</td><td>28<sup>+5</sup><sub>-8</sub></td><td>Disc inner radius</td></tr><tr><td><i>r</i><sub>c</sub> (au)</td><td>114<sup>+7</sup><sub>-5</sub></td><td>Disc peak radius</td></tr><tr><td>γ<sub>1</sub></td><td>0.8<sup>+0.3</sup><sub>-0.4</sub></td><td>Inner disc slope index</td></tr><tr><td>γ<sub>2</sub></td><td>-3.9<sup>+0.5</sup><sub>-0.6</sub></td><td>Outer disc slope index</td></tr><tr><td><i>r</i><sub>g</sub> (au)</td><td>74 ± 3</td><td>Radius of the gap</td></tr><tr><td><i>w</i><sub>g</sub> (au)</td><td>27 ± 5</td><td>Width of the gap</td></tr><tr><td>δ<sub>g</sub></td><td>1.0</td><td>Fixed fractional gap depth</td></tr><tr><td>PA (°)</td><td>61 ± 4</td><td>Disc position angle</td></tr><tr><td><i>i</i> (°)</td><td>40 ± 3</td><td>Disc inclination from face-on</td></tr><tr><td>α</td><td>2.54 ± 0.17</td><td>Disc spectral index</td></tr></tbody></table>

forcing the presence of a gap, since the width is allowed to be of negligible size, which is analogous to removing the gap.

In addition to this, we impose a lower limit for $r_{min}$ of 14 au, which is roughly the minimum radius where solids could remain in the system on stable orbits if HD 206893 B was on a circular orbit at 10 au and had a mass of ~12$M_{Jup}$ (Wisdom 1980; Morrison & Malhotra 2015). Given this surface density definition, the disc does not have a well-defined outer edge. Nevertheless, we set a fixed outer edge of 250 au. This is justified by the absence of significant emission beyond 250 au and because we do not see any sharp outer edge in the disc brightness radial profiles. Finally, note that the surface density could be parametrized differently (e.g. with two belts instead of one with a gap); however, this choice of parametrization allows to constrain the gap between the two peaks of emission better using a single parameter.

### 3.1 Results

The best-fitting parameters are presented in Table 2. First, we find that the disc dust density peaks at $r_c = 114^{+7}_{-3} \text{ au}$. Exterior to that, the disc surface density declines steeply as $r^{-3.9 \pm 0.5}$. Interior to $r_c$, we find that in order to reconcile the dust levels away from the gap at ~30 and 120 au the surface density must increase moderately as $r^{0.8 \pm 0.4}$. The disc inner edge is most likely to be around 30 au, although smaller values down to 14 au are still compatible with the data within ~2σ. We can also constrain $r_{min} < 40 \text{ au}$ at a 99.7 per cent confidence level.

In between $r_{min}$ and $r_c$, we find that the gap centre is well constrained to $74 \pm 3 \text{ au}$. This location is consistent within $3\sigma$ with the location of the gaps discovered in three other systems at mm wavelengths (Marino et al. 2018b, 2019; MacGregor et al. 2019). The gap width is constrained to $27 \pm 5 \text{ au}$. This constraint on the width can be directly associated with a planet mass, assuming the gap is truly empty (as assumed in our model to derive its width) and was cleared by a planet on a circular orbit through scattering. Roughly speaking, the width of the gap is expected to be the same as the size of the chaotic zone where mean motion resonances overlap, i.e. $w_g \approx 3a_p\mu^{2/7}$ (where $a_p$ is the planet semimajor axis and $\mu$ is the planet-star mass ratio; Wisdom 1980; Morrison & Malhotra 2015; Marino et al. 2018b). Given the estimated stellar mass of 1.3 M$_{\odot}$ (Delorme et al. 2017), the measured gap width translates to a planet mass of $0.9^{+0.8}_{-0.5}M_{Jup}$, (based on the posterior distribution of $r_g$ and $w_g$), and
with $3\sigma$ upper and lower limits of $3.5M_{\text{Jup}}$ and $0.03M_{\text{Jup}}$ (10 M$_{\oplus}$), respectively. Therefore, the putative outer planet has a lower mass than the two inner companions. A caveat in this interpretation is that we have assumed the gap is empty at its centre. If this is not the case and the gap is not fully empty at its centre, the derived width would be biased towards smaller values in order to reproduce the same equivalent width. This means that the gap could be shallower and broader in reality, and more sensitive observations are needed to assess that. Nevertheless, if the gap was truly carved by a planet, we expect the surface density of particles to reach zero near the orbit of the planet (unless there is a large population of solids in horseshoe orbits; Marino et al. 2018b). Hence, fixing the gap depth to 1 is a reasonable assumption and consistent with our interpretation.

The disc orientation is well constrained with $i = 40 \pm 3^\circ$ and $\mathrm{PA} = 61 \pm 4^\circ$. Note that this inclination is consistent with the stellar pole inclination ($30 \pm 5^\circ$; Delorme et al. 2017). We will use the disc orientation later in Section 4 as a prior to fit the orbit of HD 206893 B assuming they are co-planar. The disc spectral index is found to be $2.53 \pm 0.16$, i.e. consistent with the fluxes measured from the images and with other debris discs (MacGregor et al. 2016).

The parameters to model the SMG that we found within the primary beam are also well constrained. We find it is centred at a separation of ($-5.6$ arcsec, $-10.0$ arcsec), i.e. towards the SW, and has an integrated flux of 1.5 mJy at 0.88 mm and 0.24 mJy at 1.3 mm. It is resolved, and its size (or standard deviation) is best fit by a $0.6 \ \mathrm{arcsec} \times 0.3 \ \mathrm{arcsec}$ 2D Gaussian profile.

Using our best-fitting model, we compute residual maps at the same resolution as the images in Fig. 1, and residual radial profiles with the same resolution as in Fig. 2. We do not find any residuals stronger than $3\sigma$ within 200 au. To visualize how well the model fits the data in the visibility space, in Fig. 5 we compare the observed deprojected and azimuthally averaged visibilities with the best-fitting model (green line). Note that this procedure assumes the disc is axisymmetric; thus, before averaging, we subtract the SMG component from the observed and model visibilities. Our best-fitting model reproduces well the multiple wiggles in the real component of the visibilities, while the imaginary component is consistent with zero. In the same figure, we overlay a model without a gap for comparison (green line). The model without a gap fails to reproduce the amplitude of the wiggles beyond $50 \ \mathrm{k\lambda}$.

## 4 DISC ORIENTATION AND HD 206893 B'S ORBIT

One of the advantages of having tight constraints on the orientation of the disc is that it can be compared with the orbit of companions (which could be misaligned), or can be used as prior information when deriving the orbits of companions assuming both lie in the same plane. In the case of HD 206893, the orientation of B's orbit is consistent with the disc orientation derived here (Grandjean et al. 2019), although the orbit inclination and position angle are not constrained as well as for the disc. Hence, there could be a moderate misalignment between the two. However, as shown by Pearce & Wyatt (2014), such misalignments do not last long if the companion is much more massive than the disc. We expect that in a few secular time-scales the disc will re-orient to the orbital plane of the companion (if the misalignment is $\lesssim 30^\circ$). The secular time-scale is only about 10 Myr at 100 au (or shorter if B is more massive than $12M_{\text{Jup}}$); thus, it would be unlikely to observe a misaligned disc given the age of this system (50–700 Myr; Delorme et al. 2017).

Assuming that the disc and B are co-planar, we can refine its orbital parameters using as priors the disc inclination and position angle

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Deprojected visibilities of band 7 (left) and band 6 (right) data, assuming a disc position angle of $61^\circ$ and an inclination of $40^\circ$. The real and imaginary components of the observed visibilities are azimuthally averaged within $8 \ \mathrm{k\lambda}$ wide radial bins, and are presented as black error bars in the top and bottom panels, respectively. The error bars represent the uncertainty estimated as the standard deviation of the observed visibilities in each bin divided by the square root of the number of independent points. The continuous green and red lines represent best-fitting models with and without a gap, respectively (binned using the same procedure as for the data). For better display, we crop data points beyond $240 \ \mathrm{k\lambda}$ since they are all consistent with zero.</figcaption></figure>

(longitude of ascending node) derived in Section 3.1. To constrain B's orbital parameters, we use its astrometry reported in Grandjean et al. (2019, table 2) and in Stolker et al. (2020, table B.1), and the package ORBITIZE (Blunt et al. 2020). This package allows to recover the posterior distribution of the orbital elements using a parallel-tempering MCMC algorithm (Vousden, Farr & Mandel 2016). Additionally, we assume that the SE side of the disc is the near side, and thus we set normally distributed priors with an orbital inclination $i \sim \mathcal{N}(140^\circ, 3^\circ)$ and longitude of ascending node $\Omega \sim \mathcal{N}(61^\circ, 4^\circ)$. We find that adding the disc orientation prior breaks some of the degeneracies and reduces the uncertainties in the estimates. We find B's semimajor axis, eccentricity, and argument of pericentre to be constrained to $a_p = 11.4^{+1.1}_{-0.8}$ au, $e_p = 0.14^{+0.05}_{-0.04}$, and $\omega = 74^{+36}_{-29}^\circ$, respectively (see Fig. 6 for the posterior distribution of the most relevant parameters using $10^7$ points after convergence). These results show that B's orbit is eccentric (larger than zero with an $\sim 3\sigma$ confidence) if co-planar with the disc, and the eccentricity is not larger than 0.36 (99.7 per cent confidence). The posterior distributions of the inclination and $\Omega$ are centred on the priors and have uncertainties that are slightly smaller than the prior distributions. Therefore, there is no indication that the co-planar assumption is in any tension with the astrometric data. In Fig. 7, we show 100 orbits drawn from the posterior distribution to illustrate the orientation of the orbit. We find that the pericentre is more likely to be found in the NW half of the system, and if very eccentric ($e_p > 0.3$) B is currently close to its pericentre. Using both semimajor axis and eccentricity distributions, we find a $3\sigma$ upper limit of 22 au on B's apocentre.

Given the estimated eccentricity, we can now assess whether B could have forced an eccentricity in the disc, producing both an offset in the disc radial structure and a higher disc flux near
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Posterior distribution of the semimajor axis, eccentricity, inclination, argument of pericentre, and longitude of ascending node, derived by fitting the astrometric position of HD 206893 B using the package ORBITIZE (Blunt et al. 2020), and using Gaussian priors on the inclination and longitude of ascending node. The vertical dashed lines represent the 16th, 50th, and 84th percentiles. Contours correspond to 68, 95, and 99.7 per cent confidence regions. This plot was generated using the PYTHON module CORNER (Foreman-Mackey 2016).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** HD 206893 B's possible orbits, drawn from 100 different points of the posterior distribution derived using the package ORBITIZE (Blunt et al. 2020). The orange points with error bars represent the astrometric data points.</figcaption></figure>

apocentre (i.e. apocentre glow; Pan, Nesvold & Kuchner 2016). Since in Section 2.1.1 we found that the NE half of the disc was brighter, it is fair to assume $\omega \sim 180^\circ$ and thus $e_p \lesssim 0.3$ (i.e. we assume the planet and disc are apsidally aligned, which is valid if the companion is much more massive than the disc). This means that the forced eccentricity, $e_f = 5a_{BD}e_{BD}/(4a)$ (Murray & Dermott 1999; Mustill & Wyatt 2009), imposed by the BD on particles with semimajor axis *a* near the disc peak density ($r_c$) is $\lesssim 0.03$. Such a low forced eccentricity would not produce a detectable contrast in brightness between apocentre and pericentre ($\sim 1.03$; Pan et al. 2016). For any other value of $\omega$, we still expect low forced eccentricities ($\lesssim 0.06$), mainly due to the small ratio between B's semimajor axis and the disc distance, and thus not detectable asymmetries due to interactions with B alone.

## 5 DISCUSSION

In this section, we summarize and discuss the different dynamical constraints on companions in HD 206893 (Section 5.1), and we discuss the different potential origins of the gap (Section 5.2), disc stirring (Section 5.3), and the ubiquity of gaps in exo-Kuiper belts (Section 5.4).
## 5.1 Three companions

In this section, we try to put in context the different observational constraints on HD 206893 to conclude on some basic properties of its companions. First, by combining the different astrometric data of HD 206893 B and the derived disc orientation in this work, we constrain the orbit of B to a semimajor axis of $11.4_{-0.8}^{+1.1}$ au and an eccentricity of $0.14_{-0.04}^{+0.05}$. Its mass remains uncertain in the range $12-50 M_{Jup}$, given its uncertain age (50–700 Myr; Delorme et al. 2017), which also means that B could have truncated the disc through scattering or via secular resonances (see Section 5.2.2). Since no additional companions have been detected with direct imaging beyond ~4 au, we can use this as an upper limit on the magnitude or mass of additional companions. We use the $5\sigma$ detection limits from Milli et al. (2017) and Delorme et al. (2017) and transform these to planet masses using AMES-COND models (Baraffe et al. 2003). Assuming an age of 50 Myr (700 Myr), we can rule out the presence of planets more massive than ~$3 M_{Jup}$ (~$10 M_{Jup}$) from 4 to 160 au.

Second, this system has a PMa of $96 \pm 28$ m s⁻¹ when subtracting the proper motion measured by Gaia DR2 from the one estimated using DR2 and Hipparcos astrometry, which has a longer baseline (Kervella et al. 2019). Although its magnitude is consistent with the dynamical effect that B would have on the star, as shown by Grandjean et al. (2019) the PMa has a position angle (233 ± 12°) that is inconsistent with the location of B (PA ~ 70° in 2015 May). This means that there is likely an additional inner companion (HD 206893 C), which would also be responsible for the RV drift of $87_{-14}^{+16}$ m s⁻¹ yr⁻¹ detected by Grandjean et al. (2019). Note that B is not sufficiently massive or close to the star to explain this RV acceleration. Given the inclination of this system (40°) and the length of this trend, C should have a semimajor axis larger than 1.4 au and a mass greater than ~$10 M_{Jup}$ to produce an acceleration of ~90 m s⁻¹ yr⁻¹. This inner companion would dominate the observed stellar velocity around the centre of mass ($v_* \propto m/\sqrt{a}$, with $m$ and $a$ the mass and semimajor axis of a companion, respectively) and thus the PA of the PMa is not expected to be correlated with the position of B. Using N-body simulations with the PYTHON package REBOUND (Rein & Liu 2012; Rein et al. 2019), we confirm this and find that both direction and magnitude of the PMa can be explained by C, even in the presence of B further out.

Third, the gap in the planetesimal belt centred at 74 ± 3 au and extending 27 ± 5 au suggests the presence of a third companion (HD 206893 b) that carved this gap around its orbit. The gap's width and centre constrain the planet mass to $0.9_{-0.5}^{+0.8} M_{Jup}$ and semimajor axis to 74 ± 3 au (based on the posterior distribution of the gap width in Section 3.1). Such a planet would have remained undetected in the existing direct imaging data. The SPHERE H2 observations reached a contrast of 16 mag at that distance, whereas AMES-COND models predict a contrast of 18 (27) mag if 50 (700) Myr old. Although detecting b is beyond the current capabilities of direct imaging instruments, it will certainly be within reach of James Webb Space Telescope with Near-Infrared Camera or Mid-Infrared Instrument at longer wavelengths. Note that these constraints assume that the gap was indeed carved by a planet in situ clearing its orbit through scattering. In Section 5.2, we discuss this possibility and alternative scenarios.

In Fig. 8 and Table 3, we summarize all these constraints on companions around HD 206893. The black points with error bars represent HD 206893 B and the putative planet b if the gap was cleared through scattering (labelled as 1; see Section 5.2.1) or if it was cleared through secular interactions (labelled as 2; see Section 5.2.5). The red hatched regions represent $5\sigma$ upper limits

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Companion mass versus semimajor axis in HD 206893. The black error bars represent HD 206893 B and two possible masses and semimajor axes of the putative planet b according to two different scenarios (see Table 3). The blue line and shaded region (68, 95, and 99.7 per cent confidence levels) represent the companion mass of an additional companion (C), necessary to explain the PMa, calculated using the Gaia DR2 PMa. The black dashed line represents the mass above which the orbit of a companion would be unstable based on a five mutual Hill radii criterion, and thus companions are excluded in the grey hatched region. The orange dotted lines represent the minimum semimajor axis and minimum mass of a companion to explain the RV drift. Taking this consideration and the stability limits, the companion responsible for the RV drift must lie in the orange shaded region. The red hatched regions represent the companion masses that would have been detected with SPHERE in addition to B assuming ages of 700 Myr (top) and 50 Myr (bottom). The green hatched region represents the companion masses and semimajor axes that are ruled out by the presence of the disc.</figcaption></figure>

<table><caption>**Table 3.** Summary of the inferred properties of the three companions in the HD 206893 system based on direct imaging, RV, PMa, and disc imaging with ALMA.</caption><thead><tr><th>Designation</th><th>a<br>(au)</th><th>e</th><th>Mass<br>(M<sub>Jup</sub>)</th><th>Reference</th></tr></thead><tbody><tr><td>HD 206893 (AC)B</td><td>11.2<sup>+1.5</sup><sub>-0.9</sub></td><td>0.14<sup>+0.07</sup><sub>-0.06</sub></td><td>12–50</td><td>(1, 2, Section 4)</td></tr><tr><td>HD 206893 (A)C</td><td>1.4–4.5</td><td>–</td><td>4–100</td><td>(1, 3, Section 5.1)</td></tr><tr><td>HD 206893 (ABC)b<sup>1</sup></td><td>74<sup>+3</sup><sub>-3</sub></td><td>–</td><td>0.9<sup>+0.8</sup><sub>-0.5</sub></td><td>(1, Section 5.2.1)</td></tr><tr><td>HD 206893 (ABC)b<sup>2</sup></td><td>30–50</td><td>–</td><td>0.07–0.4</td><td>(1, Section 5.2.5)</td></tr></tbody></table>

Note. References used in this table: (1) this work; (2) Delorme et al. (2017); (3) Grandjean et al. (2019). ¹Putative planet opening the gap around its orbit through scattering. ²Putative planet opening the gap beyond its orbit through secular interactions. Its semimajor axis is expected to be at the disc inner peak and its mass equal to the total disc mass.

from direct imaging assuming ages of 50 Myr (bottom) and 700 Myr (top). The green hatched region shows the planets that are ruled out if on a circular orbit since they would push the disc inner edge beyond 40 au (inconsistent with our observations), produce a gap at 74 au much wider than observed, or carve additional gaps in the disc that are not observed. Note that in the secular interaction scenario the planet can be located in a region where the disc density is high (Pearce & Wyatt 2015). The blue line and shaded region (68, 95, and 99.7 per cent confidence levels) represent the mass of C in order to explain the magnitude of HD 206893's PMa. This mass is calculated using the Gaia DR2 PMa (after subtracting the PM from Hipparcos–Gaia DR2), equations in Kervella et al. (2019), and the disc orientation derived in this work. The confidence levels were determined using Monte Carlo simulations (bootstrapping the PMa and disc orientation). In addition, the dashed black line represents
the mass above which the system would become unstable in a short time-scale, here simply defined as an orbital spacing smaller than five mutual Hill radii assuming B is on a circular orbit and has a mass of $12M_{\text{Jup}}$. Thus, we can exclude companions in the grey hatched region. Note that the stability criterion depends on several factors that we are not exploring here and thus this upper limit should be taken with caution (see details in Smith & Lissauer 2009). Moreover, we found that B is likely on an eccentric orbit and thus the allowed planet masses would be even lower (Gladman 1993; Lazzoni et al. 2018). In the same figure, we overlay the RV constraints on C, namely its minimum semimajor axis (vertical orange line) and minimum mass (diagonal orange line) calculated to produce an RV acceleration of at least $45 \text{ m s}^{-1} \text{ yr}^{-1}$ ($3\sigma$ lower limit; Grandjean et al. 2019). The orange shaded region therefore represents a conservative range of companions that could explain the RV trend and still be stable. When combined, these constraints imply that the innermost companion should lie in the intersection between the orange wedge and the blue shaded area to explain the RV trend and PMa, i.e. have a semimajor axis of roughly $1.4-4.5 \text{ au}$ and a mass in the range $4-100M_{\text{Jup}}$. This massive inner companion would lie at a projected separation of $26-110 \text{ mas}$, and thus it could be resolved using sparse aperture masking on SPHERE or with GRAVITY at K band using molecular mapping. It is difficult to predict its exact PA since its period is not well constrained ($1.5-8.4 \text{ yr}$) and is comparable to the period over which Gaia obtained astrometric observations ($1.8 \text{ yr}$; Gaia Collaboration 2018). Hence, the intrinsic velocity vector of the star has been significantly averaged over time (Kervella et al. 2019), and thus its direction is biased. The analysis presented here illustrates how the combination of NIR direct imaging, ALMA imaging, RV, and PMa constraints can be used in combination to tightly constrain the architecture of a system.

Note that while we do not have strong constraints on the eccentricities of HD 206893 b and C, if B is truly eccentric we expect that the two additional companions will gain eccentricity through secular interactions even if they had initially zero eccentricities. Over Myr time-scales, they will exchange eccentricities and thus b and C could have eccentricities as high as B depending on their masses and semimajor axes.

## 5.2 The gap origin

In this section, we discuss potential dynamical origins for the observed gap based on different scenarios explored in the previous work and the multiple information that exists for this system, namely the BD detected at a 11 au separation and the inner companion responsible for the RV drift and PMa (see Section 5.1).

### 5.2.1 Gap clearing by planet at 74 au

As commented in Section 3.1, the gap discovered at 74 au hints at the presence of a planet on a circular orbit clearing the material near its orbit. The putative planet would need to be a gas or ice giant of mass $0.9_{-0.5}^{+0.8} M_{\text{Jup}}$ (based on the posterior distribution of the gap width), in order to create a $27 \pm 5 \text{ au}$ wide gap via scattering of particles within its chaotic zone. Such a planet, if massive enough to retain a significant atmosphere, probably formed while a gas-rich protoplanetary disc was still present and possibly before the planetesimal belt was formed. In that hypothetical case, the distribution of solids would have been already truncated around the planet's orbit. As shown in Marino et al. (2019), if the gap was inherited from a gap in the distribution of dust that then grew to

form the planetesimal belt, then the putative planet could have a lower mass of a few tens of Earth masses. In fact, several gaps are observed at those distances in protoplanetary discs (e.g. Andrews et al. 2018; Huang et al. 2018; Long et al. 2018); hence, it is possible that the gap could have been inherited from the dust distribution in HD 206893's protoplanetary disc. While growing, this planet could have also migrated through the gas-rich disc, or via planetesimal-driven migration after disc dispersal possibly widening the gap. This means that the required planet mass to open such a gap could be much lower. Therefore, the derived planet masses must be taken with caution. Nevertheless, the gap width does impose a strict upper limit of $3.5 M_{\text{Jup}}$, otherwise the gap width would be much greater than the ALMA observations show.

If the putative planet was on a slightly eccentric orbit, it would likely clear an asymmetric and eccentric gap and force an eccentricity in the disc (e.g. Pearce & Wyatt 2014). Note that this putative planet could become eccentric simply through secular interactions with B that is likely on an eccentric orbit (Section 4). This eccentricity could then explain the tentative $30 \pm 13 \text{ per cent}$ brightness asymmetry. Detailed N-body simulations that explore the level of disc asymmetry produced by a planet in the middle of the disc with different eccentricities are needed to fully assess this scenario.

While a planet with an orbit within the gap is plausible, this scenario triggers many questions since the gaps found around HD 107146, HD 92945, and HD 15115 (Marino et al. 2018b, 2019; MacGregor et al. 2019) lie all at a very similar radius ($59-74 \text{ au}$). Naively we would expect these radii to differ by more given the multiple factors that determine a planet's final orbit and the wide range of radii covered by these discs ($40-150 \text{ au}$). This is not an issue that applies only to this scenario, but rather to any dynamical mechanism that requires the presence of planets and fine-tuning of their parameters. Alternatively, the gaps could lie at a similar radius if they are formed instead as a consequence of a change in planetesimal properties (e.g. inefficient planetesimal formation, or weaker planetesimals that collisionally evolve faster) or a hotspot for efficient planet formation. These questions are difficult to answer with a limited sample of discs observed with enough sensitivity and resolution to detect these gaps (see Section 5.4). A future ALMA survey could expand this sample and constrain the ubiquity of gaps and determine their radius and width distribution.

### 5.2.2 Secular resonances between planet at 2 au and BD at 11 au

Given the two companions on orbits interior to the disc inner edge, this system is well suited to test whether secular resonances exterior to their orbits could create a gap in the disc at 74 au where particles eccentricities are excited to large values (Yelverton & Kennedy 2018). To assess this, we solve equation (13) in Yelverton & Kennedy (2018)<fnref n="5" /> to find the gap/resonance location using the estimated masses and semimajor axes of the two inner companions. We find that secular resonances cannot open a gap at such a large distance with respect to HD 206893 B (see solid blue line in Fig. 9). Since the mass ratio between B and C is $\lesssim 10$ (see Table 3), we find that the strongest secular resonance would be located at $\lesssim 50 \text{ au}$. Therefore, this scenario cannot explain the gap at 74 au.

Nevertheless, secular resonances could have instead truncated the disc inner edge. Under the influence of a single planet on an eccentric

<footnote n="5">⁵There is a typo in equation (13) in Yelverton & Kennedy (2018). In the right-hand side, the first factor should be $a_1^{1/2} a_2^2 a^{-5/2}$ instead of $(a_1/a_2)^{1/2}$. We confirmed this through private communication with Yelverton et al.</footnote>
<figure><img src="image_11.png" /><figcaption>Figure 9. Location of the strongest secular resonance exterior to the orbits of different pairs of companions: massive inner companions at 2 and 11 au (blue); BD at 11 au and hypothetical planet at 30 au (orange); and BD at 11 au and putative planet in the gap at 74 au (green). The continuous line represent the true location of the resonance, while the dashed line its location approximated by equation (2). Only mass ratios consistent with the dynamical constraints are displayed. The discontinuity on the line is due to one of the two exterior resonances becoming stronger than the other above a mass ratio $\sqrt{a_1/a_2}$, while the other becomes negligible (we only plot the strongest). $m_1$ and $m_2$ stand for the masses of the inner and outer companions, respectively.</figcaption></figure>

orbit, Pearce & Wyatt (2014) showed that the disc should be truncated at a semimajor axis (at apocentre)

$$a_{\text{in}} \approx a_{\text{p}}(1 + e_{\text{p}}) + 5R_{\text{H,Q}}, \quad (1)$$

where $R_{\text{H,Q}}$ is the Hill radius of an eccentric planet at apocentre. Under the influence of two planets, that inner edge could be further out if the strongest secular resonance is located beyond $a_{\text{in}}$. Using equation (24) in Yelverton & Kennedy (2018), we find the approximate location of the resonance to be (dashed blue line in Fig. 9)

$$r_{\text{min}}^{\text{SR}} = 29 \left( \frac{a_2}{11 \text{ au}} \right)^{11/7} \left( \frac{a_1}{2 \text{ au}} \right)^{-4/7} \left( \frac{m_2}{m_1} \right)^{2/7} \text{ au}, \quad (2)$$

where $a_1$ and $a_2$ are the semimajor axes of the companions, while $m_1$ and $m_2$ are their masses, with 1 designating the innermost companion and 2 the BD. In order to effectively truncate the disc, the width of the secular resonance would need to be large enough to excite solids from $a_{\text{in}}$ to $r_{\text{min}}^{\text{SR}}$. This width is controlled primarily by the eccentricity of HD 206893 B. Given the uncertainties on its eccentricity ($0.14^{+0.05}_{-0.04}$) and its mass (12–50$M_{\text{Jup}}$; Delorme et al. 2017), it is hard to assess whether secular resonances will deplete a region that otherwise would be stable. Assuming the best-fitting values of our fit and using equation (36) in Yelverton & Kennedy (2018), we expect a secular resonance width of ~20 au. Therefore, it is possible that the disc inner edge is not truncated by the BD, but by secular resonances. Note that HD 206893 B alone could truncate the disc out to ~30 au if it has an eccentricity of 0.2 and a mass of 50$M_{\text{Jup}}$. A more precise characterization of HD 206893 B's orbit and the disc inner edge could provide better estimates of the masses of the two inner companions.

### 5.2.3 Secular resonances between BD and planet at 30 au

Secular resonances might still have created the gap if instead the resonance at 74 au was with the BD at 11 au and an outer planet

sitting just interior to the disc inner edge at around 30 au. We find that this putative planet and the BD at 11 au would need to have a mass ratio close to unity (see solid orange line in Fig. 9). This is ruled out by direct imaging, which did not detect any companion more massive than $3M_{\text{Jup}}$ at 30 au (assuming an age of 50 Myr; Delorme et al. 2017).

### 5.2.4 Secular resonances between BD and planet at 74 au

If the gap at 74 au was indeed carved by a planet at that distance, it is worth discussing what other observables could confirm this scenario. Given the range of possible masses of the BD and the putative planet c at 74 au, a secular resonance could be located in between the orbit of b and the disc outer edge. We find that the resonance would be located between 85 and 170 au for mass ratios of b and B of $10^{-3}$ to 0.3 (solid green line in Fig. 9). This resonance would create a gap that is expected to be very wide ($\gtrsim$50 au if b has an eccentricity of $\gtrsim$0.05 according to equation 36 in Yelverton & Kennedy 2018), and thus noticeable by our observations. This means that we can already rule out that this gap is present in between 74 and 110 au, otherwise the gap would be seen to be much wider and the disc peak further out (or at 40 au instead). Therefore, imposing that the secular resonance is located beyond 120 au, we can constrain the mass ratio to be larger than 0.04 (consistent with the expected range; see Table 3). This means the putative planet would need to be at least $0.5M_{\text{Jup}}$ in mass, and thus likely a gas giant.

There could already be evidence of the presence of this secular resonance in between 120 and 170 au. Such a resonance could be responsible for shaping the distribution of solids in the outer regions, giving HD 206893's disc its observed appearance. As Figs 2 and 11 show, its outer edge is not sharp or well defined. The surface brightness declines smoothly with radius down to the noise level at 180 au. The smooth decline could be due to solids in high-eccentricity orbits, similar to a scattered disc (Booth et al. 2009; Wyatt et al. 2010; Marino et al. 2017a; Geiler et al. 2019). In contrast, HD 107146 and HD 92945 have well-defined sharp outer edges. Therefore, this smoother outer edge in HD 206893 could be due to the effect of secular resonances exciting the eccentricities of particles over a broad region in the outer edge of this disc.

A caveat in the reasoning above is that the analysis in Yelverton & Kennedy (2018) is only strictly valid for a system with only two planets and a massless disc (or with a mass much smaller than the planet masses). Taking into account the disc mass and the effect of the innermost companion is beyond the scope of this paper, but we acknowledge that it is important for the location of secular resonances and thus is crucial for constraining the mass of this putative planet b.

### 5.2.5 Secular interactions between the disc and a scattered planet

A different possibility to explain the observed gap and tentative disc asymmetry is through secular interactions between the disc and a planet of a similar mass (Pearce & Wyatt 2015). In this scenario, a planet is scattered out by a more massive one (e.g. B) on to an eccentric orbit that overlaps with the disc. Through secular interactions, the low-mass planet is able to open a broad and asymmetric gap in the disc, which could explain our observations. For this to happen, the disc must be wide and its mass must be comparable to the planet mass. Although it is hard to assess the total disc mass since observations are not sensitive to planetesimals, we use the collisional model described in Marino et al. (2017b) and estimate the total disc mass to be around ~20–120 $M_{\oplus}$ assuming the disc
has been collisionally evolving for 50–700 Myr and planetesimals have a maximum size of 100 km. Therefore, we can already rule out that the BD opened the observed gap since it is at least 40 times more massive than the disc. Moreover, the model by Pearce & Wyatt (2015) predicts that HD 206893 B would be located just interior to the gap where the disc surface density peaks (30–50 au according to Fig. 2), which is inconsistent with B's position at 11 au.

Nevertheless, B could have been the massive companion that scattered out an additional less massive planet that could have opened the gap and would now reside in between the disc inner edge and the gap. Note that this scenario proposed by Pearce & Wyatt (2015) did not take into account the influence of an inner massive companion and thus it is unclear how exactly its presence affects the secular evolution and gap opening. Simulations tailored to this system could help to assess in detail whether this scenario could be at play. The last row of Table 3 shows the mass and orbital constraints on this putative planet for this scenario.

### 5.2.6 Planet-less scenarios

Finally, we consider whether the gap could have been opened without the influence of a planet. As mentioned in Marino et al. (2019), gaps in the distribution of mm-sized grains could be due to changes in planetesimal properties within the gap (e.g. strength, maximum size, porosity). Such changes could change the size distribution and collisional evolution in a way that it produces a depletion on mm-sized grains. It is uncertain how large those changes would need to be and whether this is something plausible according to planetesimal formation models.

Another possibility is that dust-gas interactions could lead to gaps or ring-like structure (e.g. through photoelectric instabilities; Klahr & Lin 2005; Lyra & Kuchner 2013; Richert, Lyra & Kuchner 2018). As in Marino et al. (2018b, 2019), we can also rule out that for this system gas drag could be important for mm-sized grains. Using our CO mass upper limit (Section 2.2) and dust mass (Section 3.1), we find that the CO gas-to-dust mass ratio is at least $10^{-4}$ and thus even if H$_2$ was present (with a typical interstellar medium-like H$_2$/CO abundance of $10^4$) we expect a dust-to-gas mass ratio much smaller than unity. Moreover, collisional time-scales are orders of magnitude shorter than stopping times due to gas drag, and thus its effect is expected to be negligible for mm-sized grains.

A third possibility is that the gaps were inherited from the solid distribution in protoplanetary discs, where gaps in the dust distribution can be created via different methods without the intervention of any planet, for example dead zones (Pinilla et al. 2016), magnetohydrodynamic zonal flows (Flock et al. 2015), secular gravitational instability (Takahashi & Inutsuka 2014), instabilities originating from dust settling (Lorén-Aguilar & Bate 2015), dust particle growth by condensation near ice lines (Saito & Sirono 2011), and viscous ring instability driven by dust (Dullemond & Penzlin 2018). All these mechanisms can shape the distribution of solids in protoplanetary discs, which could be later inherited by planetesimals formed from those solids, although how exactly this would be inherited is uncertain given the unconstrained planetesimal formation process. If any of those models were to predict a fixed gap radius around 70 au, then this could explain the distribution of gaps observed so far in exo-Kuiper belts.

## 5.3 Disc stirring

Since the observed dust must be replenished through collisions of larger planetesimals (e.g. Wyatt & Dent 2002), their orbits must be

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Comparison between the age of HD 206893 (grey region) and stirring time-scales due to HD 206893 B (blue), b (green), and self-stirring (orange) as a function of disc semimajor axis or radius. The width of the shaded regions represents the respective uncertainties due to the age uncertainty, which could range between 50 Myr (dashed line) and 700 Myr (dotted line).</figcaption></figure>

stirred above a certain level such that relative velocities are high enough to cause destructive collisions. Planetesimals are expected to have nearly circular orbits when protoplanetary discs disperse; therefore, stirring should take place after disc dispersal and on a time-scale shorter than the age of the observed system. There are two main stirring mechanisms proposed. Planetesimals could be stirred by planets or more massive companions (e.g. Mustill & Wyatt 2009). Alternatively, big planetesimals within the disc could stir it and ignite a collisional cascade (i.e. self-stirring; Kenyon & Bromley 2008, 2010; Kennedy & Wyatt 2010; Krivov & Booth 2018).

For HD 206893, we know the system has a massive companion at 11 au that may have stirred the disc. Based on its estimated eccentricity of $0.14_{-0.04}^{+0.05}$ (Section 4), we find it could stir solids to eccentricities above 0.01 out to 150 au (equation 8 in Mustill & Wyatt 2009), unless its eccentricity is below 0.1. In Fig. 10, we compare the age of the system with the time-scale it would take to stir solids due to secular interactions with B (equation 15 in Mustill & Wyatt 2009) for a mass in the range $12-50 M_{\text{Jup}}$ (blue shaded region). Curves corresponding to a system age of 50 (700) Myr are represented with dashed (dotted) lines. We find that B could have stirred the disc out to ~100 au if the system is young, or beyond 200 au if it is rather old. Hence, unless B has a low eccentricity or the system is young, B alone could explain that the disc is stirred at least out to 150 au. Moreover, if the putative planet b is indeed at 74 au, has a mass between $0.4 M_{\text{Jup}}$ and $1.7 M_{\text{Jup}}$, and is on a mildly eccentric orbit (0.02), we find that it could have stir the disc even on a shorter time-scale than B beyond 74 au (green shaded region; the discontinuity around 74 au is due to the use of equation 16 in Mustill & Wyatt 2009, for $r < 74$ au). Note that the innermost companion would be very inefficient at stirring solids at large radii given its small semimajor axis.

It is also worth considering self-stirring to see whether both mechanisms could be at play. Using equation (33) in Krivov & Booth (2018), a maximum planetesimal diameter of 100 km, and a disc mass of 20–120 $M_{\oplus}$, we find self-stirring is too slow to excite eccentricities beyond 100 au. Considering a larger maximum planetesimal diameter of 400 km, and correspondingly a disc mass larger by a factor of ~1.3 to fit the same surface density of dust (see equation 7 in Marino et al. 2019), we find that the disc could have been self-stirred out to 150 au if older than 500 Myr (orange shaded region). Nevertheless,
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Deprojected surface brightness profiles of HD 107146 (blue line), HD 92945 (orange dashed line), and HD 206893 at 1.3 mm (green dotted line), computed by azimuthally averaging the emission. We subtracted the emission at the centre using Gaussian profiles of FWHM equal to the beam major axes. The shaded regions correspond to $1\sigma$ uncertainties. Note that the shaded regions are representative of the uncertainty over a resolution element, whose sizes are represented by coloured horizontal lines at the bottom of the figure. The dashed vertical line represents the semimajor axis of HD 206893 B and the grey region its chaotic zone.</figcaption></figure>

the time-scale for self-stirring is overall longer than stirring by B or b; therefore, we conclude that the disc was likely stirred by the companions if born unstirred (a similar conclusion was reached by Musso Barcucci et al. 2019, for the debris disc around HD 193571).

## 5.4 Are gaps common among exo-Kuiper belts?

Based on the literature, there are only six exo-Kuiper belts that so far have been observed with ALMA with enough resolution and sensitivity to detect gaps (e.g. $\gtrsim$4 beams across their width and an S/N $\gtrsim$ 10 in the deprojected radial profile). These are HD 107146, HD 92945, $\beta$ Pic, HD 15115, AU Mic (Marino et al. 2018b, 2019; Daley et al. 2019; MacGregor et al. 2019; Matrà et al. 2019), and now HD 206893. Four of these show evidence of gaps (HD 107146, HD 92945, HD 15115, and HD 206893), suggesting gaps could be common in exo-Kuiper belts, at least among wide and bright discs. This number could grow to five if we consider AU Mic's best-fitting model of a narrow inner ring in addition to a broader outer disc, as evidence of a gap in between these two components. Such ubiquity is not rare among protoplanetary discs, with the majority of discs larger than 50 au showing radial substructure in the form of gaps and rings (e.g. Andrews et al. 2018; Long et al. 2018, 2019). This means gaps could be truly ubiquitous in both large protoplanetary and planetesimal belts. Whether this type of substructure can be directly linked between the two, e.g. gaps in planetesimal belts could be inherited from the dust distribution in protoplanetary discs, is still unclear. Future ALMA observations of a large sample of debris disc could reveal the distribution of gap radii and widths, and thus allow for a statistical comparison between the two. Moreover, gaps could also be present in narrow exo-Kuiper belts as well, as it is in some protoplanetary discs (e.g. HD 169142; Pérez et al. 2019a), and this is something yet to be explored.

Although the gap locations of all four debris discs with gaps observed with ALMA lie at a similar radius, they have significant differences in their radial profiles. In Fig. 11, we compare the

surface brightness of these discs (except for HD15115 that is edge-on). HD 206893's disc seems to be the widest disc with the smallest inner edge and largest outer edge. Given the lower S/N in HD 206893 observations, it is hard to compare the gap width and depth with the other two. Nevertheless, a rough comparison indicates that HD 206893's gap is wider than HD 92945's and deeper than HD 107146's (the only gap that is well resolved with several beams across). Another strong difference is that in HD 206893, the peak in surface density is beyond the gap, while this is the opposite in HD 107146 and HD 92945. This difference could be due to collisional evolution if HD 206893 is significantly older (e.g. 700 Myr old) compared to HD 92945 and HD 107146 (e.g. 100 Myr old), in which case we expect a larger relative depletion between the region interior and exterior to the gap (e.g. Wyatt et al. 2007; Kennedy & Wyatt 2010).

In scattered light, gaps have also been detected in the distribution of small grains in six exo-Kuiper belts: HD 92945 and HD 15115 with gaps also seen at mm wavelengths (Golimowski et al. 2011; Schneider et al. 2014; Engler et al. 2019), and in HD 141569, HD 131835, HD 120326, and HD 141943, only seen so far in scattered light observations (Perrot et al. 2016; Bonnefoy et al. 2017; Feldt et al. 2017; Boccaletti et al. 2019). These gaps would be consistent with gaps in the distribution of planetesimals. However, these could also be due to gas-dust interactions, which could generate multiple rings and gaps in the distribution of µm-sized dust (Lyra & Kuchner 2013; Richert et al. 2018). In fact, two of these systems are known to have large levels of gas (HD 141569 and HD 131835; Zuckerman, Forveille & Kastner 1995; Moór et al. 2015; Kral et al. 2019). ALMA observations constraining the presence of gas and the existence of these gaps in the distribution of large dust are crucial to assess whether these gaps are also present in the distribution of planetesimals.

An important consequence of gaps being common among exo-Kuiper belts is that they indicate that a large fraction of their mass was scattered away (unless they are inherited from protoplanetary discs). As discussed in Marino et al. (2019), this material could encounter additional inner planets (e.g. at the belt inner edges) that could scatter the material even closer in until being accreted by temperate low-mass planets (e.g. Marino et al. 2018a). Such accretion could lead to volatile delivery (Kral et al. 2018; Wyatt, Kral & Sinclair 2020) and the build-up of secondary atmospheres, or instead atmospheric erosion depending on planetesimal properties. Understanding the frequency of this process is therefore important to constrain the evolution of atmospheres of close-in planets in systems with exo-Kuiper belts.

## 6 CONCLUSIONS

In this work, we have presented the first ALMA observations (at 0.88 and 1.3 mm) of the system HD 206893 to image its debris disc. This system is known to host a directly imaged companion with a mass of 12–50$M_{\text{Jup}}$ and a semimajor axis of 11 au (HD 206893 B), and likely an additional inner companion at around 2 au responsible for an RV trend and PMa (HD 206893 C). Through analysis in the image and visibility space, we have found that the disc extends from roughly 30 to 180 au, with a peak in surface brightness and density at 110 au, and a local minimum or gap at 74 au. This gap is found to be $27 \pm 5$ au wide, which if carved by a planet *in situ* through scattering can be translated to a planet mass of $0.9_{-0.5}^{+0.8}M_{\text{Jup}}$ (HD 206893 b). This gap in a debris/planetesimal disc is the fourth to be found, and is centred at a similar radius as the rest, namely around 70 au. Why these gaps seem to be located at the same radius is unclear.
In addition to studying the radial structure, we searched for asymmetries. We find a marginal evidence of an asymmetry in the disc with the NE half being ~30 per cent brighter, which cannot be explained by a single background galaxy. If real, this asymmetry could be due to dynamical interactions with b. We also searched for CO emission in the system, but we found no significant emission – still consistent with solids being volatile rich and having compositions similar to Solar system comets.

Since it has been proposed that B's spectrum is reddened by circumstellar material, we searched for dust emission at B's position. We did not find any emission, which rules out the presence of a massive dusty disc larger than 0.1 au. Moreover, we find that accreted dust from the outer debris discs would also be insufficient to cause any reddening. Therefore, if reddening is caused by dust, this is probably lifted inside and above its photosphere.

Based on the derived disc orientation, we were able to better constrain the orbit of B by assuming it is co-planar with the disc. We find B is likely on an eccentric orbit with an eccentricity of $0.14^{+0.05}_{-0.04}$ and semimajor axis of $11.4^{+1.1}_{-0.8}$ au. Given these constraints and B's estimated mass, the disc could have been truncated by B, explaining its observed inner edge of $27 \pm 5$ au. However, the exact position of the disc inner edge and the predicted truncation radius (given by B's mass and orbit) are still very uncertain; hence, it is still possible that the inner edge is farther away than expected.

We have used all available dynamical and observational constraints (RV, PMa, stability) to determine the mass and semimajor axis of C. We have found that its semimajor axis should be in the range of 1.4–4.5 au (consistent with previous work) and have a mass between $4M_{\text{Jup}}$ and $100M_{\text{Jup}}$. Therefore, C could be a massive gas giant, BD, or a low-mass star. Based on their estimated orbits and masses, secular interactions between C and B could place a secular resonance at about 30 au, near the disc inner edge. Therefore, it is plausible that the disc inner edge was truncated by this resonance.

While the gap at 74 au could have been carved by a planet in *situ*, there are other dynamical mechanisms by which a planet could carve such a gap, namely secular resonances with two inner planets and secular interactions between the disc and a planet on a highly eccentric orbit. For the former scenario to work, the putative outer planet would need to be located at 30 au and have a similar mass compared to B. However, such a planet is ruled out by direct imaging observations that did not detect additional companions beyond B with a similar mass. The latter scenario could work if the scattered planet has a mass similar to the disc (~20–120 $M_{\oplus}$), which rules out that B is responsible for the gap. If B instead scattered out a low-mass planet, such a planet could have opened the gap and today reside at ~40 au on a low-eccentricity orbit.

Since the observed dust indicates an ongoing collisional cascade and thus ongoing destructive planetesimal collisions, the disc must have been stirred in the past. This could have happened either via secular interactions with B or b, or via self-stirring. We find that both mechanisms could be efficient at stirring the disc in time-scales shorter than the age of the system, but it is likely that planet stirring by b or B dominates.

HD 206893 is a unique laboratory to study planetary dynamics and the interaction between planetesimal discs and massive companions. Future deeper ALMA observations could better constrain the dust distribution within the gap, the level of asymmetry in the disc, and the exact position of the disc inner edge. Such constraints could favour a specific scenario of the ones discussed in this paper to explain the gap, and help to estimate better the masses of the inner companions. Finally, there is growing evidence indicating that gaps could be common in exo-Kuiper belts, although the sample size

of debris discs that have been observed with ALMA with enough resolution and sensitivity is still very small. If gaps carved by planets in exo-Kuiper belts are common, it is possible that inward scattering of volatile-rich material from the belt to inner planets and subsequent accretion of volatiles is a frequent process in exoplanetary systems.

## ACKNOWLEDGEMENTS

We thank Jeff Jennings and Richard Booth for their help and discussion on how to use FRANK. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2017.1.00828.S and 2017.1.00825.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. TH acknowledges support from the European Research Council under the Horizon 2020 Framework Program via the ERC Advanced Grant Origins 83 24 28. GMK was supported by the Royal Society as a Royal Society University Research Fellow. VF's postdoctoral fellowship was supported by the Exoplanet Science Initiative at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). Finally, we would also like to thank the anonymous referee for a very constructive report that improved the clarity of this paper.

## DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author. The ALMA band 6 data are publicly available and can be queried and downloaded directly from the ALMA archive at https://almascience.nrao.edu/asax/. The band 7 data will become publicly available on 2020 August 13 at the same web address.

## REFERENCES

* Acke B. et al., 2012, A&A, 540, A125

* Andrews S. M. et al., 2018, ApJ, 869, L41

* Baraffe I., Chabrier G., Barman T. S., Allard F., Hauschildt P. H., 2003, A&A, 402, 701

* Beust H. et al., 2014, A&A, 561, A43

* Blunt S. et al., 2020, AJ, 159, 89

* Boccaletti A. et al., 2019, A&A, 625, A21

* Boley A. C., Payne M. J., Corder S., Dent W. R. F., Ford E. B., Shabram M., 2012, ApJ, 750, L21

* Bonnefoy M. et al., 2017, A&A, 597, L7

* Bonsor A., Wyatt M. C., Kral Q., Kennedy G., Shannon A., Ertel S., 2018, MNRAS, 480, 5560

* Booth M., Wyatt M. C., Morbidelli A., Moro-Martín A., Levison H. F., 2009, MNRAS, 399, 385

* Booth M. et al., 2016, MNRAS, 460, L10

* Boss A. P., 1997, Science, 276, 1836

* Boss A. P., 2003, ApJ, 599, 577

* Boss A. P., 2011, ApJ, 731, 74

* Burns J. A., Lamy P. L., Soter S., 1979, Icarus, 40, 1

* Chen C. H., Mittal T., Kuchner M., Forrest W. J., Lisse C. M., Manoj P., Sargent B. A., Watson D. M., 2014, ApJS, 211, 25

* Chiang E., Kite E., Kalas P., Graham J. R., Clampin M., 2009, ApJ, 693, 734

* Daley C. et al., 2019, ApJ, 875, 87

* Delorme P. et al., 2017, A&A, 608, A79

* Dent W. R. F. et al., 2014, Science, 343, 1490
Dullemond C., Juhasz A., Pohl A., Sereshti F., Shetty R., Peters T., Commercon B., Flock M., 2017, RADMC3D v0.41, http://www.itau.uni-heidelberg.de/dullemond/software/radmc-3d/

Dullemond C. P., Penzlin A. B. T., 2018, A&A, 609, A50

Engler N. et al., 2019, A&A, 622, A192

Faramaz V., Beust H., Augereau J.-C., Kalas P., Graham J. R., 2015, A&A, 573, A87

Faramaz V. et al., 2019, AJ, 158, 162

Feldt M. et al., 2017, A&A, 601, A7

Flock M., Ruge J. P., Dzyurkevich N., Henning T., Klahr H., Wolf S., 2015, A&A, 574, A68

Foreman-Mackey D., 2016, J. Open Source Softw., 1, 24

Fujimoto S., Ouchi M., Shibuya T., Nagai H., 2017, ApJ, 850, 83

Gaia Collaboration, 2018, A&A, 616, A1

Gaspar A., Rieke G. H., 2020, Proc. Natl Acad. Sci. USA, 117, 9712

Geiler F., Krivov A. V., Booth M., Löhne T., 2019, MNRAS, 483, 332

Gladman B., 1993, Icarus, 106, 247

Golimowski D. A. et al., 2011, AJ, 142, 30

Grandjean A. et al., 2019, A&A, 627, L9

Greaves J. S., Holland W. S., Jayawardhana R., Wyatt M. C., Dent W. R. F., 2004, MNRAS, 348, 1097

Huang J. et al., 2018, ApJ, 869, L42

Hughes A. M., Duchêne G., Matthews B. C., 2018, ARA&A, 56, 541

Jennings J., Booth R. A., Tazzari M., Rosotti G. P., Clarke C. J., 2020, MNRAS, 495, 3209

Kalas P., Graham J. R., Clampin M., 2005, Nature, 435, 1067

Kalas P. et al., 2008, Science, 322, 1345

Kennedy G. M., Wyatt M. C., 2010, MNRAS, 405, 1253

Kennedy G. M., Wyatt M. C., 2011, MNRAS, 412, 2137

Kennedy G. M. et al., 2015a, ApJS, 216, 23

Kennedy G. M. et al., 2015b, MNRAS, 449, 3121

Kenyon S. J., Bromley B. C., 2008, ApJS, 179, 451

Kenyon S. J., Bromley B. C., 2010, ApJS, 188, 242

Kervella P., Arenou F., Mignard F., Thévenin F., 2019, A&A, 623, A72

Klahr H., Lin D. N. C., 2005, ApJ, 632, 1113

Konopacky Q. M. et al., 2016, ApJ, 829, L4

Kral Q., Marino S., Wyatt M. C., Kama M., Matrà L., 2019, MNRAS, 489, 3670

Kral Q., Matrà L., Wyatt M. C., Kennedy G. M., 2017, MNRAS, 469, 521

Kral Q., Wyatt M. C., Triaud A. H. M. J., Marino S., Thebault P., Shorttle O., 2018, MNRAS, 479, 2649

Krivov A. V., Booth M., 2018, MNRAS, 479, 3300

Lagrange A. M. et al., 2012, A&A, 542, A40

Lagrange A. M. et al., 2019, A&A, 621, L8

Launhardt R. et al., 2020, A&A, 635, A162

Lazzoni C. et al., 2018, A&A, 611, A43

Lindroos L. et al., 2016, MNRAS, 462, 1192

Long F. et al., 2018, ApJ, 869, 17

Long F. et al., 2019, ApJ, 882, 49

Lorén-Aguilar P., Bate M. R., 2015, MNRAS, 453, L78

Lyra W., Kuchner M., 2013, Nature, 499, 184

MacGregor M. A. et al., 2016, ApJ, 823, 79

MacGregor M. A. et al., 2017, ApJ, 842, 8

MacGregor M. A. et al., 2019, ApJ, 877, L32

Marino S., Bonsor A., Wyatt M. C., Kral Q., 2018a, MNRAS, 479, 1651

Marino S., Flock M., Henning T., Kral Q., Matrà L., Wyatt M. C., 2020, MNRAS, 492, 4409

Marino S., Wyatt M. C., Kennedy G. M., Holland W., Matrà L., Shannon A., Ivison R. J., 2017b, MNRAS, 469, 3518

Marino S., Yelverton B., Booth M., Faramaz V., Kennedy G. M., Matrà L., Wyatt M. C., 2019, MNRAS, 484, 1257

Marino S. et al., 2016, MNRAS, 460, 2933

Marino S. et al., 2017a, MNRAS, 465, 2595

Marino S. et al., 2018b, MNRAS, 479, 5423

Marois C., Zuckerman B., Konopacky Q. M., Macintosh B., Barman T., 2010, Nature, 468, 1080

Marshall J. P. et al., 2014, A&A, 565, A15

Matrà L., Marino S., Kennedy G. M., Wyatt M. C., Öberg K. I., Wilner D. J., 2018b, ApJ, 859, 72

Matrà L., Panić O., Wyatt M. C., Dent W. R. F., 2015, MNRAS, 447, 3936

Matrà L., Wilner D. J., Öberg K. I., Andrews S. M., Loomis R. A., Wyatt M. C., Dent W. R. F., 2018a, ApJ, 853, 147

Matrà L., Wyatt M. C., Wilner D. J., Dent W. R. F., Marino S., Kennedy G. M., Milli J., 2019, AJ, 157, 135

Matrà L. et al., 2017, ApJ, 842, 9

McMullin J. P., Waters B., Schiebel D., Young W., Golap K., 2007, in Shaw R. A., Hill F., Bell D. J., eds, ASP Conf. Ser. Vol. 376, Astronomical Data Analysis Software and Systems XVI. Astron. Soc. Pac., San Francisco, p. 127

Meshkat T. et al., 2017, AJ, 154, 245

Milli J. et al., 2017, A&A, 597, L2

Moro-Martín A., Malhotra R., Bryden G., Rieke G. H., Su K. Y. L., Beichman C. A., Lawler S. M., 2010, ApJ, 717, 1123

Moro-Martín A. et al., 2007, ApJ, 658, 1312

Moro-Martín A. et al., 2015, ApJ, 801, 143

Morrison S., Malhotra R., 2015, ApJ, 799, 41

Mouillet D., Larwood J. D., Papaloizou J. C. B., Lagrange A. M., 1997, MNRAS, 292, 896

Moór A., Ábrahám P., Derekas A., Kiss C., Kiss L. L., Apai D., Grady C., Henning T., 2006, ApJ, 644, 525

Moór A. et al., 2015, ApJ, 814, 42

Moór A. et al., 2017, ApJ, 849, 123

Mumma M. J., Charnley S. B., 2011, ARA&A, 49, 471

Murray C. D., Dermott S. F., 1999, Solar System Dynamics. Cambridge Univ. Press, Cambridge

Musso Barcucci A. et al., 2019, A&A, 627, A77

Mustill A. J., Wyatt M. C., 2009, MNRAS, 399, 1403

Nederlander A. et al., 2020, ApJ, submitted

Pan M., Nesvold E. R., Kuchner M. J., 2016, ApJ, 832, 81

Pearce T. D., Wyatt M. C., 2014, MNRAS, 443, 2541

Pearce T. D., Wyatt M. C., 2015, MNRAS, 453, 3329

Perrot C. et al., 2016, A&A, 590, L7

Pérez S., Casassus S., Baruteau C., Dong R., Hales A., Cieza L., 2019a, AJ, 158, 15

Pérez S., Marino S., Casassus S., Baruteau C., Zurlo A., Flores C., Chauvin G., 2019b, MNRAS, 488, 1005

Pinilla P., Flock M., Ovelar M. d. J., Birnstiel T., 2016, A&A, 596, A81

Quillen A. C., 2006, MNRAS, 372, L14

Rameau J. et al., 2016, ApJ, 822, L29

Read M. J., Wyatt M. C., Marino S., Kennedy G. M., 2018, MNRAS, 475, 4953

Regály Z., Dencs Z., Moór A., Kovács T., 2018, MNRAS, 473, 3547

Rein H., Liu S.-F., 2012, A&A, 537, A128

Rein H. et al., 2019, MNRAS, 485, 5490

Ricci L., Carpenter J. M., Fu B., Hughes A. M., Corder S., Isella A., 2015, ApJ, 798, 124

Richard A. J. W., Lyra W., Kuchner M. J., 2018, ApJ, 856, 41

Saito E., Sirono S.-i., 2011, ApJ, 728, 20

Schneider G. et al., 2014, AJ, 148, 59

Shannon A., Mustill A. J., Wyatt M., 2015, MNRAS, 448, 684

Simpson J. M. et al., 2015, ApJ, 799, 81

Smith A. W., Lissauer J. J., 2009, Icarus, 201, 381

Stolker T. et al., 2020, A&A, 635, A182

Su K. Y. L. et al., 2017, AJ, 154, 225

Takahashi S. Z., Inutsuka S.-I., 2014, ApJ, 794, 55

van Leeuwen F., 2007, A&A, 474, 653

van Lieshout R., Dominik C., Kama M., Min M., 2014, A&A, 571, A51

Visser R., van Dishoeck E. F., Black J. H., 2009, A&A, 503, 323

Vorobyov E. I., 2013, A&A, 552, A129

Vousden W. D., Farr W. M., Mandel I., 2016, MNRAS, 455, 1919

Wilner D. J., MacGregor M. A., Andrews S. M., Hughes A. M., Matthews B., Su K., 2018, ApJ, 855, 56

Wisdom J., 1980, AJ, 85, 1122

Wyatt M. C., 2008, ARA&A, 46, 339
Wyatt M. C., Booth M., Payne M. J., Churcher L. J., 2010, MNRAS, 402, 657

Wyatt M. C., Dent W. R. F., 2002, MNRAS, 334, 589

Wyatt M. C., Greaves J. S., Dent W. R. F., Coulson I. M., 2005, ApJ, 620, 492

Wyatt M. C., Kral Q., Sinclair C. A., 2020, MNRAS, 491, 782

Wyatt M. C., Smith R., Greaves J. S., Beichman C. A., Bryden G., Lisse C. M., 2007, ApJ, 658, 569

Wyatt M. C. et al., 2012, MNRAS, 424, 1206

Yelverton B., Kennedy G. M., 2018, MNRAS, 479, 2673

Yelverton B., Kennedy G. M., Su K. Y. L., 2020, MNRAS, 495, 1943

Yelverton B., Kennedy G. M., Su K. Y. L., Wyatt M. C., 2019, MNRAS, 488, 3588

Zapata L. A., Ho P. T. P., Rodríguez L. F., 2018, MNRAS, 476, 5382

Zuckerman B., Forveille T., Kastner J. H., 1995, Nature, 373, 494

Zuckerman B., Song I., 2012, ApJ, 758, 77

Zurlo A. et al., 2016, A&A, 587, A57

## APPENDIX A: CONTINUUM IMAGING

In Fig. A1, we present the continuum clean images used to compute the radial profiles in Fig. 2.

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

**Figure A1.** Continuum clean images at 0.88 mm (12 m + ACA, left-hand panel) and 1.3 mm (right-hand panel) of HD 206893 obtained using Briggs weighting and a robust parameter of 2. Additionally, we applied a UV tapering of 0.4 arcsec to the band 7 data. The images are also corrected by the primary beam; hence, the noise increases towards the edges. The contours represent two, three and five times the image rms (12 and 4.9 µJy per beam at the centre of the band 7 and 6 images, respectively). The stellar position is marked with a white plus symbol near the centre of the image (based on Gaia DR2) and the beams are represented by white ellipses in the bottom-left corners (0.57 arcsec × 0.50 arcsec and 0.70 arcsec × 0.57 arcsec, respectively).

This paper has been typeset from a TEX/LATEX file prepared by the author.