key: cord-0311411-paz0pnjk
authors: Marino, Sebastian
title: Constraining planetesimal stirring: how sharp are debris disc edges?
date: 2021-04-05
journal: nan
DOI: 10.1093/mnras/stab771
sha: ae42ef06206d3d4639fbe6898422b8e6ee9f9735
doc_id: 311411
cord_uid: paz0pnjk

The dust production in debris discs by grinding collisions of planetesimals requires their orbits to be stirred. However, stirring levels remain largely unconstrained, and consequently the stirring mechanisms as well. This work shows how the sharpness of the outer edge of discs can be used to constrain the stirring levels. Namely, the sharper the edge is the lower the eccentricity dispersion must be. For a Rayleigh distribution of eccentricities ($e$), I find that the disc surface density near the outer edge can be parametrised as $tanh[(r_{max}-r)/l_{rm out}]$, where $r_{max}$ approximates the maximum semi-major axis and $l_{rm out}$ defines the edge smoothness. If the semi-major axis distribution has sharp edges $e_mathrm{rms}$ is roughly $1.2 l_{rm out}/r_{max}$, or $e_mathrm{rms}=0.77 l_{rm out}/r_{max}$ if semi-major axes have diffused due to self-stirring. This model is fitted to ALMA data of five wide discs: HD 107146, HD 92945, HD 206893, AU Mic and HR 8799. The results show that HD 107146, HD 92945 and AU Mic have the sharpest outer edges, corresponding to $e_mathrm{rms}$ values of $0.121pm0.05$, $0.15^{+0.07}_{-0.05}$ and $0.10pm0.02$ if their discs are self-stirred, suggesting the presence of Pluto-sized objects embedded in the disc. Although these stirring values are larger than typically assumed, the radial stirring of HD 92945 is in good agreement with its vertical stirring constrained by the disc height. HD 206893 and HR~8799, on the other hand, have smooth outer edges that are indicative of scattered discs since both systems have massive inner companions.

Planetesimal discs at tens of au analogous to the Kuiper belt are a frequent component of planetary systems (Su et al. 2006; Eiroa et al. 2013; Sibthorpe et al. 2018) . Mutual collisions within these belts grind down solids producing high dust levels, which cause debris discs (see reviews by Wyatt 2008; Krivov 2010; Hughes et al. 2018 ). This collisional cascade requires planetesimal orbits to be dynamically excited or stirred by large planetesimals or planets, increasing their eccentricities and inclinations (and thus their relative velocities) enough such that mutual collisions lead to fragmentation. However, even though the disc evolution through collisional processes has been extensively studied analytically (e.g. Dominik & Decin 2003; Wyatt et al. 2007; Löhne et al. 2008; Shannon & Wu 2011; Pan & Schlichting 2012; Geiler & Krivov 2017 ) and via numerical simulations (e.g. Krivov et al. 2006; Thébault & Augereau 2007; Gáspár et al. 2012) , the stirring levels of debris discs remain largely unknown.

Constraining the current eccentricity and inclination excitation (i.e. stirring levels) of planetesimal discs is important since it provides unique information on the dynamical history of planetary ★ E-mail: sebastian.marino.estay@gmail.com systems. Proof of this is how the stirring or dynamical structure of the Kuiper belt has been used to unveil the dynamical history of Neptune (see recent review by Morbidelli & Nesvorný 2020) . For example, the Kuiper belt's resonant population is direct evidence of Neptune's radial migration in the past (Malhotra 1993 (Malhotra , 1995 , and its migration rate is well constrained by the inclination distribution of the Kuiper belt's hot population (Nesvorný 2015) . Moreover, the low inclinations of some detached and resonant Kuiper Belt objects favours scenarios in which Neptune reached an eccentricity of ∼ 0.1 during its migration (Nesvorný et al. 2021) , possibly due to a mild dynamical instability (e.g. as in the Nice model, Tsiganis et al. 2005) . Hence stirring levels can provide a very complex picture of the dynamical history of planetary systems.

Determining the stirring level of planetesimal discs in exoplanetary systems is, however, much more challenging because individual orbits cannot be observationally determined. Instead, observations typically constrain the projected surface brightness at multiple wavelengths. Because the disc surface brightness depends on the dust size distribution, orbital elements and dust temperatures, this is a degenerate problem. Nevertheless, there are special cases in which information on the orbital excitation of disc particles can be directly obtained, especially at long wavelengths that are sensitive to the emission of large grains whose dynamics are unaffected by radiation pressure. Observations with the Atacama Large Millimeter/submillimeter Array (ALMA) have been fundamental for this, and below I summarise the special cases in which ALMA observations have allowed constraining the excitation of eccentricities and/or inclinations.

Recent ALMA observations of the edge-on disc around Pic resolved its vertical extent , and revealed a vertical distribution that is inconsistent with a single Rayleigh distribution of inclinations ( ). Its vertical distribution was better reproduced with two dynamical populations with typical inclination dispersions ( rms ) of 1 • and 9 • . Such a scenario is analogous to the classical Kuiper belt, which also has a bimodal inclination distribution (e.g. Brown 2001) . Similarly, ALMA observations of AU Mic marginally resolved its vertical extent (Daley et al. 2019 ) and constrained as well the distribution of inclinations to rms ∼ 3 • .

For moderately inclined and narrow discs, the vertical extent can be constrained by comparing the surface brightness as a function of azimuth if the resolution is comparable or smaller than the disc width (Marino et al. 2016) . This is because the non-zero vertical thickness of a disc widens the surface brightness distribution near the minor axis, hence reducing slightly its surface brightness compared to a flat disc. This principle has helped to marginally constrain the vertical thickness of HR 4796's debris ring , finding rms ∼ 3 • . Similarly, it is also possible to determine the height of a wide disc if it contains narrow features such as a gap, as done by Marino et al. (2019) for HD 92945's disc, finding rms ∼ 4 • . Note that the latter case is very similar to modelling done in protoplanetary discs with gaps to measure the degree of dust settling (e.g. Pinte et al. 2016) .

The width of debris discs provides as well a direct upper limit on the eccentricity levels. This is because for a uniform distribution of longitudes of pericentres, the difference between pericentre and apocentre distance must be equal or smaller than the disc width. This is particularly interesting when applied to narrow debris discs interacting with an eccentric planet since the final width of the disc strongly depends on the forced eccentricity induced by the planet (e.g. Pearce & Wyatt 2014) . Kennedy (2020) used these arguments to show that the narrowness of the debris disc around Fomalhaut (MacGregor et al. 2017 ) and HD 202628 ) is inconsistent with the standard secular evolution theory. The dispersion of eccentricities in those discs, constrained by the disc widths, is smaller than expected under the presence of an eccentric planet. This inconsistency could be due to planetesimals being born on eccentric and aligned orbits, or their relative velocities being damped reducing the proper eccentricities (e.g. through collisions, Nesvold et al. 2013) .

The multiwavelength analysis of disc sizes can also constrain stirring levels. By combining FIR and mm spatial information with collisional evolution models, Geiler et al. (2019) found that the FIR emission of HR8799's debris disc cannot be reproduced by a dynamically cold wide belt that best matches the mm data. They concluded that to fit both the large extention of the FIR emission and high mm flux, their model required the presence of two populations: a dynamically cold and wide disc and an excited disc with eccentricities of 0.3-0.5.

Similarly, mm observations that determine the disc emission spectral index can constrain the grain size distribution (e.g. Ricci et al. 2015; MacGregor et al. 2016; Löhne 2020) . This is important since the grain size distribution carries information about the dynamical excitation of solids (e.g. Thébault & Wu 2008; Pan & Schlichting 2012) , and the strength of solids (O'Brien & Greenberg 2003; Wyatt et al. 2011) . Therefore constraining the mm spectral indices of discs could provide additional constraints on the stirring levels. However, as pointed out by Löhne (2020) , the insufficient knowledge of the composition, porosity and optical properties of solids in debris discs means that estimates of the grain size distribution suffer from large systematic uncertainties.

Finally, while observations at optical or NIR wavelengths can also provide high fidelity characterisations of the disc surface brightness in scattered light, using such observations to constrain the dynamical excitation of the disc is even more challenging. This is because such observations trace small grains at the bottom of the collisional cascade, whose dynamics are strongly affected by radiation pressure and stellar winds (e.g. Burns et al. 1979; Plavchan et al. 2005) . Such processes increase the eccentricities of small grains once released from larger bodies, and thus their dynamical stirring does not directly trace the stirring of the larger bodies. Even their vertical distribution is affected by radiation pressure and collisions (Thébault 2009) , and thus observations that constrain the height of small grains (e.g. Olofsson et al. 2020) do not necessarily trace the vertical stirring of planetesimals. In principle, dynamical models combined with radiative transfer simulations (e.g Pawellek et al. 2019a,b) could be used to retrieve the stirring levels of the parent bodies. However, so far this has proven very difficult due to degeneracies in model parameters that set the radial distribution of small grains. Furthermore, regardless of the planetesimal's stirring levels, their bright halos that extend beyond the planetesimal discs are expected to have surface brightness distributions following a universal power law index of -3.5 (Thébault & Wu 2008) , erasing most of the dynamical information of their parent bodies. Therefore, observations at longer wavelengths that trace the distribution of larger grains are better suited to constrain the stirring levels.

This paper focuses on how the eccentricity levels can be constrained by studying how smooth or sharp the outer edge of a disc looks at long wavelengths. While the outer edge smoothness depends both in the semi-major axis distribution and level of eccentricities, the smoothness places a direct upper limit on the eccentricity of particles. In §2 I present a simple model of a stirred disc that can be approximated by an analytic expression that characterises the smoothness of the inner and outer edges due to stirring. In §3 this model is used to fit ALMA disc observations and constrain their eccentricity excitation. I discuss the results in §4 and I summarise the main conclusions in §5.

This section aims to show how the surface density radial profile of a disc changes depending on the level of orbital excitation. Specifically, I will focus on the outer edge of a disc which will become smoother for higher eccentricities and inclinations of particles in a disc as they span a larger range of radii.

The disc is composed of an ensemble of particles in Keplerian orbits, with randomised orbital parameters. Each particle has a different semi-major axis ( ) that is distributed as ( ) ∼ +1 between min and max , with being the surface density slope or power law index 1 . Eccentricities ( ) and inclinations ( ) follow a Rayleigh distribution, i.e. where 2 and 2 are the mean square eccentricity ( 2 rms ) and inclination ( 2 rms ). The choice of a Rayleigh distribution is motivated by the expected orbital distribution of ensembles of interacting planetesimals (Ida & Makino 1992) . In addition, the dispersion of eccentricities and inclinations is set to rms = 2 rms , which is the expected equilibrium when particles interact in a Keplerian potential (Ida & Makino 1992; . The rest of the orbital parameters (longitude of pericentre, longitude of ascending node and mean anomaly), are considered to be uniformly distributed such that the modelled discs are axisymmetric. Figure 1 shows an example of a Rayleigh distribution for with rms = 0.03 and 0.15. The mean is approximately 0.9 rms , the mode 0.7 rms and the standard deviation 0.5 rms .

A total of 10 7 orbits are drawn, from which the surface density is computed as a function of the projected radius (hereafter referred to simply as radius). Figure 2 shows the surface density for a disc with = 0, min = 50, max = 150 and rms = 0.03 (blue) and 0.15 (orange). The solid lines correspond to the surface density computed using 1 au wide radial bins. The dashed lines correspond to the same profiles that are additionally convolved with a Gaussian kernel/beam with a full width at half maximum (FWHM) of 15 au or 10% of max . This is to mimic the response of an interferometer such as ALMA. By comparing both solid and dashed lines, some first conclusions can be drawn on the effect of particles eccentricities. The higher is, the smoother the outer edge is as particles spend time at a wider range of radii. Note that although particles have a distribution of inclinations, the smoothing is dominated by the eccentricities. This is because the radial smoothing grows as O ( ) and O ( 2 ), thus eccentricities dominate. Another effect of higher rms is that the radius at which the density starts to drop is also shifted inwards. This is because for large eccentricities the surface density at a given radius has contributions from a larger range of semi-major axes, but for close to max this range is significantly smaller since no particles have > max . Note that the inner edge also becomes smoother with higher eccentricities, but this change is less noticeable. This is because for a fixed eccentricity, the radial span of an orbit is proportional to its semi-major axis, and thus at a smaller radius the same fractional changes become less noticeable when convolved with the same beam compared with the outer edge.

I am interested in finding a parametric model that could reproduce the smoothness of the inner and outer edges, and thus to constrain the level of eccentricity by fitting such a model to disc observations. In principle, it is also possible to fit distributions of orbital elements (e.g. MacGregor et al. 2017; Kennedy 2020); however, this is computationally expensive and thus much slower since typically 10 6 orbits must be drawn to obtain smooth models without shot noise. I identify and test two potential simple models to parametrise the edges: A model in which the edges have Gaussian profiles and a model where the edges are approximated by hyperbolic tangents. Figure 3 compares the model surface density obtained from randomised orbits and the parametric models that best approximate to the surface density. The Gaussian edges (dotted lines) produce smooth edges in the surface density profiles, but are not smooth enough to reproduce the surface density near to where the density starts to drop. On the other hand, the hyperbolic tangent model reproduces better the profile across all radii for low and high rms . Therefore, hereafter I adopt this parametrisation as default, which is defined as

where controls the surface density slope within the disc, min and max approximate min and max , and in and out control the smoothness of the inner and outer edges. Below in §2.1 I show how the ratio out / max can be used to constrain rms .

With a good analytical model in hand to reproduce the smoothness of the disc edges, it is now important to test how well rms can be constrained by fitting such a model. I anticipate that the smoothness of the disc outer edge here defined as out / max should scale linearly with rms for small values. In order to test that I simulate the orbits of multiple discs with rms = 0.03, 0.05, 0.1, 0.18, 0.3, 0.5, and fit them with the model defined by Eq. 2 leaving Σ , min , min , , in , and out as free parameters. Figure 4 shows how the retrieved smoothness defined by out / max varies as a function of rms . The different colours correspond to discs of different fractional widths, here defined as the range of semi-major axes divided by the mid Best fit values for out / max vs different input values of rms when assuming a distribution of semi-major axes that has sharp edges. The colours represent different disc fractional widths in semi-major axis (i.e. width over the mid semi-major axis), with blue, orange and green representing fractional widths of 1, 0.5 and 0.2, respectively. The dashed line represents a linear relation between out / max and rms that approximates well the results, except for rms > Δ / mid . semi-major axis. I find that the derived values for out / max correspond approximately to rms /1.2, providing a simple relation to retrieve rms when applying this model to data.

Only when rms becomes comparable to the fractional width, the linear relation starts to break and 1.2 out / max underpredicts rms . In such cases the surface density is made of particles from a wide range of radii, is strongly peaked, and it resembles a Gaussian distribution. I find that for very narrow semi-major axis distributions, the FWHM of the surface density distribution is close to 1.7 mid rms (where mid is the mid semi-major axis). This relation resembles what Kennedy (2020) found when studying eccentric and narrow debris rings; namely, that the observed disc width is set by the orbital excitation of particles in the disc. Best fit values for out / max vs different input values of rms when assuming a distribution of semi-major axes that has sharp edges and a fractional width of 1. The shape of the markers represent different slopes of the surface density distribution, with circles, squares and diamonds representing power law indices of 0, -1 and 1, respectively. The dashed line represents a linear relation between out / max and rms that approximates well the results.

I also test how well rms can be retrieved when the surface density slope is different than zero and is still left as a free parameter. Figure 5 shows out / max as a function of rms for = −1, 0, 1. I find that 1.2 out / max still approximates rms well up to values of ∼ 0.5, at which point rms would be overestimated for = −1 and slightly underestimated for = 1. Therefore I conclude that the chosen parametric model is robust in retrieving the value of rms , except when the fractional width of semi-major axes becomes comparable to rms .

A caveat in the analysis presented above is that it assumes that the distribution of semi-major axes has a sharp cut, and thus any smoothness in the surface density profile is attributed primarily to the eccentricity of particles. While this might not be the case, the scenario assumed here can be used to set an upper limit in the eccentricity of particles in a disc. This is because higher eccentricities would make the outer edge appear even smoother for any distribution of semi-major axes. Therefore, the sharp cut in semi-major axes beyond max is the best case scenario to produce a sharp edge. Any other continuous distribution of would make the outer edge even smoother. Therefore this model offers a useful tool to at least set an upper limit on rms .

As stated above, so far I have assumed that the semi-major axis distribution has sharp edges; however, there are scenarios in which the edges in the semi-major axis distribution of dust or planetesimals could become smoother with time. This could be expected in a self-stirring scenario, in which the increase in eccentricities due to close encounters leads to a diffusion in semi-major axis as well (see e.g. Fig 3 in . Assuming the pericentre distance between encounters is kept constant, particles stirred to eccentricities diffuse in semi-major axis by a length ∼ /(1 − ) (e.g. Duncan et al. 1987) . If this is the case, the observed smoothness in the surface density would be a combination of both eccentricities and the smooth semi-major axis distribution.

In order to test the effect of both effects combined, new cal- Best fit values for out / max vs different input values of rms when assuming a distribution of semi-major axes that have smooth edges. The colours represent different disc fractional widths in semi-major axis (i.e. width over the mid semi-major axis), with blue, orange and green representing fractional widths of 1, 0.5 and 0.2. The dashed line represents a linear relation between out / max and rms that approximates well the results, except for broad discs with high rms values. culations similar to those shown above in §2 are performed. Now instead of having a semi-major axis distribution ( ) following a simple power law distribution, now ( ) is defined as a power law with smooth edges (assuming the same functional form as Equation 2). Namely

where the smoothness parameters in and out are set to min,max rms /(1 − rms ) at the inner and outer edge, respectively. This additional smoothing means that for the same stirring level rms , the edges of a disc are smoother than in the case considered in §2. With this additional semi-major axis smoothing I fit the same parametric model defined by Equation 2, and find that the best fit value of out / max now corresponds to approximately 1.3 rms for rms 0.2 as shown in Figure 6 . Hence the same observed level of smoothness corresponds to a lower degree of orbital excitation in this scenario. When discussing the observed smoothness levels in discs and its connection to self-stirring in §4.1, I will consider this scenario to interpret the results.

Finally, I also find that the derived values of min and max approximate well min and max (< 10% error) for rms 0.2 when assuming smooth edges in the semi-major axis distribution, and for rms 0.5 when assuming sharp edges in the semi-major axis distribution. Therefore, min and max are also good proxies for determining the semi-major axis distribution, and 2( max − min )/( max + min ) approximates the fractional width of semi-major axes.

In order to measure how well rms can be retrieved from a noisy radial profile, I produce multiple models with different resolutions (20%, 10% and 5% of max ) by convolving with a Gaussian beam, different signal-to-noise ratios (SNR) per beam 2 (15, 30 and 60), and different rms values (0.025, 0.05, 0.1, 0.2, 0.4). I fit the parametric model convolved with the same Gaussian beam to these mock profiles using an MCMC procedure (Foreman-Mackey et al. 2013 ) and estimate rms as 1.2 out / max . I find that the fractional error depends on the beam size and SNR as follows

where is the distance in pc, beam is the beam size in arcseconds and max is the maximum semi-major axis. If rms beam / max 0.1 , then rms is not well constrained but a 99.7% confidence upper limit can be set which approximates to

Therefore, observations that have resolved well the disc emission (with a few beams across) and with SNR 20, can be used to set meaningful constraints on rms . In §3 I fit the parametric model described by Equation 2 to ALMA observations of 5 discs, and based on the derived smoothness of their outer edges I constrain their eccentricity levels considering the two models presented above, i.e. sharp and smooth semi-major axes distributions.

I wish to apply the model presented above to real data to constrain the level of stirring of some planetesimal discs. The best data to do this are ALMA observations of discs that have resolved their widths with multiple resolution elements (or beams) with SNR > 10. I identify 5 disc observations as the most promising for this (Faramaz et al. submitted) . Pic would also be a good candidate since it is very well resolved with ALMA. However, it is known to be asymmetric, warped, and possibly has multiple dynamical components , and references therein). Thus it is excluded from this sample. Figure 7 shows ALMA continuum images of the five discs analysed here.

The observed visibilities of each disc are fitted by first creating synthetic images using RADMC3D (Dullemond et al. 2017 ) of the dust continuum emission. For this, a dust surface density defined as Equation 2 is input into RADMC3D, together with a stellar model spectrum consistent with previous studies of these systems. Because HD 107146, HD 92945, HD 206893 and AU Mic are known to have a gap in between a broad disc (AU Mic's disc was described instead as double belts by Daley et al. 2019) , their models include a Gaussian gap to reproduce the surface density as best as possible. In cases where multi band data are available (HD 107146, HD 206893 and HR 8799), only one image is produced with RADMC3D and the rest obtained by scaling the disc and central emission. The disc These images are not primary beam corrected, and thus the noise levels are constant across the images. HD 107146's image (top left) is obtained from 1.14 mm data after subtracting the best fit of a background source near the disc inner edge, and using Briggs weighting (robust=0.5). HD 92945's image (top middle) is obtained from 0.86 mm data using natural weights and a uvtaper of 0. 7. HD 206893's image (top right) is obtained from 1.35 mm data using Briggs weighting (robust=2.0) and a uvtaper of 0. 4. HR 8799's image (bottom left) is obtained from 0.88 mm data after subtracting the best fit of a background source near the disc inner edge, and using natural weights and a uvtaper of 0. 8. AU Mic's image (bottom middle) is obtained from 1.35 mm data using Briggs weighting (robust=0.5). The white cross in the middle of the images represents the stellar position. The ticks in all panels are spaced by 1 . A scale bar is shown at the bottom right of each panel.

emission is scaled using a disc spectral index, mm , which I leave as a free parameter. The stellar emission in the central pixel is also scaled assuming a spectral index of 2.

I further assume a vertical mass distribution that is Gaussian (consistent with a Rayleigh distribution of inclinations, Matrà et al. 2019) , with a standard deviation that scales linearly with radius as ℎ . Given the low inclination of HD 107146, HD 206893 and HR 8799, ℎ is hardly constrained and thus I simply fix it to 0.05. For HD 92945 and AU Mic, ℎ is left as a free parameter since these two discs are seen at a higher inclination and thus ℎ can be constrained Daley et al. 2019) .

Each synthetic image is then multiplied by the corresponding antenna primary beam, which is different for different wavelengths and antenna arrays (either 12m or ACA). Model visibilities are obtained based on those images using a Fast Fourier Transform algorithm, and a bilinear interpolation to extract visibilities at the same − points as the observations (as in Marino et al. 2018a ). Additionally, I consider phase center offsets for each data set, which are applied by multiplying each model visibility by exp[−2 ( ΔRA + Δ DEC )], where ΔRA and ΔDEC are the sky offsets in right ascension and declination directions, respectively. Finally the 2 is computed for each set of parameters, and this is used to feed an Affine Invariant MCMC Ensemble sampler using (Goodman & Weare 2010; Foreman-Mackey et al. 2013 ) to estimate the posterior distribution of each parameter.

It is well known that the estimated uncertainties or weights of the observed visibilities delivered by the ALMA calibration pipeline are often off by a small factor, even after using tools such as that determine the visibility dispersion with CASA (McMullin et al. 2007 ). To correct for this I multiply the weights of each data set by a fixed value to force the reduced 2 ( 2 ) to be approximately 1 as in Marino et al. (2018a) . Because the SNR of each visibility is 1, these fixed values can be found with a null model. A subtle but important detail here is that 2 is defined as 2 divided by the number of data points 3 , which in this case is twice the number of observed visibilities since each visibility has a real and imaginary component that are independent. This was not realised in Marino et al. (2018a Marino et al. ( , 2019 and thus the reported uncertainties were overestimated by a factor √ 2. To conclude, each model is fitted with at least 10 parameters 3 Formally 2 is defined as 2 divided by the degrees of freedom. In this case this is almost equal to the number of data points which are typically ∼ 10 6 for ALMA 1h integrations and moderate time and spectral averaging. that set the system orientation, phase center offsets and surface density. Additional parameters are used to fit the disc height, spectral indices, sub-millimeter galaxies (SMGs), gaps in the discs and central unresolved flux. In appendix A I describe a few particularities in the fitting procedure of each disc. The results of the fits are presented below.

Table 1 presents the retrieved fractional widths and rms , and Figure  8 shows the posterior distribution of the model surface density profiles. The best fit values for the disc parameters are presented in Table A1. The first immediate result is that the five discs are very broad with high fractional widths. This is not surprising since the sample only includes discs that are resolved with multiple beams across their width, naturally biasing this sample. However, HD 206893, HR 8799 and AU Mic have high uncertainties in their fractional widths. These large uncertainties are due to the degeneracies between and min and max . Nevertheless the true radial extent of the discs is better constrained by the model radial profiles shown in Figure 8. These clearly show the known gaps in HD 107146, HD 92945 and HD 206893, a broad and smooth disc around HR 8799, and low level emission extending down to ∼ 10 au for AU Mic supporting the conclusions of Daley et al. (2019) . Whether this inner emission is separated from the emission further out by a gap is less clear in this system.

The second immediate result is that the five discs have outer edges that are best fit by a smooth transition, which could be indicative of a significant level of stirring as shown in §2. Table 1 presents the inferred rms (median values and lower and upper limits) when assuming a semi-major axis distribution ( ( )) with sharp edges (3rd, 4th and 5th columns) or smooth edges (6th, 7th and 8th column). The sharpest surface density outer edges are found around HD 107146, HD 92945 and AU Mic, with estimated rms values of 0.18 ± 0.01, 0.23 +0.11 −0.08 and 0.16 +0.03 −0.04 , respectively, assuming sharp edges for ( ). When considering a self-stirred case, i.e. ( ) having smooth edges, the outer edges of these three discs suggest rms values of 0.121 ± 0.005, 0.15 +0.07 −0.05 and 0.10 ± 0.02. Whether these values truly correspond to the level of stirring in these three discs is uncertain, but they offer at least an upper limit to the stirring levels in these discs.

A particular test can be done by comparing the derived values of rms and ℎ since it is expected that rms ∼ 2 rms if discs are in equilibrium, and ℎ is directly related to rms . As shown by Matrà et al. (2019) for a Rayleigh distribution of eccentricities rms = √ 2ℎ, and thus rms = √ 8ℎ if these discs are in equilibrium. For HD 92945 and AU Mic, ℎ was left as a free parameter since it can be constrained given the high inclination of these discs from face on orientation. The results show that the vertical aspect ratios are constrained to values of 0.043 ± 0.014 and 0.015 ± 0.003 for these discs, respectively. Thus, given the vertical extent of HD 92945 and AU Mic, the expected values for rms are 0.12 ± 0.04 and 0.04 ± 0.008, respectively. Thus, the condition rms ∼ 2 rms is met for HD 92945, with inferred levels of vertical and radial stirring that are consistent with being in equilibrium. This is not the case for AU Mic for which a significantly larger value of radial stirring is found, which is in a 3 tension. Three reasons are identified for why this could be the case. First, in a self-stirred disc the eccentricities grow much faster than the inclinations Krivov & Booth 2018) 

found a value for ℎ (0.027 +0.004 −0.005 ) that is 30% larger than the best fit presented here, which would solve the tension. This difference on the estimated ℎ could be due to differences in the chosen parametric model. Third, it could be that its semi-major axis distribution has an outer edges that is smoother than assumed here, and thus the true value of rms is lower and closer to 0.04.

The eccentricity values found for HD 206893 and HR 8799 are much larger and this can be seen directly from their smoother outer edges in the best fit model radial profiles in Figure 8 . These two systems host massive companions that have been directly imaged, one being likely a brown dwarf at 11 au around HD 206893 (Milli et al. 2017; Delorme et al. 2017; Grandjean et al. 2019; Ward-Duong et al. 2021 ) and four giant planets around HR 8799 (Marois et al. 2008 (Marois et al. , 2010 . The formation, evolution and interaction of these companions with the disc could explain the particularly high eccentricities or extended outer edges. In fact, as proposed by Geiler et al. (2019) , HR 8799's multiwavelength disc emission is best fit as the combination of a dynamically cold and hot populations. The latter could be a scattered disc of particles having close encounters with Table 1 . Retrieved fractional width of semi-major axes, rms and rms . The fractional width is defined as 2( max − min )/( max + min ) (2 nd column). The values for rms are estimated as 1.2 out / max when assuming a semi-major axis distribution with sharp edges (3 rd , 4 th and 5 th column), or as 0.77 out / max when assuming a semi-major axis distribution with smooth edges (6 rd , 7 th and 8 th column Figure 9 summarises the derived rms of the five discs analysed here as a function of their fractional width 4 , and under the assumption of a semi-major axis distribution with smooth edges. The grey dashed line shows the maximum rms for a circular belt, which corresponds to a case of a narrow semi-major axis distribution where the the belt FWHM is set entirely by planetesimal eccentricities (Δ = 1.7 rms , §2.1). The same figure shows the derived proper eccentricity dispersion of the narrow and eccentric belts around Fomalhaut and HD 202628 (Kennedy 2020) 5 , and upper limits of rms for the narrow belts around HR 4796, Eri, HD 170773, HD 181327 and Corvi based on their fractional widths ( rms ≤ 0.6Δ / ). This comparison illustrates how the eccentricity values of wide discs derived here are lower than the most basic upper limit based on their widths, and comparable with the upper limits of some of the narrowest discs. Future ALMA observations that resolve the width of both narrow and wide discs with multiple beams could constrain the outer edge smoothness of an even larger sample.

Note that the very smooth outer edge of HR 8799's appears to be inconsistent with the maximum rms represented by the dashed line. This apparent inconsistency is mainly due to two effects. First, the eccentricity values for HR 8799 and HD 206893 are likely overestimated. Figure 6 shows that for rms > 0.2, the smoothness is greater than the linear regression. This means that the derived eccentricity for a very smooth edge will be overestimated, and thus the values derived for HR 8799 and HD 206893 must be taken with caution. Second, the width of HR 8799's disc ( max − min ) does not correspond necessarily to its FWHM. This is especially true for values different from zero, in which case max − min approximates 4 The derived fractional width of these five discs does not necessarily correspond to the fractional width of the surface density distribution defined by its FWHM, especially for cases with different from zero and high eccentricities. The fractional width defined by 2( max − min )/( max + min ) better represents the fractional width of semi-major axes. 5 In a disc that is eccentric due to a forced eccentricity, the dispersion of eccentricities that set the relative velocities depends both on the mean and standard deviation of proper eccentricities. Since in Kennedy (2020) The rms values for HD 107146, HD 92945, HD 206893, HR 8799 and AU Mic are based on the their outer edges and assume a smooth edge in their semi-major axis distribution (i.e. self-stirred scenario). The rms values for Fomalhaut and HD 202628 correspond to the dispersion of proper eccentricities based on fitting orbital distributions fits to the observed surface brightness. The rms upper limits for HR 4796, Eri, HD 170773, HD 181327 and Corvi are set equal to half their fractional widths, i.e. rms ≤ 0.6Δ / (grey dashed line). This latter limit is only strictly valid for a circular discs. The coloured dots and ellipses represent the medians and 1 confidence intervals based on the posterior distributions. Note that the ellipse of HD 107146 is smaller than the blue dot and this is not visible. better the width of the semi-major axis distribution rather than the FWHM of the surface density. In fact, I find that and the disc width are anticorrelated in the posterior distribution. Hence the values of high rms and low fractional width correspond to high values of , for which max − min underestimates the disc FWHM of the surface density distribution. From the posterior distribution of the surface density I find that the fractional width of HR 8799 based on its FWHM is constrained between 0.7 and 1.1, i.e. larger than the value used in this figure and reported in Table 1 . Therefore the results for HR 8799's are still consistent with the limit shown by the dashed line.

As expected, for HD 107146, HD 92945 and HD 206893 the results confirm the presence of gaps centred at roughly 73 au. Despite the strong similarity in the gap radii, their widths and depths seem to be significantly different. HD 107146's gap is possibly the shallowest (∼50% empty compared with the peak surface density) and the widest. In fact, the gap FWHM converges to large values (∼ 90 au) that are comparable to the radial extent of the disc. HD 206893's gap seems to be the deepest and likely fully emptied of material, while HD 92945's gap is possibly the narrowest (28 au). HR 8799's disc does not show any evidence of a gap; nevertheless, it is important to note that it does not cover the range of radii where the other three gaps are located. In fact, HR 8799 b is on an orbit at 70 au (Wang et al. 2018 ) that coincides with the gap centers. Hence, the absence of a gap centred at ∼ 70 au in HR 8799's disc could be simply due to the three inner planets that would have cleared the disc inner regions from 10-50 au. On the other hand, AU Mic's disc does not extend beyond ∼ 50 au, but a possible gap or inflection point in the surface density profile is found at 19 au. This is consistent with the findings by Daley et al. (2019) of the need of emission at around 11 au in addition to the main belt further out.

The smoothness of the discs inner edges is less constrained as anticipated in 2. With the exception of HR 8799, the inner edges are all consistent with being sharp (i.e. in / min 1) and thus potentially truncated. Note that a sharp inner edge could still be consistent with high eccentricity particles further out if these are restricted by a minimum pericentre. In a scenario where the disc is truncated by an inner planet, the sharpness of the inner edge is closely related to the planet mass and its eccentricity (Quillen 2006; Chiang et al. 2009 Dong et al. 2020) . Hence an inner planet could be inferred by characterising how sharp a disc inner edge is (e.g. Matrà et al. 2020) . Alternatively, the observed location of the disc inner edge could be simply due to collisional evolution, i.e. it marks the radius at which the collisional timescale of the largest planetesimals is equal to the age of the system (Kennedy & Wyatt 2010; Schüppler et al. 2016; Marino et al. 2017a,b; Geiler & Krivov 2017) . Interior to this radius, the disc is strongly depleted due to collisions and is expected to have a surface density that scales with radius roughly as 7/3 (Kennedy & Wyatt 2010) . Future work is needed that compare both scenarios and assess which one resembles better the observed disc inner edges.

Finally, the surface density slopes are all consistent with a flat surface density of mm-sized grains, except for AU Mic's disc that is best fit with an increasing surface density from 10 to 40 au. Flat surface densities for small and mm-sized grains are expected in broad discs with large planetesimals that have collisional lifetimes greater than the system ages (Schüppler et al. 2016; Marino et al. 2017b; Geiler & Krivov 2017) . For standard strength properties of solids, Marino et al. (2017b) found that the surface density slope of mmsized grains (i.e. ) is 0.6 + 0.9, where is the primordial surface density slope of planetesimals. The increasing surface density slope of mm-sized grains in AU Mic could indicate that the primordial surface density slope of planetesimals increases with radius and has a slope of ∼ 1.5. Alternatively, it could be that largest planetesimals are already in collisional equilibrium and thus the surface density follows the 7/3 profile described above.

In this section I discuss the potential origins for the inferred eccentricities, its implications, and some of the caveats and limitations of the modelling approach presented in this paper.

Here I consider the possibility that the inferred levels of stirring in HD 107146, HD 92945 and AU Mic are due to dynamical excitations between large planetesimals or protoplanets within the disc. It is well established that planetesimals can stir themselves to levels in which the velocity dispersion is comparable to the escape velocity of the largest planetesimals in the disc. Once that level is reached, collisions tend to dominate over gravitational stirring, halting the stirring (e.g. Safronov 1972; Goldreich et al. 2004; Schlichting 2014) . Therefore, the inferred stirring levels can be used to find a lower limit on the escape velocity of the most massive bodies embedded in the disc. For a Rayleigh distribution, the typical relative velocities between planetesimals is √︃ 1.25 2 rms + 2 rms Kep , where Kep is the local Keplerian velocity for a circular orbit (Lissauer 1993) . The escape velocity for a spherical object of Diameter and bulk density is √︁ 2 2 /3. Equating both, an expression for the minimum size of the largest planetesimal responsible for stirring is obtained, 

where ★ is the stellar mass. Consequently, the values found for HD 107146, HD 92945 and AU Mic would imply planetesimal diameters as large as 600, 400 and 500 km, respectively (using the 5% percentiles of rms in Table 1 , max values in Table A1 , and stellar masses of 1.0, 0.9 and 0.5 , respectively). These diameters are equivalent to masses of 2%, 0.5% and 0.8% of Pluto's mass.

The next simple condition that has to be met is that the level of stirring could have been reached within the age of the system. Following the numerical results of and the more recent extention of these to a planetesimal size distribution by Krivov & Booth (2018) , the evolution of rms as a function of time due to viscous stirring in the dispersion-dominated regime is given by

where Ω is the Keplerian frequency, Σ the disc surface density, max the maximum planetesimal size, and ≈ 40 a numerical factor. The factor takes a value between 1 and 2 depending on the eccentricity levels of the largest planetesimals. A value of unity corresponds to stirrers with zero eccentricities and inclinations, and a value of 2 to a uniform level of stirring for all solids in the disc. Equation 8 is valid for a planetesimal mass distribution / ∝ − . By inverting Equation 8, an expression for the minimum disc surface density to reach the level of stirring within . Surface density vs diameter of planetesimals that are responsible for stirring the discs around HD 107146 (blue), HD 92945 (orange) and AU Mic (purple). The solid lines represent the minimum surface densities to stir the discs within the maximum age of these systems. The dashed line represents the maximum surface densities that could survive against collisions for the maximum estimated age of these systems and minimum inferred level of stirring. The blue, orange and purple shaded areas are ruled out since the discs would be stirred to high levels that are inconsistent with the data. The horizontal dotted lines represent the MMSN level at the outer edge of these discs. the age of the system is obtained 

assuming = 1.6 as in Krivov & Booth (2018) . Figure 10 shows the minimum surface density (solid lines) for HD 107146, HD 92945 and AU Mic as a function of maximum planetesimal diameter to reach rms levels of 0.11, 0.07 and 0.06 at the disc outer edges in a timescale of 200, 250 and 21 Myr, respectively. These timescales correspond roughly to their age upper limits (Miret-Roig et al. 2020; Mesa et al. 2021) . Only values above min are shown. These curves show that HD 107146 and AU Mic require the largest surface densities to explain the inferred eccentricity levels at the disc outer edges. Similar to Matrà et al. (2019) , these can be compared with the maximum surface density that could have survived collisions (assuming the same level of stirring for the largest planetesimals) that is represented with dashed lines. This condition does not add further constraints since the minimum surface density is below the maximum surface density curve for maximum planetesimal diameters above min . The horizontal dotted lines represent the Minimum Mass Solar Nebula (MMSN, that is equivalent to ( /1 au) −3/2 ⊕ au −2 , Weidenschilling 1977; Hayashi 1981) at the given disc outer radii of these three systems. HD 107146 and AU Mic require similar surface density levels, although at different radii, while HD92945 requires a surface density that is an order of magnitude lower. Nevertheless, for a maximum planetesimal diameter of ∼ 1000 km, the three surface densities are roughly consistent with the derived dust masses when assuming a size distribution with a power law index of -3.5. Therefore I conclude that it is plausible that these discs are self-stirred if the stirrers have a surface density close to a MMSN and masses similar to Pluto or larger. Less massive stirrers would require much larger surface densities, which translate to disc masses in the range 100-1000 ⊕ .

How and when those massive stirrers formed is an open question. Simulations by Kenyon & Bromley (2008) show that forming such massive bodies through a gas-free planetesimal disc at large radii is possible, but it can take more than 1 Gyr at these distances, thus they would have had to form during the gas rich protoplanetary disc phase. Note that for the high surface densities considered here, the condition of overlapping feeding zones of the largest planetesimals is readily met, i.e. there would be enough stirrers across the disc.

Finally, it is still possible that these discs are self-stirred, but their semi-major axis distributions are even smoother than considered. In this case the upper limits I derived (8th column in Table 1) can be used to set an upper limit on the planetesimal surface density as a function of the maximum planetesimal diameter. The shaded regions in 10 represent the surface densities and diameters that are ruled out. In these regions the discs would be stirred to levels that are inconsistent with the smoothness of their outer edges in timescales shorter than the age lower limits of these systems (120, 120 and 16 Myr for HD 107146, HD 92945 and AU Mic, respectively). In other words, the discs would appear smoother than observed.

As shown in §4.1, stirring by planetesimals is plausible; however, it requires very high surface densities or Pluto-sized bodies embedded in the disc. An alternative scenario is that the smoothness of the outer edges is due to interactions with planets. These could either stir the whole disc or create additional hot dynamical populations of particles in analogy to the scattered and hot populations in the Kuiper belt (see Morbidelli & Nesvorný 2020 , for a recent review). I identify multiple possibilities that could raise the eccentricities and these are discussed below.

First, a planet formed close to the initial disc inner edge could have undergone planetesimal-driven outward migration while scattering planetesimals onto high eccentricity orbits with larger semi-major axes, similar to the proposed Neptune migration into the Kuiper belt (e.g. Malhotra 1993 Malhotra , 1995 Nesvorný 2015) . This scattered population of planetesimals would be characterised by pericentres close to the disc inner edge, or further out if they become detached (as in the detached or fossilised scattered population in the Kuiper belt). This population would also have apocentres greater than the inferred disc outer edges such that their distribution smooths the surface density in the outer regions. This scenario could be consistent with HR 8799's properties, especially since a scattered population could be caused by the four known giant planets as they cleared the regions near their current orbits and possibly an even wider region if they migrated to their current orbits (Read et al. 2018; Goździewski & Migaszewski 2018 . This idea is consistent with the two population model proposed by Geiler et al. (2019) to explain Herschel FIR and ALMA mm wavelength data.

Alternatively, a massive planet in the middle of a disc could clear a gap as seen in HD 107146, HD 92945, HD 206893 and potentially AU Mic by scattering material forming a scattered component. However, N-body simulations by Marino et al. (2018a, see their Fig. 5 ) did not produce a noticeably smooth outer edge with planet masses below 90 ⊕ . Those simulations nevertheless ignored disc self-gravity, planet migration or more massive planets, and thus this scenario cannot be ruled-out yet. Note that most of the literature on scattered components has focused on the Solar System, and more general studies are needed to study how a scattered component could smooth the disc outer edge and the effect of collisions (as in Wyatt et al. 2010) . Lastly, it is also possible that a low mass planet formed near the disc outer edge could undergo planetesimaldriven inward migration, crossing and stirring the entire disc within the age of these systems (Kirsh et al. 2009 ).

Second, a planet embedded in the disc or at its inner edge could have trapped planetesimals in resonances and increase their eccentricities while it was migrating (e.g. Levison & Morbidelli 2003; Wyatt 2003) . However, for a wide disc a planet near the inner edge would hardly excite the eccentricity of particles near the outer edge via mean-motion resonances, and scattering is more likely as also concluded by Matrà et al. (2019) for the case of Pic.

Third, in analogy to the Nice model Gomes et al. 2005; Nesvorný & Morbidelli 2012) , a dynamical instability between inner planets in these systems could have raise their eccentricities and thus excited the orbits of planetesimals via secular interactions (Mustill & Wyatt 2009 ). This scenario would smooth the disc outer edges and could also produce half emptied gaps (Pearce & Wyatt 2015) . However, this scenario is less likely to explain HD 107146's axisymmetric gap and HD 206893's fully emptied gap.

Finally, secular resonances could also excite the eccentricity of solids in the disc at specific locations if one or multiple planets reside interior to the disc (Zheng et al. 2017; Sefilian et al. 2020) . The location of these resonances can shift in time as the disc losses mass due to collisional evolution or as planets migrate, thus exciting the eccentricity of planetesimals over a wide range of semi-major axis.

Therefore, I conclude that a single or multiple planets could be responsible for the inferred excitation levels via scattering, an instability and subsequent secular interactions with the disc, or sweeping secular resonances.

Throughout this paper I have assumed that the observed smoothness of the disc outer edges is due to the orbital eccentricities and that the distribution of semi-major axes has a sharp drop or at least sharper than the observed levels. If this is truly the case depends strongly on how planetesimals are formed and how this varies as a function of radius. If planetesimals at tens of au are formed via streaming instability (Youdin & Goodman 2005; Johansen et al. 2007 ) triggered in radial dust traps in protoplanetary disks (e.g. Pinilla et al. 2012 Pinilla et al. , 2020 Dullemond et al. 2018; Stammler et al. 2019) , the resulting planetesimal disc could have a smooth or sharp outer edge depending on how radially concentrated were dust particles and how those traps evolve in time (e.g. Lenz et al. 2019; Shibaike & Alibert 2020, Miller et al. in prep) . Alternatively, streaming instability could be triggered while the gas surface density drops due to photoevaporation (Throop & Bally 2005; Carrera et al. 2017; Ercolano et al. 2017) , in which case the planetesimal distribution could have a steep outer edge if the dust distribution had also a steep outer edge. However, more recent models of dust evolution in photoevaporating discs show that the dust-to-gas ratio is weakly affected due to efficient radial drift (e.g. Sellek et al. 2020) , and thus photoevaporation might not trigger the streaming instability. Therefore, dust traps might be necessary to overcome radial drift even when considering photoevaporation. In §2.1 I argue that, despite the fact that the smoothness of the surface density could be a consequence of either the semi-major axis or eccentricity distributions (or both), it can nevertheless set an upper limit on the eccentricity distribution. Similarly, by assuming eccentricities are zero the derived smoothness values constrain the maximum smoothness in the semi-major axis distribution. In this way, the smoothness could be used to test and constrain planetesimal formation models in the future.

Another strong assumption through this paper is that the distributions of eccentricities and semi-major axes are independent and that the former follows a Rayleigh distribution. There are multiple scenarios in which this is not the case. For example, in a disc of particles being scattered by an inner planet the eccentricities are strongly correlated with the semi-major axis of particles, with eccentricities approximating 1 − / , where is the common pericenter near the planet's orbit. In this case the disc outer edge should follow approximately a power law distribution with a power law index of -3 (Duncan et al. 1987; Marino et al. 2018b) , although collisional evolution tends to flatten the surface density distribution over time . In order to test what would be the result of fitting the parametric model described by Equation 2 to a scattered disc, I fit this model to several simulated discs with power law outer edges declining as −3 . I find that out / max converges to values > 0.5 for wide discs and values closer to 1 for narrower discs, thus consistent with the results for HD 206893 and HR 8799. Note that this is not proof of those discs having a scattered component; nevertheless, it shows that the results of very smooth edges are consistent with that possibility. Future work comparing scattered disc and alternative models to the data could provide more support to this hypothesis.

Another possibility is that the eccentricities do follow a Rayleigh distribution, but the level of dynamical excitation or rms varies as a function of semi-major axis. This might be the most natural outcome of self-stirring or stirring via secular interactions with a planet. In a self-stirring scenario, rms should vary as ( 1/2 Σ max ) 1/4 (see Equation 8 ). Thus if either Σ or max decreased with semi-major axis more steeply than −1/2 (e.g. Kenyon & Bromley 2008; Lenz et al. 2019; Klahr & Schreiber 2020) , rms should decrease with radius. Conversely, if stirring has reached its equilibrium (i.e. the relative velocities are roughly the same as the escape velocity of the largest planetesimals), the eccentricities could increase with radius for a constant maximum planetesimal size (see Equation 6 ). If eccentricities are instead excited via secular interactions with an inner planet, rms should decrease with semi-major axis as 1/ (e.g. Mustill & Wyatt 2009 ). In either of these cases, the estimates of rms are still valid near the outer edge of the disc.

Finally, evidence of high levels of stirring could also exist in other systems such as HD 181327, HD 61005 and HD 32297. Those systems were found to posses halos or extended components of mm-grains beyond their main belt (Marino et al. 2016; MacGregor et al. 2018) . Such halos could be indicative of scattered discs. While HD 61005 and HD 32297 are seen edge-on and thus their radial profiles are hard to disentangle directly from the images, HD 181327 is seen close to face-on. Its main belt is narrow (width of 21 au) and centered at 80 au (after correcting by the Gaia DR3 distance of 48 pc, Gaia Collaboration et al. 2020), but it displays low level emission out to ∼ 200 au (Marino et al. 2016) . If this halo was produced by particles with pericentres at the belt inner edge (70 au) and apocentres as large as 200 au, that would imply eccentricities as high as 0.5. On the other hand, HD 61005 and HD 32297 seemed to have smooth outer edges instead of an instant drop; however, MacGregor et al. (2018) found their outer edges had power law slope values of ∼ −6, significantly steeper than expected for a simple scattered disc. Further observations and modelling are necessary to assess if these halos are consistent with scattered discs or not.

In this paper, I have shown how the stirring levels of a planetesimal disc can influence its appearance ( §2). Namely, the higher the dispersion of eccentricities ( ) is, the smoother the disc surface density will be. This effect is most noticeable at the disc outer edge since features are smoothed by a length scale ∼ rms if eccentricities levels are uniform across the disc. I found that a sharp edge in the distribution of semi-major axes, translates to a smooth edge in the surface density radial profile that bears strong similarities with a hyperbolic tangent function (Equation 2).

By parametrising the disc surface density as a power law distribution with edges following hyperbolic tangents, the eccentricity levels in a disc can be retrieved. In §2.1 I found that when rms is smaller than the disc fractional width, rms is well approximated by ≈ 1.2 out / max , where out and max are parameters that determine the smoothness length and location of the disc outer edge, respectively. These relations assume, however, that the semi-major axis distribution has sharp edges, which might be unrealistic. It nevertheless represents a limit case in which the smoothness is only due to stirring. A smoother semi-major axis distribution would result in even smoother edges, therefore 1.2 out / max provides an upper limit on rms near the disc outer edge.

There are particular scenarios in which the semi-major axis distribution could become smoother. Such is the case of a selfstirred disc, where the increase in eccentricities due to close encounters with massive planetesimals should also lead to a diffusion in semi-major axis. I explored this case in §2.2, finding that rms ≈ 0.77 out / max , i.e. the necessary eccentricity levels to explain the outer edge smoothness are lower.

In §3 these findings are applied to real observations by fitting the parametric model to the ALMA data of five discs that are well resolved: HD 107146, HD 92945, HD 206893, HR 8799 and AU Mic. The five discs require out / max > 0, i.e. the shape of their outer edges is at least marginally resolved. I find that HD 107146, HD 92945 and AU Mic require lower eccentricities (i.e. have sharper outer edges) compared with HD 206893 and HR 8799. Assuming a self-stirring scenario, the former group requires rms values of 0.121 ± 0.05, 0.15 +0.07 −0.05 and 0.10 ± 0.02, respectively. Interestingly, the value of rms and scale height derived for HD 92945 are consistent with a scenario in which the radial and vertical stirring are in equilibrium. On the other hand, HD 206893 and HR 8799 require much smoother outer edges with rms values ∼ 0.5 if their semi-major axis distribution had sharp edges. Such high eccentricity levels would be more consistent with stronger interactions with massive planets. It is perhaps not a coincidence that these systems also host massive inner companions. Such companions could have migrated through planetesimal driven migration in the past, scattering a large population of planetesimals raising their eccentricities and semi-major axes.

Using the derived stirring levels, in §4.1 I explored the possibility that the discs around HD 107146, HD 92945 and AU Mic are self-stirred. I found that to be able to stir the disc to the inferred levels, planetesimals of at least 400 km in diameter should be present. Given the age of these systems, I also constrained how high should the planetesimal surface density be to reach the observed levels within the age of these systems. I found that if planetesimals are smaller than Pluto the required surface densities are significantly higher than a MMSN and vice versa. While this finding favours the presence of Pluto-size objects in the disc, the formation of such massive bodies through collisions of smaller planetesimals is too slow to explain their presence. This problem could be overcome if these massive bodies are formed already in protoplanetary discs.

Alternatively, the discs could have been stirred by planets ( §4.2). I identified four mechanisms that could stir the disc smoothing the outer edges, and conclude that scattering, secular interactions and secular resonances are the most promising to stir these discs via planet-disc interactions. While only two of these five systems have known companions (HD 206893 and HR 8799), the known gaps in the rest suggests the presence of unseen lower mass companions interacting with the discs. The presence of more massive companions around HD 206893 and HR 8799 could lead to a higher excitation of their discs, explaining their very smooth outer edges. In fact, their smoothness levels derived from the fits are consistent with typical scattered disc profiles.

Finally, I hope the work presented here serves as a future guide on how characterize how sharp or smooth are the edges of planetesimal discs. While these edges are set by both the semimajor axis and eccentricity distribution and thus it is a degenerate problem, the smoothness of the surface density provides an upper limit for both the smoothness of the semi-major axis distribution and dispersion of eccentricities. Coupling this analysis with vertically resolved observations could significantly help to lift degeneracies and test the different assumptions made in this work.

The best fit values for the disc parameters of the five systems are presented in Table A1 . Below I describe the data used in the analysis and a few particularities in the models used to fit each system.

The ALMA data of HD 107146 used in this paper is described in Marino et al. (2018a) , and corresponds to band 7 observations as part of the project 2016.1.00104.S (PI: S. Marino), and band 6 observations as part of the project 2016.1.00195.S (PI: J. Carpenter). In addition, I use new data in band 7 from the project 2019.1.00189.S (PI: S. Marino). The new data were taken in December 2019 and corresponds to a single compact antenna configuration (resolution of ∼ 1 ), thus providing marginal improvement over previous observations. The data will be described in more detail in a future paper after completion of the observations, which will include higher resolution observations to be taken once ALMA observations resume after the 2020 shutdown due to the COVID-19 pandemic.

The new band 7 data revealed that the clump identified in Marino et al. (2018a) that could have corresponded to its possible warm asteroid belt (Morales et al. 2011; , is not co-moving with HD 107146, and thus it is most likely a background galaxy. In order to account for this background source in the modelling, I add to the model images an additional non co-moving component. This additional component has a surface brightness distribution parametrised as a 2D elliptical Gaussian and with a relative position to HD 107146 that changes according to HD 107146's proper motion and epoch of observation. This background source is fit in parallel with the disc emission. I find that the background source is best fit by an offset from HD 107146 in December 2019 of −0. 01 ± 0. 03 in RA and 0. 09 ± 0. 02 in Dec. Its flux at 0.88 mm is 0.79 ± 0.07 mJy and has a spectral index of 3.2 ± 0.2, hence consistent with a background sub-mm galaxy. Its standard deviation along its major axis is 0. 64 ± 0. 05 and along its minor axis 0. 53 ± 0. 04 (PA of 94 • ± 12 • ), thus it is roughly circular.

The ALMA data of HD 92945 used in this paper is described in Marino et al. (2019) , and corresponds to band 7 observations as part of the project 2016.1.00104.S (PI: S. Marino). Because the star is marginally detected, I leave its flux as a free parameter and find a best fit value of 35 ± 15 Jy.

The ALMA data of HD 206893 used in this paper corresponds to band 7 data described in Marino et al. (2020) as part of the project 2017.1.00828.S (PI: A. Zurlo), and band 6 data described in (Nederlander et al. 2021 ) as part of the project 2017.1.00825.S (PI: A. M. Hughes). Since this system has a massive companion on an orbit near 10 au that would have cleared any debris near its orbit, the surface density is forced to zero within 14 au. This is roughly the minimum outer extent of the companion chaotic zone if on a circular orbit.

The ALMA data of HR 8799 used in this paper corresponds to band 6 data described in Booth et al. (2016) as part of the project 2012.1.00482.S (PI:A. Jordan), and band 7 data described in Faramaz et al. (submitted) as part of the project 2016.1.00907.S (PI: V. Faramaz). Channels with CO emission in both data sets were flagged. Since this system has massive companions interior to the disc that would have cleared any debris, the surface density is forced to zero within 70 au. This is roughly the semi-major axis of the outer most planet. Because HR 8799 also presents a bright background galaxy that overlays with the disc (Faramaz et al. submitted) , I include an additional component in the modelling using the same approach as for HD 107146. The background source is best fit by an offset from HR 8799 in June 2018 of −1. 28 ± 0. 05 in RA and 2. 34 ± 0. 05 in Dec. Its flux at 0.88 mm is 0.31 ± 0.02 mJy and 0.06 ± 0.02 mJy at 1.3 mm. The central flux (stellar photosphere plus any other additional emission) is also left as a free parameter, finding a best fit of 71 ± 10 Jy at 0.88 mm and 39 ± 16 Jy at 1.3mm. Table A1 . Best fit values of the disc parameters of HD 107146, HD 92945, HD 206893, AU Mic and HR 8799. The best fit values correspond to the median, with uncertainties based on the 16th and 84th percentiles of the marginalised distributions. The 2 nd column corresponds to the dust mass. The 3 rd -7 th columns correspond to the parameters that define the surface density according to Equation 2. The 8 th , 9 th and 10 th columns correspond to the gap central radius, its FWHM and its depth (with a value of 1 meaning a completely empty gap at its center), respectively. The 11 th and 12 th column show the disc inclination from a face-on orientation and its position angle. The 13 th column displays the disc vertical aspect ratio (constant across the disc), assuming a Gaussian vertical mass distribution with a standard deviation of ℎ . The final column shows the spectral index of the disc emission. This parameter was fixed to 0.05 for HD 107146, HD 206893 and HR 8799, and left as a free parameter for HD 92945 and AU Mic. Figure A1 . Radial profile of the azimuthally averaged residuals after subtracting the best fit models. The y axis represents the residuals in units of the local uncertainty after azimuthal averaging, while the x axis represents the deprojected radius. Since AU Mic's disc is edge-on, the profile is obtained simply as a cut through the disc major axis. The spacing between the dots represent the beam major axis for each image.

The ALMA data of AU Mic used in this paper corresponds to band 6 data described in Daley et al. (2019) as part of the project 2012.1.00198.S (PI: A. M. Hughes). I use a reduced data provided by L. Matrà which will be presented as part of the REASONS survey (Sepulveda et al. 2019, Matrà et al. in prep) . This data set was flagged at epochs when strong flares were present, and its three individual observations re-centred to the stellar position via visibility modelling. Because the star is strongly detected, I leave its flux as a free parameter and find a best fit value of 0.26 ± 0.02 mJy. Figure A1 shows the radial profiles of the deprojected and azimuthally averaged residual maps in units of the local uncertainty. These maps correspond dirty maps of the residuals, obtained with CASA using Briggs weighting and a robust parameter of 2.0. For HD 206893 and HR 8799 band 7 data, an additional uvtaper of 0. 4 and 0. 8 was used, respectively, to increase the sensitivity per beam. There are no residuals above or below 3 in any of the profiles, meaning the best fits are able to reproduce all the relevant radial features in the data. Since AU Mic's disc is edge on, deprojecting its disc emission is not possible without modelling. Thus the presented residual radial profile for this source is simply a cut through the disc major axis. This paper has been typeset from a T E X/L A T E X file prepared by the author.

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

Progress of Theoretical Physics Supplement

Astronomical Society of the Pacific Conference Series

Kuiper belt: formation and evolution

A12 Pan M

Evolution of the protoplanetary cloud and formation of the earth and planets

I would like to thank the anonymous referee for a very constructive review that improved the quality and clarity of this paper. I would also like to thank Luca Matrà for providing AU Mic's calibrated data set, Virginie Faramaz and Mark Booth for providing HR 8799 calibrated data sets, and Alexander Krivov for noticing an error in equation 9 that has been corrected before publication. 

The data underlying this article will be shared on reasonable request to the corresponding author. The ALMA data is publicly available and can be queried and downloaded directly from the ALMA archive at https://almascience.nrao.edu/asax/.